You may like using Statistics for Time Series, and Testing Correlation JavaScript.
One of the simplest versions of the theorem says that if is a random sample of size n (say, n larger than 30) from an infinite population, finite standard deviation , then the standardized sample mean converges to a standard normal distribution or, equivalently, the sample mean approaches a normal distribution with mean equal to the population mean and standard deviation equal to standard deviation of the population divided by the square root of sample size n. In applications of the central limit theorem to practical problems in statistical inference, however, statisticians are more interested in how closely the approximate distribution of the sample mean follows a normal distribution for finite sample sizes, than the limiting distribution itself. Sufficiently close agreement with a normal distribution allows statisticians to use normal theory for making inferences about population parameters (such as the mean ) using the sample mean, irrespective of the actual form of the parent population.
It is well known that whatever the parent population is, the standardized variable will have a distribution with a mean 0 and standard deviation 1 under random sampling. Moreover, if the parent population is normal, then it is distributed exactly as a standard normal variable for any positive integer n. The central limit theorem states the remarkable result that, even when the parent population is non-normal, the standardized variable is approximately normal if the sample size is large enough (say > 30). It is generally not possible to state conditions under which the approximation given by the central limit theorem works and what sample sizes are needed before the approximation becomes good enough. As a general guideline, statisticians have used the prescription that if the parent distribution is symmetric and relatively short-tailed, then the sample mean reaches approximate normality for smaller samples than if the parent population is skewed or long-tailed.
In this lesson, we will study the behavior of the mean of samples of different sizes drawn from a variety of parent populations. Examining sampling distributions of sample means computed from samples of different sizes drawn from a variety of distributions, allow us to gain some insight into the behavior of the sample mean under those specific conditions as well as examine the validity of the guidelines mentioned above for using the central limit theorem in practice.
Under certain conditions, in large samples, the sampling distribution of the sample mean can be approximated by a normal distribution. The sample size needed for the approximation to be adequate depends strongly on the shape of the parent distribution. Symmetry (or lack thereof) is particularly important. For a symmetric parent distribution, even if very different from the shape of a normal distribution, an adequate approximation can be obtained with small samples (e.g., 10 or 12 for the uniform distribution). For symmetric short-tailed parent distributions, the sample mean reaches approximate normality for smaller samples than if the parent population is skewed and long-tailed. In some extreme cases (e.g. binomial) samples sizes far exceeding the typical guidelines (e.g., 30) are needed for an adequate approximation. For some distributions without first and second moments (e.g., Cauchy), the central limit theorem does not hold.
Realize that fitting the "best'' line by eye is difficult, especially when there is a lot of residual variability in the data.
Know that there is a simple connection between the numerical coefficients in the regression equation and the slope and intercept of regression line.
Know that a single summary statistic like a correlation coefficient does not tell the whole story. A scatter plot is an essential complement to examining the relationship between the two variables.
Thus, when the variability that we predict (between the two groups) is much greater than the variability we don't predict (within each group) then we will conclude that our treatments produce different results.
Exponential distribution gives distribution of time between independent events occurring at a constant rate. Its density function is:
where l is the average number of events per unit of time, which is a positive number.
The mean and the variance of the random variable t (time between events) are 1/ l, and 1/l2, respectively.
Applications include probabilistic assessment of the time between arrival of patients to the emergency room of a hospital, and arrival of ships to a particular port.
Comments: Special case of both Weibull and gamma distributions.
You may like using Exponential Applet to perform your computations.
You may like using the following Lilliefors Test for Exponentially to perform the goodness-of-fit test.
Poisson process are often used, for example in quality control, reliability, insurance claim, incoming number of telephone calls, and queuing theory.
An Application: One of the most useful applications of the Poisson Process is in the field of queuing theory. In many situations where queues occur it has been shown that the number of people joining the queue in a given time period follows the Poisson model. For example, if the rate of arrivals to an emergency room is l per unit of time period (say 1 hr), then:
The mean and variance of random variable n are both l . However if the mean and variance of a random variable having equal numerical values, then it is not necessary that its distribution is a Poisson.
Applications:
and so on. In general:
You may like using Poisson Applet to perform your computations.
Replace the numerical example data with your up-to-14 pairs of Observed values & their frequencies, and then click the Calculate button. Blank boxes are not included in the calculations.
In entering your data to move from cell to cell in the data-matrix use the Tab key not arrow or enter keys.
Example: Used to generate random numbers in sampling and Monte Carlo simulation.
Comments: Special case of beta distribution.
The mass function of geometric mean of n independent uniforms [0,1] is:
P(X = x) = n x(n - 1) (Log[1/xn])(n -1) / (n - 1)!.
zL = [UL-(1-U)L] / L is said to have Tukey's symmetrical l-distribution.
You may like using Uniform Applet to perform your computations.
For more SPSS programs useful to simulation input/output analysis, visit Data Analysis Routines.
Computer programs that generate "random" numbers use an algorithm. That means if you know the algorithm and the seedvalues you can predict what numbers will result. Because you can predict the numbers they are not truly random - they are pseudorandom. For statistical purposes "good" pseudorandom numbers generators are good enough.
A FORTRAN code for a generator of uniform random numbers on [0,1]. RANECU is multiplicative linear congruential generator suitable for a 16-bit platform. It combines three simple generators, and has a period exceeding 81012.
It is constructed for more efficient use by providing for a sequence of such numbers, LEN in total, to be returned in a single call. A set of three non-zero integer seeds can be supplied, failing which a default set is employed. If supplied, these three seeds, in order, should lie in the ranges [1,32362], [1,31726] and [1,31656] respectively.
SUBROUTINE RANECU (RVEC,LEN) C Portable random number generator for 16 bit computer. C Generates a sequence of LEN pseudo-random numbers, returned in C RVEC. DIMENSION RVEC(*) SAVE ISEED1,ISEED2, ISEED3 DATA ISEED1,ISEED2,ISEED3/1234, 5678, 9876/ C Default values, used if none suppliedvia an ENTRY C call at RECUIN DO 100 I = 1,LEN K=ISEED1/206 ISEED1 = 157 * (ISEED1 - K * 206) - K * 21 IF(ISEED1.LT.0) ISEED1=ISEED1+32363 K=ISEED2/217 ISEED2 = 146 * (ISEED2 - K*217) - K* 45 IF(ISEED2.LT.O) ISEED2=ISEED2+31727 K=ISEED3/222 ISEED3 = 142 * (ISEED3 - K *222) - K * 133 IF(ISEED3.LT.0) ISEED3=ISEED3+31657 IZ=ISEED1-ISEED2 IF(IZ.GT.706)IZ = Z - 32362 IZ = 1Z+ISEED3 IF(IZ.LT.1)IZ = 1Z + 32362 RVEC(I)=REAL(IZ) * 3.0899E - 5 100 CONTINUE RETURN ENTRY RECUIN(IS1, IS2, IS3) ISEED1=IS1 ISEED2=IS2 ISEED3=IS3 RETURN ENTRY RECUUT(IS1,IS2,IS3) IS1=ISEED1 IS2=ISEED2 IS3=ISEED3 RETURN END
The Shuffling Routine in Visual Basic
The Square Histogram Method
We are given a histogram, with vertical bars having heights proportional to the probability with which
we want to produce a value indicated by the label at the base.
A simple such histogram, layed flat, might be:
The idea is to cut the bars into pieces then reassemble them into a square histogram, all heights equal, with each final bar having a lower part, as well as an upper part indicating where it came from. A single uniform random variable U can then be used to choose one of the final bars and to indicate whether to use the lower or upper part. There are many ways to do this cutting and reassembling; the simplest seems to be the Robin Hood Algorithm: Take from richest to bring the poorest up to average.
STEP 1: The original (horizontal) histogram, average "height" 20:
Take 17 from strip 'a' to bring strip 'e' up to average. Record donor and use old 'poor' level to mark lower part of donee:
Then bring 'd' up to average with donor 'b'. Record donor and use old 'poor' level to mark lower part of donee:
Then bring 'a' up to average with donor 'c'. Record donor and use old 'poor' level to mark lower part of donee:
Finally, bring 'b' up to average with donor 'c'. Record donor and use old 'poor' level to mark lower part of donee:
We now have a "squared histogram", i.e., a rectangle with 4 strips of equal area, each strip with two regions. A single uniform variate U can be used to generate a,b,c,d,e with the required probabilities, .32, .27, .26, .12 .06.
Setup: Make tables,
Generation Process:
Let j be the integer part of 1+5*U, with U uniform in (0,1). If U < T[j] return V[j], else return V[K[j]]. In many applications no V table is necessary: V[i]=i and the generating procedure becomes If U < T[j] return j, else return K[j].
For more, visit the Web site: Modeling & Simulation Resources.
References & Further Readings:
Aiello W., S. Rajagopalan, and R. Venkatesan, Design of practical and provably good random number generators, Journal of Algorithms, 29, 358-389, 1998.
Dagpunar J., Principles of Random Variate Generation, Clarendon, 1988.
Fishman G., Monte Carlo, Springer, 1996.
James, Fortran version of L'Ecuyer generator, Comput. Phys. Comm., 60, 329-344, 1990.
Knuth D., The Art of Computer Programming, Vol. 2, Addison-Wesley, 1998.
L'Ecuyer P., Efficient and portable combined random number generators, Comm. ACM, 31, 742-749, 774, 1988.
L'Ecuyer P., Uniform random number generation, Ann. Op. Res., 53, 77-120, 1994.
L'Ecuyer P., Random number generation. In Handbook on Simulation, J. Banks (ed.), Wiley, 1998.
Maurer U., A universal statistical test for random bit generators, J. Cryptology, 5, 89-105, 1992.
Sobol' I., and Y. Levitan, A pseudo-random number generator for personal computers, Computers & Mathematics with Applications, 37(4), 33-40, 1999.
Tsang W-W., A decision tree algorithm for squaring the histogram in random number generation, Ars Combinatoria, 23A, 291-301, 1987
Test for Randomness:
A. Test for independence:
Plot the xi realization vs xi+1. If there is independence, the graph will not show any distinctive patterns at all, but will be perfectly scattered.
B. Runs tests.(run-ups, run-downs):
This is a direct test of the independence assumption. There are two test statistics to consider: one based on a normal approximation and another using numerical approximations.
Test based on Normal approximation:
Suppose you have N random realizations. Let a be the total number of runs in a sequence. If the number of positive and negative runs are greater than say 20, the distribution of a is reasonably approximated by a Normal distribution with mean (2N - 1) /3 and (16N - 29) / 90. Reject the hypothesis of independence or existence of runs if | Zo| > Z(1-alpha/2) where Zo is the Z score.
C. Correlation tests:
Do the random numbers exhibit discernible correlation? Compute the sample Autcorrelation Function.
Frequency or Uniform Distribution Test:
Use Kolmogorov-Smirimov test to determine if the realizations follow a U(0,1)
References & Further Readings:
Headrick T., Fast fifth-order polynomial transforms for generating univariate and multivariate nonnormal distributions, Computational Statistics and Data Analysis, 40 (4), 685-711, 2002.
Karian Z., and E. Dudewicz, Modern Statistical Systems and GPSS Simulation, CRC Press, 1998.
Kleijnen J., and W. van Groenendaal, Simulation: A Statistical Perspective, Wiley, Chichester, 1992
Korn G., Real statistical experiments can use simulation-package software, Simulation Modelling Practice and Theory, 13(1), 39-54, 2005.
Lewis P., and E. Orav, Simulation Methodology for Statisticians, Operations Analysts, and Engineers, Wadsworth Inc., 1989
Madu Ch., and Ch-H. Kuei, Experimental Statistical Designs and Analysis in Simulation Modeling, Greenwood Publishing Group, 1993.
Pang K., Z. Yang, S. Hou, and P. Leung, Non-uniform random variate generation by the vertical strip method, European Journal of Operational Research, 142(3), 595-609, 2002.
Robert C., and G. Casella, Monte Carlo Statistical Methods, Springer, 1999.
Simulation in general is to pretend that one deals with a real thing while really working with an imitation. In operations research the imitation is a computer model of the simulated reality. A flight simulator on a PC is also a computer model of some aspects of the flight: it shows on the screen the controls and what the "pilot" (the youngster who operates it) is supposed to see from the "cockpit" (his armchair).
Why to use models? To fly a simulator is safer and cheaper than the real airplane. For precisely this reason, models are used in industry commerce and military: it is very costly, dangerous and often impossible to make experiments with real systems. Provided that models are adequate descriptions of reality (they are valid), experimenting with them can save money, suffering and even time.
When to use simulations? Systems that change with time, such as a gas station where cars come and go (called dynamic systems) and involve randomness. Nobody can guess at exactly which time the next car should arrive at the station, are good candidates for simulation. Modeling complex dynamic systems theoretically need too many simplifications and the emerging models may not be therefore valid. Simulation does not require that many simplifying assumptions, making it the only tool even in absence of randomness.
How to simulate? Suppose we are interested in a gas station. We may describe the behavior of this system graphically by plotting the number of cars in the station; the state of the system. Every time a car arrives the graph increases by one unit while a departing car causes the graph to drop one unit. This graph (called sample path), could be obtained from observation of a real station, but could also be artificially constructed. Such artificial construction and the analysis of the resulting sample path (or more sample paths in more complex cases) consists of the simulation.
Types of simulations: Discrete event. The above sample path consisted of only horizontal and vertical lines, as car arrivals and departures occurred at distinct points of time, what we refer to as events. Between two consecutive events, nothing happens - the graph is horizontal. When the number of events are finite, we call the simulation "discrete event."
In some systems the state changes all the time, not just at the time of some discrete events. For example, the water level in a reservoir with given in and outflows may change all the time. In such cases "continuous simulation" is more appropriate, although discrete event simulation can serve as an approximation.
Further consideration of discrete event simulations.
How is simulation performed? Simulations may be performed manually. Most often, however, the system model is written either as a computer program (for an example click here) or as some kind of input into simulator software.
System terminology:
State: A variable characterizing an attribute in the system such as level of stock in inventory or number of jobs waiting for processing.
Event: An occurrence at a point in time which may change the state of the system, such as arrival of a customer or start of work on a job.
Entity: An object that passes through the system, such as cars in an intersection or orders in a factory. Often an event (e.g., arrival) is associated with an entity (e.g., customer).
Queue: A queue is not only a physical queue of people, it can also be a task list, a buffer of finished goods waiting for transportation or any place where entities are waiting for something to happen for any reason.
Creating: Creating is causing an arrival of a new entity to the system at some point in time.
Scheduling: Scheduling is the act of assigning a new future event to an existing entity.
Random variable: A random variable is a quantity that is uncertain, such as interarrival time between two incoming flights or number of defective parts in a shipment.
Random variate: A random variate is an artificially generated random variable.
Distribution: A distribution is the mathematical law which governs the probabilistic features of a random variable.
A Simple Example: Building a simulation gas station with a single pump served by a single service man. Assume that arrival of cars as well their service times are random. At first identify the:
states: number of cars waiting for service and number of cars served at any moment
events: arrival of cars, start of service, end of service
entities: these are the cars
queue: the queue of cars in front of the pump, waiting for service
random realizations: interarrival times, service times
distributions: we shall assume exponential distributions for both the interarrival time and service time.
Next, specify what to do at each event. The above example would look like this: At event of entity arrival: Create next arrival. If the server is free, send entity for start of service. Otherwise it joins the queue. At event of service start: Server becomes occupied. Schedule end of service for this entity. At event of service end: Server becomes free. If any entities waiting in queue: remove first entity from the queue; send it for start of service.
Some initiation is still required, for example, the creation of the first arrival. Lastly, the above is translated into code. This is easy with an appropriate library which has subroutines for creation, scheduling, proper timing of events, queue manipulations, random variate generation and statistics collection.
How to simulate? Besides the above, the program records the number of cars in the system before and after every change, together with the length of each event.
A typical stochastic system has a large number of control parameters that can have a significant impact on the performance of the system. To establish a basic knowledge of the behavior of a system under variation of input parameter values and to estimate the relative importance of the input parameters, sensitivity analysis applies small changes to the nominal values of input parameters. For systems simulation, variations of the input parameter values cannot be made infinitely small. The sensitivity of the performance measure with respect to an input parameter is therefore defined as (partial) derivative.
Problem Formulation: Identify controllable and uncontrollable inputs. Identify constraints on the decision variables. Define measure of system performance and an objective function. Develop a preliminary model structure to interrelate the inputs and the measure of performance.
Data Collection and Analysis: Regardless of the method used to collect the data, the decision of how much to collect is a trade-off between cost and accuracy.
Simulation Model Development: Acquiring sufficient understanding of the system to develop an appropriate conceptual, logical and then simulation model is one of the most difficult tasks in simulation analysis.
Model Validation, Verification and Calibration: In general, verification focuses on the internal consistency of a model, while validation is concerned with the correspondence between the model and the reality. The term validation is applied to those processes which seek to determine whether or not a simulation is correct with respect to the "real" system. More prosaically, validation is concerned with the question "Are we building the right system?". Verification, on the other hand, seeks to answer the question "Are we building the system right?" Verification checks that the implementation of the simulation model (program) corresponds to the model. Validation checks that the model corresponds to reality. Calibration checks that the data generated by the simulation matches real (observed) data.
Validation: The process of comparing the model's output with the behavior of the phenomenon. In other words: comparing model execution to reality (physical or otherwise)
Verification: The process of comparing the computer code with the model to ensure that the code is a correct implementation of the model.
Calibration: The process of parameter estimation for a model. Calibration is a tweaking/tuning of existing parameters and usually does not involve the introduction of new ones, changing the model structure. In the context of optimization, calibration is an optimization procedure involved in system identification or during experimental design.
Input and Output Analysis: Discrete-event simulation models typically have stochastic components that mimic the probabilistic nature of the system under consideration. Successful input modeling requires a close match between the input model and the true underlying probabilistic mechanism associated with the system. The input data analysis is to model an element (e.g., arrival process, service times) in a discrete-event simulation given a data set collected on the element of interest. This stage performs intensive error checking on the input data, including external, policy, random and deterministic variables. System simulation experiment is to learn about its behavior. Careful planning, or designing, of simulation experiments is generally a great help, saving time and effort by providing efficient ways to estimate the effects of changes in the model's inputs on its outputs. Statistical experimental-design methods are mostly used in the context of simulation experiments.
Performance Evaluation and What-If Analysis: The `what-if' analysis is at the very heart of simulation models.
Sensitivity Estimation: Users must be provided with affordable techniques for sensitivity analysis if they are to understand which relationships are meaningful in complicated models.
Optimization: Traditional optimization techniques require gradient estimation. As with sensitivity analysis, the current approach for optimization requires intensive simulation to construct an approximate surface response function. Incorporating gradient estimation techniques into convergent algorithms such as Robbins-Monroe type algorithms for optimization purposes, will be considered.
Gradient Estimation Applications: There are a number of applications which measure sensitivity information, (i.e., the gradient, Hessian, etc.), Local information, Structural properties, Response surface generation, Goal-seeking problem, Optimization, What-if Problem, and Meta-modelling
Report Generating: Report generation is a critical link in the communication process between the model and the end user.
Examples are the delay {D(i), i = 1, 2, ...} of the ith customer and number of customers {Q(t), T ³ 0} in the queue at time t in an M/M/1 queue. In the first example, we have a discrete- time, continuous state, while in the second example the state is discrete and time in continuous.
The following table is a classification of various stochastic processes. The man made systems have mostly discrete state. Monte Carlo simulation deals with discrete time while in discrete even system simulation the time dimension is continuous, which is at the heart of this site.
Change in the States of the System Continuous Discrete Time Continuous Level of water
Stationary Process (strictly stationary): A stationary stochastic process is a stochastic process {X(t), t Î T} with the property that the joint distribution all vectors of h dimension remain the same for any fixed h.
First Order Stationary: A stochastic process is a first order stationary if expected of X(t) remains the same for all t.
For example in economic time series, a process is first order stationary when we remove any kinds of trend by some mechanisms such as differencing.
Second Order Stationary: A stochastic process is a second order stationary if it is first order stationary and covariance between X(t) and X(s) is function of t-s only.
Again, in economic time series, a process is second order stationary when we stabilize also its variance by some kind of transformations such as taking square root.
Clearly, a stationary process is a second order stationary, however the reverse may not hold.
In simulation output statistical analysis we are satisfied if the output is covariance stationary.
Covariance Stationary: A covariance stationary process is a stochastic process {X(t), t Î T} having finite second moments, i.e. expected of [X(t)]2 be finite.
Clearly, any stationary process with finite second moment is covariance stationary. A stationary process may have no finite moment whatsoever.
Since a Gaussian process needs a mean and covariance matrix only, it is stationary (strictly) if it is covariance stationary.
Two Contrasting Stationary Process:
Consider the following two extreme stochastic processes:
- A sequence Y0, Y1,....., of independent identically distributed, random-value sequence is a stationary process, if its common distribution has a finite variance then the process is covariance stationary.
- Let Z be a single random variable with known distribution function, and set Z0 = Z1 = ....Z. Note that in a realization of this process, the first element, Z0, may be random but after that there is no randomness. The process {Zi, i = 0, 1, 2, ..} is stationary if Z has a finite variance.
Output data in simulation fall between these two type of process. Simulation outputs are identical, and mildly correlated (how mild? It depends on e.g. in a queueing system how large is the traffic intensity r). An example could be the delay process of the customers in a queueing system.
Gather steady state simulation output requires statistical assurance that the simulation model reached the steady state. The main difficulty is to obtain independent simulation runs with exclusion of the transient period . The two technique commonly used for steady state simulation are the Method of Batch means, and the Independent Replication.
None of these two methods is superior to the other in all cases. Their performance depend on the magnitude of the traffic intensity. The other available technique is the Regenerative Method, which is mostly used for its theoretical nice properties, however it is rarely applied in actual simulation for obtaining the steady state output numerical results.
Suppose you have a regenerative simulation consisting of m cycles of size n1, n2,…nm, respectively. The cycle sums is:
yi = S xij / ni, the sum is over j=1, 2, ..,ni The overall estimate is:Estimate = Syi / S ni, the sums are over i=1, 2, ..,m The 100(1-a/2)% confidence interval using the Z-table (or T-table, for m less than, say 30), is: Estimate ± Z. S/ (n. m½) where, n = S ni /m, the sum is over i=1, 2, ..,m and the variance is:S2 = S (yi - ni . Estimate)2/(m-1), the sum is over i=1, 2, ..,m
Method of Batch Means: This method involves only one very long simulation run which is suitably subdivided into an initial transient period and n batches. Each of the batch is then treated as an independent run of the simulation experiment while no observation are made during the transient period which is treated as warm-up interval. Choosing a large batch interval size would effectively lead to independent batches and hence, independent runs of the simulation, however since number of batches are few on cannot invoke the central limit theorem to construct the needed confidence interval. On the other hand, choosing a small batch interval size would effectively lead to significant correlation between successive batches therefore cannot apply the results in constructing an accurate confidence interval.
Suppose you have n equal batches of m observations each. The means of each batch is:
meani = S xij / m, the sum is over j=1, 2, ..,m The overall estimate is: Estimate = Smeani / n, the sum is over i=1, 2, ..,nThe 100(1-a/2)% confidence interval using the Z-table (or T-table, for n less than, say 30), is: Estimate ± Z. S where the variance is: S2 = S (meani - Estimate)2/(n-1), the sum is over i=1, 2, ..,n
Method of Independent Replications: This method is the most popularly used for systems with short transient period. This method requires independent runs of the simulation experiment different initial random seeds for the simulators' random number generator. For each independent replications of the simulation run it transient period is removed. For the observed intervals after the transient period data is collected and processed for the point estimates of the performance measure and for its subsequent confidence interval.
Suppose you have n replications with of m observations each. The means of each replication is:
meani = S xij / m, the sum is over j=1, 2, ..,m The overall estimate is:Estimate = Smeani / n, the sum is over i=1, 2, ..,n The 100(1-a/2)% confidence interval using the Z-table (or T-table, for n less than, say 30), is: Estimate ± Z. Swhere the variance is: S2 = S (meani - Estimate)2/(n-1), the sum is over i=1, 2, ..,n
Further Reading:
Sherman M., and D. Goldsman, Large-sample normality of the batch-means variance estimator, Operations Research Letters, 30, 319-326, 2002.
Whitt W., The efficiency of one long run versus independent replications in steady-state simulation, Management Science, 37(6), 645-666, 1991.
Batch Means is a method of estimating the steady-state characteristic from a single-run simulation. The single run is partitioned into equal size batches large enough for estimates obtained from different batches to be approximately independent. In the method of Batch Means, it is important to ensure that the bias due to initial conditions is removed to achieve at least a covariance stationary waiting time process. An obvious remedy is to run the simulation for a period large enough to remove the effect of the initial bias. During this warm-up period, no attempt is made to record the output of the simulation. The results are thrown away. At the end of this warm-up period, the waiting time of customers are collected for analysis. The practical question is "How long should the warm-up period be?". Abate and Whitt provided a relatively simple and nice expression for the time required (tp) for an M/M/1/ queue system (with traffic intensity r) starting at the origin (empty) to reach and remain within 100p% of the steady- state limit as follows:
tp(r) = 2C(r) Ln {1/[(1-p)(1+2C(r))]}/(1-r)2
where
C(r)=[2+ r + ( r2 + 4r )½] / 4.
Some notions of tp(r) as a function of r and p, are given in following table:
Traffic
Although this result is developed for M/M/1 queues, it has already been established that it can serve as an approximation for more general; i.e., GI/G/1 queues.
Further Reading:
Abate J., and W. Whitt, Transient behavior of regular Brownian motion, Advance Applied Probability, 19, 560-631, 1987.
Chen E., and W. Kelton, Determining simulation run length with the runs test, Simulation Modelling Practice and Theory, 11, 237-250, 2003.
After deciding what method is more suitable to apply, the main question is determination of number of runs. That is, at the planning stage of a simulation investigation of the question of number of simulation runs (n) is critical.
The confidence level of simulation output drawn from a set of simulation runs depends on the size of data set. The larger the number of runs, the higher is the associated confidence. However, more simulation runs also require more effort and resources for large systems. Thus, the main goal must be in finding the smallest number of simulation runs that will provide the desirable confidence.
Pilot Studies: When the needed statistics for number of simulation runs calculation is not available from existing database, a pilot simulation is needed.
For large pilot simulation runs (n), say over 30, the simplest number of runs determinate is:
[(Za/2)2 S2] / d2
where d is the desirable margin of error (i.e., the absolute error), which is the half-length of the confidence interval with 100(1- a)% confidence interval. S2 is the variance obtained from the pilot run.
One may use the following sample size determinate for a desirable relative error D in %, which requires an estimate of the coefficient of variation (C.V. in %) from a pilot run with n over 30:
[(Za/2)2 (C.V.)2] / D2
These sample size determinates could also be used for simulation output estimation of unimodal output populations, with discrete or continuous random variables provided the pilot run size (n) is larger than (say) 30.
The aim of applying any one of the above number of runs determinates is at improving your pilot estimates at feasible costs.
You may like using the following Applet for determination of number of runs.
Further Reading:
Díaz-Emparanza I, Is a small Monte Carlo analysis a good analysis? Checking the size power and consistency of a simulation-based test, Statistical Papers, 43(4), 567-577, 2002.
Whitt W., The efficiency of one long run versus independent replications in steady-state simulation, Management Science, 37(6), 645-666, 1991.
At the planning stage of a simulation modeling the question of number of simulation runs (n) is critical. The following Java applets compute the needed Runs Size based on current avialable information ontained from a pilot simulation run, to achieve an acceptable accuracy and/or risk.
Enter the needed information, and then click the Calculate button.
The aim of applying any one of the following number of simulation runs determinates is at improving your pilot estimates at a feasible cost.
Notes: The normality condition might be relaxed for number of simulation runs over, say 30. Moreover, determination of number of simulation runs for mean could also be used for other unimodal simulation output distributions including those with discrete random variables, such as proportion, provided the pilot run is sufficiently large (say, over 30).
Pilot Runs' Size (n): Current Estimate: Current Variance Estimate: Acceptable Significant Level (a): Acceptable Absolute Error: The Required Runs' Size Is:
ACSL, APROS, ARTIFEX, Arena, AutoMod, C++SIM, CSIM, Call$im, FluidFlow, GPSS, Gepasi, JavSim, MJX, MedModel, Mesquite, Multiverse, NETWORK, OPNET Modeler, POSES++, Simulat8, Powersim, QUEST, REAL, SHIFT, SIMPLE++, SIMSCRIPT, SLAM, SMPL, SimBank, SimPlusPlus, TIERRA, Witness, SIMNON, VISSIM, and javasim.
There are several things that make an ideal simulation package. Some are properties of the package, such as support, reactivity to bug notification, interface, etc. Some are properties of the user, such as their needs, their level of expertise, etc. For these reasons asking which package is best is a sudden failure of judgment. The first question to ask is for what purpose you need the software? Is it for education, teaching, student-projects or research?
The main question is: What are the important aspects to look for in a package? The answer depends on specific applications. However some general criteria are: Input facilities, Processing that allows some programming, Optimization capability, Output facilities, Environment including training and support services, Input-output statistical data analysis capability, and certainly the Cost factor.
You must know which features are appropriate for your situation, although, this is not based on a "Yes" or "No" judgment.
For description of available simulation software, visit Simulation Software Survey.
Reference & Further Reading:
Nikoukaran J., Software selection for simulation in manufacturing: A review, Simulation Practice and Theory, 7(1), 1-14, 1999.
Without computer one cannot perform any realistic dynamic systems simulation.
SIMSCRIPT II.5 is a powerful, free-format, English-like simulation language designed to greatly simplify writing programs for simulation modelling. Programs written in SIMSCRIPT II.5 are easily read and maintained. They are accurate, efficient, and generate results which are acceptable to users. Unlike other simulation programming languages, SIMSCRIPT II.5 requires no coding in other languages. SIMSCRIPT II.5 has been fully supported for over 33 years. Contributing to the wide acceptance and success of SIMSCRIPT II.5 modelling are:
DESIGN:
A powerful worldview, consisting of Entities and Processes, provides a natural conceptual framework with which to relate real objects to the model.
PROGRAMMING:
SIMSCRIPT II.5 is a modern, free-form language with structured programming constructs and all the built-in facilities needed for model development. Model components can be programmed so they clearly reflect the organization and logic of the modeled system. The amount of program needed to model a system is typically 75% less than its FORTRAN or C counterpart.
DEBUGGER:
A well designed package of program debug facilities is provided. The required tools are available to detect errors in a complex computer program without resorting an error. Simulation status information is provided, and control is optionally transferred to a user program for additional analysis and output.
EVOLUTION:
This structure allows the model to evolve easily and naturally from simple to detailed formulation as data becomes available. Many modifications, such as the choice of set disciplines and statistics are simply specified in the Preamble.
DOCUMENTATION:
You get a powerful, English-like language supporting a modular implementation. Because each model component is readable and self-contained, the model documentation is the model listing; it is never obsolete or inaccurate.
For more information contact SIMSCRIPT
Guidelines for Running SIMSCRIPT on the VAX System
Network Access/Utilities
Connect UBE
Username:
Password:
Get $ sign
Step 1. Create and edit your source program. Type in
$EDT PROG.SIM
Step 2. Attach SIMSCRIPT. Type in
$SIMSCRIPT
Step 3. Compile the source program file. Type in
$SIMSCOMP PROG.SIM
Check for compilation errors. To locate the errors, type in
$EDT PROG.LIST
Step 4. Link the object file. Type in
$SIMLINK PROG
now you have a file containing your program in an
executable format.
Step 5. To execute the program type in
$RUN PROG
Step 6. To get your hard copy print, type
$PRINT/NAME=your own name PROG.OUT, PROG.LIS
Alternatively, use submit command to place the command procedure
in the batch job que.
Create a command file say PROG.COM containing:
$SIMSCRIPT
$SIMSCOMP PROG.SIM
$SIMLINK PROG
$RUN PROG
Then, submit
$Submit PROG.COM
An Example:
'' Solving an analytic equation arising from optimization
'' of a coherent reliability system with 3 homogeneous components
Preamble
Define V as a real 1-dimensional arrays
Define I, N as integer variables
End
Main
Open 3 for output, Name = "PROG.OUT"
Use 3 for ouTPUT
LET N=50
Reserve V(*) as N
LET V(1)=1.
For I = 1 to N-1
DO
LET U=V(I)
LET PPRAM=-9./((1+U)**2) + 9./((2.+U)**2) + 1./U**2
LET V(I+1)=V(I)+PPRAM/I
LOOP
PRINT 1 LINE WITH V(N) AND PPRAM THUS
OPTIMAL RATE IS ****.*****, DERIVATIVE IS ***.*****
END
The output
OPTIMAL RATE IS 0.76350, DERIVATIVE IS -0.00000
The modeling techniques used by system dynamics and discrete event simulations are often different at two levels: The modeler way of representing systems might be different, the underlying simulators' algorithms are also different. Each technique is well tuned to the purpose it is intended. However, one may use a discrete event approach to do system dynamics and vice versa.
Traditionally, the most important distinction is the purpose of the modeling. The discrete event approach is to find, e.g., how many resources the decision maker needs such as how many trucks, and how to arrange the resources to avoid bottlenecks, i.e., excessive of waiting lines, waiting times, or inventories. While the system dynamics approach is to prescribe for the decision making to, e.g., timely respond to any changes, and how to change the physical structure, e.g., physical shipping delay time, so that inventories, sales, production, etc.
System dynamics is the rigorous study of problems in system behavior using the principles of feedback, dynamics and simulation. In more words system dynamics is characterized by:
Some interesting and useful areas of system dynamics modeling approach are:
A modeler must consider both as complementary tools to each other. Systems dynamic to look at the high level problem and identify areas which need more detailed analysis. Then, use discrete event modeling tools to analyze (and predict) the specific areas of interest.
Faster hardware and improved software have made building complex simulations easier. Computer simulation methods can be effective for the development of theories as well as for prediction. For example, macro-economic models have been used to simulate future changes in the economy; and simulations have been used in psychology to study cognitive mechanisms.
The field of 'social simulation' seems to be following an interesting line of inquiry. As a general approach in the field, a 'world' is specified with much computational detail. Then the 'world' is simulated (using computers) to reveal some of the 'non-trivial' implications (or 'emergent properties') of the 'world'. When these 'non trivial' implications are made known (fed back) in world, apparently it constitutes some 'added values'.
Artificial Life is an interdisciplinary study enterprise aimed at understanding life-as-it-is and life-as-it-could-be, and at synthesizing life-like phenomena in chemical, electronic, software, and other artificial media. Artificial Life redefines the concepts of artificial and natural, blurring the borders between traditional disciplines and providing new media and new insights into the origin and principles of life.
Simulation allows the social scientist to experiment with ‘artificial societies' and explore the implications of theories in ways not otherwise possible.
Reference and Further Readings:
Gilbert N., and K. Troitzsch, Simulation for the Social Scientist, Open University Press, Buckingham, UK, 1999.
Sichman J., R. Conte, and N. Gilbert, (eds,), Multi-Agent Systems and Agent-Based Simulation, Berlin, Springer-Verlag, 1998.
The appearance of the network-friendly programming language, Java, and of distributed object technologies like the Common Object Request Broker Architecture (CORBA) and the Object Linking and Embedding / Component Object Model (OLE/COM) have had particularly acute effects on the state of simulation practice.
Currently, the researchers in the field of web-based simulation are interested in dealing with topics such as methodologies for web-based model development, collaborative model development over the Internet, Java-based modeling and simulation, distributed modeling and simulation using web technologies, and new applications.
Synchronization, scheduling, memory management, randomized and reactive/adaptive algorithms, partitioning and load balancing.
Synchronization in multi-user distributed simulation, virtual reality environments, HLA, and interoperability.
System modeling for parallel simulation, specification, re-use of models/code, and parallelizing existing simulations.
Language and implementation issues, models of parallel simulation, execution environments, and libraries.
Theoretical and empirical studies, prediction and analysis, cost models, benchmarks, and comparative studies.
Computer architectures, VLSI, telecommunication networks, manufacturing, dynamic systems, and biological/social systems.
Web based distributed simulation such as multimedia and real time applications, fault tolerance, implementation issues, use of Java, and CORBA.
References & Further Readings:
Bossel H., Modeling & Simulation, A. K. Peters Pub., 1994.
Delaney W., and E. Vaccari, Dynamic Models and Discrete Event Simulation, Dekker, 1989.
Fishman G., Discrete-Event Simulation: Modeling, Programming and Analysis, Springer-Verlag, Berlin, 2001.
Fishwick P., Simulation Model Design and Execution: Building Digital Worlds, Prentice-Hall, Englewood Cliffs, 1995.
Ghosh S., and T. Lee, Modeling & Asynchronous Distributed Simulation: Analyzing Complex Systems, IEEE Publications, 2000.
Gimblett R., Integrating Geographic Information Systems and Agent-Based Modeling: Techniques for Simulating Social and Ecological Processes, Oxford University Press, 2002.
Harrington J., and K. Tumay, Simulation Modeling Methods: An Interactive Guide to Results-Based Decision, McGraw-Hill, 1998.
Haas P., Stochastic Petri Net Models Modeling and Simulation, Springer Verlag, 2002.
Hill D., Object-Oriented Analysis and Simulation Modeling, Addison-Wesley, 1996.
Kouikoglou V., and Y. Phillis, Hybrid Simulation Models of Production Networks, Kluwer Pub., 2001.
Law A., and W. Kelton, Simulation Modeling and Analysis, McGraw-Hill, 2000.
Nelson B., Stochastic Modeling: Analysis & Simulation, McGraw-Hill, 1995.
Oakshott L., Business Modelling and Simulation, Pitman Publishing, London, 1997.
Pidd M., Computer Simulation in Management Science, Wiley, 1998.
Rubinstein R., and B. Melamed, Modern Simulation and Modeling, Wiley, 1998.
Severance F., System Modeling and Simulation: An Introduction, Wiley, 2001.
Van den Bosch, P. and A. Van der Klauw, Modeling, Identification & Simulation of Dynamical Systems, CRC Press, 1994.
Woods R., and K. Lawrence, Modeling and Simulation of Dynamic Systems, Prentice Hall, 1997.
DES's system performance is often measured as an expected value. Consider a system with continuous parameter v Î V Í Rn, where V is an open set. Let
J(v) = EY | v [Z (Y)]
be the steady state expected performance measure, where Y is a random vector with known probability density function (pdf), f(y; v) depends on v, and Z is the performance measure.
In discrete event systems, Monte Carlo simulation is usually needed to estimate J(v) for a given value v = v0. By the law of large numbers
J(v0) = 1/n S Z (y i)
converges to the true value, where yi, i = 1, 2,..., n are independent, identically distributed, random vector realizations of Y from f (y; v 0), and n is the number of independent replications.
We are interested in sensitivities estimation of J(v) with respect to v.
Local information: An estimate for dJ/dv is a good local measure of the effect of on performance. For example, simply knowing the sign of the derivative dJ/dv at some point v immediately gives us the direction in which v should be changed. The magnitude of dJ/d? also provides useful information in an initial design process: If dJ/dv is small, we conclude that J is not very sensitive to changes in , and hence focusing concentration on other parameters may improve performance.
Structural properties: Often sensitivity analysis provides not only a numerical value for the sample derivative, but also an expression which captures the nature of the dependence of a performance measure on the parameter v . The simplest case arises when dJ/dv can be seen to be always positive (or always negative) for any sample path; we may not be able to tell if the value of J(v) is monotonically increasing (or decreasing) in v . This information in itself is very useful in design and analysis. More generally, the form of dJ/dv can reveal interesting structural properties of the DES (e.g., monotonicity, convexity). Such properties must be exploited in order to determine optimal operating policies for some systems.
Response surface generation: Often our ultimate goal is to obtain the function J(v), i.e., a curve describing how the system responds to different values of v . Since J(v) is unknown, one alternative is to obtain estimates of J(v) for as many values of v as possible. This is clearly a prohibitively difficult task. Derivative information, however may include not only first-order but also higher derivatives which can be used to approximate J(v). If such derivative information can be easily and accurately obtained, the task of response surface generation may be accomplished as well.
Goal-seeking and What-if problems: Stochastic models typically depend upon various uncertain parameters that must be estimated from existing data sets. Statistical questions of how input parameter uncertainty propagates through the model into output parameter uncertainty is the so-called "what-if" analysis. A good answer to this question often requires sensitivity estimates. The ordinary simulation output results are the solution of a direct problem: Given the underlying pdf with a particular parameter value v , we may estimate the output function J(v). Now we pose the goal-seeking problem: given a target output value J0 of the system and a parameterized pdf family, find an input value for the parameter, which generates such an output. There are strong motivations for both problems. When v is any controllable or uncontrollable parameter [the decision maker is, for example, interested in estimating J(v) for a small change in v ], the so called "what-if" problem, which is a "direct problem" and can be solved by incorporating sensitivity information in the Taylor's expansion of J(v) in the neighborhood of v . However, when v is a controllable input, the decision maker may be interested in the goal-seeking problem: what change in the input parameter will achieve a desired change in output value J(v). Another application of goal-seeking arises when we want to adapt a model to satisfy a new equality constraint (condition) for some stochastic function. The solution to the goal-seeking problem is to estimate the derivative of the output function with respect to the input parameter for the nominal system; use this estimate in a Taylor's expansion of the output function in the neighborhood of the parameter; and finally, use Robbins-Monro (R-M) type of stochastic approximation algorithm to estimate the necessary controllable input parameter value within the desired accuracy.
Optimization: Discrete-event simulation is the primary analysis tool for designing complex systems. However, simulation must be linked with a mathematical optimization technique to be effectively used for systems design. The sensitivity dJ/dv can be used in conjunction with various optimization algorithms whose function is to gradually adjust v until a point is reached where J(v) is maximized (or minimized). If no other constraints on v are imposed, we expect dJ/dv = 0 at this point .
Step 1: How does a change in the value of a parameter vary the sample duration related to that parameter?
Step 2: How does the change in an individual sample duration reflect itself as a change in a subsequent particular sample realization?
Step 3: Finally, what is the relationship between the variation of the sample realization and its expected value?
Now we look at the general (non-regenerative) case. In this case any simulation will give a biased estimator of the gradient, as simulations are necessarily finite. If n (the length of the simulation) is large enough, this bias is negligible. However, as noted earlier, the variance of the SF sensitivity estimator increases with increase in n so, a crude SF estimator is not even approximately consistent. There are a number of ways to attack this problem. Most of the variations in an estimator comes from the score function. The variation is especially high, when all past inputs contribute to the performance and the scores from all are included. When one uses batch means, the variation is reduced by keeping the length of the batch small.
A second way is to reduce the variance of the score to such an extent that we can use simulations long enough to effectively eliminate the bias. This is the most promising approach. The variance may be reduced further by using the standard variance reduction techniques (VRT), such as importance sampling. Finally, we can simply use a large number of iid replications of the simulation.
Frequency domain simulation experiments identify the significant terms of the polynomial that approximates the relationship between the simulation output and the inputs. Clearly, the number of simulation runs required to identify the important terms by this approach is much smaller than those of the competing alternatives, and the difference becomes even more conspicuous as the number of parameters increases.
In terms of ease of implementation, PA estimators may require considerable analytical work on the part of algorithm developer, with some "customization" for each application, whereas SF has the advantage of remaining a general definable algorithm whenever it can be applied.
Perhaps the most important criterion for comparison lies in the question of accuracy of an estimator, typically measured through its variance. If an estimator is strongly consistent, its variance is gradually reduced over time and ultimately approaches to zero. The speed with which this happens may be extremely important. Since in practice, decisions normally have to be made in a limited time, an estimator whose variance decreases fast is highly desirable. In general, when PA does provide unbiased estimators, the variance of these estimators is small. PA fully exploits the structure of DES and their state dynamics by extracting the needed information from the observed sample path, whereas SF requires no knowledge of the system other than the inputs and the outputs. Therefore when using SF methods, variance reduction is necessary. The question is whether or not the variance can be reduced enough to make the SF estimator useful in all situations to which it can be applied. The answer is certainly yes. Using the standard variance reduction techniques can help, but the most dramatic variance reduction occurs using new methods of VR such as conditioning, which is shown numerically to have a mean squared error that is essentially the same as that of PA.
For more, visit the Web site: Modeling & Simulation Resources.
References & Further Readings:
Arsham H., Algorithms for Sensitivity Information in Discrete-Event Systems Simulation, Simulation Practice and Theory, 6(1), 1-22, 1998.
Fu M., and J-Q. Hu, Conditional Monte Carlo: Gradient Estimation and Optimization Applications, Kluwer Academic Publishers, 1997.
Rubinstein R., and A. Shapiro, Discrete Event Systems: Sensitivity Analysis and Stochastic Optimization by the Score Function Method, John Wiley & Sons, 1993.
Whitt W., Minimizing delays in the GI/G/1 queue, Operations Research, 32(1), 41-51, 1984.
Many man-made systems can be modeled as Discrete Event Systems (DES); examples are computer systems, communication networks, flexible manufacturing systems, production assembly lines, and traffic transportation systems. DES evolve with the occurrence of discrete events, such as the arrival of a job or the completion of a task, in contrast with continuously variable dynamic processes such as aerospace vehicles, which are primarily governed by differential equations. Owing to the complex dynamics resulting from stochastic interactions of such discrete events over time, the performance analysis and optimization of DES can be difficult tasks. At the same time, since such systems are becoming more widespread as a result of modern technological advances, it is important to have tools for analyzing and optimizing the parameters of these systems.
Analyzing complex DES often requires computer simulation. In these systems, the objective function may not be expressible as an explicit function of the input parameters; rather, it involves some performance measures of the system whose values can be found only by running the simulation model or by observing the actual system. On the other hand, due to the increasingly large size and inherent complexity of most man-made systems, purely analytical means are often insufficient for optimization. In these cases, one must resort to simulation, with its chief advantage being its generality, and its primary disadvantage being its cost in terms of time and money. Even though, in principle, some systems are analytically tractable, the analytical effort required to evaluate the solution may be so formidable that computer simulation becomes attractive. While the price for computing resources continue to dramatically decrease, one nevertheless can still obtain only a statistical estimate as opposed to an exact solution. For practical purposes, this is quite sufficient.
These man-made DES are costly, and therefore it is important to operate them as efficiently as possible. The high cost makes it necessary to find more efficient means of conducting simulation and optimizing its output. We consider optimizing an objective function with respect to a set of continuous and/or discrete controllable parameters subject to some constraints.
Usually when modelers choose a DES approach they often model the system as open loop or nearly open loop system, making the system behave as if there where no superior agent controlling the whole production/service/ process. Closing the loops should be an elemental task that simulation modeler should take care of, even if the scope does not involve doing it, there must be awareness of system behavior, particularly if there is known to be that the system if under human decision making processes/activities.
J(v) = EY | v [Z (Y)]
be the steady state expected performance measure, where Y is a random vector with known probability density function (pdf), f(y; v) depends on v, and Z is the performance measure.
In discrete event systems, Monte Carlo simulation is usually needed to estimate J(v) for a given value
v = v 0. By the law of large numbers
J(v0) = 1/n S Z (y i)
converges to the true value, where yi, i = 1, 2,..., n are independent, identically distributed, random vector realizations of Y from f (y; v 0), and n is the number of independent replications.
The aim is to optimize J(v) with respect to v .
We shall group the optimization techniques for simulation into seven broad categories; namely, Deterministic Search, Pattern Search, Probabilistic Search, Evolutionary Techniques, Stochastic Approximation, Gradient Surface, and some Mixtures of the these techniques;
Deterministic search techniques include heuristic search, complete enumeration, and random search techniques.
The heuristic search technique is probably most commonly used in optimizing response surfaces. It is also the least sophisticated scheme mathematically, and it can be thought of as an intuitive and experimental approach. The analyst determines the starting point and stopping rule based on previous experience with the system. After setting the input parameters (factors) to levels that appear reasonable, the analyst makes a simulation run with the factors set at those levels and computes the value of the response function. If it appears to be a maximum (minimum) to the analyst, the experiment is stopped. Otherwise the analyst changes parameter settings and makes another run. This process continues until the analyst believes that the output has been optimized. Suffice it to say that, if the analyst is not intimately familiar with the process being simulated, this procedure can turn into a blind search and can expend an inordinate amount of time and computer resources without producing results commensurate with input. The heuristic search can be ineffective and inefficient in the hand of a novice.
The random search technique resembles the complete enumeration technique except that one selects a set of inputs at random. The simulated results based on the set that yields the maximum (minimum) value of the response function is taken to be the optimal point. This procedure reduces the number of simulation runs required to yield an 'optimal' result; however, there is no guarantee that the point found is actually the optimal point. Of course, the more points selected, the more likely the analyst is to achieve the true optimum. Note that the requirement that each factor assumes only a finite number of values is not a requirement in this scheme. Replications can be made on the treatment combinations selected, to increase the confidence in the optimal point. Which strategy is better, replicating a few points or looking at a single observation on more points, depends on the problem.
Although it is desirable for search procedures to be efficient over a wide range of response surfaces, no current procedure can effectively overcome non-unimodality (surfaces having more than one local maximum or minimum). An obvious way to find the global optimal would be to evaluate all the local optima. One technique that is used when non-unimodality is known to exist, is called the "Las Vegas" technique. This search procedure estimates the distribution of the local optima by plotting the estimated J( v ) for each local search against its corresponding search number. Those local searches that produce a response greater than any previous response are then identified and a curve is fitted to the data. This curve is then used to project the "estimated incremental" response that will be achieved by one more search. The search continues until the value of the estimated improvement in the search is less than the cost of completing one additional search.
It should be noted that a well-designed experiment requires a sufficient number of replications so that the average response can be treated as a deterministic number for search comparisons. Otherwise, since replications are expensive, it becomes necessary to effectively utilize the number of simulation runs. Although each simulation is at a different setting of the controllable variables, one can use smoothing techniques such as exponential smoothing to reduce the required number of replications.
Two directions are defined to be conjugate whenever the cross-product terms are all zero. The conjugate direction technique tries to find a set of n dimensions that describes the surface such that each direction is conjugate to all others.
Using the above result, the technique attempts to find two search optima and replace the nth dimension of the quadratic surface by the direction specified by the two optimal points. Successively replacing the original dimension yields a new set of n dimensions in which, if the original surface is quadratic, all directions are conjugate to each other and appropriate for n single variable searches. While this search procedure appears to be very simple, we should point out that the selection of appropriate step sizes is most critical. The step size selection is more critical for this search technique because - during axis rotation - the step size does not remain invariant in all dimensions. As the rotation takes place, the best step size changes, and becomes difficult to estimate.
Although quadratic functions are sometimes used, one assumes that performance is linearly related to the change in the controllable variables for small changes. Assume that a good approximation is a linear form. The basis of the linear steepest ascent is that each controllable variable is changed in proportion to the magnitude of its slope. When each controllable variable is changed by a small amount, it is analogous to determining the gradient at a point. For a surface containing N controllable variables, this requires N points around the point of interest. When the problem is not an n-dimensional elliptical surface, the parallel-tangent points are extracted from bitangents and inflection points of occluding contours. Parallel tangent points are points on the occluding contour where the tangent is parallel to a given bitangent or the tangent at an inflection point.
A variation of the random search technique determines the maximum of the objective function by analyzing the distribution of J(v) in the bounded sub-region. In this variation, the random data are fitted to an asymptotic extreme-value distribution, and J* is estimated with a confidence statement. Unfortunately, these techniques cannot determine the location of J* , which can be as important as the J value itself. Some techniques calculate the mean value and the standard deviation of J(v) from the random data as they are collected. Assuming that J is distributed normally in the feasible region., the first trial, that yields a J-value two standard deviations within the mean value, is taken as a near-optimum solution.
SA as an optimization technique was first introduced to solve problems in discrete optimization, mainly combinatorial optimization. Subsequently, this technique has been successfully applied to solve optimization problems over the space of continuous decision variables. SA is a simulation optimization technique that allows random ascent moves in order to escape the local minima, but a price is paid in terms of a large increase in the computational time required. It can be proven that the technique will find an approximated optimum. The annealing schedule might require a long time to reach a true optimum.
GT differ from traditional optimization procedures in that GT work with a coding of the decision parameter set, not the parameters themselves; GT search a population of points, not a single point; GT use objective function information, not derivatives or other auxiliary knowledge; and finally, GT use probabilistic transition rules, not deterministic rules. GT are probabilistic search optimizing techniques that do not require mathematical knowledge of the response surface of the system, which they are optimizing. They borrow the paradigms of genetic evolution, specifically selection, crossover, and mutation.
Selection: The current points in the space are ranked in terms of their fitness by their respective response values. A probability is assigned to each point that is proportional to its fitness, and parents (a mating pair) are randomly selected.
Crossover: The new point, or offspring, is chosen, based on some combination of the genetics of the two parents.
Mutation: The location of offspring is also susceptible to mutation, a process, which occurs with probability p, by which a offspring is replaced randomly by a new offspring location.
A generalized GT generates p new offspring at once and kills off all of the parents. This modification is important in the simulation environment. GT are well suited for qualitative or policy decision optimization such as selecting the best queuing disciplines or network topologies. They can be used to help determine the design of the system and its operation. For applications of GT to inventory systems, job-shop, and computer time-sharing problems. GT do not have certain shortcomings of other optimization techniques, and they will usually result in better calculated optima than those found with the traditionally techniques. They can search a response surface with many local optima and find (with a high probability) the approximate global optimum. One may use GT to find an area of potential interest, and then resort to other techniques to find the optimum. Recently, several classical GT principles have been challenged.
Differential Evolution: Differential Evolution (DE) is a genetic type of algorithm for solving continuous stochastic function optimization. The basic idea is to use vector differences for perturbing the vector population. DE adds the weighted difference between two population vectors to a third vector. This way, no separate probability distribution has to be used, which makes the scheme completely self-organizing.
References & Further Readings:
Choi D.-H., Cooperative mutation based evolutionary programming for continuous function optimization, Operations Research Letters, 30, 195-201, 2002.
Reeves C., and J. Rowe, Genetic Algorithms: Principles and Perspectives, Kluwer, 2002.
Saviotti P., (Ed.), Applied Evolutionary Economics: New Empirical Methods and Simulation Techniques, Edward Elgar Pub., 2002.
Wilson W., Simulating Ecological and Evolutionary Systems in C, Cambridge University Press, 2000.
Interest was renewed in the R-M technique as a means of optimization, with the development of the perturbation analysis, score function (known also as likelihood ratio method), and frequency domain estimates of derivatives. Optimization for simulated systems based on the R-M technique is known as a "single-run" technique. These procedures optimize a simulation model in a single run simulation with a run length comparable to that required for a single iteration step in the other methods. This is achieved essentially be observing the sample values of the objective function and, based on these observations, updating the values of the controllable parameters while the simulation is running, that is, without restarting the simulation. This observing-updating sequence is done repeatedly, leading to an estimate of the optimum at the end of a single-run simulation. Besides having the potential of large computational savings, this technique can be a powerful tool in real-time optimization and control, where observations are taken as the system is evolving in time.
The gradient surface method (GSM) combines the virtue of RSM with that of the single- run, gradient estimation techniques such as Perturbation Analysis, and Score Function techniques. A single simulation experiment with little extra work yields N + 1 pieces of information; i.e., one response point and N components of the gradient. This is in contrast to crude simulation, where only one piece of information, the response value, is obtained per experiment. Thus by taking advantage of the computational efficiency of single-run gradient estimators. In general, N-fold fewer experiments will be needed to fit a global surface compared to the RSM. At each step, instead of using Robbins-Monro techniques to locate the next point locally, we determine a candidate for the next point globally, based on the current global fit to the performance surface.
The GSM approach has the following advantages; The technique can quickly get to the vicinity of the optimal solution because its orientation is global [23, 39]. Thus, it produces satisfying solutions quickly; Like RSM, it uses all accumulated information; And, in addition, it uses gradient surface fitting, rather than direct performance response-surface fitting via single-run gradient estimators. This significantly reduces the computational efforts compared with RSM. Similar to RSM, GSM is less sensitive to estimation error and local optimality; And, finally, it is an on-line technique, the technique may be implemented while the system is running.
A typical optimization scheme involves two phases: a Search Phase and an Iteration Phase. Most results in analytic computational complexity assume that good initial approximations are available, and deal with the iteration phase only. If enough time is spent in the initial search phase, we can reduce the time needed in the iteration phase. The literature contains papers giving conditions for the convergence of a process; a process has to be more than convergent in order to be computationally interesting. It is essential that we be able to limit the cost of computation. In this sense, GSM can be thought of as helping the search phase and as an aid to limit the cost of computation. One can adopt standard or simple devices for issues such as stopping rules.
For on-line optimization, one may use a new design in GSM called 'single direction' design. Since for on-line optimization it may not be advisable or feasible to disturb the system, random design usually is not suitable.
Simulation has numerous advantages over other approaches to performance and dependability evaluation; most notably, its modelling power and flexibility. For some models, however, a potential problem is the excessive simulation effort (time) required to achieve the desired accuracy. In particular, simulation of models involving rare events, such as those used for the evaluation of communications and highly-dependable systems, is often not feasible using standard techniques. In recent years, there have been significant theoretical and practical advances towards the development of efficient simulation techniques for the evaluation of these systems.
Methodologies include: Techniques based on importance sampling, The "restart" method, and Hybrid analytic/simulation techniques among newly devised approaches.
General comparisons among different techniques in terms of bias, variance, and computational complexity are not possible. However, a few studies rely on real computer simulations to compare different techniques in terms of accuracy and number of iterations. Total computational effort for reduction in both the bias andvariance of the estimate depends on the computational budget allocated for a simulation optimization. No single technique works effectively and/or efficiently in all cases.
The simplest technique is the random selection of some points in the search region for estimating the performance measure. In this technique, one usually fixes the number of simulation runs and takes the smallest (or largest) estimated performance measure as the optimum. This technique is useful in combination with other techniques to create a multi-start technique for global optimization. The most effective technique to overcome local optimality for discrete optimization is the Tabu Search technique. In general, the probabilistic search techniques, as a class, offer several advantages over other optimization techniques based on gradients. In the random search technique, the objective function can be non-smooth or even have discontinuities. The search program is simple to implement on a computer, and it often shows good convergence characteristics in noisy environments. More importantly, it can offer the global solution in a multi-modal problem, if the technique is employed in the global sense. Convergence proofs under various conditions are given in.
The Hooke-Jeeves search technique works well for unconstrained problems with less than 20 variables; pattern search techniques are more effective for constrained problems. Genetic techniques are most robust and can produce near-best solutions for larger problems. The pattern search technique is most suitable for small size problems with no constraint, and it requires fewer iterations than the genetic techniques. The most promising techniques are the stochastic approximation, simultaneous perturbation, and the gradient surface methods. Stochastic approximation techniques using perturbation analysis, score function, or simultaneous perturbation gradient estimators, optimize a simulation model in a single simulation run. They do so by observing the sample values of the objective function, and based on these observations, the stochastic approximation techniques update the values of the controllable parameters while the simulation is running and without restarting the simulation. This observing-updating sequence, done repeatedly, leads to an estimate of the optimum at the end of a single-run simulation. Besides having the potential of large savings in computational effort in the simulation environment, this technique can be a powerful tool in real-time optimization and control, where observations are taken as the system is evolving over time.
Response surface methods have a slow convergence rate, which makes them expensive. The gradient surface method combines the advantages of the response surface methods (RSM) and efficiency of the gradient estimation techniques, such as infinitesimal perturbation analysis, score function, simultaneous perturbation analysis, and frequency domain technique. In the gradient surface method (GSM) the gradient is estimated, and the performance gradient surface is estimated from observations at various points, similar to the RSM. Zero points of the successively approximating gradient surface are then taken as the estimates of the optimal solution. GSM is characterized by several attractive features: it is a single run technique and more efficient than RSM; at each iteration step, it uses the information from all of the data points rather than just the local gradient; it tries to capture the global features of the gradient surface and thereby quickly arrive in the vicinity of the optimal solution, but close to the optimum, they take many iterations to converge to stationary points. Search techniques are therefore more suitable as a second phase. The main interest is to figure out how to allocate the total available computational budget across the successive iterations.
For when the decision variable is qualitative, such as finding the best system configuration, a random or permutation test is proposed. This technique starts with the selection of an appropriate test statistic, such as the absolute difference between the mean responses under two scenarios. The test value is computed for the original data set. The data are shuffled (using a different seed); the test statistic is computed for the shuffled data; and the value is compared to the value of the test statistic for the original, un-shuffled data. If the statistics for the shuffled data are greater than or equal to the actual statistic for the original data, then a counter c, is incremented by 1. The process is repeated for any desired m number of times. The final step is to compute (c+1)/(m+1), which is the significant level of the test. The null hypothesis is rejected if this significance level is less than or equal to the specified rejection level for the test. There are several important aspects to this nonparametric test. First, it enables the user to select the statistic. Second, assumptions such as normality or equality of variances made for the t-test, ranking-and-selection, and multiple-comparison procedures, are no longer needed. A generalization is the well-known bootstrap technique.
What Must Be Done
For more, visit the Web site: Modeling & Simulation Resources.
References & Further Readings:
Arsham H., Techniques for Monte Carlo Optimizing, Monte Carlo Methods and Applications, 4(3), 181-230, 1998.
Arsham H., Stochastic Optimization of Discrete Event Systems Simulation, Microelectronics and Reliability, 36(10), 1357-1368, 1996.
Fu M., and J-Q. Hu, Conditional Monte Carlo: Gradient Estimation and Optimization Applications, Kluwer Academic Publishers, 1997.
Rollans S. and D. McLeish, Estimating the optimum of a stochastic system using simulation, Journal of Statistical Computation and Simulation, 72, 357 - 377, 2002.
Rubinstein R., and A. Shapiro, Discrete Event Systems: Sensitivity Analysis and Stochastic Optimization by the Score Function Method, John Wiley & Sons, 1993.
In many simulation applications such as systems analysis and design applications, the decision maker may not be interested in optimization but wishes to achieve a certain value for J(v), say J0. This is the goal-seeking problem: given a target output value J0 of the performance and a parameterized pdf family, one must find an input value for the parameter, which generates such an output.
J(v) = J(v0) + dv.J' (v0) + (dv)2 J' ' (v0) / 2 + ...,
where dv = v-v0 and the primes denote derivatives. This metamodel approximates J(v) for small dv. To estimate J(v) in the neighborhood of v0 by a linear function, we need to estimate the nominal J(v) and its first derivative. Traditionally, this derivative is estimated by crude Monte Carlo; i.e., finite difference which requires rerunning the simulation model. Methods which yield enhanced efficiency and accuracy in estimating, at little additional computational (Not simulation) cost, are presented in this site. The Score Function method of estimating the first derivative is:
J'(v) = EY|v [ Z(y(v)) .S]
where S=f'(y; v) / f(y; v)=d Lnf(y; v) / dv is the Score function and differentiations is with respect to v, provided that, f'(y; v) exist, and f(y; v) is positive for all v in V.
The Score function approach can be extended in estimating the second and higher order of derivatives. For example, an estimate for the second derivative based on the Score Function method is:
J'' (v0) = E Z(yi).H(yi ; v0),
where,
H(yi; v0) = f ' ' (yi; v0) / f (yi; v0).
Where S and H = S' + S2 are the score and information functions, respectively, widely used in statistics literature, such as in the construction of Cramer-Rao bounds. By having gradient and Hessian in our disposal, we are able to construct a second order local metamodel using the Taylor's series.
An Illustrative Numerical Example: For most complex reliability systems, the performance measures such as mean time to failure (MTTF) are not available in analytical form. We resort to Monte Carlo Simulation (MCS) to estimate MTTF function from a family of single-parameter density functions of the components life with specific value for the parameter. The purpose of this section is to solve the inverse problem, which deals with the calculation of the components' life parameters (such as MTTF) of a homogeneous subsystem, given a desired target MTTF for the system. A stochastic approximation algorithm is used to estimate the necessary controllable input parameter within a desired range of accuracy. The potential effectiveness is demonstrated by simulating a reliability system with a known analytical solution.
Consider the coherent reliability sub-system with four components component 1, and 2 are in series, and component 3 and 4 also in series, however these two series of components are in parallel, as illustrated in the following Figure. All components are working independently and are homogeneous; i.e., manufactured by an identical process, components having independent random lifetimes Y1, Y2, Y3, and Y4, which are distributed exponentially with rates v = v0 = 0.5.
The system lifetime is
Z (Y1,Y2,Y3,Y4; v0) = max [min (Y3,Y4), min (Y1,Y2)].It is readily can be shown that the theoretical expected lifetime of this sub-system is
J(v0) = 3/(4 v0). The underlying pdf for this system is:f(y; v) = v4exp(-v Syi),
Applying the Score function method, we have:
S(y) = f ' (y; v) / f(y; v) = 4/v - S yi,
and
H(y) = f ' ' (y; v) / f(y; v) = [v2 (S yi)2 - 8v (S yi) + 12] / v2,
The estimated average lifetime and its derivative for the nominal system with v = v0 = 0.5, are:
J(v0 ) = S max [min (Y3, j,Y4, j), min (Y1, j, Y2, j)] / n,
J'(v0 ) = S max [min (Y3, j,Y4, j), min (Y1, j, Y2, j)] . S(Yi,j) / n,
and
J"(v0 ) = S max[min (Y3, j,Y4, j), min (Y1, j, Y2, j)] . H(Yi, j)/n,
respectively, where Yi, j is the jth observation for the ith component (i = 1, 2, 3, 4). We have performed a Monte Carlo experiment for this system by generating n = 10000 independent replications using SIMSCRIPT II.5 random number streams 1 through 4 to generate exponential random variables Y1, Y2, Y3, Y4 , respectively, on a VAX system. The estimated performance is J(0.5) = 1.5024, with a standard error of 0.0348. The first and second derivatives estimates are -3.0933 and 12.1177 with standard errors of 0.1126 and 1.3321, respectively.
The response surface approximation in the neighborhood of v = 0.5 is:
J(v) = 1.5024 + (v - 0.5) (-3.0933) + (v - 0.5)2 (12.1177)/2 =
A numerical comparison based on exact and the approximation by this metamodel reveals that the largest absolute error is only 0.33% for any v in the range of [0.40, 0.60]. This error could be reduced by either more accurate estimates of the derivatives and/or using a higher order Taylor expansion. A comparison of the errors indicates that the errors are smaller and more stable in the direction of increasing v. This behavior is partly due to the fact that lifetimes are exponentially distributed with variance 1/v. Therefore, increasing v causes less variance than the nominal system (with v = 0.50).
For a given J(v) the estimated dv, using the first order approximation is:
dv = [J(v)-J(v0)] / J' (v0),
provided that the denominator does not vanish for all v0 in set V.
The Goal-seeker Module: The goal-seeking problem can be solved as a simulation problem. By this approach, we are able to apply variance reduction techniques (VRT) used in the simulation literature. Specifically, the solution to the goal-seeking problem is the unique solution of the stochastic equation J(v) - J0 = 0. The problem is to solve this stochastic equation by a suitable experimental design, to ensure convergence. The following is a Robbins- Monro (R-M) type technique.
vj+1 = vj+dj [J0 - J(vj)] / J'(vj),
where dj is any divergent sequence of positive numbers. Under this conditions, dv =J0 - J(vj) converges to approach zero while dampening the effect of the simulation random errors. These conditions are satisfied, for example, by the harmonic sequence dj = 1/j. With this choice, the rate of reduction of di is very high initially but may reduce to very small steps as we approach the root. Therefore, a better choice is, for example dj = 9 / (9 + j). This technique involves placing experiment i+1 according to the outcome of experiment i immediately preceding it, as is depicted in the following Figure:
Under these not unreasonable conditions, this algorithm will converge in mean square; moreover, it is an almost sure convergence. Finally, as in Newton's root-finding method, it is impossible to assert that the method converges for just any initial v = v0, even though J'(v) may satisfy the Lipschits condition over set V. Indeed, if the initial value v0 is sufficiently close to the solution, which is usually the case, then this algorithm requires only a few iterations to obtain a solution with very high accuracy.
An application of the goal-seeker module arises when we want to adapt a model to satisfy a new equality constraint (condition) for some stochastic function. The proposed technique can also be used to solve integral equations by embedding the Importance Sampling techniques within a Monte Carlo sampling.
One may extend the proposed methodology to the inverse problems with two or more unknown parameters design by considering two or more relevant outputs to ensure uniqueness. By this generalization we could construct a linear (or even nonlinear) system of stochastic equations to be solved simultaneously by a multidimensional version of the proposed algorithm. The simulation design is more involved for problems with more than a few parameters.
For more, visit the Web site: Modeling & Simulation Resources.
References and Further Readings:
Arsham H., The Use of Simulation in Discrete Event Dynamic Systems Design, Journal of Systems Science, 31(5), 563-573, 2000.
Arsham H., Input Parameters to Achieve Target Performance in Stochastic Systems: A Simulation-based Approach, Inverse Problems in Engineering, 7(4), 363-384, 1999.
Arsham H., Goal Seeking Problem in Discrete Event Systems Simulation, Microelectronics and Reliability, 37(3), 391-395, 1997.
Batmaz I., and S. Tunali, Small response surface designs for metamodel estimation, European Journal of Operational Research, 145(3), 455-470, 2003.
Ibidapo-Obe O., O. Asaolu, and A. Badiru, A New Method for the Numerical Solution of Simultaneous Nonlinear Equations, Applied Mathematics and Computation, 125(1), 133-140, 2002.
Lamb J., and R. Cheng, Optimal allocation of runs in a simulation metamodel with several independent variables, Operations Research Letters, 30(3), 189-194, 2002.
Simpson T., J. Poplinski, P. Koch, and J. Allen, Metamodels for Computer-based Engineering Design: Survey and Recommendations, Engineering with Computers, 17(2), 129-150, 2001.
Tsai C-Sh., Evaluation and optimisation of integrated manufacturing system operations using Taguch's experiment design in computer simulation, Computers And Industrial Engineering, 43(3), 591-604, 2002.
Simulation continues to be the primary method by which system analysts obtain information about analysis of complex stochastic systems. In almost all simulation models, an expectedvalue can express the system's performance. Consider a system with continuous parameter v Î V, where V is the feasible region. Let
J(v) = EY | v [Z (Y)]
be the steady state expected performance measure, where Y is a random vector with known probability density function (pdf), f(y; v) depends on v, and Z is the performance measure.
In discrete event systems, Monte Carlo simulation is usually needed to estimate J(v) for a given value v. By the law of large numbers
J(v) = 1/n S Z (y i)
where yi, i = 1, 2,..., n are independent, identically distributed, random vector realizations of Y from f (y; v ), and n is the number of independent replications. This is an unbiased estimator for J(v) and converges to J(v) by law of large numbers.
There are strong motivations for estimating the expected performance measure J(v) for a small change in v to v + dv, that is to solve the so-called "what if" problem.
The simulationist must meet managerial demands to consider model validation and cope with uncertainty in the estimation of v. Adaptation of a model to new environments also requires an adjustment in v.
An obvious solution to the "what if" problem is the Crude Monte Carlo (CMC) method, which estimates J(v+ dv) for each v separately by rerunning the system for each v+ dv. Therefore costs in CPU time can be prohibitive The use of simulation as a tool to design complex computer stochastic systems is often inhibited by cost. Extensive simulation is needed to estimate performance measures for changes in the input parameters. As as an alternative, what people are basically doing in practice is to plot results of a few simulation runs and use a simple linear interpolation/extrapolation.
In this section we consider the "What-if" analysis problem by extending the information obtained from a single run at the nominal value of parameter v to the closed neighborhood. We also present the use of results from runs at two or more points over the intervening interval. We refer to the former as extrapolation and the latter as interpolation by simulation. The results are obtained by some computational cost as opposed to simulation cost. Therefore, the proposed techniques are for estimating a performance measure at multiple settings from a simulation at a nominal value.
J(v+ dv) = EY | v + dv[Z (Y)] = EY | v [Z (Y).W]
where the likelihood ratio W is:
W = f(y; v+ dv) / f(y; v)
adjusts the sample path, provided f(y; v) does not vanish. Notice that by this change of probability space, we are using the common realization as J(v).
The generated random vector y is roughly representative of Y, with f(v). Each of these random observations, could also hypothetically came from f(v+ dv). W weights the observations according to this phenomenon.
Therefore, the "What-if" estimate is:
J(v) = 1/n S Z (y i).Wi
which is based on only one sample path of the system with parameter v and the simulation for the system with v+ dv is not required.
Unfortunately LR produces a larger variance compared with CMC. However, since E(W)=1, the following variance reduction techniques (VRT) may improve the estimate.
J(v) = S Z (y i).Wi / S Wi
S(y) = d Ln f(y; v) / dv
We consider the exponential (approximation) model for J(v+ dv) in a first derivative neighborhood of v by:
J(v+ dv) = E [ Z(y). exp[dvS(y)] / E[exp(dS(y))]
Now we are able to estimate J(v+ dv) based on n independent replications as follows:
J(v+ dv) = S Ai / S Bi
where,
Ai= Z(yi).exp[dv.S(yi)]
and
Bi= exp[dv. S(yi)].
J(v + dv) = J(v) + dv.J' (v) + ...,
where the prime denotes derivative. This metamodel approximates J(v + dv)) for small dv. For this estimate, we need to estimate the nominal J(v) and its first derivative. Traditionally, this derivative is estimated by crude Monte Carlo; i.e., finite difference, which requires rerunning the simulation model. Methods which yield enhanced efficiency and accuracy in estimating, at little additional cost, are of great value.
There are few ways to obtain efficiently the derivatives of the output with respect to an input parameter as presented earlier on this site. The most straightforward method is the Score Function (SF). The SF approach is the major method for estimating the performance measure and its derivative, while observing only a single sample path from the underlying system. The basic idea of SF is that the derivative of the performance function, J'(v), is expressed as expectation with respect to the same distribution as the performance measure itself.
Therefore, for example, using the estimated values of J(v) and its derivative J'(v), the estimated J(v + dv) is:
J(v + dv) = J(v) + dv J'(v0)
with variance
Var[J(v+dv) ] = Var[J(v)] + (dv)2 Var[J'(v)] + 2dv Cov[J(v), J'(v)].
This variation is needed for constructing a confidence interval for the perturbed estimate.
J(v) = EY | v [Z (Y)] = f . EY | v [Z (Y)] +
Similar to the Likelihood Ratio approach, this can be written as:
J(v) = EY | v [Z (Y)] = f . EY | v1 [Z (Y).W1] +
where the likelihood ratios W1 and W2 are W1 = f(y; v) / f(y; v1) and W2 = f(y; v) / f(y; v2), respectively.
One obvious choice is f = f(y; v1) / [f(y; v1)+f(y; v2)]. This method can easily extended to k-point interpolation.
For 2-point interpolation, if we let f to be constant within the interval [0, 1], then the linear interpolated "what-if" estimated value is:
J(v) = f . J1 + (1-f) . J2
where the two estimates on the RHS of are two independent Likelihood Ratio extrapolations using the two end-points.
We define f* as the f in this convex combination with the minimum error in the estimate. That is, it minimizes
Var[J(v)] = f2Var[J1] + (1-f)2Var[J2)]
By the first order necessary and sufficient conditions, the optimal f is:
f* = Var[J2] / {Var[J1] + Var[J2]}
Thus, the "best linear" interpolation for any point in interval [v1, v2] is:
J(v) = f* . J1 + (1- f*) . J2
which is the optimal interpolation in the sense of having minimum variance.
Yi=(1/v) Ln (1/Ui)
where Ln is the natural logarithm and Ui is a random number distributed Uniformly [0,1]. In the case of perturbed v, the counterpart realization using the same Ui is
Yi=[1/(v+dv)] Ln (1/Ui).
Clearly, this single run approach is limited, since the inverse transformation is not always available in closed form. The following Figure illustrates the Perturbation Analysis Method:
Since the Perturbation Analysis Approach has this serious limitation, for this reason, we presented some techniques for estimating performance for several scenarios using a single-sample path, such as the Likelihood Ratio method, which is illustrated in the following Figure.
Research Topics: Items for further research include:
i) to introduce efficient variance reduction and bias reduction techniques with a view to improving the accuracy of the existing and the proposed methods;
ii) to incorporate the result of this study in a random search optimization technique. In this approach one can generate a number of points in the feasible region uniformly distributed on the surface of a hyper-sphere each stage the value of the performance measure is with a specified radius centered at a starting point. At estimated at the center (as a nominal value). Perturbation analysis is used to estimate the performance measure at the sequence of points on the hyper-sphere. The best point (depending whether the problem is max or min) is used as the center of a smaller hyper- sphere. Iterating in this fashion one can capture the optimal solution within a hyper-sphere with a specified small enough radius. Clearly, this approach could be considered as a sequential self-adaptive optimization technique.
iii) to estimate the sensitivities i.e. the gradient, Hessian, etc. of J(v) can be approximated using finite difference. For example the first derivative can be obtained in a single run using the Likelihood Ratio method as follows:
J'(v) = [S(Zi Wi - Zi)]/(n dv),
J'(v) = [S(Zi Wi - Zi)]/(dv. SWi),
or
J(v) = SZi[Wi+ - Wi-] / { 2 S [(1+ dv)Wi+ - (1-dv)Wi-]}
the sums are over all i, i =1, 2, 3,.., n, where
Wi+ = f(y:v+ v)/f(y; v), and Wi- = f(y; v- v)/f(y; v).
The last two estimators may induce some variance reductions.
iv) Other interpolation techniques are also possible. The most promising one is based on Kriging. This technique gives more weight to `neighboring' realizations, and is widely used in geo-statistics.
Other items for further research include some experimentation on large and complex systems such as a large Jacksonian network with routing that includes feedback loops in order to study the efficiency of the presented technique.
For more, visit the Web site: Modeling & Simulation Resources.
References & Further Readings:
Arsham H., Performance Extrapolation in Discrete-event Systems Simulation, Journal of Systems Science, 27(9), 863-869, 1996.
Arsham H., A Simulation Technique for Estimation in Perturbed Stochastic Activity Networks, Simulation, 58(8), 258-267, 1992.
Arsham H., Perturbation Analysis in Discrete-Event Simulation, Modelling and Simulation, 11(1), 21-28, 1991.
Arsham H., What-if Analysis in Computer Simulation Models: A Comparative Survey with Some Extensions, Mathematical and Computer Modelling, 13(1), 101-106, 1990.
Arsham H., Feuerverger, A., McLeish, D., Kreimer J. and Rubinstein R., Sensitivity analysis and the what-if problem in simulation analysis, Mathematical and Computer Modelling, 12(1), 193-219, 1989.
PDF Version
The Copyright Statement: The fair use, according to the 1996 Fair Use Guidelines for Educational Multimedia, of materials presented on this Web site is permitted for non-commercial and classroom
purposes only.
This site may be mirrored intact (including these notices), on any server
with public access. All files are available at http://home.ubalt.edu/ntsbarsh/Business-stat for mirroring.
Kindly e-mail me your comments, suggestions, and concerns. Thank you.
Professor Hossein Arsham
This site was launched on 2/11/1995, and its intellectual materials have been thoroughly revised on a yearly basis. The current version is the 9th Edition. All external links are checked once a month.
Back to:
Dr Arsham's Home Page
EOF: Ó 1995-2015.
In the corporate world we know that having a well drilled sales force is a key aspect of achieving successful sales targets. If your sales team works remotely then this is critically important to the success of your sales strategy. Having a team of competent and adaptable staff will give your sales team the ability to:
Sales training is critically important to the success of product launches, consistent sales performances and meeting corporate targets. There are many different methods of training a workforce but the most important method for your business is the one that gives the best results.
Sales simulation training is based on David Kolb’s Experiential Learning Model. It is a learning method which is simple to understand and is backed by scientific studies. The basic theory behind Kolb’s Experiential Model is described as simply as learning through experience.
Learning through experience with sales simulation is exactly the type of training that will turn lackluster sales employees into a well drilled sales force in the shortest time possible.
With that in mind, these are my 5 killer sales simulation features which make sales simulation training the smart choice for your business.
Traditional forms of training often involve a sales training presentation in a classroom environment where the sales team are trained through repetition and simple visual observation teaching.
This type of training does very little to stimulate the learning centers of the brain, the sales training pitch could be very useful, but unless the trainees are stimulated then the pitch only lasts a short moment in the short term memory before being forgotten.
Practice makes perfect and by putting each sales team through a set of scenarios they can learn in an interactive and stimulating manner. People learn best by making mistakes, not by having their hand held throughout a boring presentation.
The difference between traditional sales training methods and newer sales simulation is that the latter forces your sales team to think for themselves. Sales simulation lets your staff experience the consequences of their mistakes instead of a trainer spelling it out. There is no right or wrong answer in these scenarios, your sales team will get to experience their mistakes and identify areas of improvement for themselves, but most of all, they’ll remember it for future real life sales opportunities.
Your sales staff are a team of individuals, you will likely have a diverse spectrum of staff which have grown within the business over a number of years or have recently joined the team. This will eventually develop into wider inconsistencies in training. The training that new staff receive is more likely fairly evolved and modern. In comparison, some longer reigning staff who may be stuck in their ways but still lead with experience over newer staff will have been trained in early versions of your sales strategy. You can bring consistency to your sales training by using eLearning simulation scenarios which are easy to deploy even through mobile phones.
Sales simulation uses custom built scenarios to help drill consistency into a team of individuals. Over time the team will have developed bad habits and their way of doing things, but by using a range of scenarios on the team as a whole they can begin to identify areas of improvement for each other and follow a planned sales strategy more closely.
Bringing your sales team together as a unit is extremely important not just for sales consistency but for identifying areas of improvement in your sales strategy. You need to be able to trust your sales performance and analysis data, but you can only do that if you can trust that your sales team is working consistently.
During the simulation scenario all data is recorded for future analysis and use. In order to further develop your team you will have a plethora of data which will allow you to identify and give accurate feedback on weak points in your sales team’s performance.
You need your team to work towards the same goals, they can best achieve this by working together as a unit. Sales simulation creates an openness amongst your team and shows them that they have the same goals to work towards and how they should work together to achieve them. It is common for more experienced staff to become proud and unwilling to help others, but by drilling the team to work together and admit but work on their mistakes they can improve their team performance targets.
"To learn from their experience, teams must create a conversational space where members can reflect on and talk about their experience together." - David Kolb
Forging a sales team is not about getting your employees to do what you want, it is about creating a culture of success. Simply telling an employee what they did wrong is not an effective method of improvement, but helping them to experience their mistakes through a simulation scenarios is. You need your team to be a unit that shares burdens and strives towards meeting targets, sales simulation is a key part of forging an an effective unit such as this.
Staff who have been great assets for many years can sometimes become stubborn in their ways. What they have been doing for the past years has worked well for them, but this can be a sticking point when implementing new products and new strategies.
With these experienced staff members you need to show them through sales simulation scenarios that there are other ways to improve and hit even better performance targets. Simply telling them what they should be doing is not enough to inspire and motivate them to learn something new. But if you can show them through scenarios how easily they could improve then you may find that a new product launch strategy is more easily adopted.
Some staff members will have only been in a few different roles within the company in their entire career, this can create conflict and a lack of empathy between other roles and departments. Sales simulation can allow you to place those staff into roles they’ve never been in before. This will help them gain experience and understanding of what the staff they often liaise with are required to do on a regular basis. Doing so will develop empathy between the staff and greater understanding, thus helping communication between departments.
Traditional classroom and presentation based training methods are mostly a one-way street. There is no room for evaluating staff competencies in the way that a sales simulation scenario can. Identifying staff competencies is an integral part of understanding what it is your staff need help with.
Through the use of sales simulation training, your management team can evaluate your sales team’s understanding and their capability to execute on your company's sales process.
Your staff will engage with virtual customers and be given different options for responses. There are no wrong answers, the purpose of these scenarios is to allow your staff to go back and correct their mistakes, learning by doing.
Technology is monitoring and recording every step of the process which can help identify your strong and weak teams and the same for aspects of your sales strategy.
Key leaders in your team will become apparent in ways you may not have noticed before. This will give you opportunities to plan and rearrange sales teams to give balance to your sales force.
To recap, these are just some of the advantages of sales simulation over traditional sales training:
With the advent of eLearning, sales simulation training has become more accessible and affordable than ever. eLearning scenarios allow you to simulate online customers easily and create a range of different scenarios for employee to train with. You can even provide eLearning training to your team using their mobile phone, allowing them to brush up on competencies and improve their sales performance whenever they want.