scipy.optimize.leastsq возвращает параметры наилучшего предположения, а не новое наилучшее соответствие

Я хочу подогнать лоренцев пик к набору данных x и y, данные в порядке. Другие программы, такие как OriginLab, идеально подходят для этого, но я хотел автоматизировать настройку с помощью Python, поэтому у меня есть приведенный ниже код, основанный на http://mesa.ac.nz/?page_id=1800

Проблема, с которой я сталкиваюсь, заключается в том, что scipy.optimize.leastsq возвращает как наилучшее соответствие тем же начальным параметрам предположения, которые я ему передал, по сути ничего не делая. Вот код.

#x, y are the arrays with the x,y  axes respectively
#defining funcitons
def lorentzian(x,p):
  return p[2]*(p[0]**2)/(( x - (p[1]) )**2 + p[0]**2)

def residuals(p,y,x):
  err = y - lorentzian(x,p)
  return err      

p = [0.055, wv[midIdx], y[midIdx-minIdx]]   
pbest = leastsq(residuals, p, args=(y, x), full_output=1)
best_parameters = pbest[0]
print p
print pbest

p — начальные предположения, а best_parameters — возвращаемые параметры «наилучшего соответствия» из наименьшего квадрата, но они всегда одинаковы.

это то, что возвращает full_output=1 (длинные числовые массивы были сокращены, но по-прежнему репрезентативны)

    [0.055, 855.50732, 1327.0]
    (array([  5.50000000e-02,   8.55507324e+02,   1.32700000e+03]), 
    None, {'qtf':array([ 62.05192947,  69.98033905,  57.90628052]), 
    'nfev': 4, 
    'fjac': array([[-0.,  0.,  0., 0.,  0.,  0.,  0.,],
    [ 0., -0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
    0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
    0.,  0.],
    [ 0.,  0., -0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
    0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
    0.,  0.]]), 
    'fvec': array([  62.05192947,   69.98033905,   
    53.41218567,   45.49879837,   49.58242035,   36.66483688,
    34.74443436,   50.82238007,   34.89669037]), 
    'ipvt': array([1, 2, 3])},  
    'The cosine of the angle between func(x) and any column of the\n  Jacobian 
    is at most 0.000000 in absolute value', 4)

кто-нибудь может увидеть, что не так?


python scipy data-fitting
person Eoin Murray    schedule 18.09.2012    source источник
comment
О чем говорит информация, возвращенная в full_output=1? Может есть намек...   -  person seberg    schedule 18.09.2012
comment
обновлен, чтобы добавить, и удалил то, что я проверил, чтобы быть кодом, избыточным для проблемы   -  person Eoin Murray    schedule 18.09.2012

Ответы (1)


arrow_upward
5
arrow_downward

Быстрый поиск в Google намекает на проблему с данными с одинарной точностью (другие ваши программы почти наверняка также преобразуются в двойную точность, хотя это явная проблема и с scipy, см. также это отчет об ошибке). Если вы посмотрите на свой результат full_output=1, вы увидите, что якобиан везде приближается к нулю. Таким образом, указание якобиана явно может помочь (хотя даже в этом случае вы можете захотеть выполнить преобразование вверх, потому что минимальная точность для относительной ошибки, которую вы можете получить с одинарной точностью, очень ограничена).

Решение: самое простое и численно лучшее решение (конечно, предоставление реального якобиана также является бонусом) — просто привести ваши x и y данные к двойной точности (например, x = x.astype(np.float64) подойдет).

Я бы не советовал этого, но вы также можете исправить это, установив аргумент ключевого слова epsfcn (а также, возможно, аргументы ключевого слова допуска) вручную, что-то вместе с epsfcn=np.finfo(np.float32).eps. Кажется, это в какой-то мере решает проблему, но (поскольку большинство расчетов выполняется со скалярами, а скаляры не вызывают приведение вверх в ваших вычислениях) вычисления выполняются в float32, и потеря точности кажется довольно большой, по крайней мере, когда не предоставление Dfunc.

person seberg    schedule 18.09.2012