KUDOS https://www.fuw.edu.pl/~zarnecki/SAED/saed23.html¶
Bonamente Chap. 22 - Statistics and Analysis of Scientific Data - Springer https://link.springer.com/book/10.1007/978-981-19-0365-6
A simple explorer based on $\chi^2$¶
In [28]:
import math
import numpy as np
from matplotlib import pyplot as plt
from uncertainties import ufloat
from scipy import integrate
In [112]:
!pip install uncertainties
Collecting uncertainties Downloading uncertainties-3.2.2-py3-none-any.whl.metadata (6.9 kB) Downloading uncertainties-3.2.2-py3-none-any.whl (58 kB) Installing collected packages: uncertainties Successfully installed uncertainties-3.2.2
In [156]:
# Read background fit output
#infile = "http://hep.fuw.edu.pl/u/zarnecki/saed23/11_homework_data.dat"
### KUDOS https://www.fuw.edu.pl/~zarnecki/SAED/saed23.html
infile = "./11_homework_data.dat"
ttab = np.loadtxt(infile)
Nevt=ttab.size
In [18]:
ttab.shape
Out[18]:
(1000,)
In [19]:
print(ttab[0:12])
[4.8860e+00 3.0000e-03 4.2780e+00 4.4680e+00 1.3860e+00 6.4130e+00 1.1416e+01 2.5940e+00 2.7370e+00 4.4900e-01 8.1800e-01 1.0990e+00]
$$ \frac{N_{dec}}{ t_{dec}} \cdot \exp(-x/t_{dec}) + N_{bg}$$
In [25]:
def myExponentialPlusConts(x,par):
'''
Background is assumed to be a constant equal to par[2]
'''
# decay time
tdec = par[0]
# Decay normalisation
ndec = par[1]
# Background normalisation
nbg = par[2]
# Exponential probability distribution with flat background
val = ndec/tdec * np.exp(-x/tdec) + nbg
return val
In [56]:
def ChiSquareTestStatistics(xvec,yvec,errvec,myfun,par):
fval = myfun(xvec,par)
dy = (yvec-fval)/errvec
chi2 = np.sum(dy*dy)
return chi2
In [60]:
ChiSquareTestStatistics(centerbins,c,np.sqrt(c),myExponentialPlusConts,[2.492,450,2.74])
Out[60]:
np.float64(22.861215201942166)
In [62]:
Nb=30
bins=np.linspace(0,15,Nb+1)
c, b, p = plt.hist(ttab,bins=bins)
centerbins=(b[:-1]+b[1:])/2.
binsize=bins[1]-bins[0]
plt.xlabel('Time $\\mu s$')
plt.ylabel('Counts/'+str(binsize)+'$\\mu s^{-1}$')
xs=np.linspace(0,15,100)
ys=myExponentialPlusConts(xs,[2.492,450,2.74])
plt.plot(xs,ys,'r')
plt.yscale('log')
ChiSquareTestStatistics(centerbins,c,np.sqrt(c),myExponentialPlusConts,[2.492,450,2.74])
Out[62]:
np.float64(22.861215201942166)
In [157]:
ufloat(3,0.12)
Out[157]:
3.0+/-0.12
In [159]:
# Main loop
Nmc = 1000
# Desired ranges of parameter variation (resolution of the exploration)
step = np.array([0.20,50,1.])
Npars = step.shape[0]
# container for the history of the parameters explored
parhist=np.empty(shape=(Nmc,Npars))
parametersInitialValues=[7.2,722,7.12]
par = np.array(parametersInitialValues)
chi2old = ChiSquareTestStatistics(centerbins,c,np.sqrt(c),myExponentialPlusConts,par)
chiHistory=[chi2old] #as a list so that later I can add to it with +=[] insted of verbose numpy
for imc in range(Nmc):
#displace the parameter by a random jump
u=np.random.uniform(-1.,1.,size=Npars)
dpar = step*u
chi2new = ChiSquareTestStatistics(centerbins,c,np.sqrt(c),myExponentialPlusConts,par+dpar)
# probability you want to keep this parameter choice
prob = np.exp(-0.5*(chi2new-chi2old))
r = np.random.uniform()
# ACCEPT or REJECT the point based on its comparison with a random [0,1] uniform number
if r < prob:
# probability of the new set of parameters is high, we keep it and we update the reference chi2 value
par += dpar
chi2old = chi2new
# if the prob was low par has not been updated and the algorithm remains in the old point
chiHistory+=[chi2new]
parhist[imc] = par
In [166]:
np.array(chiHistory)[1:0]
Out[166]:
array([], dtype=float64)
In [224]:
def info_histogram(axs,r,c,_data,start=None, end=None):
data=_data[start:end]
_u=ufloat(np.mean(data), np.std(data) )
axs[r,c].hist(data,label=_u,density=True,histtype='step')
axs[r,c].legend()
In [178]:
fig, axs = plt.subplots(2,2,layout='constrained',figsize=np.array([11,11])/1.4)
info_histogram(axs,0,0,parhist[:,0])
info_histogram(axs,0,0,parhist[:,0],start=200)
axs[0,0].set_xlabel('$\\tau$')
info_histogram(axs,0,1,parhist[:,1])
info_histogram(axs,0,1,parhist[:,1],start=200)
axs[0,1].set_xlabel('$N_\\mu$')
info_histogram(axs,1,0,parhist[:,2])
info_histogram(axs,1,0,parhist[:,2],start=200)
axs[1,0].set_xlabel('$N_B$')
info_histogram(axs,1,1,np.array(chiHistory)/Nb)
info_histogram(axs,1,1,np.array(chiHistory)/Nb,start=200)
axs[1,1].set_xlabel('$\\chi^2$')
Out[178]:
Text(0.5, 0, '$\\chi^2$')
In [175]:
fig, axs = plt.subplots(4,1,layout='constrained',figsize=np.array([11,22])/1.4)
def info_evolution(axs,r,_data,start=None,end=None,offset=0):
data=_data[start:end]
_Nmc=data.shape[0]
xs = np.linspace(1.,_Nmc,_Nmc)
xs=xs+offset
_u=ufloat(np.mean(data), np.std(data) )
axs[r].plot(xs,data,label=_u)
axs[r].legend()
info_evolution(axs,0,parhist[:,0])
info_evolution(axs,0,parhist[:,0],start=200,offset=200)
info_evolution(axs,0,parhist[:,0],start=400,offset=400)
axs[0].set_xlabel('$\\tau$')
info_evolution(axs,1,parhist[:,1])
axs[1].set_xlabel('$N_\\mu$')
info_evolution(axs,2,parhist[:,2])
axs[2].set_xlabel('$N_B$')
info_evolution(axs,3,np.array(chiHistory)/Nb)
info_evolution(axs,3,np.array(chiHistory)/Nb,start=200,offset=200)
info_evolution(axs,3,np.array(chiHistory)/Nb,start=400,offset=400)
axs[3].set_xlabel('$\\chi^2$')
Out[175]:
Text(0.5, 0, '$\\chi^2$')
Progress of the overall fit¶
In [208]:
last100=parhist[100:200]
xs=np.linspace(0,15,100)
for i in range(len(last100)):
ys=myExponentialPlusConts(xs,last100[i])
plt.plot(xs,ys,'g',linewidth=0.05)
plt.errorbar(centerbins,c,np.sqrt(c))
plt.yscale('log')
In [209]:
last100=parhist[200:300]
xs=np.linspace(0,15,100)
for i in range(len(last100)):
ys=myExponentialPlusConts(xs,last100[i])
plt.plot(xs,ys,'g',linewidth=0.05)
plt.errorbar(centerbins,c,np.sqrt(c))
plt.yscale('log')
In [210]:
last100=parhist[300:400]
xs=np.linspace(0,15,100)
for i in range(len(last100)):
ys=myExponentialPlusConts(xs,last100[i])
plt.plot(xs,ys,'g',linewidth=0.05)
plt.errorbar(centerbins,c,np.sqrt(c))
plt.yscale('log')
In [207]:
last100=parhist[-100:]
xs=np.linspace(0,15,100)
for i in range(len(last100)):
ys=myExponentialPlusConts(xs,last100[i])
plt.plot(xs,ys,'g',linewidth=0.05)
plt.errorbar(centerbins,c,np.sqrt(c))
plt.yscale('log')
Priors and constraints¶
Some numbers can become negative, where they have no meaning, because they are "counts".¶
Let us review the Metropolis-Hastings recipe, recalling Bayes
$$ \mathcal{P}(A|B)\mathcal{P}(B)=\mathcal{P}(B|A)\mathcal{P}(A) $$
The acceptance probability for $\theta^\prime$ from a current state $\theta$ is \begin{equation} \alpha(\theta \to \theta^\prime) = \min\left(\frac{\pi(\theta^\prime)\cdot q(\theta|\theta^\prime)}{\pi(\theta)\cdot q(\theta^\prime | \theta)} ,1\right) \label{eq:one} \tag{1} \end{equation}
where
- $q$ is the "proposal" distribution, that is the one proposing the jump
dpar - $\pi$ is desired stationaty probability of the chain, i.e. the probability densitiy function for the parameters to be determined based in the observed data Z, that is to say $\pi(\theta)$ is $\pi(\theta|Z)$.
Due to Bayes we have for a generic $\theta_1$ $$ \pi(\theta_1) \equiv P(\theta_1|Z) = \frac{p(\theta_1)\cdot P(Z|\theta_1) }{p(Z)} $$ which entails from eq.$\eqref{eq:one}$
\begin{equation} \alpha(\theta \to \theta^\prime) = \min\left(\ \frac{p(\theta^\prime)\cdot P(Z|\theta^\prime)\cdot q(\theta|\theta^\prime)}{p(\theta) \cdot P(Z|\theta) \cdot q(\theta^\prime | \theta)} ,1\right) \label{eq:two} \tag{2} \end{equation}
We recognize that
- $P(Z|\theta)$ is by definition $\mathcal{L}(\theta)$ for any $\theta$
All the quantities are known as we choose $p$ and $q$. In particular it is worth to notice that a the denominated we have
- $p(\theta^\prime)\cdot q(\theta^\prime|\theta)$ that is the probabiity of proposing the given $\theta^\prime $ from the present point $\theta$.
For uniform priors $p$ and $q$ this reduces to
$$\alpha(\theta \to \theta^\prime) = \min\left(
\frac{\mathcal{L(\theta^\prime)}}{\mathcal{L}(\theta)} \
,1\right) \tag{3} $$
as in the first implementation.
In [213]:
all(np.array([1,2,3])>0)
Out[213]:
True
In [214]:
# Main loop
Nmc = 1000
# Desired ranges of parameter variation (resolution of the exploration)
step = np.array([0.20,50,1.])
Npars = step.shape[0]
def prior(v:np.array):
if all(v>0):# all par > 0
return 1
else:
return 0
# container for the history of the parameters explored
parhist=np.empty(shape=(Nmc,Npars))
parametersInitialValues=[7.2,722,7.12]
par = np.array(parametersInitialValues)
chi2old = ChiSquareTestStatistics(centerbins,c,np.sqrt(c),myExponentialPlusConts,par)
chiHistory=[chi2old] #as a list so that later I can add to it with +=[] insted of verbose numpy
for imc in range(Nmc):
#displace the parameter by a random jump
u=np.random.uniform(-1.,1.,size=Npars)
dpar = step*u
newpar=par+dpar
chi2new = ChiSquareTestStatistics(centerbins,c,np.sqrt(c),myExponentialPlusConts,newpar)
# 2*LogL = - chi^2
# L = exp (-chi^2/2)
# probability you want to keep this parameter choice
LikelihoodRatio= np.exp(-0.5*(chi2new-chi2old))
PriorsRatio=prior(newpar)/prior(par)
prob = min(1, LikelihoodRatio*PriorsRatio)
r = np.random.uniform()
# ACCEPT or REJECT the point based on its comparison with a random [0,1] uniform number
if r < prob:
# probability of the new set of parameters is high, we keep it and we update the reference chi2 value
par += dpar
chi2old = chi2new
# if the prob was low par has not been updated and the algorithm remains in the old point
chiHistory+=[chi2new]
parhist[imc] = par
In [215]:
fig, axs = plt.subplots(4,1,layout='constrained',figsize=np.array([11,22])/1.4)
def info_evolution(axs,r,_data,start=None,end=None,offset=0):
data=_data[start:end]
_Nmc=data.shape[0]
xs = np.linspace(1.,_Nmc,_Nmc)
xs=xs+offset
_u=ufloat(np.mean(data), np.std(data) )
axs[r].plot(xs,data,label=_u)
axs[r].legend()
info_evolution(axs,0,parhist[:,0])
info_evolution(axs,0,parhist[:,0],start=200,offset=200)
info_evolution(axs,0,parhist[:,0],start=400,offset=400)
axs[0].set_xlabel('$\\tau$')
info_evolution(axs,1,parhist[:,1])
axs[1].set_xlabel('$N_\\mu$')
info_evolution(axs,2,parhist[:,2])
axs[2].set_xlabel('$N_B$')
info_evolution(axs,3,np.array(chiHistory)/Nb)
info_evolution(axs,3,np.array(chiHistory)/Nb,start=200,offset=200)
info_evolution(axs,3,np.array(chiHistory)/Nb,start=400,offset=400)
axs[3].set_xlabel('$\\chi^2$')
Out[215]:
Text(0.5, 0, '$\\chi^2$')
In [225]:
fig, axs = plt.subplots(2,2,layout='constrained',figsize=np.array([11,11])/1.4)
info_histogram(axs,0,0,parhist[:,0])
info_histogram(axs,0,0,parhist[:,0],start=200)
info_histogram(axs,0,0,parhist[:,0],start=900)
axs[0,0].set_xlabel('$\\tau$')
info_histogram(axs,0,1,parhist[:,1])
info_histogram(axs,0,1,parhist[:,1],start=200)
axs[0,1].set_xlabel('$N_\\mu$')
info_histogram(axs,1,0,parhist[:,2])
info_histogram(axs,1,0,parhist[:,2],start=200)
axs[1,0].set_xlabel('$N_B$')
info_histogram(axs,1,1,np.array(chiHistory)/Nb)
info_histogram(axs,1,1,np.array(chiHistory)/Nb,start=200)
axs[1,1].set_xlabel('$\\chi^2$')
Out[225]:
Text(0.5, 0, '$\\chi^2$')
Progress of the overall fit¶
In [216]:
last100=parhist[100:200]
xs=np.linspace(0,15,100)
for i in range(len(last100)):
ys=myExponentialPlusConts(xs,last100[i])
plt.plot(xs,ys,'g',linewidth=0.05)
plt.errorbar(centerbins,c,np.sqrt(c))
plt.yscale('log')
In [218]:
last100=parhist[200:300]
xs=np.linspace(0,15,100)
for i in range(len(last100)):
ys=myExponentialPlusConts(xs,last100[i])
plt.plot(xs,ys,'g',linewidth=0.05)
plt.errorbar(centerbins,c,np.sqrt(c))
plt.yscale('log')
In [219]:
last100=parhist[300:400]
xs=np.linspace(0,15,100)
for i in range(len(last100)):
ys=myExponentialPlusConts(xs,last100[i])
plt.plot(xs,ys,'g',linewidth=0.05)
plt.errorbar(centerbins,c,np.sqrt(c))
plt.yscale('log')
In [207]:
last100=parhist[-100:]
xs=np.linspace(0,15,100)
for i in range(len(last100)):
ys=myExponentialPlusConts(xs,last100[i])
plt.plot(xs,ys,'g',linewidth=0.05)
plt.errorbar(centerbins,c,np.sqrt(c))
plt.yscale('log')
Hit or miss à la MH¶
Note that the Metropolis-Hasting formula also provides a sort of "hit or miss" procedure. Adapting eq.($\ref{1}$) for flat $q$ the M-H procedure becomes the "hit or miss" for the aymptotyc pdf $\pi$
In [5]:
# generate in the hyper-cube [0,1]^n
try:
del pior
except:
def prior(x:np.array):
if any(x>1):
return 0
if any(x<0):
return 0
return 1
In [7]:
prior(np.array([1,2,2]))
Out[7]:
0
In [48]:
def pdf(x):
return 0.2+np.exp(-3*x)
In [49]:
integral=integrate.quad(pdf, 0, 1)#,args=(lam,))
In [50]:
def npdf(x):
return pdf(x)/integral[0]
In [55]:
# Main loop
Nmc = int(10e4)
# Desired ranges of parameter variation (resolution of the exploration)
step = np.array([0.1])
Npars = step.shape[0]
# container for the history of the parameters explored
parhist=np.empty(shape=(Nmc,Npars))
parametersInitialValues=[0.5]
par = np.array(parametersInitialValues)
for imc in range(Nmc):
#displace the parameter by a random jump
u=np.random.uniform(-1.,1.,size=Npars)
dpar = step*u
newpar=par+dpar
priorRatio=prior(newpar)/prior(par)
pdfRatio=pdf(newpar)/pdf(par)
prob = min(1, pdfRatio*priorRatio)
r = np.random.uniform()
# ACCEPT or REJECT the point based on its comparison with a random [0,1] uniform number
if r < prob:
# probability of the new set of parameters is high, we keep it and we update the reference chi2 value
par += dpar
# if the prob was low par has not been updated and the algorithm remains in the old point
parhist[imc] = par
In [56]:
xs=np.linspace(0,1,1000)
ys=npdf(xs)
plt.plot(xs,ys)
plt.hist(parhist,density=True)
Out[56]:
(array([1.98082937, 1.61670559, 1.31178567, 1.04696838, 0.90325899,
0.78895153, 0.67864432, 0.61954046, 0.54383552, 0.51013332]),
array([2.50225989e-05, 1.00018492e-01, 2.00011960e-01, 3.00005429e-01,
3.99998898e-01, 4.99992367e-01, 5.99985836e-01, 6.99979305e-01,
7.99972774e-01, 8.99966243e-01, 9.99959712e-01]),
<BarContainer object of 10 artists>)
however it s does not work much if the scanning step is too small
In [57]:
# Main loop
Nmc = int(10e4)
# Desired ranges of parameter variation (resolution of the exploration)
step = np.array([0.01])
Npars = step.shape[0]
# container for the history of the parameters explored
parhist=np.empty(shape=(Nmc,Npars))
parametersInitialValues=[0.5]
par = np.array(parametersInitialValues)
for imc in range(Nmc):
#displace the parameter by a random jump
u=np.random.uniform(-1.,1.,size=Npars)
dpar = step*u
newpar=par+dpar
priorRatio=prior(newpar)/prior(par)
pdfRatio=pdf(newpar)/pdf(par)
prob = min(1, pdfRatio*priorRatio)
r = np.random.uniform()
# ACCEPT or REJECT the point based on its comparison with a random [0,1] uniform number
if r < prob:
# probability of the new set of parameters is high, we keep it and we update the reference chi2 value
par += dpar
# if the prob was low par has not been updated and the algorithm remains in the old point
parhist[imc] = par
In [58]:
xs=np.linspace(0,1,1000)
ys=npdf(xs)
plt.plot(xs,ys)
plt.hist(parhist,density=True)
Out[58]:
(array([2.61725276, 2.10360316, 1.70806496, 0.87928492, 0.57855587,
0.4369422 , 0.37633635, 0.43654216, 0.42604115, 0.43834233]),
array([7.02051412e-06, 9.99973629e-02, 1.99987705e-01, 2.99978048e-01,
3.99968390e-01, 4.99958732e-01, 5.99949075e-01, 6.99939417e-01,
7.99929760e-01, 8.99920102e-01, 9.99910444e-01]),
<BarContainer object of 10 artists>)
A much longer run starts to give sensible results
In [64]:
# Main loop
Nmc = int(10e5)
# Desired ranges of parameter variation (resolution of the exploration)
step = np.array([0.01])
Npars = step.shape[0]
# container for the history of the parameters explored
parhist=np.empty(shape=(Nmc,Npars))
parametersInitialValues=[0.5]
par = np.array(parametersInitialValues)
for imc in range(Nmc):
#displace the parameter by a random jump
u=np.random.uniform(-1.,1.,size=Npars)
dpar = step*u
newpar=par+dpar
priorRatio=prior(newpar)/prior(par)
pdfRatio=pdf(newpar)/pdf(par)
prob = min(1, pdfRatio*priorRatio)
r = np.random.uniform()
# ACCEPT or REJECT the point based on its comparison with a random [0,1] uniform number
if r < prob:
# probability of the new set of parameters is high, we keep it and we update the reference chi2 value
par += dpar
# if the prob was low par has not been updated and the algorithm remains in the old point
parhist[imc] = par
In [60]:
xs=np.linspace(0,1,1000)
ys=npdf(xs)
plt.plot(xs,ys)
plt.hist(parhist,density=True)
Out[60]:
(array([2.22292289, 1.85699748, 1.46043161, 1.1761074 , 0.99488472,
0.69805033, 0.54653809, 0.42048622, 0.32728484, 0.29644439]),
array([1.48691670e-07, 9.99986690e-02, 1.99997189e-01, 2.99995710e-01,
3.99994230e-01, 4.99992750e-01, 5.99991271e-01, 6.99989791e-01,
7.99988311e-01, 8.99986832e-01, 9.99985352e-01]),
<BarContainer object of 10 artists>)
In [65]:
xs=np.linspace(0,1,1000)
ys=npdf(xs)
plt.plot(xs,ys)
plt.hist(parhist,density=True)
Out[65]:
(array([1.9261059 , 1.63796502, 1.45475446, 1.20246368, 0.84385259,
0.6197119 , 0.59968184, 0.60019184, 0.57479176, 0.54051166]),
array([2.96101880e-07, 9.99999897e-02, 1.99999683e-01, 2.99999377e-01,
3.99999071e-01, 4.99998764e-01, 5.99998458e-01, 6.99998151e-01,
7.99997845e-01, 8.99997539e-01, 9.99997232e-01]),
<BarContainer object of 10 artists>)
In [63]:
xs=np.linspace(0,1,1000)
ys=npdf(xs)
plt.plot(xs,ys)
plt.hist(parhist,density=True)
Out[63]:
(array([2.11779343, 1.60334259, 1.233732 , 1.00940163, 0.93053151,
0.79326128, 0.65791106, 0.60525098, 0.54569088, 0.50310081]),
array([8.91492965e-07, 1.00000730e-01, 2.00000568e-01, 3.00000406e-01,
4.00000244e-01, 5.00000083e-01, 5.99999921e-01, 6.99999759e-01,
7.99999597e-01, 8.99999436e-01, 9.99999274e-01]),
<BarContainer object of 10 artists>)
In [ ]: