[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")
This was off-topic. So something about business: isn't the (Win)BUGS
author working with a R port?
cheers, jari oksanen
--
Jari Oksanen -- Dept Biology, Univ Oulu, 90014 Oulu, Finland
email jari.oksanen at oulu.fi, homepage http://cc.oulu.fi/~jarioksa/
More information about the R-help
mailing list