3
$\begingroup$

I generate my data with this model:

$y=ax+b+\nu$

where $\nu$ is a random value from a random variable which follows a normal distribution with mean equal to zero and standard deviation equal to $\sigma$. The following R script generates my data.

set.seed(42)
a_t=1 #true value for a
b_t=1 #true value for b
sigma_t=1 #true value for sigma
x=1:10
noise=rnorm(length(x),0,sigma_t)
y=a_t*x+b_t+noise
plot(x,y)
abline(a_t,b_t)

enter image description here

Now I would like to "recover" from $x$ and $y$ the values of $a$, $b$ and $\sigma$ using a Bayesian approach. I will follow the grid approximation1,2

  1. Define the grid. This means you decide how many points to use in estimating the posterior, and then you make a list of the parameter values on the grid.
  2. Compute the value of the prior at each parameter value on the grid.
  3. Compute the likelihood at each parameter value.
  4. Compute the unstandardized posterior at each parameter value, by multiplying the prior by the likelihood. In this step, we are approximating the posterior density by a discrete probability distribution on the grid.
  5. Standardize the posterior, by dividing each value by the sum of all values. Finally, using the R command sample, take a random sample with replacement from the discrete distribution.
# define the grid
a=seq(0,2,by=.05)
b=seq(0,2,by=.05)
sigma=seq(.1,2,by=.01)

pars <- expand.grid(a=a,b=b,sigma=sigma)

# define the priors
pars$a_prior <- dunif(pars$a,min=0,max=2)
pars$b_prior <- dunif(pars$b,min=0,max=2)
pars$sigma_prior <- dunif(pars$sigma,min=0,max=2)
pars$prior <- pars$a_prior * pars$b_prior * pars$sigma_prior

# compute the likelihood
sz <- nrow(pars)
for(i in 1:sz){
    likelihoods <- dnorm(y,pars$a[i]*x+pars$b[i],pars$sigma[i])
    pars$likelihood[i] <- prod(likelihoods)
    print(i/sz)
}

# compute the unnormalized posterior
pars$unnormalized_posterior <- pars$likelihood*pars$prior
# normalize the posterior
pars$posterior <- pars$unnormalized_posterior/sum(pars$unnormalized_posterior)

# get the sample from the posterior
sample_indices <- sample(1:nrow(pars), size=1e4, replace = TRUE, prob=pars$posterior)
pars_sample <- pars[sample_indices,c('a','b','sigma')]

# draw a 2d histogram for a and b
library(ggplot2)
p <- ggplot(pars_sample,aes(a,b))
p <- p+stat_bin2d(bins=10)
plot(p)

At the end of the procedure I have a three columns matrix (or the pars_sample dataframe in the R script) that is a sample from the posterior distribution $p(a,b,\sigma|x,y)$; now I would like to think at $\sigma$ as a nuisance parameter (is that meaningful?) and I would like to "marginalize it away", i.e. I would like to get

$p(a,b|x,y)=\int p(a,b,\sigma|x,y) d \sigma$

Hogg et al.3 wrote:

This marginalization can be performed by direct numerical integration (think gridding up the nuisance parameters, evaluating the posterior probability at every grid point, and summing)

but I have no idea on how to proceed in R.

enter image description here

1MCELREATH, Richard. Statistical rethinking: A Bayesian course with examples in R and Stan. CRC press, 2020. Page 40.

2ALBERT, Jim. Bayesian Computation with R. New York, NY : Springer New York, 2009. Page 27.

3HOGG, David W., BOVY, Jo and LANG, Dustin. Data analysis recipes: Fitting a model to data. arXiv:1008.4686

$\endgroup$

1 Answer 1

4
$\begingroup$

It is as simple as dropping the $\sigma$ column from your pars_sample dataframe.

The three column dataframe is a sample from the joint distribution of all three parameters. Dropping $\sigma$ leaves a sample from the joint distribution of the two remaining parameters. The procedure is to generate a, b, and $\sigma$ while accounting for influence all three values, but then ignore $\sigma$ to compute marginal distribution of a and b.

This is different from the procedure you quote from Hogg et al. Their method is a modification of your prior step. Instead of sampling from the full joint posterior you would reduce the three-way joint posterior to a two-way (a and b) marginal posterior and sample from that.

The two methods are both valid. You would choose the one that is more computationally convenient in terms of desired precision, computation time, complexity, etc. From where you are it looks like the easiest thing is to generate samples of all three parameters, drop the $\sigma$, and then compute the a and b marginal distribution.

$\endgroup$
2
  • $\begingroup$ So the figure i.stack.imgur.com/yWe5N.png in the question is already a 2D histogram from $\int p(a,b,\sigma|x,y) d \sigma$, isn't it? $\endgroup$ Nov 15, 2020 at 15:17
  • 1
    $\begingroup$ Yes, that would be it. Sampling from that density gives you the marginal distribution of a and b. $\endgroup$ Nov 15, 2020 at 15:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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