scipy.optimize.curve_fit не может установить сдвинутую перекошенную гауссовую кривую

Question

Я пытаюсь подогнать искаженную и сдвинутую гауссовую кривую, используя функцию scipy curve_fit, но я обнаружил, что при определенных условиях установка довольно плохая, часто давая мне приблизиться или ровно прямую.

Код, приведенный ниже, получен из документации curve_fit. Предоставляемый код представляет собой произвольный набор данных для тестирования, но достаточно хорошо отображает проблему.

import numpy as np 
from scipy.optimize import curve_fit 
import matplotlib.pyplot as plt 
import math as math 
import scipy.special as sp 

#def func(x, a, b, c): 
# return a*np.exp(-b*x) + c 

def func(x, sigmag, mu, alpha, c,a): 
    #normal distribution 
    normpdf = (1/(sigmag*np.sqrt(2*math.pi)))*np.exp(-(np.power((x-mu),2)/(2*np.power(sigmag,2)))) 
    normcdf = (0.5*(1+sp.erf((alpha*((x-mu)/sigmag))/(np.sqrt(2))))) 
    return 2*a*normpdf*normcdf + c 

x = np.linspace(0,100,100) 
y = func(x, 10,30, 0,0,1) 
yn = y + 0.001*np.random.normal(size=len(x)) 

popt, pcov = curve_fit(func, x, yn,) #p0=(9,35,0,9,1)) 

y_fit= func(x,popt[0],popt[1],popt[2],popt[3],popt[4]) 

plt.plot(x,yn) 
plt.plot(x,y_fit) 

Проблема, кажется, всплывают, когда я смещаться гауссовой слишком далеко от нуля (используя mu). Я попытался дать начальные значения, даже те, которые идентичны моей исходной функции, но это не решает проблему. Для значения mu=10, curve_fit работает отлично, но если я использую mu>=30, он больше не подходит для данных.

Answer 1

Давая отправную точку для минимизации, часто творит чудеса. Попробуйте дать минимизатор некоторую информацию о положении максимума и ширину кривой:

popt, pcov = curve_fit(func, x, yn, p0=(1./np.std(yn), np.argmax(yn) ,0,0,1)) 

меняющегося эту единственную строку в коде с sigma=10 и mu=50 производит

Answer 2

Вы можете позвонить по номеру curve_fit со случайным начальным предположением и выбрать параметры с минимальной ошибкой.

import numpy as np 
from scipy.optimize import curve_fit 
import matplotlib.pyplot as plt 
import math as math 
import scipy.special as sp 

def func(x, sigmag, mu, alpha, c,a): 
    #normal distribution 
    normpdf = (1/(sigmag*np.sqrt(2*math.pi)))*np.exp(-(np.power((x-mu),2)/(2*np.power(sigmag,2)))) 
    normcdf = (0.5*(1+sp.erf((alpha*((x-mu)/sigmag))/(np.sqrt(2))))) 
    return 2*a*normpdf*normcdf + c 

x = np.linspace(0,100,100) 
y = func(x, 10,30, 0,0,1) 
yn = y + 0.001*np.random.normal(size=len(x)) 

results = [] 
for i in xrange(50): 
    p = np.random.randn(5)*10 
    try: 
     popt, pcov = curve_fit(func, x, yn, p) 
    except: 
     pass 
    err = np.sum(np.abs(func(x, *popt) - yn)) 
    results.append((err, popt)) 
    if err < 0.1: 
     break 

err, popt = min(results, key=lambda x:x[0]) 
y_fit= func(x, *popt) 

plt.plot(x,yn) 
plt.plot(x,y_fit) 
print len(results)