Help with Bayesian inference in OpenBugs

Question

I have a task that involved Bayesian inference and could use some pointers and hints. I've already got some parts figured out but others remain blurred. Also, my OpenBUGS abilities are frankly limited so advice how to concretely proceed with my problem would not hurt.

So the task at hand is following. I am supposed to model the proportion of bicycles in areas with and without bicycle roads. I have 10 and 8 observations respectively. Since the variable is essentially binomial (the vehicle observed is either a bicycle or something else) a binomial model is my choice for the likelihood function. For a prior a beta distribution that is skewed towards other vehicles (for the "no bicycle roads" category even more so) is a valid pick. (it's only natural to assume that the proportion of bicycles is smaller than the proportion of other modes of transportation so there's no need to go uninformative prior. In the area that lacks bicycle roads that should be far more pronounced still) So far so good. I've gathered that the OpenBUGS - code that (hopefully) calculates the posterior is following:

model{
x ~ dbin(px,n)
y ~ dbin(py,m)
px ~ dbeta(a,b)
py ~ dbeta(a,b)
}
list(N=,x=c(data),y=c(data))

Next I am supposed to simulate from the posterior distribution a difference of the expected values x and y (ie. posterior means) and draw a histogram from it. This is where I stumble. I don't know how to concretely do that.

Also, what does it mean for a prior distribution to be independent of parameters θy and θx? That is specifically asked. Isn't a prior distribution by definition independent of the population parameters?

Any hints how to proceed with the task and the coding would be appreciated. I could also use R to do this task but I think it's easier to do with OpenBugs.

Answer 1

########################################################################################
#These values for n, m, x, y,a, and b are just randomly chosen for illustration.  These values will come from your data
n = 100
m = 100
x = rbinom(n,1,.7)
y = rbinom(m,1,.3)
a = .2
b = .1

########################################################################################
#R Code 
#You need to pick values for your priors, i.e., a = ? and b = ?

#Set B to be any large number you want
B = 1000
diff = rep(NA,B)

for(i in 1:B){

#Obtain posterior samples from x and y
posterior.x.samples = rbeta(1000,a+sum(x),b+n-sum(x))
posterior.y.samples = rbeta(1000,a+sum(y),b+m-sum(y))

#Calculate posterior expectations
Ex = mean(posterior.x.samples)
Ey = mean(posterior.y.samples)

#Take the difference in posterior means
diff[i] = Ex-Ey

}

hist(diff,xlab="Difference in Posterior Expecttaions",prob=T,
     main="Histogram")