Generating random numbers

I want to write a Matlab code for generating “m” numbers from N(0,1) using following algorithm. My code should be a general code for “m”, “mu” & “sigma”.

1.Generate u1, u2 from UNIF(0,1), then set y=tan(pi*(u1-1/2))

2.If u2 <= (sqrt(e)/2)*(1+y^2)*e^((-y^2)/2) then set x=y, otherwise go to step 1

3.Repeat 1-2 until you generate m numbers.

this is my code:

function [randnormal]=atakan(a,b,m) randnormal=[];

count=1;

while (count<=m) R = normrnd(a,b);

    u1=unifrnd(0,1);
    u2=unifrnd(0,1);

    y=tan(pi*(u1-1/2));

    if (u2<=((sqrt(exp(1))/2)*(1+y^2)*(exp(1)^(-y^2/2))))
        randnormal=[R;y];

    end

     count=count+1;
end

[Jan]

You forgot to ask a question.

Please format the code correctly using the information provided at the "Markup help" link.

Answer 1

[Atakan]

function x = atakan(mu,sigmasqroot,m)

count=0; x=[]; y=[];

while (count<m)

    u1=rand;
    u2=rand;

    y=tan(pi*(u1-(1/2)));
    test=((sqrt(exp(1))/2)*(1+(y^2))*((exp(1)^((-y^2)/2))));

    if (u2<=test)  
        z=sqrt(sigmasqroot)*(y+mu); % if x~N(0,1) then sigma*(x+mu)~N(mu,sigma^2) then set x=y;
        x=[x; z];  
        count=count+1;
    end

end

[the cyclist]

Just so you know, "growing" your array x in this way is extremely inefficient, computationally, because MATLAB will need to continually reallocate memory for the ever-larger array. The method in my solution, in which I preallocate the memory before the while loop, will be stupendously faster for large values of m.

Answer 2

[the cyclist]

You were overwriting your random numbers, rather than storing all "m" of them, and you were not tracking the counting of successes properly. Does this work better? (I did not check any other part of your algorithm.)

function [randnormal]=atakan(a,b,m)

randnormal=zeros(2,m);

count=1;

while (count<=m)

      R = normrnd(a,b);

      u1=unifrnd(0,1);
      u2=unifrnd(0,1);
      y=tan(pi*(u1-1/2));
      if (u2<=((sqrt(exp(1))/2)*(1+y^2)*(exp(1)^(-y^2/2))))
          randnormal(:,count)=[R;y];
          count=count+1;
      end

end

end

[Atakan]

Thanks again. I have solved it by another way.

[James Tursa]

What way? Can you post your method so others can see and not leave this thread dangling?