Simulating Normal Distribution with Polar Method

The polar method allows for quick generation of normally distributed random variables. For each trial, generate r^2 as exponentially distributed with λ = 1/2, and θ uniformly distributed between 0 and 2π. Within each trial, in fact 2 independent copies of the normally distributed random variables are generated: Z1 = sqrt(r^2)*cos(θ) and Z2 = sqrt(r^2)*sin(θ), both with μ = 0 and σ = 1. For normally distributed random variables with different parameters of μ and σ, simply use the transformation σ*Z+μ.

n = input('Input the number of copies: ');
u = input('Input the mean: ');
s = input('Input the standard deviation: ');
results = zeros(1,n);
for i = 1:2:n
    rsq = -2*log(rand);
    theta = 2*pi*rand;
    results(i) = s*sqrt(rsq)*cos(theta)+u;
    results(i+1) = s*sqrt(rsq)*sin(theta)+u;
end
disp 'Calculated mean and standard deviation using the polar method:';
[mean(results), std(results)]

The program at the end calculates the μ and σ values of the n copies of the random variables created as a result of this simulation, stored in the array named "results."