Posterior of $\text{Normal}(\theta,1)$ with a Cauchy prior distribution

Question

If $X \sim N(\theta,1)$ with Cauchy as robust prior

$$\pi(\theta) = \frac{1}{\pi(1+\theta^2)} \qquad -\infty < \theta < \infty$$

What will be the posterior distribution when Cauchy is $(-2 < \theta <2)$. I tried when $C(0,1)$

(2) How to describe a rejection sampler to simulate samples of $\theta$ from $\pi$($\theta$|x), with $\pi$($\theta$) as the proposal distribution and same how to do it in R

So please help to figure out.

Answer 1

What do you mean by "find?" I can tell you $\pi(\theta \mid x)$ is proportional to $$ f(x \mid \theta) \pi(\theta) \propto \exp\left[-\frac{(\theta - x)^2}{2} - \log(1 + \theta^2) \right], $$ but I don't recognize this density. You can use this fact in a number of techniques to sample from the posterior (e.g. accept-reject, importance sampling with resampling, MCMC, etc.).

By the way, your prior is a $\mathcal{C}(0,1)$ prior, which has support $-\infty < \theta < \infty$.

Edit:

You seem to be interested in accept-reject sampling. Say you're targeting $\pi(\theta \mid x) \propto f(x \mid \theta) \pi(\theta)$, but you want to simulate from the prior. First, check that $$ \frac{f(x \mid \theta) \pi(\theta)}{\pi(\theta)} \le\sup_{\theta} f(x \mid \theta) = (2\pi)^{-1/2}. $$ So you want to do the following steps over and over again:

1. Sample $\theta^i \sim \pi(\theta)$
2. Calculate $f(x \mid \theta^i) (2 \pi)^{1/2}$
3. Draw $U \sim \text{Uniform}(0,1]$
4. If $U < f(x \mid \theta^i) (2 \pi)^{1/2}$ accept the sample, otherwise throw it away and start again at step 1.

You can also sample from any other $g(\theta)$ that has support on the reals. You just have to find a $c$ (hopefully as small as possible) such that $\pi(\theta) f(x \mid \theta) \le c g(\theta)$. It works because:

\begin{align*} &P\left( \theta \le t \bigg\rvert U \le \frac{\pi(\theta) f(x \mid \theta)}{c g(\theta) } \right) \\ &= \frac{P\left( \theta \le t , U \le \frac{\pi(\theta) f(x \mid \theta)}{c g(\theta) } \right)}{P\left(U \le \frac{\pi(\theta) f(x \mid \theta)}{c g(\theta) } \right)} \\ &= \frac{\int_{-\infty}^t \int_0^{ \frac{\pi(\theta) f(x \mid \theta)}{c g(\theta) } }g(\theta)1 \text{d}u \text{d}\theta }{ \int_{-\infty}^{\infty} \int_0^{ \frac{\pi(\theta) f(x \mid \theta)}{c g(\theta) } }g(\theta)1 \text{d}u \text{d}\theta } \\ &= \frac{\int_{-\infty}^t g(\theta)\frac{\pi(\theta) f(x \mid \theta)}{c g(\theta) } \text{d}\theta }{ \int_{-\infty}^{\infty} g(\theta) \frac{\pi(\theta) f(x \mid \theta)}{c g(\theta) } \text{d}\theta } \\ &= \frac{\int_{-\infty}^t \pi(\theta) f(x \mid \theta) \text{d}\theta }{ \int_{-\infty}^{\infty} \pi(\theta) f(x \mid \theta)\text{d}\theta } \\ &= \pi(\theta \le t \mid x). \end{align*} So these accepted draws $\theta^1,\theta^2, \ldots, \theta^n$ have as their distribution the posterior of interest, and we never had to figure out the normalizing constant. We did have to figure out an upper bound, though. In my opinion, that's the downside of this method compared with other methods where that step isn't necessary.

Answer 2

The answer by Taylor is excellent, and already correctly gives the posterior kernel, which gives an intractable integral. If your goal is to find the posterior density (as opposed to sampling from the posterior distribution) then you are effectively just looking for the constant-of-integration:

$$H(x) \equiv \frac{1}{\pi} \int \limits_{-\infty}^\infty \frac{\exp(-\tfrac{1}{2} (\theta - x)^2)}{(1+\theta^2)} d\theta.$$

In this case acceptance-sampling is an inefficient method of numerically estimating the integral, since it discards a lot of generated values. A more efficient method is to calculate numerically using a large number $m$ of quantiles of the standard normal distribution:

$$\hat{H}_m(x) \equiv \sqrt{\frac{2}{\pi}} \sum_{k=1}^m \frac{1}{1+(x+z_{(k)})^2} \quad \quad \quad z_{(k)} = \Phi^{-1}(\tfrac{2k-1}{2m}).$$

It is simple to establish that $\hat{H}_m \rightarrow H$ as $m \rightarrow \infty$, which gives you the approximate posterior:

$$\pi(\theta|x) \approx \frac{1}{\pi} \cdot \frac{\exp(-\tfrac{1}{2} (\theta - x)^2)}{(1+\theta^2) \hat{H}_m(x)} \quad \quad \text{for all } \theta \in \mathbb{R}.$$

---

Simulating the required constant: This estimate of the definite integral can be implemented in R with a vectorised function of x using the code below. The integral is symmetric about zero so we plot it (on a logarithmic scale) for positive values of $x$.

#Create function to estimate definite integral
H <- function(x, m) { P     <- ((1:m) - 1/2)/m;
                      Z     <- qnorm(P, 0, 1);
                      HHAT  <- rep(0, length(x));
                      for (i in 1:length(x)) {
                        DENOM   <- 1 + (x[i]+Z)^2;
                        HHAT[i] <- sqrt(2/pi)*sum(1/DENOM); }
                      HHAT; }

#Plot this function for values of x
library(ggplot2);
library(scales);

theme_update(plot.title    = element_text(size = 15, hjust = 0.5),
         plot.subtitle = element_text(size = 10, hjust = 0.5),
         axis.title.x  = element_text(size = 10, hjust = 0.5),
         axis.title.y  = element_text(size = 10, vjust = 0.5),
         plot.margin   = unit(c(1, 1, 1, 1), "cm"));

m <- 10^6;
DATA <- data.frame(xx = (0:100),
                   HH = H((0:100), m));
ggplot(data = DATA, aes(x = xx, y = HH)) +
  geom_line(size = 1.2, colour = 'red') +
  scale_y_log10(breaks = trans_breaks('log10', function(x) 10^x),
                labels = trans_format('log10', math_format(10^.x))) +
  ggtitle('Plot of estimated value of integral') + 
  labs(subtitle = paste0('(Quantile sampling with ', 
                         format(m, scientific = FALSE, big.mark = ','), 
                         ' values)')) + 
  xlab('x') + ylab('hat(H)');