Scipy optimize Curve_fit дает различные графики для тех же параметров при подборе пользовательской функции - PullRequest
Новая
       91

Scipy optimize Curve_fit дает различные графики для тех же параметров при подборе пользовательской функции

0 голосов
/ 27 октября 2019

У меня проблема с подгонкой пользовательской функции с использованием scipy.optimize в Python, и я не знаю, почему это происходит. Я генерирую данные по центрированному и нормализованному биномиальному распределению (гауссова кривая), а затем подгоняю кривую. Ожидаемый результат показан на рисунке, когда я излагаю свою функцию над подобранными данными. Но когда я делаю примерку, она терпит неудачу.

Я убежден, что это вещь питоническая, потому что она должна давать параметр a = 1 (именно так я его определяю), и это дает его, но тогда подгонкаплохо (см. рисунок). Однако, если я изменяю сигму на 0,65 * сигма в: 

p_halfg, p_halfg_cov = optimize.curve_fit(lambda x, a:piecewise_half_gauss(x, a, sigma = 0.65*sigma_fit), x, y, p0=[1])

, это дает почти идеальное соответствие (то есть 5/3, как предсказывает математика). Эти припадки должны быть одинаковыми, а они нет! Я даю больше комментариев ниже. Не могли бы вы сказать мне, что происходит и где может быть проблема?

График с a = 1 и sigma = sigma_fit

График с sigma = 0.65* sigma_fit

Я генерирую данные из нормализованного биномиального распределения (я могу предоставить свой код, но значения теперь важнее). Это распределение с N = 10 и p = 0,5, и я центрирую его и беру только правую часть кривой. Затем я подгоняю его к своей функции полугаусса, которая должна иметь то же распределение, что и бином, если ее параметр a = 1 (а сигма равна сигме распределения, sqrt (np (1-p))),Во-первых, проблема в том, что он не соответствует данным, показанным на рисунке, несмотря на получение правильного значения параметра a.

Обратите внимание на странные вещи ... если я установлю sigma = 3 * sigma_fit, я получу= 1/3 и очень плохая посадка (недооценка). Если я установлю его на 0,2 * sigma_fit, я получу также плохую посадку и a = 1 / 0,2 = 5 (завышение). И так далее. Почему? (между прочим, a = 1 / sigma, поэтому процедура подгонки должна работать). 

import numpy as np
import matplotlib.pyplot as plt

import math
pi = math.pi

import scipy.optimize as optimize

# define my function
sigma_fit = 1
def piecewise_half_gauss(x, a, sigma=sigma_fit):
    """Half of normal distribution curve, defined as gaussian centered at 0 with constant value of preexponential factor for x < 0
    Arguments: x values as ndarray whose numbers MUST be float type (use linspace or np.arange(start, end, step, dtype=float),
    a as a parameter of width of the distribution,
    sigma being the deviation, second moment
    Returns: Half gaussian curve

    Ex:
    >>> piecewise_half_gauss(5., 1)
    array(0.04839414)

    >>> x = np.linspace(0,10,11)
    ... piecewise_half_gauss(x, 2, 3)
    array([0.06649038, 0.06557329, 0.0628972 , 0.05867755, 0.05324133,
       0.04698531, 0.04032845, 0.03366645, 0.02733501, 0.02158627,
       0.01657952])

    >>> piecewise_half_gauss(np.arange(0,11,1, dtype=float), 1, 2.4)
    array([1.66225950e-01, 1.52405153e-01, 1.17463281e-01, 7.61037856e-02,
       4.14488078e-02, 1.89766470e-02, 7.30345854e-03, 2.36286717e-03,
       6.42616248e-04, 1.46914868e-04, 2.82345875e-05])

    """
    return np.piecewise(x, [x >= 0, x < 0],
                        [lambda x: np.exp(-x ** 2 / (2 * ((a * sigma) ** 2))) / (np.sqrt(2 * pi) * sigma * a),
                         lambda x: 1 / (np.sqrt(2 * pi) * sigma)])


# Create normalized data for binomial distribution Bin(N,p)
n = 10
p = 0.5

x = np.array([0., 1., 2., 3., 4., 5.])
y = np.array([0.25231325, 0.20657662, 0.11337165, 0.0417071 , 0.01028484,
       0.00170007])

# Get the estimate for sigma parameter
sigma_fit = (n*p*(1-p))**0.5

# Get fitting parameters
p_halfg, p_halfg_cov = optimize.curve_fit(lambda x, a:piecewise_half_gauss(x, a, sigma = sigma_fit), x, y, p0=[1])
print(sigma_fit, p_halfg, p_halfg_cov)


## Plot the result
# unpack fitting parameters
a = np.float64(p_halfg)

# unpack uncertainties in fitting parameters from diagonal of covariance matrix
#da = [np.sqrt(p_halfg_cov[j,j]) for j in range(p_halfg.size)] # if we fit more parameters
da = np.float64(np.sqrt(p_halfg_cov[0]))

# create fitting function from fitted parameters
f_fit = np.linspace(0, 10, 50)
y_fit = piecewise_half_gauss(f_fit, a)

# Create figure window to plot data
fig = plt.figure(1, figsize=(10,10))

plt.scatter(x, y, color = 'r', label = 'Original points')
plt.plot(f_fit, y_fit, label = 'Fit')
plt.xlabel('My x values')
plt.ylabel('My y values')

plt.text(5.8, .25, 'a = {0:0.5f}$\pm${1:0.6f}'.format(a, da))
plt.legend()

Однако, если я нанесу его вручную, он точно подойдет! 

plt.scatter(x, y, c = 'r', label = 'Original points')
plt.plot(np.linspace(0,5,50), piecewise_half_gauss(np.linspace(0,5,50), 1, sigma_fit), label = 'Fit')
plt.legend()
...