Precision of Metropolis Monte Carlo sampling of a PDF

Question

I want to sample from the probability distribution function (PDF)

$$p(t)=3e^{-t}\big[\frac{1}{2}+\frac{1}{2}e^{-8t}\cos({40t})\big]$$

for $0<t<1$. The function looks like this:

PDF

The PDF can not be analytically inverted, so I tried to use Metropolis Monte Carlo sampling with the following python code:

pdf = lambda t: 3*np.exp(-t)*(0.5+0.5*np.exp(-8*t)*np.cos(40*t))
lw  = 0
up  = 1
sig = 0.1
ini_t   = 0.01
new_t   = ini_t
samples = []
for i in range(100000):
    cand_t = scipy.stats.truncnorm.rvs((lw-new_t)/sig, (up-new_t)/sig, loc=new_t, scale=sig)
    if np.random.uniform(0, 1) < pdf(cand_t)/pdf(new_t):
        new_t = cand_t
    samples.append(new_t)
burn_in = 10000
samples = np.array(samples[burn_in:])

The result looks like the green plot of the following image:

Simulated PDF

As it is obvious, the performance of sampling is rather poor. How can I improve it?

Answer 1

Your program looks fine except that there's too much indentation before samples.append(new_t). That line shouldn't be conditional on the if statement; a sample is generated regardless of whether you accept or reject.

By the way, you're using ”performance“ in an unorthodox way. In this context it's usually used to mean speed; judging from the image, your problem rather seems to be correctness.