1
$\begingroup$

I have been searching for some answers about the need to use monte carlo methods to sampling distributions.
When comparing with the standard random sampling - in the computational sense, we have a table and take n random samples - I can't see any major advantage easily.
Thinking a bit more, all that I can think of is related to the central limit theorem (our samples will be normally distributed and will not maintain the original distribution of the data). Is my assumption correct?
Can you provide more details into this? interesting sources for me to read?

Thanks

$\endgroup$
2
$\begingroup$

We may be interested in means and standard deviations of some function of $\theta$,. When comparing two or more populations we may be interested in the posterior distributionof $|\theta_{1} −\theta_{2}|, \theta_{1}/\theta_{1}$, or max$({\theta_{1}, . . . , \theta_{m}})$, all of which are functions of more than one parameter. Obtaining exact values for these posterior quantities can be difficult or impossible, but if we can generate random sample values of the parameters from their posterior distributions, then all of these posterior quantities of interest can be approximated to an arbitrary degree of precision using the Monte Carlo method.

Let $\theta$ be a parameter of interest and let $y_{1}, . . . ,y_{n}$ be the numerical values of a sample from a distribution $p(y_{1}, . . . ,y_{n})$. Suppose we could sample some number $S$ of independent, random $\theta$-values from the posterior distribution $p(\theta|y_{1}, . . . ,y_{n})$: $$\theta^{(1)}, . . . , \theta^{(s)} \sim p(\theta|y_{1}, . . . ,y_{n})$$.

Then the empirical distribution of the samples $\theta^{(1)}, . . . , \theta^{(s)}$ would approximate $p(\theta|y_{1}, . . . ,y_{n})$, with the approximation improving with increasing $S$.The empirical distribution of $\theta^{(1)}, . . . , \theta^{(s)} $is known as a Monte Carlo approximation to $p(\theta|y_{1}, . . . ,y_{n})$.

Additionally, let $g(\theta)$ be (just about) any function of $\theta$. The law of large numbers says that if $\theta^{(1)}, . . . , \theta^{(s)} $ are i.i.d. samples from $p(\theta|y_{1}, . . . ,y_{n})$, then as $S\rightarrow \infty$ $$\frac{1}{S}\sum_{s=11}^{S}g(\theta) ^{(s)}\rightarrow E(g(\theta)|y_{1}, . . . ,y_{n}) $$

Example :

Suppose we are interested in the posterior distribution of some computable function $g(\theta)$ of $\theta$. In the binomial model, for example, we are sometimes interested in the log odds :

$$log odds (\theta) = log \frac{\theta}{1-\theta}=\gamma$$

this quantity of interest can be computed using a Monte Carlo approach:

sample $\theta^{(1)} \sim p(\theta|y_{1}, . . . ,y_{n})$, compute $\gamma^{(1)} = g(\theta^{(1)})$

sample $\theta^{(2)} \sim p(\theta|y_{1}, . . . ,y_{n})$, compute $\gamma^{(2)} = g(\theta^{(2)})$ .

.

.

sample $\theta^{(s)} \sim p(\theta|y_{1}, . . . ,y_{n})$, compute $\gamma^{(s)} = g(\theta^{(s)})$

Then the sequence $\{\gamma^{(s)}, . . . , \gamma^{(s)}\}$ constitutes $S$ independent samples from $p(\gamma|y_{1}, . . . ,y_{n})$.

Let $\theta$ be a population proportion . Using a binomial sampling model and a uniform prior distribution the posterior distribution of $\theta$ will be beta distributed say with beta(45, 420). Using the Monte Carlo algorithm described above, we can obtain samples of the log-odds from the posterior distribution of $\theta$ . In R, the Monte Carlo algorithm involves only a few commands:

set.seed(1)
a<-1 ; b<-1 # Parameter for the prior 

# Prior for theta 
theta.prior.sim<-rbeta(10000,a,b) # Beta(1,1) is also uniform

# Prior of the function of interest log odds (gamma)
gamma.prior.sim<- log( theta.prior.sim/(1-theta.prior.sim) )

# Posterior for theta 
theta.post.sim<-rbeta(10000,45,420)

# Posterior for log odds (gamma) 
gamma.post.sim<- log( theta.post.sim/(1-theta.post.sim) )

# Plot the prior and the posterior for gammma 
par(mar=c(3,3,1,1),mgp=c(1.75,.75,0))
par(mfrow=c(2,3))
par(cex=.8)
par(mfrow=c(1,2),mar=c(3,3,1,1), mgp=c(1.75,.75,.0))
plot(density(gamma.prior.sim,adj=2),xlim=c(-5,5),main="", 
 xlab=expression(gamma),
 ylab=expression(italic(p(gamma))),col="gray")
plot(density(gamma.post.sim,adj=2),xlim=c(-5,5),main="",xlab=expression(gamma),
 ylab=expression(paste(italic("p("),gamma,"|",y[1],"...",y[n],")",
                       sep="")) )
lines(density(gamma.prior.sim,adj=2),col="gray")

For further reading I recommend A First Course in Bayesian Statistical Methods and Doing Bayesian Data Analysis: A Tutorial Introduction with R and BUGS

$\endgroup$
2
$\begingroup$

Whenever you can do simple random sampling, you would. In fact, when there's just a simple distribution involved, you can often just calculate what you are interested in without doing random sampling. You use Markov chain Monte Carlo (MCMC) techniques when that is not possible.

One example would be a Bayesian hierarchical model. Very rarely do you know the distribution of the posterior. It will not usually be some known standard distribution (so even if you could derive the distribution, a sampler likely does not exist) and because you lack the normalizing constant of the distribution, but MCMC sampling is possible.

$\endgroup$

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