mcupy

This is a python wrapper for mcupy

Installation

This library relies on

1. : the core part for describing the graph structure of bayesian networks
   
2. : for performing cosmology calculation (this library is actually for astronomer
   
3. : parsing expression
   

Be sure to puth all above three libs and this lib in the same parent directory.

And it relies on

After above libs are installed, just execute go into pymcutil and run

./setup.sh build
./setup.sh install # as root if needed

##Example Check .
This is an example given in thie book section 8.2.

First let's import necessary packages

import sys
from mcupy.graph import *
from mcupy.nodes import *
from mcupy.utils import *
try:
    import pydot
except(ImportError):
    import pydot_ng as pydot

Create a graph object, which is used to hold nodes.

g=Graph()

Create some nodes

A=UniformNode(0.001,1-1e-5).withTag("A")
B=UniformNode(0.001,1-1e-5).withTag("B")
mu=UniformNode(.001,100-1e-5).withTag("mu")
sigma=UniformNode(.001,100-1e-5).withTag("sigma")

And some more nodes

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

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.

mA=g.getMonitor(A)
mB=g.getMonitor(B)
mSigma=g.getMonitor(sigma)
mMu=g.getMonitor(mu)

We need a variable to hold the results

result=[]

Then we perform the sampling for 1000 time for burning

for i in log_progress(range(1000)):
    g.sample()    

Then we perform 30000 sampling and record the results

for i in log_progress(range(30000)):
    g.sample()
    result.append([mA.get(),mB.get(),mMu.get(),mSigma.get()])

Then we plot the results.

%matplotlib inline
import seaborn
import scipy
result=scipy.array(result)
seaborn.jointplot(result[:,0],result[:,1],kind='hex')

seaborn.jointplot(result[:,0],result[:,2],kind='hex')

seaborn.jointplot(result[:,0],result[:,3],kind='hex')

seaborn.jointplot(result[:,1],result[:,2],kind='hex')

seaborn.jointplot(result[:,1],result[:,3],kind='hex')

seaborn.jointplot(result[:,2],result[:,3],kind='hex')