In [3]:
!date
import matplotlib.pyplot as plt, numpy as np, seaborn as sns, pandas as pd
import pymc
%matplotlib inline
sns.set_context("poster")
sns.set_style('whitegrid')

Wed Dec 31 10:01:45 PST 2014


In that example I would say that there were three fitting parameters, amp, size and ps. Let us call the number of parameters being examined N. Now let us call the number of samples to be drawn, P (1e4 in this case). I have observed that the @deterministic function gauss is called roughly N x P times.

I would like to know the reason why it is N x P?

Chris answered this, but I would add that the number of times it needs to be called is dependent on the step method being used. The default in PyMC2 is to use a Metropois step to each Stochastic node, but you can change this, for example with MDL.use_step_method(pymc.AdaptiveMetropolis, [MDL.amp, MDL.size, MDL.ps]).

Is there an attribute inside MDL to find out how many times gauss has been called?

As far as I know, there is not, but you can roll your own with a global var. Here is an example

In [7]:
x = np.ones(10)
f = np.ones(10)
f_error = np.ones(10)

# define the model/function to be fitted.
def model(x, f):
amp = pymc.Uniform('amp', 0.05, 0.4, value= 0.15)
size = pymc.Uniform('size', 0.5, 2.5, value= 1.0)
ps = pymc.Normal('ps', 0.13, 40, value=0.15)

@pymc.deterministic(plot=False)
def gauss(x=x, amp=amp, size=size, ps=ps):
global calls_to_gauss
calls_to_gauss += 1

e = -1*(np.pi**2*size*x/(3600.*180.))**2/(4.*np.log(2.))
return amp*np.exp(e)+ps
y = pymc.Normal('y', mu=gauss, tau=1.0/f_error**2, value=f, observed=True)
return locals()

calls_to_gauss = 0
MDL = pymc.MCMC(model(x,f))
MDL.sample(1e4)
print '\ncalls to gauss:', calls_to_gauss

 [-----------------100%-----------------] 10000 of 10000 complete in 2.3 sec
calls to gauss: 21400

In [9]:
calls_to_gauss = 0
MDL = pymc.MCMC(model(x,f))

# change step method, change number of calls to gauss (and acceptance rate...)

MDL.sample(1e4)
print '\ncalls to gauss:', calls_to_gauss

 [-----------------100%-----------------] 10000 of 10000 complete in 1.6 sec
calls to gauss: 5265

In [10]:
calls_to_gauss = 0
MDL = pymc.MCMC(model(x,f))

# change step method, change number of calls to gauss (and acceptance rate...)

MDL.sample(1e4)
print '\ncalls to gauss:', calls_to_gauss

 [-----------------100%-----------------] 10000 of 10000 complete in 4.8 sec
calls to gauss: 23997

In [ ]: