1
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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.

$\endgroup$
27
  • $\begingroup$ Can you please show the proportional in detail? $\endgroup$ – shuvam agrawal Dec 4 '18 at 14:23
  • $\begingroup$ What does this mean "By the way, your prior is a C(0,1) prior, which has support −∞<θ<∞." $\endgroup$ – shuvam agrawal Dec 4 '18 at 14:24
  • $\begingroup$ @shuvamagrawal regarding your first comment, multiply the two densities together, and drop any factor that doesn't depend on $\theta$. Regarding your second, I just mentioned that because your phrasing in the question sounded to me like you thought these were two different problems, when they're not. Are you interested in using accept-reject sampling? I just noticed the tag $\endgroup$ – Taylor Dec 4 '18 at 14:39
  • $\begingroup$ Yes I am interested in using accept-reject sampling $\endgroup$ – shuvam agrawal Dec 4 '18 at 15:00
  • 1
    $\begingroup$ @shuvamagrawal please post a separate question for other questions $\endgroup$ – Taylor Dec 9 '18 at 16:09
2
$\begingroup$

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)');

enter image description here

$\endgroup$
3
  • $\begingroup$ Thanks @ Ben. one thing where how to code the rejection sampler in R, and use it to generate 10, 000 samples from the posterior distribution. with using r function rcauchy to simulate from $\pi$($\theta$) $\endgroup$ – shuvam agrawal Dec 6 '18 at 13:43
  • $\begingroup$ Can you please help with that $\endgroup$ – shuvam agrawal Dec 6 '18 at 13:45
  • $\begingroup$ @shuvam agrawal: My answer does not involve rejection sampling, so I have not coded this. If you would like to ask a question about coding rejection sampling, please ask it as a new question, rather than in comments. $\endgroup$ – Ben Dec 6 '18 at 23:55

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.