0

Im looking forward to minimize this function below to estimate the parameters of the normal distribution

Function image

My code looks like this:

import numpy as np
from scipy import stats
from scipy.optimize import minimize
x = [1,2,3,4,5]
def oro(theta, x):
    norma = 0 
    b = 1
    u = theta[0]
    o = theta[1]
    x = np.array(x)
    x0 = 0
    f0 = -(((1/(o*(2*3.14)**(0.5)))*(2.718)**-(((x0-u)**2)/(2*(o**2))))**b)**-1
    for i in range(x.size):
        f = (1/(o*(2*3.14)**(0.5)))*(2.718)**-(((x[i]-u)**2)/(2*(o**2)))**b
        norma += f0*f
    return norma
theta_init = [0, 1]
res = minimize(oro, theta_init, args=x)
res

But in the end I get this:

<ipython-input-81-ee81472a023a>:8: RuntimeWarning: divide by zero encountered in double_scalars
  f0 = -(((1/(o*(2*3.14)**(0.5)))*(2.718)**-(((x0-u)**2)/(2*(o**2))))**b)**-1
<ipython-input-81-ee81472a023a>:11: RuntimeWarning: invalid value encountered in double_scalars
  norma += f0*f
<ipython-input-81-ee81472a023a>:8: RuntimeWarning: divide by zero encountered in double_scalars
  f0 = -(((1/(o*(2*3.14)**(0.5)))*(2.718)**-(((x0-u)**2)/(2*(o**2))))**b)**-1
<ipython-input-81-ee81472a023a>:11: RuntimeWarning: invalid value encountered in double_scalars
  norma += f0*f
<ipython-input-81-ee81472a023a>:8: RuntimeWarning: divide by zero encountered in double_scalars
  f0 = -(((1/(o*(2*3.14)**(0.5)))*(2.718)**-(((x0-u)**2)/(2*(o**2))))**b)**-1
<ipython-input-81-ee81472a023a>:11: RuntimeWarning: invalid value encountered in double_scalars
  norma += f0*f
      fun: nan
 hess_inv: array([[9.57096191e+02, 2.41349815e+01],
       [2.41349815e+01, 8.33412317e-01]])
      jac: array([nan, nan])
  message: 'Desired error not necessarily achieved due to precision loss.'
     nfev: 357
      nit: 4
     njev: 119
   status: 2
  success: False
        x: array([165623.69347712,   1751.95100725])

Tell me, please, what am I doing wrong?

Update after 1 answer(added bounds). I get less errors but still unsuccessful:

<ipython-input-271-b51d0c455468>:8: RuntimeWarning: divide by zero encountered in double_scalars
  f0 = -(((1/(std*(2*np.pi)**(0.5)))*(np.exp(1))**-(((x0-mean)**2)/(2*(std**2))))**b)**-1
<ipython-input-271-b51d0c455468>:11: RuntimeWarning: invalid value encountered in double_scalars
  norma += f0*f
      fun: nan
 hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
      jac: array([-0.00012861,  0.00018581])
  message: 'ABNORMAL_TERMINATION_IN_LNSRCH'
     nfev: 75
      nit: 2
     njev: 25
   status: 2
  success: False
        x: array([250.13040562, 343.06899721])

2 Answers2

0

minimize does not know that the second term in your theta vector (I assume is the standard deviation) needs to be positive.

Try adding bounds:

res = minimize(oro, theta_init, args=x, bounds=((None, None), (1e-4, None))
Tarifazo
  • 4,118
  • 1
  • 9
  • 22
0

You can breakdown your complicated function to identify offending expression, add asserts to limit its value to a correct one.

Add bounds to limit the variables to a correct value range.

Code

import numpy as np
from scipy import stats
from scipy.optimize import minimize
x = [1,2,3,4,5]
def oro(theta, x):
    norma = 0 
    b = 1
    u = theta[0]
    o = theta[1]
    x = np.array(x)
    x0 = 0

    assert o > 0

    t = (o*(2*3.14)**(0.5))
    assert t > 0

    s = (2*(o**2))
    assert s > 0

    z = (((x0-u)**2)/s)

    m = (((1/t)*(2.718)**-z)**b)
    assert m != 0

    # f0 = -(((1/(o*(2*3.14)**(0.5)))*(2.718)**-(((x0-u)**2)/(2*(o**2))))**b)**-1
    f0 = -m**-1

    for i in range(x.size):
        f = (1/(o*(2*3.14)**(0.5)))*(2.718)**-(((x[i]-u)**2)/(2*(o**2)))**b
        norma += f0*f
    return norma
theta_init = [0, 1]
res = minimize(oro, theta_init, args=x, bounds=((1, 10), (1, 10)))
# res = minimize(oro, theta_init, args=x)
print(res)

Output

      fun: -1.932543956399183e+16
 hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>      
      jac: array([-9.65382468e+16,  1.44838865e+18])
  message: 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     nfev: 9
      nit: 2
     njev: 3
   status: 0
  success: True
        x: array([10.,  1.])
ferdy
  • 4,396
  • 2
  • 4
  • 16