!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
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
calls_to_gauss = 0
MDL = pymc.MCMC(model(x,f))
# change step method, change number of calls to gauss (and acceptance rate...)
MDL.use_step_method(pymc.AdaptiveMetropolis, [MDL.amp, MDL.size, MDL.ps])
MDL.sample(1e4)
print '\ncalls to gauss:', calls_to_gauss
[-----------------100%-----------------] 10000 of 10000 complete in 1.6 sec calls to gauss: 5265
calls_to_gauss = 0
MDL = pymc.MCMC(model(x,f))
# change step method, change number of calls to gauss (and acceptance rate...)
MDL.use_step_method(pymc.AdaptiveMetropolis, [MDL.amp, MDL.size, MDL.ps])
MDL.use_step_method(pymc.AdaptiveMetropolis, [MDL.amp, MDL.size])
MDL.use_step_method(pymc.AdaptiveMetropolis, [MDL.amp, MDL.ps])
MDL.use_step_method(pymc.AdaptiveMetropolis, [MDL.size, MDL.ps])
MDL.sample(1e4)
print '\ncalls to gauss:', calls_to_gauss
[-----------------100%-----------------] 10000 of 10000 complete in 4.8 sec calls to gauss: 23997