Experiments.py

import numpy as np
from pymc import *
import matplotlib.pyplot as plt
from matplotlib.pyplot import scatter,show
H,R=4,0
x=DiscreteUniform('X',lower=1,upper=4,size=4)


@potential
def L(H=4,R=0,x=x):
    p=0
    for i in range(1,4):
        if x[i]==x[0]:
            p+=H
        else:
            p+=R      
    return p 

fig,axes=plt.subplots(nrows=2,ncols=1)
              
N,t=5000,1
model=MCMC([x,L])
model.sample(N,thin=t)   
print '!'    
#u=x.trace()[:,0]+np.random.normal(0,0.1,N/t)
#v=x.trace()[:,1]+np.random.normal(0,0.1,N/t)
#scatter(u,v)
#show()
z=model.trace('X')[:]
s=0.0
print
for n in z:
    if n[0]==n[1]:
        s+=1
print s/len(z)
print np.exp(H)/(np.exp(H)+3*np.exp(R))
u=z[:,0]+np.random.normal(0,.1,len(z[:,0]))
v=z[:,1]+np.random.normal(0,.1,len(z[:,1]))
axes[0].scatter(u,v)
Hval,Rval=5,1
model.sample(N,thin=t)       

z=model.trace('X')[:]
u=z[:,0]+np.random.normal(0,.1,len(z[:,0]))
v=z[:,1]+np.random.normal(0,.1,len(z[:,1]))
axes[1].scatter(u,v)

show()



#import pydot
#graph=pydot.Dot(graph_type='graph')
#graph.set_graphviz_executables({'dot':'C:\\Program Files (x86)\Graphviz2.38\\bin\dot.exe'})

#import scipy.misc
#pymc.graph.graph(model,path='Desktop',name='my_graph',format='png',prog='C:\Program Files (x86)\Graphviz2.38\\bin\dot.exe')
#plt.figure(figsize=(8,8))
#plt.imshow(scipy.misc.imread('my_graph.png'))
#plt.close()