Simulation and Modeling

Monte Carlo simulation is a method that solves a problem by repeated random sampling. A random experiment is performed many times using random numbers drawn from the probability distributions of the inputs, and the average of the results estimates the quantity of interest. It does not use the passing of time, so it is suited to static problems and to problems with uncertainty, such as inventory, risk and numerical integration.

Steps

1. Build the model and identify the random inputs.
2. Define the probability distribution of each input and the cumulative distribution.
3. Assign ranges of random numbers to each outcome according to the cumulative probability.
4. Generate random numbers and map them to outcomes.
5. Calculate the result of each trial.
6. Repeat many times and average the results.

Example: daily demand of newspapers

Random numbers 45, 83, 12, 67, 91 give demands 20, 30, 10, 20, 30, so the average demand over five days is $(20+30+10+20+30)/5=22$. With thousands of trials the average approaches the expected demand $10(0.2)+20(0.5)+30(0.3)=21$.

Another example: estimating $\pi$ by throwing random points in a unit square; the fraction that falls inside the quarter circle $x^2+y^2\le1$ approaches $\pi/4$.

The accuracy improves as the number of trials grows, with the error falling as $1/\sqrt n$.

Monte Carlo simulation is a technique that uses repeated random sampling to estimate the behaviour of a system or the value of a quantity that is hard to compute directly. Random numbers are used to draw values from the probability distributions of the uncertain inputs, the model is evaluated for each sample, and the results are averaged. It is generally used for static models (time plays no role) such as risk analysis, inventory, queues without a clock, and numerical integration.

Steps

1. Define the problem and the random variables.
2. Determine the probability distribution and the cumulative distribution of each variable.
3. Assign random number intervals to the outcomes.
4. Generate random numbers and read off the corresponding outcomes.
5. Compute the result for each trial and repeat for many trials.
6. Average the results to obtain the estimate.

Example: daily sales of a shop

Random numbers for 5 days: 52, 06, 83, 29, 71 give sales 7, 5, 8, 6, 7, so the average is $33/5=6.6$ units. The exact expected value is $5(0.1)+6(0.3)+7(0.4)+8(0.2)=6.7$. With more trials the simulated average comes closer to 6.7.

Classic example: estimating $\pi$ by choosing random points $(x,y)$ in the unit square; the fraction of points with $x^2+y^2\le1$ approaches $\pi/4$.

The error of the estimate decreases in proportion to $1/\sqrt n$.

Monte Carlo simulation is a method that uses repeated random sampling from probability distributions to estimate the behaviour of a system or the value of a quantity. It is a static simulation (time does not play a role) and is used for problems that are random or too complex to solve analytically, such as numerical integration, risk analysis and estimating $\pi$.

Steps

1. Define the problem and the random input variables.
2. Obtain the probability distribution of each input variable and build its cumulative distribution.
3. Assign ranges of random numbers to each outcome, according to the cumulative probability.
4. Generate random numbers and map them to outcomes.
5. Compute the result of each trial and repeat many times.
6. Average the results to obtain the estimate.

Example: daily newspaper demand. Demand and its probability are: 10 papers (0.2), 20 (0.5), 30 (0.3). The cumulative probabilities and random number ranges are:

Using two-digit random numbers 15, 73, 48, 91, 02 gives demands 10, 30, 20, 30, 10. The average demand is $(10+30+20+30+10)/5=20$ papers per day. The theoretical expected demand is $0.2(10)+0.5(20)+0.3(30)=21$, and the estimate approaches it as the number of trials increases.

Advantages: simple, handles random inputs and complex systems. Disadvantage: needs many trials for accuracy, and results are only approximate.

Queuing system

A queuing system consists of customers (jobs, packets, patients) that arrive, wait in a queue if the server(s) are busy, receive service, and leave. Its components are the arrival pattern (e.g. Poisson with rate $\lambda$), the service mechanism (service-time distribution, rate $\mu$, number of servers), the queue capacity (finite or infinite), the calling population (finite or infinite) and the queue discipline. It is described in Kendall notation $A/B/c/K/N/D$, e.g. M/M/1 (Poisson arrivals, exponential service, one server). Examples: bank counters, CPU job scheduling, a call centre.

Queuing disciplines

The queue discipline is the rule that decides which waiting customer is served next.

1. FIFO (FCFS): first in, first out; customers are served in order of arrival (bank, ticket counter).
2. LIFO (LCFS): the last to arrive is served first (stack of items, some storage systems).
3. SIRO (random selection): the next customer is chosen at random.
4. Priority discipline: customers with higher priority are served first, either pre-emptive (service of a lower-priority customer is interrupted) or non-pre-emptive (current service finishes first); e.g. emergency patients.
5. SPT / SJF (shortest processing time first): customers with the shortest service time go first, which reduces the average wait (CPU scheduling).
6. Round robin / processor sharing: each customer gets a fixed time slice in turn.
7. Balking, reneging and jockeying are customer behaviours: refusing to join a long queue, leaving after waiting, or switching to another queue.

Performance measures

For an M/M/1 queue with utilisation $\rho = \lambda/\mu < 1$:

• Server utilisation: the fraction of time the server is busy, $\rho = \lambda/\mu$.
• Average number in the system: $L = \dfrac{\rho}{1-\rho} = \dfrac{\lambda}{\mu-\lambda}$.
• Average number in the queue: $L_q = \dfrac{\rho^2}{1-\rho} = \dfrac{\lambda^2}{\mu(\mu-\lambda)}$.
• Average time in the system: $W = \dfrac{1}{\mu-\lambda}$.
• Average waiting time in the queue: $W_q = \dfrac{\lambda}{\mu(\mu-\lambda)}$.
• Probability of $n$ customers in the system: $P_n = (1-\rho)\rho^n$; $P_0 = 1-\rho$ is the probability that the system is idle.
• Probability that a customer must wait, the maximum queue length, throughput and the proportion of customers lost (for finite capacity).

These are related by Little's law, $L = \lambda W$ and $L_q = \lambda W_q$.

Difference between chi-square and KS tests for uniformity

KS test

$H_0$: the numbers are uniformly distributed on $[0,1]$; $H_1$: they are not. $N=6$, $\alpha=0.05$, critical value $D_\alpha = 0.565$.

Rank the numbers: $R_{(1)}=0.11,\ 0.45,\ 0.54,\ 0.68,\ 0.73,\ R_{(6)}=0.98$.

$D^+ = \max\{i/N - R_{(i)}\}$ and $D^- = \max\{R_{(i)} - (i-1)/N\}$.

$D^+ = 0.1033$, $D^- = 0.2833$, so

$D=\max(D^+,D^-)=0.2833$

Since $D = 0.2833 < D_\alpha = 0.565$, we fail to reject $H_0$. There is no evidence against uniformity, so the numbers can be accepted as uniformly distributed at the 5% level.

