2
$\begingroup$

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?

$\endgroup$
3
  • 1
    $\begingroup$ Your pdf is not a pdf. Its integral from $0$ to $1$ is only about $0.6559173858$ $\endgroup$ Nov 28, 2019 at 13:13
  • $\begingroup$ @RobertIsrael Yeah. But that is just scaling. The main problem is the functionality. $\endgroup$ Nov 28, 2019 at 13:26
  • 1
    $\begingroup$ 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. $\endgroup$ Nov 28, 2019 at 13:28

1 Answer 1

4
$\begingroup$

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.

$\endgroup$
2
  • 2
    $\begingroup$ ... one of the downsides in a language where whitespace has semantic significance ... $\endgroup$ Nov 28, 2019 at 13:43
  • $\begingroup$ @joriki Thanks. I applied your comments. Now that is so much better. How can I improve it, especially at the area close to zero? $\endgroup$ Nov 28, 2019 at 16:57

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .