Stagnant: a changepoint problem and an illustration of how NOT to do MCMC!)

Carlin, Gelfand and Smith (1992) analyse data from Bacon and Watts (1971) concerning a changepoint in a linear regression.

Note the repeated x's.

We assume a model with two straight lines that meet at a certain changepoint xk --- this is slightly different from the model of Carlin, Gelfand and Smith (1992) who do not constrain the two straight lines to cross at the changepoint. We assume

Yi   ~   Normal(μi, τ)
μi   =   α + βJ[i] (xi - xk)   J[i]=1 if i <= k    J[i]=2 if i > k

giving E(Y) = α at the changepoint, with gradient β1 before, and gradient β2 after the changepoint. We give independent "noninformative'' priors to α, β1, β2 and τ.

Note: alpha is E(Y) at the changepoint, so will be highly correlated with k. This may be a very poor parameterisation.

Note way of constructing a uniform prior on the integer k, and making the regression
parameter depend on a random changepoint.
model
{
for( i in 1 : N ) {
Y[i] ~ dnorm(mu[i],tau)
mu[i] <- alpha + beta[J[i]] * (x[i] - x[k])
J[i] <- 1 + step(i - k - 0.5)
punif[i] <- 1/N
}
tau ~ dgamma(0.001,0.001)
alpha ~ dnorm(0.0,1.0E-6)
for( j in 1 : 2 ) {
beta[j] ~ dnorm(0.0,1.0E-6)
}
k ~ dcat(punif[])
sigma <- 1 / sqrt(tau)
}

Data
``` list(Y = c(1.12, 1.12, 0.99, 1.03, 0.92, 0.90, 0.81, 0.83, 0.65, 0.67, 0.60, 0.59, 0.51, 0.44, 0.43,0.43, 0.33, 0.30, 0.25, 0.24, 0.13, -0.01, -0.13, -0.14, -0.30, -0.33, -0.46, -0.43, -0.65),x = c(-1.39, -1.39, -1.08, -1.08, -0.94, -0.80, -0.63, -0.63, -0.25, -0.25, -0.12, -0.12, 0.01, 0.11, 0.11, 0.11, 0.25, 0.25, 0.34, 0.34, 0.44, 0.59, 0.70, 0.70, 0.85, 0.85, 0.99, 0.99, 1.19),N = 29) ```
Inits for chain 1
``` list(alpha = 0.2, beta = c(-0.45, -1.0), tau = 5, k = 16) ```

Inits for chain 2
``` list(alpha = 0.6, beta = c(-0.45, -1.0), tau = 5, k = 8) ```

Traces of two chains shows complete dependence on starting values

Results are hopeless - no mixing at all. Note: alpha is E(Y) at the changepoint, so will be highly correlated with k. This may be a very poor parameterisation.

TRY USING CONTINUOUS PARAMETERISATION

model
{
for(i in 1 : N) {
Y[i] ~ dnorm(mu[i], tau)
mu[i] <- alpha + beta[J[i]] * (x[i] - x.change)
J[i] <- 1 + step(x[i] - x.change)
}
tau ~ dgamma(0.001, 0.001)
alpha ~ dnorm(0.0,1.0E-6)
for(j in 1 : 2) {
beta[j] ~ dnorm(0.0,1.0E-6)
}
sigma <- 1 / sqrt(tau)
#x.change ~ dunif(-1.3,1.1)
x.change ~ dunif(x[5],x[26])
}
Data
``` list(Y = c(1.12, 1.12, 0.99, 1.03, 0.92, 0.90, 0.81, 0.83, 0.65, 0.67, 0.60, 0.59, 0.51, 0.44, 0.43,0.43, 0.33, 0.30, 0.25, 0.24, 0.13, -0.01, -0.13, -0.14, -0.30, -0.33, -0.46, -0.43, -0.65),x = c(-1.39, -1.39, -1.08, -1.08, -0.94, -0.80, -0.63, -0.63, -0.25, -0.25, -0.12, -0.12, 0.01, 0.11, 0.11, 0.11, 0.25, 0.25, 0.34, 0.34, 0.44, 0.59, 0.70, 0.70, 0.85, 0.85, 0.99, 0.99, 1.19),N = 29) ```
Inits for chain 1
``` list(alpha = 0.47, beta = c(-0.45, -1.0), tau = 5, x.change = 0.5) ```
Inits for chain 2
``` list(alpha = 0.47, beta = c(-0.45, -1.0), tau = 5, x.change = -0.5) ```
Results