BRN
- 107
- 10
- TL;DR
- [Python] Minimization likelihood function in python.
Hallo at all!
I'm learning statistic in python and I have a problem to show you.
I have this parametric function:
$$P(S|t, \gamma, \beta)=\langle s(t) \rangle
\left( \frac{\gamma-\beta}{\gamma\langle
s(t) \rangle -\beta}\right)^2\left( 1- \frac{\gamma-\beta}{\gamma\langle
s(t) \rangle -\beta}\right)^{S-1}$$
$$\langle s(t) \rangle= e^{(\gamma-\beta)t}$$
and i want to estimate ##\beta## and ##\gamma## parameters at ##t=8## and ##t=10## by fitting experimental data with MLE method.
My code is this:
[CODE lang="python" title="My code"]
# parameters
beta_initial = 0.2
Gamma2_initial = 0.55
N = 2000 #osservations number
#theoretical function p(S|t,gamma,beta) for t=8 (rho8) and t=10 (rho10)
def rhos2(Gamma2, beta):
ms_t8 = scipy.exp((Gamma2 - beta)*8)
rho2_8 = scipy.zeros(N)
for n in range(0,N):
ratio8 = (Gamma2 - beta)/(Gamma2*ms_t8 - beta)
rho2_8[n] = ms_t8*((ratio8)**2)*(1 - ratio8)**(n-1)
rho2_8[0] = 0.
ms_t10 = scipy.exp((Gamma2 - beta)*10)
rho2_10=scipy.zeros(N)
for n in range(0,N):
ratio10 = (Gamma2 - beta)/(Gamma2*ms_t10 - beta)
rho2_10[n] = ms_t10*((ratio10)**2)*(1 - ratio10)**(n-1)
rho2_10[0] = 0.
return rho2_8, rho2_10
# Likelihood function (day8size, day10size experimental data)
def likelihood2(Gamma2, beta):
rho2_8, rho2_10 = rhos2(Gamma2, beta)
likelihood2 = 0.
for i in range(len(day8size)):
likelihood2 += -scipy.log(rho2_8[day8size])
for i in range(len(day10size)):
likelihood2 += -scipy.log(rho2_10[day10size])
return likelihood2
# Max of Likelihood (minimize -L)
bestfit2 = scipy.optimize.minimize(likelihood2, x0 = Gamma2_initial, args = (beta_initial,),
method = "Nelder-Mead")
#Gamma2_best2, beta_best = bestfit2
#print(Gamma2_best2, likelihood2(Gamma2_best2, beta_best))
#rho2_8, rho2_10 = rhos2(Gamma2_best2, beta_best)
[/CODE]
Now, I dont'n know if my scipy.optimize.minimize is loaded correctly and I don't know how to extract my best ##\gamma## and ##\beta## values from the fit.
Can someone help me?
Thank you!
I'm learning statistic in python and I have a problem to show you.
I have this parametric function:
$$P(S|t, \gamma, \beta)=\langle s(t) \rangle
\left( \frac{\gamma-\beta}{\gamma\langle
s(t) \rangle -\beta}\right)^2\left( 1- \frac{\gamma-\beta}{\gamma\langle
s(t) \rangle -\beta}\right)^{S-1}$$
$$\langle s(t) \rangle= e^{(\gamma-\beta)t}$$
and i want to estimate ##\beta## and ##\gamma## parameters at ##t=8## and ##t=10## by fitting experimental data with MLE method.
My code is this:
[CODE lang="python" title="My code"]
# parameters
beta_initial = 0.2
Gamma2_initial = 0.55
N = 2000 #osservations number
#theoretical function p(S|t,gamma,beta) for t=8 (rho8) and t=10 (rho10)
def rhos2(Gamma2, beta):
ms_t8 = scipy.exp((Gamma2 - beta)*8)
rho2_8 = scipy.zeros(N)
for n in range(0,N):
ratio8 = (Gamma2 - beta)/(Gamma2*ms_t8 - beta)
rho2_8[n] = ms_t8*((ratio8)**2)*(1 - ratio8)**(n-1)
rho2_8[0] = 0.
ms_t10 = scipy.exp((Gamma2 - beta)*10)
rho2_10=scipy.zeros(N)
for n in range(0,N):
ratio10 = (Gamma2 - beta)/(Gamma2*ms_t10 - beta)
rho2_10[n] = ms_t10*((ratio10)**2)*(1 - ratio10)**(n-1)
rho2_10[0] = 0.
return rho2_8, rho2_10
# Likelihood function (day8size, day10size experimental data)
def likelihood2(Gamma2, beta):
rho2_8, rho2_10 = rhos2(Gamma2, beta)
likelihood2 = 0.
for i in range(len(day8size)):
likelihood2 += -scipy.log(rho2_8[day8size])
for i in range(len(day10size)):
likelihood2 += -scipy.log(rho2_10[day10size])
return likelihood2
# Max of Likelihood (minimize -L)
bestfit2 = scipy.optimize.minimize(likelihood2, x0 = Gamma2_initial, args = (beta_initial,),
method = "Nelder-Mead")
#Gamma2_best2, beta_best = bestfit2
#print(Gamma2_best2, likelihood2(Gamma2_best2, beta_best))
#rho2_8, rho2_10 = rhos2(Gamma2_best2, beta_best)
[/CODE]
Now, I dont'n know if my scipy.optimize.minimize is loaded correctly and I don't know how to extract my best ##\gamma## and ##\beta## values from the fit.
Can someone help me?
Thank you!