Skip to content

Is there a way to not use an explicit loop?

3 messages · Juancarlos Laguardia, Victor Hernando Cervantes Botero, Dan Davison

Juancarlos Laguardia
# 
I have a problem in where i generate m independent draws from a  
binomial distribution,
say

draw1 = rbinom( m , size.a, prob.a )


then I need to use each draw to generate a beta distribution.  So,  
like using a beta prior, binomial likelihood, and obtain beta  
posterior, m many times.  I have not found out a way to vectorize  
draws from a beta distribution, so I have an explicit for loop within  
my code



for( i in 1: m ) {

beta.post = rbeta( 10000, draw1[i] + prior.constant  ,  prior.constant  
+ size.a  - draw1[i] )

beta.post.mean[i] = mean(beta.post)
beta.post.median[i] = median(beta.post)

etc.. for other info

}

Is there a way to vectorize draws from an beta distribution?

UC Slug
Victor Hernando Cervantes Botero
# 
Hi,

you might try this:

set.seed(100)

m         <- 10
size.a    <- 10
prob.a    <- 0.3
prior.constant = 0

draw1 = rbinom( m , size.a, prob.a )

beta.draws <- function(draw, size.a, prior.constant, n) {
	rbeta(n, prior.constant + draw, prior.constant + size.a - draw)
}

bdraws <- sapply(draw1, beta.draws, size.a = size.a, prior.constant =
prior.constant, n = 10000)
beta.post <- apply(bdraws, 2, function(x) c(post.mean = mean(x),
post.median = median(x)) )
beta.post
                 [,1]      [,2]      [,3]       [,4]      [,5]
[,6]      [,7]      [,8]
post.mean   0.2017118 0.1996809 0.2991173 0.10069613 0.3001924
0.2991149 0.4033310 0.2003104
post.median 0.1804893 0.1791630 0.2845427 0.07505278 0.2858155
0.2844503 0.3961419 0.1790511
                 [,9]     [,10]
post.mean   0.3013020 0.1990232
post.median 0.2886199 0.1786447



best


V?ctor H Cervantes




2008/9/17 Juancarlos Laguardia <brassman785 at gmail.com>:
Dan Davison
# 
Both shape parameters of rbeta can be vectors; for

x <- rbeta(n, shape1, shape2)

x[i] ~ Beta(shape1[i], shape2[i])

so

bbsim <- function(m=1000, num.post.draws=1e4, size.a=100, prob.a=.27, prior.count=1) {
    data.count <- rbinom(m, size.a, prob.a)
    shape1 <- rep(prior.count + data.count, each=num.post.draws)
    shape2 <- rep(prior.count + size.a - data.count, each=num.post.draws)
    matrix(rbeta(m * num.post.draws, shape1, shape2), num.post.draws, m)
}

Then you can do

beta.draws <- bbsim()
means <- apply(beta.draws, 2, mean)
medians <- apply(beta.draws, 2, median)
etc

Dan
On Wed, Sep 17, 2008 at 11:56:36AM -0700, Juancarlos Laguardia wrote: