generating random numbers from binomial distribution

To generate Bin(4, 3/4), use the following table:

x 0 1 2 3 4

P(X=x) 1/256 3/64 54/256 27/64 81/256

I want to write the MATLAB code to generate 5 random numbers from this distribution. I must use the following algorithm:

Algorithm:

1.Generate u from UNIF(0,1).

2.If P(X <= j-1) <= u <= P(X <= j) then set X=j for j=1,2,3,4

I have written the following code. But I think this code is a general code for generating from Binomial. How can I specialize it to my algorithm?

function X = dene(n,p,N)

X = zeros(1,N); % Generate the uniform random numbers: % N variates of n trials.

U = rand(N,n); % Loop over the rows, finding the number % less than p

for i = 1:N

ind = find(U(i,:) <= p);

X(i) = length(ind);

end

end

Answer 1

[the cyclist]

OK. I was right about cumsum, but there was another error. For each trial, when you generate the uniform random, you do not need to generate a vector of 5 of those values. You should just generate 1 value, to then compare against the cumulative value of p. I've fixed it up here. I also added a little plot to illustrate the distribution. (Hope I did not just do your homework for you.)

X = zeros(1,N); % Generate the uniform random numbers: % N variates of n trials.

U = rand(N,1); % Loop over the rows, finding the number % less than p

for i = 1:N

ind = find(U(i) <= cumsum(p));

X(i) = length(ind);

end

h = hist(X,unique(X))

count = fliplr(h)/sum(h)

figure
bar(unique(X),count)

[Atakan]

Ok Mr. Cyclist. Thank you very much. This wasn't exactly a homework. It is just a part of my project.

[Atakan]

The generated numbers are only 0 and 1. Is it normal?

[the cyclist]

When you say "generated numbers", I assume you are talking about the output "X". I do not get only 0s and 1s, so I guess we are calling the function differently. Here is what I do. First, the code as I have written above is inside the function dene.m, where the first line of the function is:

function X = dene(n,p,N)

Then, in the base workspace, I define the inputs:

n = 5; % The function does not actually use this input anymore!

N = 10000;

p = [1/256 3/64 54/256 27/64 81/256];

Then, I call dene() with the inputs:

X = dene(n,p,N);

What I get for output is a series of integers, ranging from 1 to 5. The frequency distribution is (approximately) the probability p that you input. (The ordering is reversed relative to p, which I did not quite figure out, I must admit.)

[Atakan]

I think you are right. If I assign only one value to p the output contains only 1 and 0. Also I wonder about that our code involves the 2. algorithm which is:

If P(X <= j-1) <= u <= P(X <= j) then set X=j for j=1,2,3,4

Answer 2

[the cyclist]

I have not looked at your algorithm in detail, but here's a guess. Where you have p inside your loop, I think you might want cumsum( p ) instead. Then you are checking to see if your randomly generated value is less that the cumulative probability distribution, which I think is what you want.

[Atakan]

Thank you for your suggestion. But still i cannot combine it with the algorithm.

[the cyclist]

OK. I'll try to take a closer look. Maybe you could add a little more details about why you think the results are incorrect. How many trials are you running?