[R] getting started on Bayesian analysis

Jari Oksanen jarioksa at sun3.oulu.fi
Wed Sep 15 09:40:52 CEST 2004

```On Wed, 2004-09-15 at 03:27, HALL, MARK E wrote:
>   I've found
>
> Bayesian Methods: A Social and Behavioral Sciences Approach
> by Jeff Gill
>
> useful as an introduction.  The examples are written in R and S with generalized scripts for doing
> a variety of problems.  (Though I never got change-point analysis to successfully in R.)
>
Change point analysis? I haven't seen the book, but I read lecture
handouts of one Bayesian course over here in Finland (Antti Penttinen,
JyvÃ€skylÃ€), and translated his example to R during one (rare) warm
summer day in a garden. So do you mean this (binary case):

> source("/mnt/flash/cb.update.R")
> cb.update
function (y, A=1, B=1, C=1, D=1, N=1200, burnin=200)
{
n <- length(y)
lambda <- numeric(N)
mu <- numeric(N)
k <- numeric(N)
lambda[1] <- A/(A+B)
mu[1] <- C/(C+D)
k[1] <- n/2
sn <- sum(y)

for (i in 2:N) {
kold <- k[i-1]
sk <- sum(y[1:kold])
lambda[i] <- rbeta(1, A+sk, B + kold - sk)
mu[i] <- rbeta(1, C + sn - sk, D + n - sn + sk - kold  )
knew <- sample(n-1, 1)
sknew <- sum(y[1:knew])
r <- (sknew - sk) *
(log(lambda[i])-log(mu[i]))-(knew-kold)*(lambda[i]-mu[i])
if(min(0,r) > log(runif(1))) k[i] <- knew
else k[i] <- k[i-1]
}
out <- cbind(lambda, mu, k)
out[(burnin+1):N, ]
}
> y <- c(rbinom(60, 1, 0.8), rbinom(40, 1, 0.3))
> uh <- cb.update(y, N=5200)
> colMeans(uh)
lambda         mu          k
0.8189303  0.4169367 59.0770000
> mean(y[1:60])
[1] 0.7833333
> mean(y[41:100])
[1] 0.45
> plot(density(uh[,1]))
> plot(density(uh[,2]))
> plot(table(uh[,3]), type="h")