# MCMC algorithm going wrong [closed]

Given this integral

$$$$\int_0^\infty \chi_{[1,2]}(x)\Gamma(C,x)\left|\cos(R x)\right| \, dx$$$$

where $$\chi_{[1,2]}=\begin{cases}1, & x \in [1,2] \\ 0, & x\not \in [1,2]\end{cases}$$

I need to estimate its value using MCMC with Metropolis algorithm and acceptance.

I wrote the sampler as it follows:

def sampler(N):
step=1
old=1.5
i=0
points=[]
while i<N:
points.append(old)
ant=func(old)
inov=norm.rvs(0,step)
new=old+inov
prox=func(new)
alphah=min(1,prox/ant)
if alphah>np.random.uniform():
old=new
i+=1
points.append(old)
return points


sampler(N) generates $$N$$ points by using the Metropolis Algorithm with N(0,step) as kernel and starting in 1.5.

The function func(x) evaluates $$\chi_{[1,2]}(x)\Gamma(C,x)\left|\cos(R x) \right|$$.

After generating $$N$$ points $$x_0,x_1,\ldots,x_N$$ the estimative is $$$$\frac{1}{N}\sum_{i=1}^N \chi_{[1,2]}(x_i)\Gamma(C,x_i) \left|\cos(R x_i)\right|$$$$

The problem is that doing like that I find approximately 0.28 while I need to find 0.26

My algorithm is wrong? Or there is something that I'm not seeing?

EDIT: Implementation of func and the Monte Carlo algorithm mc

def func(x):
c=1.398031095
d=1.9298378
if x<=2 and x>=1:
return gamma.pdf(c,x)*abs(np.cos(d*x))
return 0

def mc(N):
points=sampler(N)
i=0
g=0
while i <N:
g+=func(points[i])
i+=1
return g/N



## closed as off-topic by mkt, mdewey, Glen_bMay 12 at 12:28

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – mkt, mdewey, Glen_b
If this question can be reworded to fit the rules in the help center, please edit the question.

• It would help if you provided the implemention of func, too. Right now I can't find anything wrong with this. Perhaps you're not running it for long enough – Taylor May 10 at 19:54
• I will post the implemention of func and monte carlo – Ettore Moura May 10 at 19:55
• Add the self study tag. – Michael Chernick May 10 at 20:47
• I found the error, I wrote on func gamma.pdf(c,x). But the correct function is gamma.pdf(x,c) – Ettore Moura May 11 at 3:19

Your main confusion is that you appear to not be correctly identifying the difference between the functional and the target distribution. I would say your target distribution appears to be a simple uniform distribution, while the functional is the fancy guy. Confusingly, your functional involves a gamma density.

MCMC gives you an estimate of an expectation. Your expectation is this:

$$\int_0^\infty \chi_{[1,2]}(x)\Gamma(C,x)\left|\cos(R x)\right| \, dx = E_X[h(x)]$$ where $$h(x) = \frac{1}{\Gamma(c)}x^{c-1}e^{-x} \left|\cos(R x)\right|$$ is the function that your func evaluates, and $$X \sim \text{Uniform}(1,2]$$. $$h$$ is the functional, and $$\text{Uniform}(1,2]$$ is the target distribution.

Metropolis-Hastings wouldn't be my first choice for this problem, but you seem to be using random walk Metropolis-Hastings. I think this because you are proposing samples like this: $$X^*\mid X_{t-1} \sim \text{Normal}(X_{t-1}, 1).$$ You call $$X^*$$ new, and $$1$$ step (this name confuses me, but that's not that important).

You can write your target distribution as $$\pi(x) = \chi_{[1,2]}(x)$$, and this makes your acceptance ratio $$a(x_{t-1},x^*) = \min(1, \pi(x^*)/\pi(x_{t-1})).$$ Note that you are incorrectly calculating this with min(1,prox/ant). You are using the $$h$$ function here, when you shouldn't be (if you're targeting this uniform distribution). Instead, use this code (yours with minimal changes)

from scipy.stats import gamma,norm
import numpy as np
import matplotlib.pyplot as plt

def sampler(N):
step=1
old=1.5
i=0
points=[]
while i<N:
points.append(old)
ant= 1 if 1 < old <= 2 else 0
inov=norm.rvs(0,step)
new=old+inov
prox = 1 if 1 < new <= 2 else 0
alphah=min(1,prox/ant)

if alphah>np.random.uniform():
old=new
i+=1
points.append(old)
return points


Your mc function appears to be correct. It is calculating $$N^{-1} \sum_{i=1}^{N} h(x_t)$$ where $$x_t$$ are all your samples (which you are calling points).

Finally, yes, your answer is valid as well. You are likely reversing the argument order in the gamma.pdf function.

• How do you know that the target distribution is the Uniform(1,2) ? – Ettore Moura May 12 at 14:24
• It doesn’t have to be. It could also be a gamma distribution. In that case the functional would only be the trigonometric part and the indicator function multiplied together. If you decided to go this route, you would have to change the code in different ways. – Taylor May 12 at 14:28
• But taking the gamma as target distribution the integral must be the same right? – Ettore Moura May 12 at 14:34
• As long as you change the functional accordingly. – Taylor May 12 at 14:54
• Known that $g \propto \chi_{[1,2]}h(x)$ the integral of g is calculated the way I calculated in my post right? And even knowing this proportionality the target function can be the uniform? – Ettore Moura May 12 at 15:02

I believe the following should work for you. My question is why do you have new=old+inov in your sampler? I replaced that with new=inov. Also, you should expect to see some deviations. I didn't really dig into this problem much but I do not think there is much of a difference between 0.26 and 0.28. Is it given that 0.26 is the target value (mean) for this distribution?

import numpy as np
from scipy.stats import gamma, norm

def func(x):
c=1.398031095
d=1.9298378
if x<=2 and x>=1:
return gamma.pdf(c,x)*abs(np.cos(d*x))
return 0

def mc(N):
points=sampler(N)
i=0
g=0
while i <N:
g+=func(points[i])
i+=1
return g/N

def sampler(N):
step=1
old=1.5
i=0
points=[]
while i<N:
points.append(old)
ant=func(old)
inov=norm.rvs(0,step)
new= inov #old+inov
prox=func(new)
alphah=min(1,prox/ant)
if alphah>np.random.uniform():
old=new
i+=1
points.append(old)
return points

out = sampler(10000)
out = np.array(out)

results = np.array(list(map(func, out)))
print(np.mean(results))

• Yes it is given that we should aim for 0.26 – Ettore Moura May 10 at 22:03
• Doing what you did doesn't chance the markov process? – Ettore Moura May 11 at 0:28