Static mathematical model

A mathematical model represents a system by mathematical relations (equations, inequalities, functions) among its variables. A static model describes the system at a single point in time, or in equilibrium, and has no time variable; the relationships do not change with time and the output does not depend on previous states. A dynamic model, in contrast, shows how the system evolves with time.

Examples:

• Break-even model: profit $P = (s - v)x - F$, where $s$ is selling price per unit, $v$ variable cost per unit, $F$ fixed cost, $x$ units. Break-even point is $x = F/(s-v)$.
• Ohm's law $V = IR$.
• Area/volume formulas, or a linear programming model for product mix: maximise $Z = 3x_1 + 5x_2$ subject to $x_1 + 2x_2 \le 10$, $x_1,x_2 \ge 0$.
• Economic order quantity $Q = \sqrt{2DS/H}$.

Stochastic vs deterministic activities

Merits of system simulation

• Allows study of systems that are too costly, dangerous or impossible to experiment with directly (aircraft, nuclear plants, large networks).
• Experiments are fast and cheap: time can be compressed or expanded, and "what-if" alternatives tested without disturbing the real system.
• Handles complex systems with randomness and many interacting parts that analytic models cannot solve.
• Gives insight into which variables matter and how the system behaves (sensitivity analysis).
• Provides a controlled environment where conditions can be repeated exactly.
• Useful for training (flight simulators) and for design and planning before the system is built.
• Can be used where no mathematical solution exists.

Demerits

• Building a model can be expensive and time-consuming, and needs skilled people.
• Results are only as good as the model and data; wrong assumptions give misleading answers.
• Simulation gives estimates, not exact or optimal answers; results of stochastic models vary with each run and need many runs and statistical analysis.
• It does not optimise; it only evaluates the alternatives tried.
• Results can be hard to interpret, and may give a false sense of precision.
• Validation against the real system can be difficult.
• Simulation is sometimes used when a simpler analytic solution would do.

Markov chain

A Markov chain is a stochastic process $\{X_n\}$ with a finite or countable set of states in which the probability of moving to the next state depends only on the present state and not on the earlier history (the Markov property):

$P(X_{n+1}=j \mid X_n=i, X_{n-1},\dots,X_0) = P(X_{n+1}=j \mid X_n=i) = p_{ij}$

The $p_{ij}$ are transition probabilities, collected in the transition matrix $P=[p_{ij}]$, where every row sums to 1. If the state probability row vector at step $n$ is $\pi^{(n)}$, then $\pi^{(n+1)}=\pi^{(n)}P$. For a regular chain, $\pi^{(n)}$ approaches a steady-state vector $\pi$ with $\pi = \pi P$ and $\sum\pi_i = 1$.

Example: weather model

States: S = sunny, R = rainy. Suppose

$P=\begin{pmatrix}0.7 & 0.3\\ 0.4 & 0.6\end{pmatrix}$

