f(y) = c e^(-y^4)(1+|y|)^3

where

*c*is a normalizing constant and

*y*is defined on the whole real line.

Use of the probability integral transform (section 1.10.8) is not feasible in this setting, given the complexity of inverting the cumulative density function.

The Metropolis--Hastings algorithm is a Markov Chain Monte Carlo (MCMC) method for obtaining samples from a probability distribution. The intuition behind this algorithm is that it chooses proposal probabilities so that after the process has converged we are generating draws from the desired distribution. A further discussion can be found on page 610 of Section 11.3 of the Evans and Rosenthal text, or on page 25 of Gelman et al.

We find the acceptance probability

*a(x, y)*in terms of two densities,

*f(y)*and

*q(x,y)*(a proposal density, in our example, normal with specified mean and unit variance) so that

a(x,y) = min(1, [c*f(y)*q(y,x)]/[c*f(x)*q(x,y)]

= min(1, [e^(-y^4+x^4)(1+|y|)^3]/[(1+|x|)^3]

where we omit some algebra.

Begin by picking an arbitrary value for

*X_1*. The Metropolis--Hastings algorithm then proceeds by computing the value

*X_{n+1}*as follows:

1. Generateyfrom a Normal(X_n, 1).

2. Compute a(x, y) as above.

3. With probability a(x, y), letX_{n+1} = y

(i.e., use proposal value).

Otherwise, with probability 1-a(x, y), letX_{n+1} = X_n

(i.e., keep previous value).

The code allows for a burn-in period, which we set at 50,000 iterations, and a desired sample from the target distribution, which we make 5,000. To reduce auto-correlation, we take only every twentieth variate. In a later entry, we'll compare the resulting variates with the true distribution.

**SAS**

The SAS code is fairly straightforward.

data mh;

burnin = 50000; numvals = 5000; thin = 20;

x = normal(0);

do i = 1 to (burnin + (numvals * thin));

y = normal(0) + x;

switchprob = min(1, exp(-y**4 + x**4) *

(1 + abs(y))**3 * (1 + abs(x))**(-3));

if uniform(0) lt switchprob then x = y;

* if we don't change x, the previous value is kept--

no code needed;

if (i gt burnin) and mod(i-burnin, thin) = 0 then output;

* only start saving if we're past the burn-n period,

then thin out;

end;

run;

**R**

In R we first define a function to compute a(x,y):

alphafun = function(x, y) {

return(exp(-y^4+x^4)*(1+abs(y))^3*

(1+abs(x))^-3)

}

Then we proceed as in the SAS example.

numvals = 5000; burnin = 50000; thin = 20

res = numeric(numvals)

i = 1

xn = rnorm(1) # arbitrary value to start

for (i in 1:(burnin + numvals * thin)) {

propy = rnorm(1, xn, 1)

alpha = min(1, alphafun(xn, propy))

xn = sample(c(propy, xn), 1, prob=c(alpha,1-alpha))

if ((i > burnin) & ((i-burnin) %% thin == 0))

res[(i-burnin)/thin] = xn

}

The resulting draws from the distribution are available in the

`res`vector.

## No comments:

Post a Comment