Minimization likelihood function with parameters

[BRN]

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!

 

[pasmith]

It looks like you're trying to optimize over two variables. So your likelihood function needs to take a 1D array of length 2 as an argument, and return a float. Your likelihood2 takes two float arguments. Thus you need to pass a lambda to minimize (this is allows the arguments to likelihood2 to retain their descriptive names). Similarly the initial guess needs to be a 1D array of length 2, in the same order as the arugments to likelihood2.

import scipy.optimize.minimize
import numpy as np

# Likelihood function (day8size, day10size experimental data)
def likelihood2(Gamma2, beta):
    rho2_8, rho2_10 = rhos2(Gamma2, beta)
  
    # Best not to use the function's name as a variable in its body.
    result = 0.
    for i in range(len(day8size)):
        result += -scipy.log(rho2_8[day8size[i]])
      
    for i in  range(len(day10size)):
        result += -scipy.log(rho2_10[day10size[i]])
      
    return result

bestfit2 = scipy.optimize.minimize(
   lambda x, *args: likelihood2(x[0],x[1]), 
   # I would prefer to have likelihood2 calculate the actual likelihood, 
   # and then you can multiply it by -1
   # as part of this lambda.
   x0 = np.array([Gamma2_initial, beta_initial]),
   method = "Nelder-Mead"
)

# If it didn't work, tell someone.
if not bestfit2.success:
   raise RuntimeError('Optimisation failed with status {0:%d}: {1:%s}'.format(bestfit2.status, bestfit2.message)

# If it did, you can get the results as follows:
Gamma2_best = bestfit2.x[0]
beta_best = bestfit2.x[1]

All of this is explained in the documentation.

 

[BRN]

Hallo pasmith!
Thank you for help me.

I tried to modify my code with your way but it still doesn't work ...
Now the error is:
[CODE title="Errorr"]
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-29-96302e7f6f19> in <module>
9 lambda x, *args: likelihood2(x[0], x[1]),
10 x0 = np.array([[Gamma2_initial], [beta_initial]]),
---> 11 method = "Nelder-Mead"
12 )
13

~/.local/lib/python3.7/site-packages/scipy/optimize/_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
586 callback=callback, **options)
587 elif meth == 'nelder-mead':
--> 588 return _minimize_neldermead(fun, x0, args, callback, **options)
589 elif meth == 'powell':
590 return _minimize_powell(fun, x0, args, callback, **options)

~/.local/lib/python3.7/site-packages/scipy/optimize/optimize.py in _minimize_neldermead(func, x0, args, callback, maxiter, maxfev, disp, return_all, initial_simplex, xatol, fatol, adaptive, **unknown_options)
583
584 for k in range(N + 1):
--> 585 fsim[k] = func(sim[k])
586
587 ind = numpy.argsort(fsim)

TypeError: float() argument must be a string or a number, not 'function'
[/CODE]

It is not possible to pass a function like argument...

Any idea?

 

[pasmith]

What are the types of Gamma2_initial and beta_initial? They need to be floats.
What type does likelihood2 return? It should also be a float.

 

[BRN]

I'm sorry for the delay in answering you, but I had to study a lot for an exam.

in fact I had made a mistake in defining my function. but now i settled. However, minimization still doesn't work...

Thisi is the error:
[CODE title="output"]/home/brando/.local/lib/python3.7/site-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in double_scalars
import sys
/home/brando/.local/lib/python3.7/site-packages/ipykernel_launcher.py:14: RuntimeWarning: invalid value encountered in double_scalars

/home/brando/.local/lib/python3.7/site-packages/ipykernel_launcher.py:8: RuntimeWarning: overflow encountered in double_scalars

/home/brando/.local/lib/python3.7/site-packages/ipykernel_launcher.py:15: RuntimeWarning: overflow encountered in double_scalars
from ipykernel import kernelapp as app
/home/brando/.local/lib/python3.7/site-packages/ipykernel_launcher.py:11: RuntimeWarning: invalid value encountered in cdouble_scalars
# This is added back by InteractiveShellApp.init_path()
/home/brando/.local/lib/python3.7/site-packages/scipy/optimize/optimize.py:613: ComplexWarning: Casting complex values to real discards the imaginary part
fsim[-1] = fxr
/home/brando/.local/lib/python3.7/site-packages/scipy/optimize/optimize.py:636: ComplexWarning: Casting complex values to real discards the imaginary part
fsim[-1] = fxcc
/home/brando/.local/lib/python3.7/site-packages/scipy/optimize/optimize.py:596: RuntimeWarning: invalid value encountered in subtract
numpy.max(numpy.abs(fsim[0] - fsim[1:])) <= fatol):

---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-9-83a044d47d73> in <module>
14 if not bestfit2.success:
15 raise RuntimeError('Optimisation failed with status {0:%d}: {1:%s}'.format(bestfit2.status,
---> 16 bestfit2.message))
17
18 # If it did, you can get the results as follows:

ValueError: Invalid format specifier[/CODE]

The two datasets (day8size.dat and day10size.dat) contain only values in range [1 - 1000] and [1 - 2000] respectively.

I dont'n know the correct initial beta and gamma values. So I tried to change them in range [0 - 1], but it wasn't enough...

 

[pasmith]

My mistake.

'Optimisation failed with status {0:%d}: {1:%s}'

should be

'Optimisation failed with status {0:d}: {1:s}'

However that line should not have been executed unless something went wrong with the actual optimization.

 

[BRN]

Anyway, thanks for the help!