Posterior distribution of normally distributed variable within interval

Question

This might be a rather simple question on Bayesian statistics.

Consider the case of inference on properties of a normal distribution $N(µ,σ²)$.

If we assume the standard non-informative prior on σ² $1/σ²$, the resulting posterior distribution $$ p(µ,σ²|y) ∝ (\frac{1}{σ²}) ^{\frac{N+1}{2}} exp(-\frac{Ns^{2}_{y}}{2σ²})\frac{1}{σ²}exp(-\frac{(y-\bar{y})²}{2σ²/N}) $$

is said to have a closed form because it factors in the product of an inverted Gamma distribution for σ² and a conditional normal distribution for µ. This means that we can directly draw from the posterior distribution of σ² via drawing from the inverted gamma distribution, and then draw from the posterior distribution of µ via drawing from the normal distribution.

This can be done with the following R code (plus relying on some summary statistics):

S <- sum((data - mean(data))^2) 
n <- length(data) 
Sigma2 <- S/rchisq(5000,n-1) 
mu <- rnorm(5000,mean=mean(data),sd=sqrt(sigma2)/sqrt(n))

There is not even a need for normalizing.

So far so good.

Now consider the case of a different form of the likelihood function:

$$ p(y|µ,σ²) = (\frac{1}{2πσ²})^{N/2} exp(-\frac{Ns^{2}_{y}}{2σ²}) exp(-\frac{(y-\bar{y}+0.5)²}{2σ²/N}) - (\frac{1}{2πσ²})^{N/2} exp(-\frac{Ns^{2}_{y}}{2σ²}) exp(-\frac{(y-\bar{y}-0.5)²}{2σ²/N}) $$

that is, the probability that a true value falls within the rounding interval that results in the observation.

The posterior density in now proportional to the difference of two productis of inverted Gamma and conditional normal distributions. Can I still use the same approach for drawing from the posterior distribution, or is that impossible because the resulting posterior distribution is not a product of well-known distributions anymore?

Answer 1

Let us assume the unobserved observations $y_i$ fall into intervals $[a_i, b_i)$, and that the observed observations (sorry, couldn't resist) are the pairs of $[a_i, b_i)$. In that case, the likelihood function for a single observation is:

$l(\mu, \sigma^2 | a, b) = \frac{1}{\sigma} \int_a^b \exp\{-\frac{1}{2\sigma^2}(y-\mu)^2\}\text{d}y$

You can't use the same approach as before, unfortunately. However, that doesn't mean there aren't any approaches which work.

One such approach is MCMC, as you can explicitly calculate the likelihood function w/o much trouble in, for example, R. I've attached some code that implements a simple sampler:

# Generate underlying data
y <- rnorm(1000)
a <- floor(y*10)/10     # Intervals are of the form (0.1, 0.2), (0.2, 0.3), ...
b <- ceiling(y*10)/10   # We ignore the computer's ability to generate, e.g., 1.0 exactly.

# Now for the algorithm
log_lf <- function(mu, sigma, a, b) {
  sum(log(pnorm(b, mu, sigma) - pnorm(a, mu, sigma)))
}

mu_current <- log_sigma_current <- 0
log_lf_current <- log_lf(mu_current, exp(log_sigma_current), a, b)

mu_mcmc <- sigma_mcmc <- rep(0,1000)

for (i in 1:1000) {  # No thinning or any attempt to optimize
  mu_proposed <- mu_current + rnorm(1,0,0.1)
  log_sigma_proposed <- log_sigma_current + rnorm(1,0,0.1)
  log_lf_proposed <- log_lf(mu_proposed, exp(log_sigma_proposed), a, b)

  # With uniform priors on mu and log(sigma), no priors to take into account
  # With symmetric proposal distributions, don't need to take the proposal
  # distribution into account when calculating alpha (it cancels)
  alpha <- exp(log_lf_proposed - log_lf_current)
  if (runif(1) < alpha) {
    mu_mcmc[i] <- mu_proposed
    sigma_mcmc[i] <- exp(log_sigma_proposed)
    mu_current <- mu_proposed
    log_sigma_current <- log_sigma_proposed
    log_lf_current <- log_lf_proposed
  } else {
    mu_mcmc[i] <- mu_current
    sigma_mcmc[i] <- exp(log_sigma_current)
  }
}

and the results...

> summary(mu_mcmc)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-0.068810 -0.007038  0.022150  0.015570  0.034990  0.106300 
> summary(sigma_mcmc)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.9539  0.9915  1.0070  1.0070  1.0220  1.0790 

Pretty good, but, then again, the sample size was large (1000) and the rounding not that great (intervals of length 0.1).

Answer 2

You can use the same sampling strategy. Draw a parameter (say, $\mu$) from the marginal distribution $p(\mu|y)$, and then draw a value of $\sigma^2$ from the conditional distribution $p(\sigma^2|\mu,y)$. This exploits the identity that $p(\mu, \sigma|y)=p(\mu|\sigma^2,y)p(\sigma^2|y)=p(\sigma^2|\mu,y)p(\mu|y)$.