# Demo entry 3541727

pypmc

Submitted by anonymous on Jan 11, 2016 at 16:40
Language: Python 3. Code size: 5.2 kB.

```import cmath
import math
import sys

central = {'mtms'   : 160.0,
'mwos'   : 80.0,
'lambda' : 0.22,
'ackm'   : 0.812,
'rho'    : 0.117,
'eta'    : 0.353}

class EffectiveTheory:
"""A Generic Effective Theory"""
def __init__(self, model):
if model == 'Standard Model':
self.name = model
self.g_sl_ll = 1.16637*10**-5
self.c_sl_ll = [[1,0,0],
[0,1,0],
[0,0,1]]
self.c_b_to_u_cb_d = [0.03,0.31]
else:
print('unknown model')
sys.exit(1)

def VubOverVcb(et):
return abs(et.c_sl_ll[0][2]/et.c_sl_ll[1][2])

def Gamma(et):
return (180.0/math.pi)*cmath.phase(-(et.c_sl_ll[0][0]*et.c_sl_ll[0][2].conjugate())/(et.c_sl_ll[1][0]*et.c_sl_ll[1][2].conjugate()))

def Sin2Beta(et):
return math.sin(2*cmath.phase(-(et.c_sl_ll[1][0]*et.c_sl_ll[1][2].conjugate())/(et.c_sl_ll[2][0]*et.c_sl_ll[2][2].conjugate())))

def Alpha(et):
return (180.0/math.pi)*cmath.phase(-(et.c_sl_ll[2][0]*et.c_sl_ll[2][2].conjugate())/(et.c_sl_ll[0][0]*et.c_sl_ll[0][2].conjugate()))

evaluate = {
'vubovcb'       : VubOverVcb,
'gamma'         : Gamma,
'sintwobeta'    : Sin2Beta,
'alpha'         : Alpha,
}

class StandardModel:
"""The Standard Model Class"""
def __init__(self, central):
self.para = central
self.ckm  = [[1,0,0],[0,1,0],[0,0,1]]
self.set_parameter([])

def set_parameter(self,fitparlist):
for parlist in fitparlist:
self.para[parlist[0]] = parlist[1]
self._set_ckm()

def _set_ckm(self):
lam = self.para['lambda']
rme = self.para['rho']-1j*self.para['eta']
omre = 1-self.para['rho']-1j*self.para['eta']
lt2 = lam*lam
lt3 = lt2*lam
alt2 = lt2*self.para['ackm']
alt3 = lt3*self.para['ackm']
self.vckm = [[1-lt2/2  , lam    , alt3*rme],
[-lam     , 1-lt2/2, alt2    ],
[alt3*omre,-alt2   , 1       ]]

def write_effective_theory(self, et):
et.c_sl_ll = self.vckm

et = EffectiveTheory('Standard Model')
model = StandardModel(central)

def observable_print(list, et):
model.set_parameter(list)
model.write_effective_theory(et)
print("input: ", list)
print(" --- observables start ---")
print("vub/vcb", evaluate['vubovcb'](et))
print("gamma"  , evaluate['gamma'](et))
print("s2beta" , evaluate['sintwobeta'](et))
print("alpha"  , evaluate['alpha'](et))
print(" --- observables  end  ---")

def equation(x, s, u, et, rlist):
for xvalue in x:
if (xvalue < -1) or (xvalue > 1):
return -1000.0
model.set_parameter([['rho',x[0]],['eta',x[1]]])
model.write_effective_theory(et)
t = [evaluate['vubovcb'](et),evaluate['gamma'](et),(evaluate['sintwobeta'](et)),evaluate['alpha'](et)]
total = 0
#    for i in range(len(t)):
for i in rlist:
#        print("ll : ", i, t[i], s[i], - ((t[i] - s[i])**2)/(u[i] * 2))
total += - ((t[i] - s[i])**2)/(u[i] * 2)

# # add proper mean and sigma of measurements
#mu = [0.0998,68,0.682,91.1]
mu = [0.0998,61.5,0.762,85.4]
sigma = [0.0122,7,0.064,3.8]

log_target = lambda x: equation(x, mu, sigma, et, [3])

# # test observables
# observable_print([['rho',1.0],['eta',0.01]],et)
# observable_print([['rho',0.2],['eta',0.5]],et)
# observable_print([['rho',0.2],['eta',0.4]],et)
# observable_print([['rho',0.1],['eta',0.4]],et)
# observable_print([['rho',0.2],['eta',0.3]],et)

# testx = [0.1,0.1]
# print("log likelihood: ", log_target(testx))

# testx = [0.1,0.2]
# print("log likelihood: ", log_target(testx))

#
# here comes the pypmc stuff to be commented out
#

import numpy as np
import pypmc

# define a proposal
prop_dof   = 1.
prop_sigma = np.array([[0.1 , 0.  ]
,[0.  , 0.1]])

prop = pypmc.density.student_t.LocalStudentT(prop_sigma, prop_dof)

start = np.array([.1, .1])

# define the markov chain object

# run burn-in
mc.run(10**5)

# delete burn-in from history
mc.clear()

# run 400,000 steps adapting the proposal every 2000 steps
# hereby save the accept count which is returned by mc.run
accept_count = 0
for i in range(200):
accept_count += mc.run(2000)

# extract a reference to the history of all visited points
values = mc.history[:]
accept_rate = float(accept_count) / len(values)
print("The chain accepted %4.2f%% of the proposed points" % (accept_rate * 100) )

# plot the result
try:
import matplotlib.pyplot as plt
except ImportError:
print('For plotting "matplotlib" needs to be installed')
exit(1)

plt.hexbin(values[:,0], values[:,1], gridsize = 50, cmap='gray_r')
# plt.axis([-1, 1, -1, 1])
cb = plt.colorbar()
cb.set_label('counts')
plt.show()
```

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.