R help - Aug 2012 - How to generate a matrix of Beta or Binomial distribution

Hi folks,

I have a question about how to efficiently produce random numbers from Beta
and Binomial distributions. 

For Beta distribution, suppose we have two shape vectors shape1 and shape2.
I hope to generate a 10000 x 2 matrix X whose i th rwo is a sample from
reta(2,shape1[i]mshape2[i]). Of course this can be done via loops:

for(i in 1:10000)
{
X[i,]=rbeta(2,shape1[i],shape2[i])
}

However, the above code is time consuming. It would be better to directly
generate X without any loops. I tried the following code

X<- rbeta(2,shape1,shape2)

but it looks not right. So I was wondering if there is an R code doing this.

For Binomial distribution, I have a similar question.  I hope to generate a
10000 x n matrix Y whose i th row is a sample from rbinom(n,2,p[i]), where p
is a vector of binomial probabilities. We can also do it by loops

for(i in 1:10000)
{
Y[i,]=rbinom(n,2,p[i])
}

But it would be nice to do it without using any loops. Is this possible in
R?

Many thanks for great help and suggestions!
Jeff









--
View this message in context:
http://r.789695.n4.nabble.com/How-to-generate-a-matrix-of-Beta-or-Binomial-distribution-tp4641422.html
Sent from the R help mailing list archive at Nabble.com.

Hello,

I'm glad it helped.
Inline.
Em 27-08-2012 19:00, Zuofeng Shang escreveu:
> Hi Rui,
>
> Many thanks for your excellent code! That works very well. But I still 
> have a kind of naive question about the set.seed(1)
>
> Consider shape1=c(1,3) and shape2=c(2,4)
>
> Using R I got that
>
> > set.seed(1)
> > rbeta(2,c(1,3),c(2,4))
> [1] 0.5803922 0.4679171
>
> > set.seed(1)
> > c(rbeta(1,1,2),rbeta(1,3,4))
> [1] 0.5803922 0.4679171
>
> They have the same outputs. So rbeta(2,shape1,shape2) is the correct 
> solution which is equivalent to operating rbeta() componentwise to 
> shape1 and shape2. This is exactly the solution I wanted. Thanks a lot!
>
> But I also found the following:
>
> set.seed(1)
> > rbeta(1,1,2)
> [1] 0.5803922
>
> > set.seed(1)
> > rbeta(1,3,4)
> [1] 0.3016368
>
> This is pretty strange. Why is the second one 0.3016368 instead of 
> being 0.4679171? Probably I was wrong somewhere here.
No, you are not wrong, and no, it is not strange. You are restarting the 
random number generator so your _first_ random beta number cannot be 
eqaul to the _second_ one. Just look:

set.seed(1)
rbeta(1,3,4) # 0.3016368

set.seed(1)
rbeta(2,3,4) # 0.3016368 0.4679171

The second is, as expected, what you had. When yo change the beta 
parameters, the numbers must also change (obvious!)

Rui Barradas
>
> The binomial code works very well. Thanks for this very much!
>
> Have a nice day.
> Jeff
>
> ? 8/27/2012 1:25 PM, Rui Barradas ??:
>> Hello,
>>
>> Try the following.
>>
>> set.seed(1)
>> X <- rbeta(2*length(shape1), rep(shape1, each = 2), rep(shape2, each
>> = 2))
>> X <- matrix(X, ncol = 2, byrow = TRUE)
>>
>> #-- binomial
>>
>> n <- 10
>> p <- runif(10000)
>>
>> Y1 <- matrix(nrow = 10000, ncol = 10)
>>
>> set.seed(6292)
>> for(i in 1:10000){
>>       Y1[i, ] <- rbinom(n, 2, p[i])
>> }
>>
>> set.seed(6292)
>> Y2 <- rbinom(n*length(p), 2, rep(p, each = 10))
>> Y2 <- matrix(Y2, ncol = n, byrow = TRUE)
>>
>> identical(Y1, Y2) #TRUE
>>
>> The trick is not to use argument recycling.
>>
>> Hope this helps,
>>
>> Rui Barradas
>>
>> Em 27-08-2012 15:49, JeffND escreveu:
>>> Hi folks,
>>>
>>> I have a question about how to efficiently produce random numbers 
>>> from Beta
>>> and Binomial distributions.
>>>
>>> For Beta distribution, suppose we have two shape vectors shape1 and
>>> shape2.
>>> I hope to generate a 10000 x 2 matrix X whose i th rwo is a sample
from
>>> reta(2,shape1[i]mshape2[i]). Of course this can be done via loops:
>>>
>>> for(i in 1:10000)
>>> {
>>> X[i,]=rbeta(2,shape1[i],shape2[i])
>>> }
>>>
>>> However, the above code is time consuming. It would be better to 
>>> directly
>>> generate X without any loops. I tried the following code
>>>
>>> X<- rbeta(2,shape1,shape2)
>>>
>>> but it looks not right. So I was wondering if there is an R code 
>>> doing this.
>>>
>>> For Binomial distribution, I have a similar question.  I hope to 
>>> generate a
>>> 10000 x n matrix Y whose i th row is a sample from
rbinom(n,2,p[i]),
>>> where p
>>> is a vector of binomial probabilities. We can also do it by loops
>>>
>>> for(i in 1:10000)
>>> {
>>> Y[i,]=rbinom(n,2,p[i])
>>> }
>>>
>>> But it would be nice to do it without using any loops. Is this 
>>> possible in
>>> R?
>>>
>>> Many thanks for great help and suggestions!
>>> Jeff
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>> -- 
>>> View this message in context: 
>>>
http://r.789695.n4.nabble.com/How-to-generate-a-matrix-of-Beta-or-Binomial-distribution-tp4641422.html
>>> Sent from the R help mailing list archive at Nabble.com.
>>>
>>> ______________________________________________
>>> R-help at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide 
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>
>

