# Dealing with missing data: 
# x1...xn ~ Normal(mu, 1), but values over 2 are censored! 
# We use a flat prior on mu, but could also use a normal prior. 

data <- c(1, 1.4, 1.9, 0.3, 1.5, 1.7) 
# and then there are two observations censored by 2. 
N <- 1000

# method 1, using the censored likelihood and Random Walk MH: 
loglik <- function(mu) {sum(dnorm(data, mu, 1, log=T)) + 
  2*pnorm(2, mu, 1, lower.tail=F, log.p=T)}
result1 <- rep(0, N)
for (i in 2:N) {
  prop <- result1[i-1] + rnorm(1,0,0.1)
  accept <- exp(loglik(prop)-loglik(result1[i-1]))
  if (runif(1)<accept) result1[i] <- prop
  else result1[i] <- result1[i-1]
}

# method 2, using augmented data and Gibbs sampling: 
result2 <- matrix(c(0, 3, 3), N, 3, byrow=T)
nFullData <- length(data) + 2
for (i in 2:N) {
  augmentedData <- c(data, result2[i-1,2:3])
  result2[i,1] <- rnorm(1, mean(augmentedData), sd(augmentedData)/sqrt(nFullData))
  while((result2[i,2] <- rnorm(1, result2[i,1], 1))<=2) TRUE
  while((result2[i,3] <- rnorm(1, result2[i,1], 1))<=2) TRUE
}

par(mfrow=c(2,2))
plot(result1, type="l")
for (i in 1:3) plot(result2[,i], type="l")

#######################################
## Slice sampling
#######################################

    #install.packages("mcsm")
    #library(mcsm)
    data(challenger)

    Nsim <- 1000
    x <- challenger[, 2]
    y <- challenger[, 1]
    n <- length(x)
    a <- b <- rep(0, Nsim)
    
    plot(x,y)
    
    #  m <- mean(x)
    #  x <- x - m 
    
    res <- glm(y ~ x, binomial()) 
    print(res$coef)
    a[1] <- res$coef[1]
    b[1] <- res$coef[2]
    
    for (t in 2:Nsim) {
      uni <- runif(n) * exp(y * (a[t - 1] + b[t - 1] * x))/(1 + 
                                    exp(a[t - 1] + b[t - 1] * x))
      
      z <- log(uni/(1-uni))
      mina <- max(( z-b[t-1]*x)[y==1])
      maxa <- min((-z-b[t-1]*x)[y==0])
      a[t] <- runif(1, mina, maxa)
      
      minb <- -Inf
      maxb <- Inf
      if (sum((y==1)&(x>0))>0) 
        minb <- max(minb, max((( z-a[t])/x)[(y==1) & (x>0)]))
      if (sum((y==0)&(x<0))>0)
        minb <- max(minb, max(((-z-a[t])/x)[(y==0) & (x<0)]))
      if (sum((y==1)&(x<0))>0)
        maxb <- min(maxb, min((( z-a[t])/x)[(y==1) & (x<0)]))
      if (sum((y==0)&(x>0))>0) 
        maxb <- min(maxb, min(((-z-a[t])/x)[(y==0) & (x>0)])) 
      b[t] <- runif(1, minb, maxb)
    }    
    plot(a, type="l")
    plot(b, type="l")
    plot(a,b)
    
    
    # data.frame(a=a-b*m, b=b)
    
    # data.frame(a=a,b=b)