Task: Write a script to create a set of N uniform random numbers and normal random numbers. Calculate the mean and standard deviations of each set of numbers. To be specific, let N range from 10 to 500 in steps of 10. For each of these steps, calculate the mean and standard deviation. Plot these values in two different plots and show that more values give better statistics. (What are the true values of these statisitics if you take a very large sample?)
All computer-generated "random" numbers are actually "pseudo-random", meaning that they have the character of random numbers but are actually produced as a cycle. For the same starting point (seed value), the same "random" numbers will be produced. This also means that there are a finite number of values in the sequence. After a (relatively large) set of numbers from the sequence, the pattern will repeat.
A set of pseudo-random values from a uniform distribution in the open interval from 0 to 1 (not including these values) is created with U=rand(r,c);. This will produce a set of numbers with r rows and c columns. Be aware that the command U=rand(n); produces a matrix of values with n rows and n columns. To get a list of numbers, be sure to set either r or c to 1.
A set of pseudo-random values from a standard normal distribution is obtained with N=randn(r,c);, where again r and c give the number of rows and columns desired. The standard normal distribution has a mean of zero and a standard deviation of 1.
Because the random numbers being produced are predictable from a given starting place, it is often important to choose different starting places (a "seed" value) to ensure that a different sequence of numbers is produced.
The command rng(SD) provides a seed (the value of SD should be a non-negative integer) to the random number generator. You can get the current seed with the command seed=rng;. It is possible to randomize the random number generators with rng('shuffle') which uses some form of the current time as a seed, so every start is from a different seed value.
Monte Carlo simulations are used to gather statistics from a problem with random starting points or random parameter values. The "random walk" simulation is one example of such a calculation.
N=1000; % total number of steps
x=zero(N,1); y=zero(N,1); % location in a plane
% draw the random numbers representing step lengths in each direction
R= 2*(rand(N,2)-.5); % convert to the range (-1 to 1)
for i=2:N
% take steps of random size in random directions
x(i)=x(i-1)+R(i,1); y(i)=y(i-1)+R(i,2);
end
disp(['Final location: (' num2str(x(N)) ',' num2str(y(N)) ')'];
disp(['Average location: (' num2str(mean(x)) ',' num2str(mean(y)) ')'];
disp(['Area covered: (' num2str(std(x)) ',' num2str(std(y)) ')'];
figure
plot(x,y,'LineWidth',1.5)
title('Drunken Walk')
Given the starting location at the origin, the drunk should stay near the
origin. However, the area covered may be interestingly large.
There are many other problems of this kind that are analyzed with this sorts of randomized calculations.
Flow chart for the task:
%%% set up values of N %%% set up arrays to hold the mean and std for to two types of random number %%% loop through values of N %%% generate the random numbers %%% calculate the mean and std for each set, saving the values %%% plot the mean and std for each type of random number %%%
Solution script for the task:
%%% set up values of N
N=[10:10:500];nV=length(N);
%%% set up arrays to hold the mean and std for to two types of random number
mU=zeros(nV,1);sU=zeros(nV,1);
mN=zeros(nV,1);sN=zeros(nV,1);
%%% loop through values of N
for i=1:nV
num=N(i);
%%% generate the random numbers
U=rand(num,1);G=randn(num,1);
%%% calculate the mean and std for each set, saving the values
mU(i)=mean(U);sU(i)=std(U);
mG(i)=mean(G);sG(i)=std(G);
end
figure
%%% plot the mean and std for each type of random number
subplot(2,2,1)
plot(N,mU,'k',N,.5*ones(nV,1),'k--')
title('Mean of Uniform')
subplot(2,2,3)
plot(N,sU,'k',N,(1/sqrt(12))*ones(nV,1),'k--')
title('Std of Uniform')
subplot(2,2,2)
plot(N,mG,'k',N,zeros(nV,1),'k--')
title('Mean of Normal')
subplot(2,2,4)
plot(N,sG,'k',N,ones(nV,1),'k--')
title('Std of Normal')