On Mon, Aug 27, 2012 at 9:49 AM, JeffND <Zuofeng.Shang.5 at nd.edu>
wrote:
> Hi folks,
>
> I have a question about how to efficiently produce random numbers from Beta
> and Binomial distributions.
>
> For Beta distribution, suppose we have two shape vectors shape1 and shape2.
> I hope to generate a 10000 x 2 matrix X whose i th rwo is a sample from
> reta(2,shape1[i]mshape2[i]). Of course this can be done via loops:
>
> for(i in 1:10000)
> {
> X[i,]=rbeta(2,shape1[i],shape2[i])
> }
>
> However, the above code is time consuming. It would be better to directly
> generate X without any loops. I tried the following code
>
> X<- rbeta(2,shape1,shape2)
No, not quite: this only makes 2 random variates, while before you
made 2*10,000 = 20k. You can make them all in one fell swoop like
this:

X <- rbeta(20000, shape1, shape2)

Now you might ask, "shape1" and "shape2" are only 10k long:
how does
this work? R will recycle shape1 and shape2 as needed, so you get a
call as if you actually used

rbeta(20000, c(shape1, shape1), c(shape2, shape2))

or more generally

rbeta(20000, rep(shape1, length.out = 20000), rep(shape2, length.out = 20000))

If I understand your code correctly, you want each row to use the same
shape1/shape2 parameters, so we will need to reshape the result of
rbeta() into a matrix: it turns out that this is actually pretty easy
to do for your case:

matrix(rbeta(20000, shape1, shape2), ncol = 2)

R will get the correct number of rows by default (it does the 20000
elements / 2 columns = 10000 rows division internally) and because R
uses column-major ordering, you get the parameter arrangement by happy
coincidence. (Basically, R fills up the first column first, then the
second; this aligns nicely with the recycling behavior in this
problem, so we get back to shape1[1] as soon as we start filling the
second column)

This is a little bit of happy hackery, and you might instead prefer
something more explicit like this:

matrix(rbeta(20000, rep(shape1, each = 2), rep(shape2, each = 2)),
ncol = 2, byrow = TRUE)

which will repeat shape1 and shape2 elementwise two times (so c(1,2,3)
becomes c(1,1,2,2,3,3)) and then fills rowwise. It's a little clearer,
at the expense of verbosity.
>
> but it looks not right. So I was wondering if there is an R code doing
this.
>
> For Binomial distribution, I have a similar question.  I hope to generate a
> 10000 x n matrix Y whose i th row is a sample from rbinom(n,2,p[i]), where
p
> is a vector of binomial probabilities. We can also do it by loops
>
> for(i in 1:10000)
> {
> Y[i,]=rbinom(n,2,p[i])
> }
This should follow pretty easily from the above.

Hope this helps,
Michael
>
> But it would be nice to do it without using any loops. Is this possible in
> R?
>
> Many thanks for great help and suggestions!
> Jeff
>
>
>
>
>
>
>
>
>
> --
> View this message in context:
http://r.789695.n4.nabble.com/How-to-generate-a-matrix-of-Beta-or-Binomial-distribution-tp4641422.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

Dear Michael,

Thanks very much for your explicit explanation! That makes much sense.  

Best wishes,
Jeff



--
View this message in context:
http://r.789695.n4.nabble.com/How-to-generate-a-matrix-of-Beta-or-Binomial-distribution-tp4641422p4641472.html
Sent from the R help mailing list archive at Nabble.com.