Some comments on the solution of project 2

2a) When \(\alpha=1\), we have a Poisson process and the number of arrivals \(n\) on the interval \([0,t]\) is then known to be a sufficient statistic for the intensity \(\lambda\), that is, given \(\lambda\), the arrival times \(y_1,y_2,\dots,y_n\) contain no further information about \(\lambda\). In this special case, the inference we would make just based on \(n\) should therefore be the the same as the inference we make from \(y_1,y_2,\dots,y_n\). This is indeed the case since \(n\sim \operatorname{Poisson}(\lambda t)\), so that the likelihood based on \(n\) only is \(L(\lambda)=\frac{(\lambda t)^n e^{-\lambda t}}{n!}\) which is identical to the likelihood based on \(y_1,y_2,\dots,y_n\) (if you have derived this correctly) up to a constant of proportionality.

2e) This task was maybe a bit tricky (but worth studying in more detail since the idea behind it is applicable in many situations). It is worth noting that we no longer have the symmetry \(Q(\alpha',\lambda'|\alpha,\lambda,\mathbf{y})=Q(\alpha,\lambda|\alpha',\lambda',\mathbf{y})\). Instead, the joint proposal (including normalising constant of the density of \(\lambda'\) conditional on \(\alpha'\)) becomes \[ Q(\alpha',\lambda'|\alpha,\lambda,\mathbf{y})\propto e^{-\frac1{2\sigma^2}(\alpha'-\alpha)^2}\frac{(y_{i+1}-y_0)^{(n+1)\alpha'+1}}{\Gamma(\alpha'(n+1)+1)}\lambda'^{(n+1)\alpha'}e^{-\lambda'(y_{i+1}-y_0)} \] This leads to (after several but not all terms in the proposal cancels against corresponding terms in the target and because of symmetry of the random walk part of the proposal) \[ \frac{\pi(\alpha',\lambda'|\mathbf{y})Q(\alpha,\lambda|\alpha',\lambda',\mathbf{y})}{\pi(\alpha,\lambda|\mathbf{y})Q(\alpha',\lambda'|\alpha,\lambda,\mathbf{y})} =e^{-(\alpha'-\alpha)} (\Gamma(\alpha')/\Gamma(\alpha))^{-(n+1)} \left(\prod_{i=1}^{n+1}(y_i-y_{i-1})\right)^{\alpha'-\alpha} (y_{i+1}-y_0)^{-(n+1)(\alpha'-\alpha)} \frac{\Gamma(\alpha'(n+1)+1)}{\Gamma(\alpha(n+1)+1)} \]

Taking the log leads to the expression in following R code. Tuning \(\sigma\) to optimal values of about .7 and 1.1, respectively, the effective size with the block update is about 3 times that of the sampler based on separate Gibbs updates for \(\lambda\) and \(\alpha\).

library(coda)
gibbs <- function(y, t, method="gibbs", sigma=.1, mcmc=1e+4) {
  n <- length(y)
  y <- c(NA,y,NA)
  
  chain <- matrix(NA, mcmc, 4)
  accepted <- 0
  
  lambda <- 1
  alpha <- 1
  for (i in 1:mcmc) {
    y[1] <- y[2] - qgamma(runif(1, min=pgamma(y[2], alpha, lambda)), alpha, lambda)
    y[n+2] <- y[n+1] + qgamma(runif(1, min=pgamma(t-y[n+1], alpha, lambda)), alpha, lambda)
    switch(method,
           "gibbs"={
             lambda <- rgamma(1, (n+1)*alpha + 1, y[n+2] - y[1])
             alphap <- rnorm(1, alpha, sigma)
             if (alphap<=0)
               accept <- 0
             else
               accept <- min(1, exp(
                 -(n+1)*(lgamma(alphap) - lgamma(alpha)) 
                 + (alphap - alpha)*((n+1)*log(lambda) - 1 + sum(log(diff(y))))))
             if (runif(1)<accept) {
               alpha <- alphap
               accepted <- accepted + 1               
             }
           },
           "block"={
             alphap <- rnorm(1, alpha, sigma)
             if (alphap<=0)
               accept <- 0
             else {
               lambdap <- rgamma(1, (n+1)*alphap + 1, y[n+2] - y[1])
               accept <- min(1, exp(
                 (alphap - alpha)*(- 1 + sum(log(diff(y))) - (n+1)*log(y[n+2] - y[1]))
                 - (n+1)*(lgamma(alphap) - lgamma(alpha))
                 + lgamma(alphap*(n + 1) + 1) - lgamma(alpha*(n + 1) + 1)
               ))
             }
             if (runif(1)<accept) {
               lambda <- lambdap
               alpha <- alphap
               accepted <- accepted + 1               
             }
           }
    )
    chain[i, ] <- c(alpha, lambda, y[1], y[n+2])
  }
  cat("Acceptance rate =", accepted/mcmc,"\n")
  colnames(chain) <- c("alpha", "lambda", "y0", "ynp1")
  mcmc(chain)
}
y <- c(5.14, 10.15, 18.4, 31.72, 35.4, 36.5, 37.95, 40.82, 64.84, 
       96.8)
t <- 100
chain <- gibbs(y, t, sigma=.7, mcmc=10000, method="gibbs")
summary(chain)
chain <- gibbs(y, t, sigma=1.1, mcmc=10000, method="block")
summary(chain)
pairs(unclass(chain), pch=".")

sigma <- seq(.1, 2, len=20)
effsize <- sapply(sigma, function(sigma) effectiveSize(gibbs(y, t, sigma=sigma, mcmc=50000, method="block")))
effsizegibbs <- sapply(sigma, function(sigma) effectiveSize(gibbs(y, t, sigma=sigma, mcmc=50000, method="gibbs")))
plot(sigma, effsize[1,])
points(sigma, effsizegibbs[1,], col="red")