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)
    return total


# # 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)

# choose a bad initialization
start = np.array([.1, .1])
 
# define the markov chain object
mc = pypmc.sampler.markov_chain.AdaptiveMarkovChain(log_target, prop, start)

# 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)
   mc.adapt()
 
# 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.

Delete this entry (admin only).