одиночный скаляр scipy optimize минимизировать - PullRequest
Новая
       89

одиночный скаляр scipy optimize минимизировать

0 голосов
/ 05 мая 2020

У меня проблема с использованием функции scipy minimize(), и я не совсем понимаю оптимизацию, чтобы asp что здесь не так ..

У меня есть функция, которая вызывает scipy.optimize.minimize(). Он отлично работает и предоставляет мне именно те выходные данные, которые мне нужны, когда x0 является массивом размером> 1, но когда x0 равно 1, он терпит неудачу. Документация говорит, что x0 должно быть np.ndarray размером (n,), но не указывает, что он должен быть> 1, поэтому я предположил, что это будет нормально. Меньшая версия моего кода, вызывающая функцию с оптимальным значением: 

import numpy as np
from scipy.optimize import minimize

def to_freq(*arrays):
    # Better version of `convert_to_freq()`
    out = []
    for a in arrays:
        converted = np.array([(x + i / len(a)) / (max(a)+1) for i, x in enumerate(a, start=1)])
        out.append(converted)
    return out

def likelihood(x, x_freq, expected, x_max):
    # Better version, supports vectorisation
    a = 2 * x * np.log(x_freq / expected) 
    b = 2 * (x_max - x) * np.log((1 - x_freq) / (1 - expected))
    return a + b

def objective(x0, labels, a, b):
    R = x0[labels=='R'].item()

    a_c, b_c = np.cumsum(a), np.cumsum(b)
    a_f, b_f = to_freq(a_c, b_c)

    # Get the expected values for signals and noises
    exp_a = ((1 - R) * b_f + R)[:-1]
    exp_b = b_f[:-1]

    # Compute the gsquared using the dual process model parameters
    #   Still getting runtime warnings about division. Function only works with numpy, so can't use math.
    a_lrat = likelihood(x=a_c[:-1], x_freq=a_f[:-1], expected=exp_a, x_max=a_c.max())
    b_lrat = likelihood(x=b_c[:-1], x_freq=b_f[:-1], expected=exp_b, x_max=b_c.max())

    return sum(a_lrat + b_lrat)

# Observations
a = [508,224,172,135,119,63]
b = [102,161,288,472,492,308]
x0 = np.array([0.520274590415736]) # Optimal value for variable
labels = np.array(['R'])

# Gives correct iotimized value of 163.27525607890783
objective(x0, labels, a, b)

И теперь случайная инициализация x0 для случаев, когда оптимальное значение неизвестно: 

x0 = np.random.uniform(-.5,0.5, len(labels)) # random initialization

# Without method='nelder-mead' occasionally gives correct value of fun, but frequently fails
opt = minimize(fun=objective, x0=x0, args=(labels, a, b), tol=1e-4)
print(opt)

Неудачный Результат оптимизации следующий: 

      fun: nan
 hess_inv: array([[1]])
      jac: array([nan])
  message: 'Desired error not necessarily achieved due to precision loss.'
     nfev: 336
      nit: 1
     njev: 112
   status: 2
  success: False
        x: array([1034.74])

Но если я продолжу запускать это и случайным образом устанавливаю начальное значение, он иногда выдаст хороший результат: 

      fun: 163.27525607888913
 hess_inv: array([[4.14149525e-05]])
      jac: array([-1.90734863e-05])
  message: 'Optimization terminated successfully.'
     nfev: 27
      nit: 7
     njev: 9
   status: 0
  success: True
        x: array([0.52027462])

Если я укажу method='nelder-mead' ( решение, возможно, не связанной проблемы ) в вызове minimize() в моей более крупной функции, он также фактически дает мне ожидаемый результат: 

 final_simplex: (array([[0.52026029],
       [0.52031204]]), array([163.27525856, 163.27527298]))
           fun: 163.2752585612531
       message: 'Optimization terminated successfully.'
          nfev: 32
           nit: 16
        status: 0
       success: True
             x: array([0.52026029])

Я действительно не понять, какой подход будет наилучшим для реализации этого, поскольку я очень неопытен в оптимизации.

[Сноска]: алгоритм минимизации иногда пытается значения, несовместимые с моей функцией (например, <0 или> 1) и вызов np.log() заканчивается предупреждением, но обычно я просто подавляю его, так как кажется, что он работает независимо ...

...