Home > Bayesian computation > Learning about a correlation – part II

## Learning about a correlation – part II

Earlier we considered the problem of learning about a correlation and illustrated a normal approximation.  Here we illustrate rejection sampling to obtain a simulated sample from the exact posterior.

The main task is to find a suitable proposal density p.  We’ll let p be a t density with mean $\mu$, variance $\sigma^2$ and small degrees of freedom $\nu$.

How to choose the parameters?  First, since the posterior is approximately N(0.50, 0.009), we’ll let $\mu = 0.50$ and $\sigma^2 = 0.009$.  Also, we’ll choose the degrees of freedom $\nu = 4$.  (It is a good strategy to have a proposal density with “flat” tails which we achieve by choosing a small number of degrees of freedom.)

Next, we have to find the constant $c$ so that $g(\rho | y) \le c p(\rho)$ for all $\rho$.  On the log scale, we want to find $c$ such that

$\log g(\rho|y) \le \log c + \log p(\rho)$.

I compute a small function that computes $\log g(\rho | y) - \log p(\rho)$:

dfunction=function(rho,data)
cor.sampling(rho,data)-dmt(rho,.5,.009,4,log=TRUE)

dfunction=function(rho,data)
cor.sampling(rho,data)-dmt(rho,.5,.009,4,log=TRUE)

I plot this function over the range of values of $\rho$ from -1 to 1 — I find that

log c = -42.9927.

To see if I have found a suitable covering function, I plot the exact posterior $g(\rho | y)$ and the covering density $c p(\rho)$.  It looks pretty good.

Now we can implement rejection sampling.   We can use a LearnBayes function rejectsampling.

R=rejectsampling(cor.sampling,list(m=.5,var=.009,df=4), dmax=-42.9927,10000,dz)

The inputs are (1) the exact log posterior function, (2) the parameters of the t proposal density (mean, variance, df), (3) the value of log c, (4) the number of values we propose, and (5) the data used in the log posterior.

Here is a density estimate of the simulated draws we get — this is very close to the exact posterior density.

plot(density(R),col=”orange”,lwd=3)