2

I need to generate samples from a pdf given by $\frac{f_Z(z)\cdot 1_{Z \in B}}{P(Z \in B)}$ where $Z \in \mathbb{R}^d$ is a normal random vector with independent components. $Z \in B$ is a set that is hard to characterize and may be defined implicitly. However, it can be assumed that it is a closed set with no holes.

After evaluating my options, the most convenient approach would be Metropolis MCMC. I am aware that the result involves correlated instead of independent samples. My objective is not to compute any mean, variance, or any other statistics. I just need to generate samples of $Z$ with pdf above. However, if the chain mixes well, the samples I get should be good enough to be roughly independent and construct histograms of $Z$, etc.

The MCMC algorithm according to Wikipedia says that if a proposed state is rejected, then the current state is set as the old state. Does this also mean that I have to "save" the old state a second time (the first time when it got accepted and the second time when the proposed state is rejected) in the list of samples that I plan to use after sampling? Or can I simply discard these "duplicates" from the list of samples that I have?

user1237300
  • 347
  • 2
  • 9
  • See here: https://stats.stackexchange.com/questions/123113/what-to-do-when-rejecting-a-proposed-point-in-mcmc?rq=1 – jcz Aug 08 '18 at 00:58

2 Answers2

2

Yes, you do have to save the old state a second time, or however many times you reject a proposal.

As you say, Metropolis-Hastings is a Markov chain Monte Carlo (MCMC) method. Given some distribution of interest $\pi$ (in your case $\pi(z)\propto f(z)\mathbf{1}_{z\in B}$), you simulate a Markov chain that is designed to have $\pi$ as its invariant distribution. That is, $\pi$ is the limiting distribution of the chain. Eventually, if you run the chain long enough, the draws in the chain are statistically indistinguishable from draws sampled from $\pi$.

A feature of the MH algorithm is that sometimes you reject proposals, so the chain will contain duplicates. And the possibility of duplicates is part of the algorithm. A sampling scheme where you generate a MH chain and then discard any duplicates will not necessarily have the desired statistical properties. $\pi$ may not be the invariant distribution any longer.

jcz
  • 1,345
  • 6
  • 14
  • thank you @bamts . I changed my code and it seemed that I got better results. But is it odd to use the MCMC samples to construct a histogram? I.e. is it weird to use correlated samples to construct a histogram? Do you have any references you can direct me to which discuss this issue? thanks! – user1237300 Aug 08 '18 at 03:05
  • It is not weird. People do it all the time, and it seems to work in practice. That said, I'm not sure what the theoretical justification for the practice *actually is*. Let's find out: https://stats.stackexchange.com/questions/361237/in-exactly-what-sense-do-mcmc-draws-approximate-the-target – jcz Aug 08 '18 at 06:39
  • I'm looking at the link now. Thanks so much for pointing it out! – user1237300 Aug 08 '18 at 20:08
0

Usually when you are doing MCMC sampling you generate a list of sampled values from the algorithm. If the proposed state is rejected then you set $X_t = X_{t-1}$, where $X_{t-1}$ is the current state of the MCMC algorithm. Alternately, you could have 2 (associated) lists first one w/the values and the second one with the counts of the number of times you stay at the state. However, you usually want to be accepting new states with some "reasonably" high probability (search for articles by Rosenthal for more on this) so usually better off just duplicating the same value multiple times in your list.

Lucas Roberts
  • 3,819
  • 16
  • 45
  • so duplicates are "not bad"? Even if my goal is to construct a histogram of the samples instead of computing means, variances, or other statistics? – user1237300 Aug 07 '18 at 23:44
  • 1
    The duplicates are not "bad" they are samples from the stationary distribution. If you do not the samples are not from the distribution you want. – Lucas Roberts Aug 08 '18 at 14:05