Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Bernoulli random variables

Subject: Bernoulli random variables

From: Bert Ji

Date: 28 Mar, 2009 19:21:01

Message: 1 of 11

Hello, I'm fairly new to programming with Matlab and was wondering if somebody could assist me write a small code concerning Bernoulli random variables.

I need to code the following:

i = 1...n
j = 1...n

X_i and Y_ij are independent bernoulli random variables with
P[X_i = 1] = p
P[Y_ij = 1] = q

then

Z_i = X_i + (1-X_i) * (1 - Product_Sequence[for i ~= j] { 1-X_j*Y_ji})

How can I create Bernoulli random variables with specified probabilities p and q?

Also could you please advise me a nice way of vectorizing this program so it will be the most efficient?

Thanks!
Bert

Subject: Bernoulli random variables

From: Nasser Abbasi

Date: 29 Mar, 2009 01:55:21

Message: 2 of 11


"Bert Ji" <bert.ji@gmail.com> wrote in message
news:gqltat$q1g$1@fred.mathworks.com...

>
> How can I create Bernoulli random variables with specified probabilities p
> and q?
>

If you have student version which comes with statistics toolbox, then you
can use binopdf, from help:

"The Bernoulli distribution is a special case of the binomial distribution,
with n = 1."


--Nasser

Subject: Bernoulli random variables

From: Bert Ji

Date: 29 Mar, 2009 11:16:01

Message: 3 of 11

I've been able to make the bernoulli random variables using the rand function. However, I am having trouble trying to completely vectorize my program. It is awfully slow when I try to do 1 million simulations.

Could you please have a look and try to tell me where I can change things to vectorize the code?

        X = bernoulli(p(i),n(i),simulations);
        XX = bernoulli(p(i),n(i),simulations);
        
        for c=1:simulations
            Y = bernoulli(q,n(i),n(i));
            XXX = XX(:,c);
            M = 1-XXX(:,ones(n(i),1)).*Y;
            M = full(spdiags(ones(n(i)),0,M));
            Z = X(:,c) + (1-X(:,c)).*(1-prod(M)');
            N(i,c) = sum(Z);
        end
        toc

Subject: Bernoulli random variables

From: Bruno Luong

Date: 29 Mar, 2009 11:34:01

Message: 4 of 11

"Bert Ji" <bert.ji@gmail.com> wrote in message <gqnl9h$9s7$1@fred.mathworks.com>...
> I've been able to make the bernoulli random variables using the rand function. However, I am having trouble trying to completely vectorize my program. It is awfully slow when I try to do 1 million simulations.
>
> Could you please have a look and try to tell me where I can change things to vectorize the code?
>

You should describe what bernouilli() returns. Hope this avoid us for making wild guessing.

Bruno

Subject: Bernoulli random variables

From: Bert Ji

Date: 29 Mar, 2009 12:01:03

Message: 5 of 11

bernoulli(p,n,m) returns a matrix of size n x m with bernoulli random variables of probability p.

Subject: Bernoulli random variables

From: Bruno Luong

Date: 29 Mar, 2009 13:22:02

Message: 6 of 11

"Bert Ji" <bert.ji@gmail.com> wrote in message <gqnl9h$9s7$1@fred.mathworks.com>...
> I've been able to make the bernoulli random variables using the rand function. However, I am having trouble trying to completely vectorize my program. It is awfully slow when I try to do 1 million simulations.
>
> Could you please have a look and try to tell me where I can change things to vectorize the code?
>
> X = bernoulli(p(i),n(i),simulations);
> XX = bernoulli(p(i),n(i),simulations);
>
> for c=1:simulations
> Y = bernoulli(q,n(i),n(i));
> XXX = XX(:,c);
> M = 1-XXX(:,ones(n(i),1)).*Y;
> M = full(spdiags(ones(n(i)),0,M));
> Z = X(:,c) + (1-X(:,c)).*(1-prod(M)');
> N(i,c) = sum(Z);
> end
> toc

Two obvious changes are:
- Preallocate N
- Use M(1:n(i):end)=1 to set the diagonal of M rather than M = full(spdiags(ones(n(i)),0,M));

Bruno

Subject: Bernoulli random variables

From: Bert Ji

Date: 29 Mar, 2009 13:44:02

Message: 7 of 11

Do you think there is something I can do with respect to the for loop for c=1:simulations?
I think that's where the bottle neck is.

Also, when you mean preallocate N, do you mean M?

> Two obvious changes are:
> - Preallocate N
> - Use M(1:n(i):end)=1 to set the diagonal of M rather than M = full(spdiags(ones(n(i)),0,M));
>
> Bruno

Subject: Bernoulli random variables

From: Bruno Luong

Date: 29 Mar, 2009 13:58:02

Message: 8 of 11

"Bert Ji" <bert.ji@gmail.com> wrote in message <gqntv2$l4s$1@fred.mathworks.com>...
> Do you think there is something I can do with respect to the for loop for c=1:simulations?
> I think that's where the bottle neck is.

You might take a look at the profiler tool. When I run the below: it takes about 0.4s on my computer, and 93% of the spending time is for calling bernouilli to generate Y !!! So the answer is no, as long as you need to generate Y at each simulation.

p=0.2
q=1-p;
n=100
m=1000

X = bernoulli(p,n,m); % n x m
XX = bernoulli(p,n,m); % n x m
N2=zeros(1,m);

profile on
tic
for c=1:m
    Y = bernoulli(q,n,n); % n x n
    XXX = XX(:,c); % n x 1
    M = 1-XXX(:,ones(n,1)).*Y; % n x n
    M(1:n+1:end) = 1;
    Z = X(:,c) + (1-X(:,c)).*(1-all(M,1)'); % n x 1
    N2(c) = sum(Z);
end
toc
profile off
profile viewer

% Bruno

Subject: Bernoulli random variables

From: dpb

Date: 29 Mar, 2009 14:14:52

Message: 9 of 11

Bert Ji wrote:
> Do you think there is something I can do with respect to the for loop for c=1:simulations?
> I think that's where the bottle neck is.

As Bruno says, profile code before trying to optimize what isn't a problem.


> Also, when you mean preallocate N, do you mean M?

No, he meant N

Again, it may not help much given the bulk of the time is undoubtedly
being spent in Bernoulli() but you're setting each element in an
initially unallocated array which is a recipe for a copy operation each
time. Set N=zeros(size,size) before beginning the loop

--

Subject: Bernoulli random variables

From: Bert Ji

Date: 29 Mar, 2009 14:27:01

Message: 10 of 11

I see, I have allocated N earlier in the code.

Also, if it's the bernoulli function that is taking most of the time.. I guess at this point I can't do much more optimization right?

Bert

Subject: Bernoulli random variables

From: Bruno Luong

Date: 29 Mar, 2009 14:46:01

Message: 11 of 11

"Bert Ji" <bert.ji@gmail.com> wrote in message <gqo0fl$ofg$1@fred.mathworks.com>...
> I see, I have allocated N earlier in the code.
>
> Also, if it's the bernoulli function that is taking most of the time.. I guess at this point I can't do much more optimization right?

No, unless if you do something more sophisticated than

randber = rand(n,m)<=p;

in your bernouilli function

Bruno

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us