In [ ]:
import sys
from mcupy.graph import *
from mcupy.nodes import *
from mcupy.utils import *
from mcupy.core import ensemble_type
try:
import pydot
except(ImportError):
import pydot_ng as pydot
Create a graph object, which is used to hold nodes.
In [ ]:
g=Graph()
Create some nodes
In [ ]:
A=FixedUniformNode(0.001,1-1e-5).withTag("A")
B=FixedUniformNode(0.001,1-1e-5).withTag("B")
mu=FixedUniformNode(.001,100-1e-5).withTag("mu")
sigma=FixedUniformNode(.001,100-1e-5).withTag("sigma")
And some more nodes
In [ ]:
for l in open('eff.txt'):
e1,nrec1,ninj1=l.split()
e1=float(e1)
nrec1=float(nrec1)
ninj1=float(ninj1)
E=C_(e1).inGroup("E")
ninj=C_(ninj1).inGroup("ninj")
eff=((B-A)*PhiNode((E-mu)/sigma)+A).inGroup("eff")
nrec=BinNode(eff,ninj).withObservedValue(nrec1).inGroup("nrec")
g.addNode(nrec)
Then let us check the topology of graph
In [ ]:
display_graph(g)
It's correct.
Then we'd like to perform several sampling and record the values.
Before sampling, we need to decide which variables we need to monitor.
In [ ]:
mA=g.getMonitor(A)
mB=g.getMonitor(B)
mSigma=g.getMonitor(sigma)
mMu=g.getMonitor(mu)
We need a variable to hold the results
In [ ]:
result=[]
Then we perform the sampling for 1000 time for burning
In [ ]:
for i in log_progress(range(1000)):
g.sample()
Then we perform 30000 sampling and record the results
In [ ]:
for i in log_progress(range(30000)):
g.sample()
result.append([mA.get(),mB.get(),mMu.get(),mSigma.get()])
Then we plot the results.
In [ ]:
%matplotlib inline
import seaborn
import scipy
result=scipy.array(result)
seaborn.jointplot(result[:,0],result[:,1],kind='hex')
In [ ]:
seaborn.jointplot(result[:,0],result[:,2],kind='hex')
In [ ]:
seaborn.jointplot(result[:,0],result[:,3],kind='hex')
In [ ]:
seaborn.jointplot(result[:,1],result[:,2],kind='hex')
In [ ]:
seaborn.jointplot(result[:,1],result[:,3],kind='hex')
In [ ]:
seaborn.jointplot(result[:,2],result[:,3],kind='hex')
Now, mcupy has also implemented the ensemble-based sampler for graph, which is much more faster. To use it, first declare a data structure to store the ensemble as:
In [ ]:
em=ensemble_type()
then, iteratively call the graph.ensemble_sample(em)
In [ ]:
result=[]
for i in log_progress(range(100000)):
g.ensemble_sample(em)
result.append([em[0][0],em[0][1],em[0][2],em[0][3]])
result=scipy.array(result)
In [ ]:
seaborn.jointplot(result[:,1],result[:,0],kind='hex')
In [ ]:
seaborn.jointplot(result[:,1],result[:,2],kind='hex')
In [ ]: