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
def Model(W_bM = W_bM, n = n):
    #Load parameters
    Timelist = {'StartTime':186,'StopTime':217} 
    Parameterlist = {'n' : n,'W_b@Multiplier' : W_bM} 
    #Run Model instance  
    ModelInstance.model_run(Parameterlist, Timelist)

    #Read output section
    S031_Q_Mod = ReadOutput('IUWS.River_quantity.Dynamic.Simul.1.out.txt','.S031.Q').as_matrix()[:-1]
    S031_d_Mod = ReadOutput('IUWS.River_quantity.Dynamic.Simul.1.out.txt','.S031.d').as_matrix()[:-1]
    S008_d_Mod = ReadOutput('IUWS.River_quantity.Dynamic.Simul.1.out.txt','.S008.d').as_matrix()[:-1]
    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')

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?

Edited: Aditional Information

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))

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
def Model(a = a, b = b, c = c, x=x):
    global EvalNum
    EvalNum += 1
    return a + b*x - c * x**2

# Likelihood function
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])
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.

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

This is as expected. Remember, in Metropolis-Hastings sampling, variables are proposed one-at-a-time, so the expected number of evaluations is always multiples of n. Your model is not particularly large, so I assume you have a large sample size?

You might look at PyMC3, which has gradient based samplers that propose values in k-space, and have much higher acceptance rates.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.