This page has been updated -- New version of the program package

Introduction

Matlab files discussed in this section: simbinom.m, simdiscr.m, simexp.m, simgeom.m, simpareto.m, ranwalk.m, brownian.m, poissonti.m, poissonjp.m, ranwalk2d.m, ranwalk3d.m, poisson2d.m, poisson3d.m, bm3plot.m
Files for downloading: gzipped tar-archive

The presentation starts from the random number generators in Matlab and a discussion on generating the basic probability distributions. Two fundamental mechanisms are random walks in discrete time and Poisson processes in continuous time. Such processes are prototypes for simple simulation algorithms based on independence. The extensions to simulation in 2 and 3 dimensions are straightforward.

Probability distributions

Uniform distribution
"Pick a number randomly in the interval [0,1]"
```
rand(1,1); rand(n,m);  
```
Zero-one distribution
"Generate sequence of m random digital numbers, 1 with prob p, 0 with prob 1-p"
```
(rand(1,m)<=p);
```
Binomial distribution
"Perform n trials where the success probability each time is p, let x be the number of successes"
```
x=sum(rand(n,m)<=p);  % x is vector of m random numbers from Bin(n,p) 
```

Normal distribution
Gaussian, or normally distributed, random numbers are generated using the Matlab command RANDN.

randn(1,m);           %vector of length m of random numbers from N(0,1)
randn(n,m);           %nxm-matrix with N(0,1) entries
mu+sigma.*randn(1,m); %m numbers from the N(mu,sigma^2)-distribution

Exponential distribution
"An event occurs randomly in time with intensity lambda, i.e. the probability the event occurs between time t and t+h is lambda h, for small h. Let t be the waiting time until the event occurs."
```
t=-log(rand)/lambda;        % random number from Exp(lambda)
t=-log(rand(1,m)./lambda);  % m random numbers from Exp(lambda)
```
Geometric distribution
"Peform trials where the success probability each time is p, let x be the number of failures before the first success"
```
x=floor(-log(rand(1,m))./(-log(1-p)));  % x is vector of m random numbers from Ge(p) 
```

Poisson distribution
Matlab Statistics Toolbox has built in command POISSRND for generating random numbers from the Poisson distribution.

poissrnd(lambda); 
poissrnd(lambda,n,m); % nxm-matrix of random numbers from the Po(lambda) distribution

If not available use e.g.

arrival=cumsum(-log(rand(1,5./lambda)));
n=length(find(arrival<=lambda));

General discrete distribution
Suppose probabilities p=[p(1) ... p(n)] for the values [1:n] are given, sum(p)=1 and the components p(j) are nonnegative. To generate a random sample of size m from this distribution imagine that the interval (0,1) is divided into intervals with the lengths p(1),...,p(n). Generate a uniform number rand, if this number falls in the jth interval give the discrete distribution the value j. Repeat m times.
```
coin=rand(1,m);
cumprob=[0 cumsum(p)];
sample=zeros(1,m);
for j=1:n
  ind=find((coin>cumprob(j)) & (coin<=cumprob(j+1)));
  sample(ind)=j;
end
```
Pareto distribution
Density function: f(x)=alpha/(1+x)^(1+alpha), x>0

This is one of the simplest examples of a distribution with so called heavy tails. Simply generate a uniform sample and apply the inverse cdf:
```
sample = (1-rand(1, npoints)).^(-1/alpha)-1;
```

Random walks

Symmetric random walk, n steps

p=0.5;  
y=[0 cumsum(2.*(rand(1,n-1)<=p)-1)]; 
plot((0:n-1),y);

Random step size random walk
Take any zero mean distribution for stepsize, e.g. uniform.
```
x=rand(1,n)-1/2;
y=[0 cumsum(x)-1)]; 
plot((0:n-1),y); 
```
Brownian motion
This is the continuous analog of symmetric random walk, each increment y(s+t)-y(s) is Gaussian with distribution N(0,t^2) and increments over disjoint intervals are independent. It is typically simulated as an approximating random walk in discrete time.
```
n=1000;
sigma=1;
y=[0 cumsum(sigma.*randn(1,n))]; % standard Brownian motion
plot((0:n),y);
```

Poisson processes

Fixed number of steps

% n simulated interarrival times from Exp(a) 
interarr=-log(rand(1,n+1))./a;  
stairs(cumsum(interarr), 0:n);

Fixed interval of time

Fix a time interval [0 tmax]. Generate random events with intensity lambda and let n be the total number of events in the interval. This number is Poisson distributed with mean value lambda*tmax. The successive times between arrivals is Exp(lambda).

tmax=100;                  % simulation interval 
lambda=1;                   % arrival intensity 
arrtime=-log(rand)/lambda;  % Poisson arrivals
i=1;                  
  while (min(arrtime(i,:))<=tmax)
    arrtime = [arrtime; arrtime(i, :)-log(rand)/lambda];  
    i=i+1;
  end
n=length(arrtime);           % arrival times t_1,...t_n
stairs(arrtime, 0:n-1);

Higher-dimensional versions

2-dimensional random walk
Plots the points (u(k),v(k)), k=1:n, in the (u,v)-plane, where (u(k)) and (v(k)) are one-dimensional random walks. The sample code plots four trajectories of the same walk, four different colors.
```
n=100000;
colorstr=['b' 'r' 'g' 'y'];          
for k=1:4
  z=2.*(rand(2,n)<0.5)-1; 
  x=[zeros(1,2); cumsum(z')];
  col=colorstr(k);
  plot(x(:,1),x(:,2),col);
  hold on
end
grid
```

3-dimensional random walk

Same as above, in 3-dimensional space.

p=0.5;
n=10000;
colorstr=['b' 'r' 'g' 'y'];          
for k=1:4
  z=2.*(rand(3,n)<=p)-1; 
  x=[zeros(1,3); cumsum(z')];
  col=colorstr(k);
  plot3(x(:,1),x(:,2),x(:,3),col);
  hold on
end
grid

3D Brownian motion

npoints = 5000;
dt = 1;
bm = cumsum([zeros(1, 3); dt^0.5*randn(npoints-1, 3)]);
figure(1);
plot3(bm(:, 1), bm(:, 2), bm(:, 3), 'k');
pcol = (bm-repmat(min(bm), npoints, 1))./ ...
         repmat(max(bm)-min(bm), npoints, 1);
hold on;
scatter3(bm(:, 1), bm(:, 2), bm(:, 3), ...
         10, pcol, 'filled');
grid on;
hold off;

Poisson points in 2D and 3D space
This is the generic model for placing points in space randomly and independently. In any fixed set of space there will be a Poisson distributed number of points with intensity proportional to the volume of that set. The number of points in any two disjoint sets are independent.
```
% Poisson points in the unit cube,
% intensity lambda 
lambda=100;
nmb=poissrnd(lambda)
x=rand(1,nmb);
y=rand(1,nmb);
z=rand(1,nmb);
grid
scatter3(x,y,z,5,5.*rand(1,nmb));
```