# Pymc 2 MCMC Sample argument vs Model evaluations

I have been successfully using Pymc 2.0 for inference in a few fast models. However, I'm trying to set up my first real case inference. I'm using a relatively "heavy" likelihood (it takes ~2 min per evaluation).

I noticed that the number of samples specified in the MCMC sampler do not match with the number of actual model evaluations. I do not know if this is a known behavior of pymc. Or I might be doing something very wrong.

I will try to explain it further:

I'm calling an external model with a deterministic decorator:

# Model deterministic part
@pm.deterministic
def Model(W_bM = W_bM, n = n):
Timelist = {'StartTime':186,'StopTime':217}
Parameterlist = {'n' : n,'W_b@Multiplier' : W_bM}
#Run Model instance
ModelInstance.model_run(Parameterlist, Timelist)

Output = np.array([S031_Q_Mod, S031_d_Mod, S008_d_Mod])

return Output


Then I set up the sampler:

m = pm.MCMC([n,W_bM, sigma_q031, sigma_d031, sigma_d008, Log_Gaussian_multivariate, Model], db='pickle', dbname='Model_DommelRiver_DepthandFlowCalibration.pickle')
m.sample(1000)


In the external function called by the Model @Deterministic, I have a routine in which I run a model with a number of specified parameters and save a few output files (for later use). Thus I can see how many times this function was evaluated.

In a few runs I have performed, It saw that my function was called for 1500-1800 times instead of the 1000 specified. This is quite scary, since I can never be sure how long will the sampler take.

I'm using just a Metropolis-Hastings sampler, which to my understanding, evaluates the likelihood*prior, compares with the previous and then accepts-rejects the new proposal. This would produce 1000 samples for which there will be n repited elements (due to rejections).

But why is the sampler running more times than the specified? Is this perhaps a wrong connection of the external model? Or does pymc 2.0 behave like this?

I wanted to check whether this would happen with a model which doesn't call my external function. I implemented a simple quadratic function as @deterministic and I placed a counter inside (so you can see how many times was this evaluated).

import pymc as pm
import numpy as np

# Regression with error Synthetic data
x = np.linspace(0,10,100)
y = 30 + 35*x - 2.4 * x**2 + np.random.normal(0, 6, len(x))

# PRIOR PARAMETERS
a = pm.Uniform('a', 0, 50)
b = pm.Uniform('b', 10, 60)
c = pm.Uniform('c', 0, 5)

# Parameters error model
sigma = pm.Uniform('sigma', 0, 20)

# Model
EvalNum = 0
@pm.deterministic
def Model(a = a, b = b, c = c, x=x):
global EvalNum
EvalNum += 1
return a + b*x - c * x**2

# Likelihood function
@pm.stochastic(observed=True)
def Log_Gaussian_iid_univariate(value=y, model=Model, std_dev=sigma):
return np.sum(-np.log(std_dev) - (value-model)**2 / (2*(std_dev**2)))

m = pm.MCMC([a,b,c, sigma, Log_Gaussian_iid_univariate, Model])
m.sample(10000)
print '\n Number of Evaluations = {}'.format(EvalNum)


The counter shows that the model was evaluated 36980 times. When the m.sample argument was specifying 10000.

Is this mismatch representing rejected points? If rejection happens, I understand we should keep the previous parameter vector in the trace, which would count as a sample. I don't see why m.sample() should evaluate the function more than the specified trace length.

• You can at least shave off the "-0.5*np.log(2)" term in the likelihood function – Vladislavs Dovgalecs Feb 20 '17 at 17:46
• Yes, thanks. I can remove as well the np.log(np.pi) since is another constant. – Argantonio65 Feb 21 '17 at 9:37