# Precision of Metropolis Monte Carlo sampling of a PDF

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. 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?

• Your pdf is not a pdf. Its integral from $0$ to $1$ is only about $0.6559173858$ Nov 28, 2019 at 13:13
• @RobertIsrael Yeah. But that is just scaling. The main problem is the functionality. Nov 28, 2019 at 13:26
• Recheck your plot of P. Wolfram Alpha shows that the integral of P(t) from 0 to 1 is about 0.96, and it's plot of the function is decently close to what your simulation results show. Nov 28, 2019 at 13:28

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.