[R] getting started on Bayesian analysis

[Jari Oksanen]

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/