4
$\begingroup$

This might be a dumb question. I'm working on 2.11 in BDA (page 2 in this pdf http://www.uio.no/studier/emner/matnat/math/STK4021/h16/undervisningsmateriale/presentation_jariek.pdf). The problem is to find the posterior $\theta|y$ in the following model using a grid approximation: $$y|\theta\sim\textsf{Cauchy}(\theta,1),\ \theta\sim\textsf{Unif}(0,100).$$ On pages 3-4, it gives code for how to do this, along with how to sample from this approximate posterior (which I'll reproduce here, slightly altered):

dens <- function (y, th){
    dens0 <- NULL
    for (i in 1:length(th))
        dens0 <- c(dens0, prod (dcauchy (y, th[i], 1)))
    dens0
}
y <- c(43, 44, 45, 46, 46.5, 47.5)
step <- .01
theta <- seq(0, 100, step)
dens.unnorm <- dens(y,theta)
dens.norm <- dens.unnorm/(step*sum(dens.unnorm))
plot (theta, dens.norm, ylim=c(0,1.1*max(dens.norm)),
      type="l", xlab="theta", ylab="normalized density",
      xaxs="i", yaxs="i", cex=2, col="red", xlim=c(30,70))
thetas <- sample (theta, 1000, dens.unnorm, replace=TRUE)
hist(thetas, xlab="theta", yaxt="n", cex=2, add=TRUE, freq=F)

This produces the following plot, where the red curve is the approximate posterior and the histogram is comprised of 1000 samples from this approximate posterior: enter image description here

What I don't understand is the line

dens.norm <- dens.unnorm/(step*sum(dens.unnorm))

where instead of just normalizing the sampled points by the sum of the sampled points, they multiply this sum by the step size of the grid (thus the density sums to 1/step instead of 1). In other sources I've found (for example McElreath's Statistical Rethinking pages 39-40), when doing a grid approximation they just normalize by the sum of the sampled points.

The resulting posterior approximation is definitely correct, as I've verified it in Stan:

data 
{
  int<lower=0> N;
  real y[N];
}
parameters 
{
  real theta;
}
model 
{
  theta ~ uniform(0, 100);
  y ~ cauchy(theta,1);
}

which produces these samples: enter image description here

I've tried changing the line in question to

dens.norm <- dens.unnorm/(sum(dens.unnorm))

but this produces a density that's not properly normalized (even though it now sums to 1, it's the red that's barely visible at the bottom): enter image description here

Anyone have any insight into what's going on here?

$\endgroup$
2
$\begingroup$

Change the step to something larger, e.g. step <- 0.5 and then plot it as a bar plot (i.e. plot(..., type = "h")), now look at the resulting plot.

bar plot of the grid data

The plot shows "samples" taken on a regular grid from the continuous density. The probability densities are unnormalized, to normalize them you would need to divide them by the total area under the curve, $\int f(x | \theta)\, \pi(\theta) \,d \theta$.

Look again at the above plot, do you have a curve and the area beneath it? No, you have just a set of points, so to calculate the area under the curve you are missing everything between the vertical bars. So instead $\int f(x | \theta) \,\pi(\theta) \,d \theta$, you do $\sum_{i \in G} f(x | \theta_i) \,\pi(\theta_i)$ where $G$ is just a small, finite subset of infinite number of the possible $\theta$'s. This is not the same.

If you had samples from this distribution and used histogram to estimate the empirical density, then you would split the samples in some number of boxes and then divide them by the sizes of the boxes and the total count, to get probability densities, i.e. "probabilities per foot". In the grid case you don't have samples scattered randomly, yet each of the "samples" represents a point in some neighborhood, so in fact you are approximating the continuous density with a "histogram", but with fixed points. So you are approximating the continuous curve with something looking more or less like the curve below.

enter image description here

To normalize the "histogram" you need to divide the "counts" by areas of the neighborhoods the points represent. Instead of dividing by the "masses" of the points, you divide by the areas of the rectangles they represent.

$\endgroup$
  • $\begingroup$ Ah ok so it's basically just approximating the integral using a Reimann sum, thanks! $\endgroup$ – aleshing Dec 21 '17 at 19:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.