(rows give today's weather: Sunny, Rainy). If today is sunny, the chance of rain the day after tomorrow is found from

$P^2=\begin{pmatrix}0.7(0.7)+0.3(0.4) & 0.7(0.3)+0.3(0.6)\\ \cdots\end{pmatrix}$

so first row of $P^2$ is $(0.61,\ 0.39)$; the probability of rain in two days is 0.39.

Steady state: $\pi_1 = 0.7\pi_1 + 0.4\pi_2$ gives $0.3\pi_1 = 0.4\pi_2$, and with $\pi_1+\pi_2=1$:

$\pi_1=\frac{4}{7}\approx 0.571,\qquad \pi_2=\frac{3}{7}\approx 0.429$

In the long run, about 57% of days are sunny and 43% rainy, whatever today's weather. Markov chains are used in queueing models, PageRank, inventory and reliability problems.

Arrival pattern

The arrival pattern describes how entities (customers, jobs, packets) arrive at a system: how many arrive and when. It is specified by the inter-arrival time distribution or by the number of arrivals per unit time. Arrivals may be deterministic (at fixed intervals) or random, in single units or in batches, and may be influenced by the state of the system (balking, reneging). The most common random model is the Poisson process: arrivals occur independently at a constant average rate $\lambda$, the number of arrivals in time $t$ is

$P(N(t)=n)=\frac{e^{-\lambda t}(\lambda t)^n}{n!}$

and the inter-arrival times are exponentially distributed with mean $1/\lambda$. This is a stationary process because $\lambda$ does not change with time.

Non-stationary Poisson process

In many real systems the arrival rate changes with time (rush hours in a bank, peak traffic on a network, lunch hours in a restaurant). A non-stationary (non-homogeneous) Poisson process has an arrival rate that is a function of time, $\lambda(t)$. Properties:

• Arrivals in non-overlapping intervals are independent.
• The expected number of arrivals in $(0,T)$ is $\Lambda(T)=\int_0^T \lambda(t)\,dt$, and the number of arrivals in $(a,b)$ is Poisson with mean $\int_a^b\lambda(t)\,dt$.
• It cannot be simulated by simply using a constant mean inter-arrival time. A common method is thinning: generate arrivals at the maximum rate $\lambda^* \ge \lambda(t)$ and accept an arrival at time $t$ with probability $\lambda(t)/\lambda^*$. Another is to treat the rate as piecewise constant over intervals.

Validation vs calibration

Verification is different from both: it checks that the model is implemented correctly in the program (debugging).

How to perform validation

1. Face validity: Ask experts and users of the real system whether the model and its results appear reasonable.
2. Validate the model assumptions: Check structural assumptions (e.g. the system is correctly simplified) and data assumptions (collect data, fit distributions, test with goodness-of-fit tests such as chi-square or KS).
3. Validate input-output transformations: Run the model with the real-world input data and compare the output with the real system's output, using statistical tests such as the t-test on means or confidence intervals; use historical data and, if possible, data not used for calibration.
4. Sensitivity analysis: Change the inputs and check that the outputs change in the way the real system would.
5. Extreme-condition tests and Turing test: Check behaviour under extreme values, and see whether experts can distinguish model output from real output.
6. Iterate: If the model is not accurate enough, calibrate it and validate again until the model is credible.

The mixed congruential method generates numbers with

$X_{i+1}=(aX_i+c)\bmod m,\qquad R_i=\frac{X_i}{m}$

Here $a=17$, $c=43$, $m=100$ and the seed is $X_0=27$.

The seed 27 reappears at $i=4$, so the sequence repeats with a period of 4:

$27,\ 2,\ 77,\ 52,\ 27,\ 2,\ 77,\ 52,\dots$

The random numbers generated are $0.02,\ 0.77,\ 0.52,\ 0.27,\dots$ and then they repeat. The short period shows why the parameters must be chosen carefully; a full-period generator needs suitable $a$, $c$ and $m$.

Replication of runs means running the same simulation model several times, each time with the same model and input parameters but with a different, independent stream of random numbers (different seeds). Each run is one replication and gives one independent observation of the output measure, for example the mean waiting time $\bar X_j$ of run $j$.

Why it is necessary

1. A single run is only one sample. The output of a stochastic simulation is itself a random variable. One run may be unusually good or bad by chance, so it cannot be trusted as the true answer.
2. Independent observations are obtained. The observations inside one run are autocorrelated (for example successive waiting times in a queue), so the ordinary variance formula does not apply to them. The replication means $\bar X_1,\dots,\bar X_n$ are independent and identically distributed, so standard statistics can be used.
3. Confidence intervals can be built. With $n$ replications, the point estimate is $\bar X=\frac1n\sum \bar X_j$ and the sample variance is $s^2=\frac1{n-1}\sum(\bar X_j-\bar X)^2$. The confidence interval is

$\bar X\pm t_{\alpha/2,\,n-1}\frac{s}{\sqrt n}$

4. Precision improves. The width of the interval shrinks as $n$ increases, so the accuracy of the estimate can be controlled.
5. Initial bias can be handled by deleting the warm-up period in each replication.

Inverse transform method. Non-uniform random variates are generated from uniform random numbers $R\sim U(0,1)$ by inverting the cumulative distribution function (CDF).

Steps

1. Find the CDF $F(x)$ of the required distribution.
2. Generate a uniform random number $R$ in $(0,1)$.
3. Set $F(X)=R$.
4. Solve for $X$: $X=F^{-1}(R)$.

Why it works. Since $F$ is increasing and $R$ is uniform, $P(X\le x)=P(F^{-1}(R)\le x)=P(R\le F(x))=F(x)$. So $X$ has the CDF $F$.

Example: exponential distribution with mean $1/\lambda$. The CDF is $F(x)=1-e^{-\lambda x}$ for $x\ge0$. Put $F(X)=R$:

$1-e^{-\lambda X}=R\ \Rightarrow\ X=-\frac1\lambda\ln(1-R)$

Since $1-R$ is also uniform, $X=-\frac1\lambda\ln R$ can be used. If the mean inter-arrival time is 5 min ($\lambda=0.2$) and $R=0.5$, then $X=-5\ln0.5=3.47$ min.

Discrete case. Build the cumulative probability table and return the value $x_i$ for which $F(x_{i-1})<R\le F(x_i)$.

Other examples: uniform $U(a,b)$ gives $X=a+(b-a)R$, and Weibull gives $X=\alpha[-\ln(1-R)]^{1/\beta}$.

Limitation. The method needs a CDF that can be inverted in closed form. For distributions such as the normal, other methods (convolution, acceptance-rejection) are used.

Assumptions. The inspection time of part B is not given, so it is taken as $20\pm5$ min. Each part type has its own inspector. The time unit is 1 minute. 10% of the parts are type B, and each inspector rejects 10% of what is inspected.

Block diagram

        GENERATE 10                (one part every 10 min)
              |
        TRANSFER .100,,INSB        (10% -> INSB, 90% -> next block)
          |                                   |
   SEIZE INSA                          INSB: SEIZE INSPB
   ADVANCE 6,2                               ADVANCE 20,5
   RELEASE INSA                              RELEASE INSPB
   TRANSFER .100,GOOD,REJ                    TRANSFER .100,GOOD,REJ
          |  90%             10%                    | 90%          10%
          v                   v                     v               v
        GOOD: TERMINATE 1   REJ: TERMINATE 1      (same blocks)

GPSS program

        GENERATE   10
        TRANSFER   .100,,INSB
        SEIZE      INSA
        ADVANCE    6,2
        RELEASE    INSA
        TRANSFER   .100,GOOD,REJ
INSB    SEIZE      INSPB
        ADVANCE    20,5
        RELEASE    INSPB
        TRANSFER   .100,GOOD,REJ
REJ     TERMINATE  1
GOOD    TERMINATE  1
        START      100

In TRANSFER .100,A,B, a fraction 0.1 of transactions go to B and the rest go to A. In the first TRANSFER the empty first location means the next block, so 90% of the parts (type A) go to inspector A and 10% (type B) jump to INSB. After inspection 10% of parts go to REJ and 90% to GOOD. Both end blocks decrement the termination counter, so the run stops after 100 parts. The output gives the utilisation of each inspector and the number of rejected and accepted parts.

1. System and its environment

A system is a collection of entities (people, machines, etc.) that act and interact together towards the accomplishment of some logical end. Example: a bank with customers, tellers, and queue.

Components of a system: entities, attributes, activities, state, events, and endogenous and exogenous variables.

The environment is everything outside the system that affects it or is affected by it, but is not part of it. A change that occurs outside the system is said to be in the environment. For the bank, the arrival rate of customers, the state of the economy, or the bank's rules from the central bank form the environment. The boundary between system and environment depends on the purpose of the study. Activities in the environment produce exogenous events (e.g. customer arrival), while activities within the system produce endogenous events (e.g. completion of service). A system closed to the environment is a closed system; one that interacts with it is an open system.

2. Simulation run statistics

These are the performance measures collected during a simulation run in order to judge the behaviour of the system. They are of two kinds:

• Observation-based (tally) statistics, based on individual observations, for example the average waiting time $\bar W=\frac1n\sum W_i$, the maximum waiting time, and the number of customers served.
• Time-persistent statistics, based on the time a state variable holds a value, for example average queue length $\bar L=\frac1T\int_0^T L(t)\,dt$ and server utilisation $\rho=\frac{\text{busy time}}{T}$.

In GPSS, the standard output reports facility utilisation, number of entries, average time per transaction, queue maximum and average contents, and total and zero entries. Because a run gives a sample, statistics from several replications are combined with confidence intervals.

Queuing system. A queuing system consists of customers (jobs) who arrive, wait in a queue if the server is busy, are served by one or more servers, and then leave. Its components are the arrival process, the queue (with its capacity and discipline such as FIFO), the service mechanism, and the calling population. Examples are a bank counter, a computer CPU scheduling jobs, or a supermarket checkout.

Kendall's notation. A queue is described by $A/B/c/K/N/D$, where

• $A$ = inter-arrival time distribution (M = Markovian/exponential, D = deterministic, G = general),
• $B$ = service time distribution (M, D, G),
• $c$ = number of parallel servers,
• $K$ = system capacity (default infinity),
• $N$ = size of the calling population (default infinity),
• $D$ = queue discipline (default FIFO).

For example, M/M/1 is a single server queue with Poisson arrivals and exponential service times, and M/M/2/10 has two servers and a capacity of 10.

Performance measures of a single-server queue (M/M/1). Let $\lambda$ be the arrival rate and $\mu$ the service rate.

These are related by Little's law, $L=\lambda W$ and $L_q=\lambda W_q$, and $W=W_q+1/\mu$. The probability that the server is idle is $1-\rho$.

Stability. The utilisation $\rho=\lambda/\mu$ determines whether the system is stable. The system reaches a steady state only if $\rho<1$, that is $\lambda<\mu$, so the server can on average serve faster than customers arrive. Then $L$, $L_q$, $W$ and $W_q$ are finite. If $\rho\ge1$ the arrivals exceed the capacity of the server, the queue grows without bound, and $L$ and $W\to\infty$ (the formulas above become negative or infinite). As $\rho\to1$ the waiting measures grow very rapidly. For example, with $\lambda=4$ and $\mu=5$, $\rho=0.8$, $L=4$, $W=1$ time unit; if $\mu$ drops to 4.1 then $\rho=0.976$ and $L\approx40$.

True random numbers are produced from unpredictable physical processes such as radioactive decay, atmospheric noise, thermal noise or the timing of keystrokes. They cannot be reproduced and have no period, but they are slow to produce and not repeatable, which makes debugging a simulation difficult.

Pseudo-random numbers are produced by a deterministic algorithm (for example a congruential generator $X_{i+1}=(aX_i+c)\bmod m$) from an initial seed. They appear random and pass statistical tests, but the sequence is fully determined by the seed.

Properties of a good generator

1. Uniformly distributed on $[0,1]$.
2. Statistically independent (no correlation between successive numbers).
3. Reproducible, so a run can be repeated with the same seed.
4. Long period (before the sequence repeats).
5. Fast and needing little memory.
6. Portable and giving the same sequence on different computers.

Kolmogorov-Smirnov test. $H_0$: the numbers are uniformly distributed on $[0,1]$. $N=6$, $\alpha=0.05$.

Sort the numbers: $R_{(i)}$ = 0.25, 0.50, 0.58, 0.64, 0.72, 0.90.

$D^+=\max\left(\tfrac iN-R_{(i)}\right)=0.1133,\qquad D^-=\max\left(R_{(i)}-\tfrac{i-1}N\right)=0.3333$

$D=\max(D^+,D^-)=0.3333$

The critical value for $N=6$ and $\alpha=0.05$ from the K-S table is $D_\alpha=0.521$.

Since $D=0.333<0.521$, the null hypothesis cannot be rejected. The numbers can be treated as uniformly distributed on $[0,1]$ (with so few numbers the test has little power).

Dynamic mathematical model. A mathematical model describes a system by equations. A dynamic model is one in which the variables change with time, so the state of the system at time $t$ depends on its earlier states. It is usually written with differential equations (continuous systems) or difference equations (discrete-time systems), and its solution is a time history of the state.

Example. The growth of a population $P(t)$ with rate constant $k$:

$\frac{dP}{dt}=kP\ \Rightarrow\ P(t)=P_0e^{kt}$

Another example is a mass-spring-damper system, $m\ddot x+c\dot x+kx=F(t)$, whose displacement varies with time, or the water level in a tank, $A\frac{dh}{dt}=q_{in}-q_{out}$. A bank queue whose length changes over the day is also dynamic.

Static model. A static model represents a system at a single point in time, or in equilibrium, so time does not appear. Example: the relation $F=ma$ at an instant, the balance sheet of a firm, or the cost of an order $C=\text{fixed}+\text{unit price}\times Q$.

Phases (steps) in a simulation study

1. Problem formulation. State clearly the problem, the objectives and the questions to be answered.
2. System definition and model conceptualisation. Define the system boundary, entities, attributes, activities and the performance measures. Start with a simple model and add detail later.
3. Data collection. Collect data about the real system (inter-arrival and service times) and fit probability distributions. This runs in parallel with model building.
4. Model translation. Code the model in a general-purpose language (C, Java, Python) or a simulation language (GPSS, Arena, SIMULA).
5. Verification. Check that the program works as intended, that is, that the model is built right (debugging, tracing).
6. Validation. Check that the model represents the real system with sufficient accuracy, by comparing its output with real data.
7. Experimental design. Decide the run length, warm-up period, number of replications and the alternatives to be simulated.
8. Production runs and analysis. Run the simulation and estimate performance measures with confidence intervals.
9. More runs? Decide from the analysis whether more experiments are needed; if so, return to step 7.
10. Documentation and reporting. Document the model, program and results, and present the conclusions.
11. Implementation of the findings.

The steps are iterative; the analyst may go back to earlier phases when verification or validation fails.

Discrete event simulation (DES). In DES the state of the system changes only at distinct, separate points of time, called events, and stays constant between events. Examples are the arrival of a customer or the end of a service. The simulation clock does not advance in fixed steps; it jumps from one event time to the next.

Main components:

• System state (e.g. number in queue, server busy or idle).
• Simulation clock.
• Event list, holding future events ordered by time.
• Statistical counters to collect results.
• Initialisation, timing, and event routines, and a report generator.

Working: remove the imminent event from the list, advance the clock to its time, execute the event routine to update the state and schedule new events, update the statistics, and repeat until a stopping condition is met. A queuing system is the typical DES model.

Poisson arrival pattern. Arrivals follow a Poisson process with rate $\lambda$ when they occur one at a time, independently, and at a constant average rate. The number of arrivals $N$ in an interval of length $t$ has the Poisson distribution

$P(N=n)=\frac{e^{-\lambda t}(\lambda t)^n}{n!},\qquad n=0,1,2,\dots$

with mean and variance both $\lambda t$. Equivalently, the inter-arrival times are independent and exponentially distributed with mean $1/\lambda$, $f(t)=\lambda e^{-\lambda t}$. The process has no memory, so a future arrival does not depend on the past. In simulation inter-arrival times are generated by $T=-\frac1\lambda\ln R$. For example, if $\lambda=6$ customers/hour, the mean gap is 10 min.

Verification. The process of checking that the computer program (the implemented model) correctly represents the conceptual model, that is, "building the model right". It is done by debugging, code walk-throughs, tracing, checking with simple inputs whose results are known, and checking that outputs are reasonable.

Calibration. The process of adjusting the model parameters (and re-checking the model) by comparing the model's output with the real system's data, repeatedly, until the model matches the real system closely enough. It is an iterative loop: compare, revise, compare again.

Validation. The process of determining whether the model is an accurate representation of the real system for the intended purpose, that is, "building the right model". The output of the model is compared with actual system data using statistical tests, face validity (experts review), or Turing-type tests.

Accreditation. The official certification by an authority or the user that a model or simulation is acceptable for a specific purpose. It is a decision based on the evidence from verification and validation and is usually used in military and large government projects.

Order of use: verification and validation give the evidence, calibration improves the model, and accreditation is the final formal acceptance.

The multiplicative congruential method uses

$X_{i+1}=aX_i \bmod m,\qquad R_i=\frac{X_i}{m}$

Here $X_0=7$, $a=11$, $m=16$.

The seed 7 recurs at $i=4$, so the period is 4 and the sequence is

$7,\ 13,\ 15,\ 5,\ 7,\ 13,\ 15,\ 5,\dots$

The generated random numbers are $0.8125,\ 0.9375,\ 0.3125,\ 0.4375$, and then they repeat.

Estimation methods in simulation. A simulation run is a sample experiment on a stochastic model, so its output is a random variable. The true performance measure (for example the mean waiting time $\theta$) is unknown and must be estimated from the observations $Y_1,\dots,Y_n$, for example by the sample mean $\hat\theta=\bar Y=\frac1n\sum Y_i$.

Why estimation is used

1. Output is random. Different runs (different random numbers) give different results, so the output of one run is not the exact answer.
2. Only a finite run is possible. The simulation cannot be run forever, so the true long-run (steady-state) value can only be approximated.
3. To measure the precision. Estimation gives not only a point estimate but also an interval estimate $\bar Y\pm t_{\alpha/2,n-1}\,s/\sqrt n$, which shows how reliable the result is.
4. To compare alternatives. Decisions among alternative designs need estimates of the difference in performance with a confidence interval.
5. To determine the run length or the number of replications needed for a given accuracy.
6. To remove bias and correlation. Observations in a run are correlated and the initial conditions bias them, so special methods are needed (replication, batch means, regenerative method, deleting the warm-up period).

Types of estimate: point estimates (mean, variance, proportion) and interval estimates (confidence intervals). For terminating simulations the method of replications is used; for steady-state simulations batch means or the regenerative method is common.

Initial bias (the initialisation or start-up problem) arises when a steady-state simulation starts in an artificial state, usually "empty and idle" (no customers, all servers free). Early in the run the system is not representative of its long-run behaviour, so the output observations in this transient period are systematically lower or higher than the steady-state values. The estimate of the steady-state mean is then biased.

Importance of eliminating it

1. Accurate estimates. The bias is removed, so the sample mean is closer to the true steady-state value. For example, a queue started empty underestimates the average waiting time.
2. Valid decisions. Wrong estimates may lead to wrong conclusions about the capacity or performance of the system.
3. Reliable confidence intervals. The interval is centred around the correct value, so its coverage is as stated.
4. Better comparison of alternatives, since each is judged on its steady-state performance.
5. Shorter runs. Removing the transient part avoids the need for extremely long runs to dilute the bias.

Methods of elimination

• Deleting (warm-up period). Discard the first $d$ observations and use $\bar Y=\frac{1}{n-d}\sum_{i=d+1}^{n}Y_i$. The value of $d$ is found by plotting the output (Welch's method, moving averages).
• Intelligent initialisation. Start the system in a typical or steady-state condition, using data from the real system or a pilot run.
• Long runs, so the effect of the start-up becomes small.
• Use replications, deleting the warm-up from each.

Care is needed: deleting too little leaves the bias, and deleting too much wastes data and raises the variance.

Model. Workers arrive every $6\pm2$ min. There are two clerks (a storage of capacity 2) who take $8\pm2$ min per requisition. The requisitions then go to one storekeeper who takes $6\pm3$ min. The simulation runs for 100 requisitions. Time unit is 1 minute.

Block diagram

GENERATE 6,2                (worker/requisition arrival)
    |
QUEUE CLERKQ                (wait for a clerk)
    |
ENTER CLERKS                (take one of the 2 clerks)
    |
DEPART CLERKQ
    |
ADVANCE 8,2                 (clerk processes the requisition)
    |
LEAVE CLERKS                (free the clerk)
    |
QUEUE STOREQ                (wait for storekeeper)
    |
SEIZE STOREKPR
    |
DEPART STOREQ
    |
ADVANCE 6,3                 (storekeeper fills it)
    |
RELEASE STOREKPR
    |
TERMINATE 1                 (100 requisitions)

GPSS program

CLERKS   STORAGE   2
         GENERATE  6,2
         QUEUE     CLERKQ
         ENTER     CLERKS
         DEPART    CLERKQ
         ADVANCE   8,2
         LEAVE     CLERKS
         QUEUE     STOREQ
         SEIZE     STOREKPR
         DEPART    STOREQ
         ADVANCE   6,3
         RELEASE   STOREKPR
         TERMINATE 1
         START     100

The STORAGE 2 statement defines the two clerks as a multiple server, so up to two requisitions can be processed together. The storekeeper is a single server, a facility, served by SEIZE/RELEASE. The mean arrival interval is 6 min, the mean clerk service is 8 min with two clerks (utilisation about $8/(2\times6)=0.67$), and the mean storekeeper service is 6 min (utilisation about $6/6=1$), so the storekeeper is the likely bottleneck and the queue STOREQ is expected to grow. The output gives the utilisation of the clerks and the storekeeper and the queue statistics.

a. Digital analog simulator (DAS)

A digital analog simulator is a digital computer program that imitates an analog computer. The system is described by block diagram elements of the analog computer (integrators, summers, multipliers, coefficient potentiometers, function generators), and the program solves the equations numerically. The user writes the problem in terms of these blocks, and the program does the numerical integration step by step.

Advantages: high accuracy and repeatability, no hardware patching, no drift or noise, easy storage and modification, and the user keeps the familiar analog block-diagram viewpoint. Disadvantage: it is slower than a real analog machine, since the integrations are done sequentially. Examples: MIDAS, CSMP, DSL/90 and, in modern form, Simulink. They are used for continuous system simulation such as control systems and mechanical systems.

b. Simulation tools

Simulation can be written in three kinds of tools:

• General-purpose languages (C, C++, Java, Python). They are flexible and cheap, but all routines (event list, random numbers, statistics) must be written by the user.
• Simulation languages (GPSS, SIMSCRIPT, SIMULA, GASP). They provide the clock, event list, random number generators, statistics and report generation. GPSS is block-oriented and is suited to queuing models.
• Simulation packages (simulators) (Arena, AnyLogic, Simul8, MATLAB/Simulink, NS-2/NS-3, SimPy). They offer a graphical interface, animation, input distribution fitting and output analysis, so models are built faster with little programming.

Selection depends on the problem, cost, ease of use, speed, flexibility, animation facilities, statistical support and the skill of the user.

Transaction in GPSS. A transaction is the dynamic entity that moves through the blocks of a GPSS model. It represents a unit of traffic in the real system, such as a customer, a part, a message or a job. Transactions are created by GENERATE blocks, move from block to block, may wait in queues or be delayed by ADVANCE, and are removed by TERMINATE. Each transaction has a number, a priority and parameters that can store data, and its movement is controlled by the simulation clock and the current/future events chains.

Facility in GPSS. A facility is a single-server resource that can be used by only one transaction at a time. It is occupied by a SEIZE block and freed by a RELEASE block. If a transaction arrives when the facility is busy, it waits (in the delay chain) until the facility is released. GPSS keeps statistics for each facility: number of entries, utilisation (fraction of time busy), and average holding time per transaction. A real example is a barber, a teller or a machine.

Barbershop model. Customers arrive every $15\pm3$ min (uniform from 12 to 18). Joey cuts hair in $18\pm2$ min (16 to 20). The run is 9 hours, i.e. $9\times60=540$ minutes. Time unit is 1 minute.

Block diagram

GENERATE 15,3        customer arrives
     |
QUEUE WAITLINE       join waiting line
     |
SEIZE JOEY           occupy the barber (facility)
     |
DEPART WAITLINE      leave the line
     |
ADVANCE 18,2         haircut
     |
RELEASE JOEY         barber free
     |
TERMINATE            customer leaves

GENERATE 540         timer segment
     |
TERMINATE 1          ends the run after 540 min

GPSS program

         GENERATE   15,3
         QUEUE      WAITLINE
         SEIZE      JOEY
         DEPART     WAITLINE
         ADVANCE    18,2
         RELEASE    JOEY
         TERMINATE
         GENERATE   540
         TERMINATE  1
         START      1

The timer transaction is created once at time 540 and its TERMINATE 1 reduces the termination counter from 1 to 0, ending the simulation. About $540/15=36$ customers arrive. Since the mean service time (18 min) is larger than the mean arrival gap (15 min), the utilisation is close to 100% and the queue keeps growing, so only about $540/18=30$ haircuts are completed. The output reports the utilisation of JOEY and the queue statistics of WAITLINE.

Why the accuracy of an analog computer is low

An analog computer represents variables by continuous physical quantities such as voltages, and solves problems by electrical circuits. Its accuracy is limited, typically to about 0.1% to 1% (three or four significant digits), because:

1. Components (resistors, capacitors, amplifiers) have tolerances and change with temperature and ageing, causing drift.
2. Noise and interference affect the signals.
3. Operational amplifiers have finite gain and offsets, and integrators leak.
4. The voltage can be read only with the precision of the measuring instrument (meter or plotter).
5. The voltage range is limited, and values outside it saturate the amplifiers.
6. Scaling is needed, which restricts the resolution of small values.

A digital computer, by contrast, represents numbers in binary with as many digits as wanted, so its accuracy is limited only by word length.

Analog computer with example

Main components: summers, inverters, integrators, coefficient potentiometers, multipliers and function generators, patched together by a patch-board according to the problem. Time integration is done in real time by operational amplifier integrators.

Example: solve the equation of a spring-mass-damper system $\ddot x=-2\dot x-4x+f(t)$, or the simple decay $\frac{dy}{dt}=-ky$, $y(0)=y_0$. Rewrite as $y=-\int ky\,dt$. An integrator with initial condition $y_0$ produces $-y$; a potentiometer set to $k$ scales the signal and an inverter restores the sign. The output of the integrator is fed back through the potentiometer to its input. The voltage at the integrator output then follows $y(t)=y_0e^{-kt}$ and can be shown on an oscilloscope or recorder.

Analog versus digital computer

Poker test. The poker test checks the independence of random numbers by treating each group of digits like a poker hand and comparing the observed frequency of each pattern with the frequency expected from truly independent random digits. For four-digit numbers the patterns and their probabilities are:

Procedure

1. $H_0$: the numbers are independent; $H_1$: they are not.
2. For $N$ numbers, expected frequency of each class is $E_i=N\,p_i$.
3. Count the observed frequency $O_i$ of each pattern.
4. Compute the statistic

$\chi_0^2=\sum_i\frac{(O_i-E_i)^2}{E_i}$

5. If $\chi_0^2>\chi^2_{\alpha,\,k-1}$, reject $H_0$.

Calculation ($N=1000$)

$\chi_0^2=0.875+0.391+14.815+20.250+36.000=72.33$

The critical value is $\chi^2_{0.05}=9.49$ (as given). Since $72.33>9.49$, $H_0$ is rejected. The numbers are not independent.

Why a confidence interval is needed

The output of a stochastic simulation is random: each run gives a different value. A single number such as the sample mean of waiting times is only a point estimate and says nothing about how close it is to the true value. A confidence interval gives a range that contains the true value with a stated probability $(1-\alpha)$, so it measures the precision and reliability of the result. It is also needed to decide whether the run is long enough, and to compare alternative designs (if the intervals do not overlap, or the interval for the difference excludes zero, there is a real difference).

Establishing a confidence interval (method of replications)

1. Run the model $n$ times independently, using different random number streams, each with the same length and the warm-up removed.
2. Let $\bar X_j$ be the mean output of replication $j$. These are independent and identically distributed.
3. Compute the overall mean and the sample variance:

$\bar X=\frac1n\sum_{j=1}^n\bar X_j,\qquad s^2=\frac{1}{n-1}\sum_{j=1}^n(\bar X_j-\bar X)^2$

4. The $100(1-\alpha)\%$ confidence interval is

$\bar X\pm t_{\alpha/2,\,n-1}\frac{s}{\sqrt n}$

where $t_{\alpha/2,n-1}$ is taken from the Student t table (the normal value $z_{\alpha/2}$ for large $n$).

For a single long run, the batch means method is used instead: divide the run into batches large enough to be nearly independent and apply the same formula to the batch means. For example, with $n=10$, $\bar X=12$, $s=2$ and $t_{0.025,9}=2.262$, the 95% interval is $12\pm2.262\times2/\sqrt{10}=12\pm1.43$, that is (10.57, 13.43).

The multiplicative congruential method uses

$X_{i+1}=aX_i \bmod m,\qquad R_i=\frac{X_i}{m}$

For three-digit integers take $m=10^3=1000$, with $a=7$ and $X_0=22$.

The ten three-digit random integers are 154, 078, 546, 822, 754, 278, 946, 622, 354, 478, and the corresponding random variables are 0.154, 0.078, 0.546, 0.822, 0.754, 0.278, 0.946, 0.622, 0.354, 0.478.

Are the two statements the same? No.

• "Building a model right" is verification. It asks whether the model has been implemented correctly: does the computer program do what the conceptual model specifies, free of programming and logical errors? It compares the program with the conceptual model.
• "Building the right model" is validation. It asks whether the conceptual model is the correct representation of the real system for the purpose of the study, with the needed accuracy. It compares the model's output with the real system.

A model can be built right but still be the wrong model (a perfectly coded model of an incorrect assumption), or be the right model but wrongly coded. Both are needed, so they are different and complementary.

Importance of V&V

1. Gives confidence that the results are credible, so decision makers will accept and use them.
2. Prevents wrong decisions that could be costly or dangerous, because a faulty model gives misleading results.
3. Detects errors early, when they are cheap to correct.
4. Ensures that the model is fit for its intended purpose and the assumptions are reasonable.
5. Supports accreditation of the model and documents its limits.

Techniques: for verification, code walk-through, debugging, tracing and testing with known cases; for validation, face validity, comparison with historical data, statistical tests and sensitivity analysis. V&V is carried out throughout the project, not only at the end.

Many real systems are combined (hybrid), for example a continuous process such as steel production with discrete batches, and both types of model are needed.

A Markov chain is a stochastic process $\{X_n, n = 0,1,2,\dots\}$ with a finite or countable set of states in which the future depends only on the present state and not on how the process reached it. This is the Markov (memoryless) property:

$P(X_{n+1}=j \mid X_n=i, X_{n-1}=i_{n-1},\dots,X_0=i_0) = P(X_{n+1}=j \mid X_n=i) = p_{ij}$

Main elements

• State space: the set of values the process can take (for example Sunny, Rainy).
• Transition probability $p_{ij}$: probability of moving from state $i$ to state $j$ in one step. Each row of the transition matrix $P=[p_{ij}]$ sums to 1.
• Initial distribution $\pi^{(0)}$: probability of starting in each state.
• n-step probabilities: $\pi^{(n)} = \pi^{(0)}P^n$.
• Steady-state distribution: for an irreducible, aperiodic chain, $\pi = \pi P$ with $\sum \pi_i = 1$.

Example (weather model). Let the states be S (sunny) and R (rainy), with

$P = \begin{bmatrix} 0.8 & 0.2 \\ 0.4 & 0.6 \end{bmatrix}$

If today is sunny, $\pi^{(0)}=(1,0)$. Tomorrow: $\pi^{(1)}=(0.8,0.2)$. The day after:

$\pi^{(2)} = (0.8\times0.8+0.2\times0.4,\; 0.8\times0.2+0.2\times0.6) = (0.72,\,0.28)$

For the steady state solve $\pi_S = 0.8\pi_S + 0.4\pi_R$ with $\pi_S+\pi_R=1$. This gives $0.2\pi_S = 0.4\pi_R$, so $\pi_S = 2\pi_R$, and therefore $\pi_S = 2/3$, $\pi_R = 1/3$. In the long run about two thirds of days are sunny, whatever today's weather.

In simulation, Markov chains are used to model state-dependent random behaviour such as machine states (working/failed), queue lengths and customer brand switching.

A queueing system consists of customers who arrive, possibly wait in a queue, receive service from one or more servers and then leave. Its basic characteristics are as follows.

1. Calling population (source): the pool from which customers come. It may be infinite (arrival rate unaffected by the number in the system) or finite (arrivals fall as more customers are in the system, for example a fixed set of machines needing repair).

2. Arrival process: describes how customers arrive. It is usually specified by the inter-arrival time distribution. Arrivals are often Poisson, so inter-arrival times are exponential with mean $1/\lambda$. Arrivals may come singly or in batches, and may be scheduled or random. A customer may balk (refuse to join) or renege (leave after waiting).

3. Service mechanism: the distribution of service time, with mean $1/\mu$ (exponential, constant or general). Service may be single or batch, and its rate may depend on queue length.

4. Number of servers: single server or multiple servers in parallel, or in series (a network of queues).

5. System capacity: the maximum number of customers allowed in the system. It may be finite or infinite. When the system is full, arrivals are lost.

6. Queue discipline: the rule for choosing the next customer to serve, such as FIFO, LIFO, SIRO (random) or priority.

7. Queue behaviour and layout: a single queue or several queues, and whether jumping between queues is allowed.

These characteristics are summarised by Kendall's notation $A/B/c/K/N/D$. Performance measures include utilisation $\rho=\lambda/(c\mu)$, average waiting time and average queue length.

A physical model is a physical representation of a system, built from real objects, which behaves like the system under study. A dynamic physical model is one whose physical behaviour changes with time. It reproduces the time-varying behaviour of the real system, usually at a different scale.

Characteristics

• It is a tangible, scaled-down (or scaled-up) replica of the real system.
• Its state changes over time, so time is an essential variable.
• It is used to observe, test and predict system behaviour before building or modifying the real system.
• It needs similarity between model and system (geometric, kinematic and dynamic similarity).

Example 1: wind tunnel model of an aircraft. A scaled aircraft model is placed in a wind tunnel and air is blown past it at controlled speeds. Lift, drag and turbulence vary with air speed and angle of attack, which change over time. Engineers measure these and use them to predict the behaviour of the full-size aircraft. Without the airflow the model would be a static model.

Example 2: scaled ship or dam model. A scaled ship hull is towed through a water tank to study wave resistance. A river or dam model with flowing water is used to study flood behaviour.

Example 3: flight simulator. A cockpit on a motion platform reproduces the dynamic response of an aeroplane to pilot inputs.

Advantages: experiments are safer and cheaper than on the real system, conditions can be repeated and controlled, and designs can be tested before construction.

Limitations: a physical model is costly to build and modify, scaling effects may make results inexact, and it is not as flexible as a mathematical or computer model.

By contrast, a static physical model (for example an architectural scale model of a building) does not change with time.

(a) Hypothesis testing

Hypothesis testing is a statistical procedure for deciding, using sample data, whether a claim about a population is supported. In simulation it is used to test random number generators (uniformity, independence) and to validate models by comparing model output with real-system data.

Steps

1. State the null hypothesis $H_0$ (for example the numbers are uniform on [0,1]) and the alternative $H_1$.
2. Choose the level of significance $\alpha$, the probability of rejecting a true $H_0$ (Type I error). A typical value is 0.05.
3. Choose a test statistic (Z, t, $\chi^2$ or K-S) and compute it from the sample.
4. Find the critical value from tables for $\alpha$ and the degrees of freedom.
5. Decision: reject $H_0$ if the statistic falls in the critical region (for example $|Z_0|>Z_{\alpha/2}$); otherwise do not reject it.

A Type II error ($\beta$) is accepting a false $H_0$. Failing to reject $H_0$ does not prove it, it only means there is no evidence against it.

(b) Stationary Poisson process

A Poisson process counts random events occurring over time. It is stationary (homogeneous) when the arrival rate $\lambda$ is constant. The properties are:

• Arrivals occur one at a time.
• Numbers of arrivals in non-overlapping intervals are independent.
• The distribution of arrivals in an interval depends only on its length, not on where it starts.

The number of arrivals $N(t)$ in time $t$ has the Poisson distribution

$P(N(t)=n)=\frac{(\lambda t)^n e^{-\lambda t}}{n!},\qquad n=0,1,2,\dots$

with mean and variance both $\lambda t$. The inter-arrival times are independent and exponential with mean $1/\lambda$. A stationary process is not suitable where the rate varies, for example with rush hours, which needs a non-stationary process. Stationary Poisson arrivals are used in queueing models such as M/M/1.

Why a model of a system is built

A model is a representation of a system. Experimenting on the real system is often impossible, costly, dangerous or slow, or the system may not exist yet. A model is built to:

• understand how the system works and how its components interact;
• predict performance under different conditions and test "what-if" cases without disturbing the real system;
• compare alternative designs or policies at low cost and risk;
• study the system faster (time compression) or more slowly (time expansion) than reality;
• find bottlenecks, and train people or communicate the design.

Static model

A static model represents a system at a particular point in time, or one in which time plays no role. The relationships do not change with time and there is no time variable. Examples are a Monte Carlo model of the profit from a product with random sales and cost (one-off calculation), a balance sheet, or the equation $\text{Profit} = \text{Revenue} - \text{Cost}$.

Static versus dynamic mathematical models

Dynamic simulation models can be further classified as discrete-event or continuous, and deterministic or stochastic.

Storage in GPSS

A storage is a GPSS entity that models a resource which can serve many transactions at the same time, such as a parking lot, a warehouse or a group of identical machines. A facility serves only one transaction at a time. Its capacity is set by the STORAGE definition:

CAP  STORAGE  5

Blocks associated with storage

• ENTER A,B: a transaction takes $B$ units of storage $A$ (default 1). If space is not available it waits.
• LEAVE A,B: the transaction returns $B$ units to storage $A$, so another waiting transaction can enter.
• STORAGE (control statement): defines the capacity of the storage.
• SAVAIL / SUNAVAIL: make the storage available or unavailable.
• Standard numerical attributes: S$name (current content), R$name (remaining capacity), SA$name (average content), SC$name (entry count), SR$name (utilisation), SM$name (maximum content).

Machine and inspector problem (using a facility)

Parts are produced every 5 minutes (mean). The inspector takes $10 \pm 3$ minutes, that is uniform between 7 and 13. About 10% of parts are rejected. The inspector is a single server, so we use a facility. Since an inspection takes longer than 5 minutes on average, a queue will form.

        GENERATE   5            ; a part every 5 minutes
        QUEUE      INSQ         ; join inspector queue
        SEIZE      INSP         ; capture inspector
        DEPART     INSQ         ; leave queue
        ADVANCE    10,3         ; inspect 10 +/- 3 minutes
        RELEASE    INSP         ; free inspector
        TRANSFER   .100,,REJ    ; 10% go to REJ, 90% continue
        TERMINATE  1            ; accepted part
REJ     TERMINATE  1            ; rejected part
        START      500          ; run until 500 parts processed

Working

• GENERATE 5 creates one part (transaction) per 5 minutes.
• QUEUE/DEPART collect the waiting time statistics.
• SEIZE/RELEASE model the single inspector. A part waits if the inspector is busy.
• ADVANCE 10,3 holds the part for a time uniformly distributed in 7 to 13 minutes.
• TRANSFER .100,,REJ sends 10% of parts to label REJ and the remaining 90% to the next block.
• Both paths end in TERMINATE 1, each decrementing the termination counter by 1. START 500 stops the run when 500 parts have been processed.

At the end, the output gives the utilisation of INSP, the queue statistics and the number of parts accepted and rejected. The mean service time (10) exceeds the inter-arrival time (5), so the queue keeps growing and the inspector utilisation is close to 100%.

Two main properties of random numbers

1. Uniformity: every number in the interval [0,1) is equally likely, so $R_i \sim U(0,1)$. If the interval is divided into $n$ equal classes, each class should contain $N/n$ observations.
2. Independence: the probability of observing a value does not depend on previous values, so there is no correlation or pattern in the sequence.

Test for autocorrelation

We test numbers starting at $i=3$ with lag $m=4$ (3rd, 7th, 11th, ...) in a sequence of $N=30$ numbers.

$H_0: \rho_{im}=0$ (independent), $H_1: \rho_{im}\neq 0$ (auto-correlated). $\alpha=0.05$, $Z_{0.025}=1.96$.

Step 1: find $M$ (the largest integer such that $i+(M+1)m \le N$):

$3 + (M+1)\cdot 4 \le 30 \Rightarrow M+1 \le 6.75 \Rightarrow M=5$

The numbers used are the 3rd, 7th, 11th, 15th, 19th, 23rd and 27th: 0.23, 0.64, 0.99, 0.91, 0.75, 0.05, 0.19.

Step 2: estimate the autocorrelation

$\hat\rho_{im} = \frac{1}{M+1}\left[\sum_{k=0}^{M} R_{i+km}R_{i+(k+1)m}\right] - 0.25$

$\hat\rho_{3,4} = \frac{2.4112}{6} - 0.25 = 0.4019 - 0.25 = 0.1519$

Step 3: standard deviation

$\hat\sigma_{\rho} = \frac{\sqrt{13M+7}}{12(M+1)} = \frac{\sqrt{72}}{72} = 0.1179$

Step 4: test statistic

$Z_0 = \frac{\hat\rho_{im}}{\hat\sigma_\rho} = \frac{0.1519}{0.1179} = 1.29$

Decision: $|Z_0| = 1.29 < Z_{0.025} = 1.96$, so $H_0$ is not rejected. There is no evidence of autocorrelation, and the 3rd, 7th, 11th, ... numbers can be treated as independent.

A Markov chain is a stochastic process in which the probability of moving to the next state depends only on the present state, not on the earlier history (the Markov property):

$P(X_{n+1}=j \mid X_n=i, X_{n-1},\dots,X_0) = P(X_{n+1}=j \mid X_n=i)=p_{ij}$

The chain is described by its states and a transition matrix $P=[p_{ij}]$, where $p_{ij}\ge 0$ and every row sums to 1. If $\pi^{(0)}$ is the initial distribution, the distribution after $n$ steps is $\pi^{(n)}=\pi^{(0)}P^n$. For a regular chain, the long-run (steady-state) distribution satisfies $\pi=\pi P$ with $\sum\pi_i=1$.

Example: market share of two brands. Customers buy brand A or brand B each month. 90% of A buyers stay with A, and 20% of B buyers switch to A:

$P=\begin{bmatrix}0.9&0.1\\0.2&0.8\end{bmatrix}$

If the market is now 60% A and 40% B, $\pi^{(0)}=(0.6,0.4)$. Next month:

$\pi^{(1)}=(0.6\times0.9+0.4\times0.2,\; 0.6\times0.1+0.4\times0.8)=(0.62,\,0.38)$

Steady state: $\pi_A = 0.9\pi_A+0.2\pi_B$ gives $0.1\pi_A=0.2\pi_B$, so $\pi_A=2\pi_B$. With $\pi_A+\pi_B=1$ we get $\pi_A=2/3$ and $\pi_B=1/3$. In the long run brand A holds about 66.7% of the market.

Markov chains are used in queueing models, weather prediction, inventory and reliability models and text generation.