Original link:http://tecdat.cn/?p=23019 

The Monte Carlo method uses random numbers to generate samples from the probability distribution P(x) and from that distribution to estimate the expected value, which is often too complex to be evaluated by precise methods. In Bayesian inference, P (x) is usually defined as a joint posterior distribution over a set of random variables. However, obtaining independent samples from this distribution is not easy, depending on the dimension of the sample space. Therefore, we need to resort to a more complex Monte Carlo method to help simplify the problem; For example, materiality sampling, rejection sampling, Gibbs sampling, and Metropolis Hastings sampling. These methods usually involve sampling from the recommended density Q(x) instead of P(x).

In importance sampling, we generate samples from Q(x) and introduce weights to account for sampling from incorrect distributions. We then adjusted the importance of each point in the estimator that we needed to evaluate. In rejection sampling, we extract a point from the proposed distribution Q(x) and calculate the ratio P(x)/Q(x). Then we take a random number U from the distribution of U(0,1); ifWe accept this point x, or we reject it and go back to Q(x) and extract another point. Gibbs sampling is a method of sampling from a distribution of at least two dimensions. Here, the proposed distribution Q(x) is defined as a conditional distribution of the joint distribution P(x). We simulate the posterior sample of P(x) by iteratively sampling from the posterior conditions, while setting the other variables at their current values.

Although both materiality sampling and rejection sampling require Q(x) to be similar to P(x) (it is difficult to create such a density in a higher dimensional problem), Gibbs sampling is difficult to apply when a conditional posteriori has no known form. This assumption can be relaxed in the more general Metropolis-Hastings algorithm, where the candidate samples are probabilistically accepted or rejected. This algorithm can accommodate both symmetric and asymmetric proposed distributions. The algorithm can be described as follows

Initialize the



extracting

To calculate

fromExtract the

if 设

Otherwise, set

The end of the

Gibbs sampling is a special case of Metropolis Hastings. It involves a proposal that is always accepted (there is always a Metropolis-Hastings ratio of 1).

We apply the Metropolis Hastings algorithm to estimate the variance components of the regression coefficients in the standard G-BLUP model.

For the G-Blup model.

Among them.A vector representing a phenotype and a matrix of genotypes.Is the vector of the labeling effect,Is the vector of model residuals, the residuals are normally distributed, the mean is 0, and the variance isand.

Taking the rest of the parameters into account,The conditional posterior density of is:

This is an inverse chi-square distribution.

Let’s say we need to makeA priori is as uninformative as possible. One option is Settings.And use rejection sampling to estimate; However, setting S0=0 May cause the algorithm to get stuck at 0. Therefore, we need a priori that can replace the inverse chi-square distribution, and can be very flexible. For this purpose, we recommend using a beta distribution. Since the posterior obtained is not an appropriate distribution, the Metropolis Hastings algorithm will be obtainedA good choice for a posterior sample.

Here we’re going to takeAs our proposal distributes Q. So.

So our target distribution isAnd the normal likelihood ofThe product of beta priors. Since the beta distribution is between 0 and 1, we use variablesTo replace β a priori, where MAX is a guaranteed greater thanNumbers like this.

Where α1 and α2 are the shape parameters of β distribution, and their mean values are determined byIs given.

We calculated our acceptance rate according to the above algorithm steps, as shown below.

And then we take a random number u out of the uniform distribution, if, the sample points are acceptedOtherwise, we reject the point and retain the current value, iterating again until convergence.

Metropolis Hastings algorithm

MetropolisHastings=function(p, ...)

 chain\[1\]=x
  for (i in 1:nIter) {
      y\[i\] <-(SS+S0)/rchisq(df=DF,n=1,...)
  
    logp.old\[i\]=-(p/2)\*log(chai) - (SS/(2\*chain) + (shape1-1)*(log(chain\[i\]/(MAX)))+(shape2-1)*(log(1-(chain\[i\]/(MAX))

    logp.new\[i\]=-(p/2)\*log(y\[i\]) - (SS/(2\*y\[i\])) + (shape1-1)*(log(y\[i\]/(MAX)))+(shape2-1)*(log(1-(y\[i\]/(MAX))
    chain\[i+1\] = ifelse (runif(1)<AP\[i\] , y\[i\], chain\[i\],...)

Gibbs sampler

gibbs=function(p,...)
b = rnorm(p,0,sqrt(varb),...)
 for (i in 1:Iter) {
      chain\[i\] <-(S+S0)/rchisq(df=DF,n=1,...)

Drawing a graph

plot = function(out1,out2)
plot(density(chain1),xlim=xlim)
lines(density(chain2),xlim=xlim)
abline(v=varb,col="red",lwd=3)

Set the parameters

Run the Gibbs sampler

##################
out1=gibbs(p=sample.small,...)
out2=gibbs(p=sample.large,...)

Run Metropolis Hastings in different situations

Small sample size, prior

out.mh=mh(p=sample.small,nIter=nIter,varb=varb,shape1=shape.flat,shape2=shape.flat, MAX=MAX)

 

The sample size is small and the β value of shape 1 parameter is large

p=sample.small
nIter
varb
shape.skew\[1\]
shape.skew\[2\]
 MAX

plot(out.mh, out.gs_1)

 

The sample size is small and the β value of shape 1 parameter is large

MetropolisHastings(p)

makeplot(out.mh, out.gs_1)
## Summary of Chain for MH: ## Min.1st Q.Median Mean 3rd Q.max. ## 0.2097 0.2436 0.2524 0.2698 0.2978 0.4658

Small sample size, same shape parameter of β (large)

plot(out.mh, out1)

Large sample size, prior

plot(out.mh, out2)

Large sample size, shape 1 parameter β

plot(out.mh, out2)

Large sample size, β values for large shape 2 parameters

plot(out.mh, out_2)

With large sample size, the shape parameters of β are the same (large)

plot(out.mh, out2)

reference

  1. Gelman, Andrew, et al. Bayesian data analysis. Vol. 2. London: Chapman & Hall/CRC, 2014.

The most popular insight

1. Metroplis-in-Gibbs sampling and MCMC running with R language

2. The Bayesian model of MCMC sampling for Stan probabilistic programming in R language

3.R language implements Metropolis — Hastings algorithm and Gibbs sampling in MCMC

Bayesian analysis Markov chain Monte Carlo method (MCMC) sampling

5. Block Gibbs sampling Bayesian multivariate linear regression in R language

6. Simple Bayesian linear regression simulation analysis of Gibbs sampling in R language

7.R language uses RCPP to accelerate Metropolis-Hastings sampling to estimate the parameters of Bayesian logistic regression model

8.R uses Metropolis-Hasting sampling algorithm for logistic regression

9. HAR-RV model based on mixed data sampling (MIDAS) regression in R language predicts GDP growth