Задача подбора гауссиана с отрицательной амплитудой с помощью astropy.modeling

Question

Я пытаюсь подогнать гауссиан к дискретному потенциалу, используя пакет astropy.modeling. Хотя я назначаю отрицательную амплитуду гауссову, он возвращает нулевой гауссиан, то есть везде с нулевой амплитудой:

plot(c[0], pot_x, label='Discrete potential')
plot(c[0], g(c[0]), label='Gaussian fit') 
legend() 

У меня есть следующие строки кода для выполнения подгонки:

g_init = models.Gaussian1D(amplitude=-1., mean=0, stddev=1.)
fit_g = fitting.LevMarLSQFitter()
g = fit_g(g_init, c[0], pot_x)

Где

c[0] = array([13.31381488, 13.31944489, 13.32507491, 13.33070493, 13.33633494,
   13.34196496, 13.34759498, 13.35322499, 13.35885501, 13.36448503,
   13.37011504, 13.37574506, 13.38137507, 13.38700509, 13.39263511,
   13.39826512, 13.40389514, 13.40952516, 13.41515517, 13.42078519])

pot_x = array([ -1.72620157,  -3.71811187,  -6.01282809,  -6.98874144,
    -8.36645166, -14.31787771, -23.3688849 , -26.14679496,
   -18.85970983, -10.73888697,  -7.10763373,  -5.81176637,
    -5.44146953,  -5.37165105,  -4.6454408 ,  -2.90307138,
    -1.66250349,  -1.66096343,  -1.8188269 ,  -1.41980552])

У кого-нибудь есть идеи, в чем может быть проблема?

Решено : Мне просто нужно было назначить среднее значение в диапазоне домена, например 13,35.

Answer 1

Поскольку я не знаком с Astropy, я использовал scipy. Код ниже обеспечивает следующий результат:

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

x = np.asarray([13.31381488, 13.31944489, 13.32507491, 13.33070493, 13.33633494,
   13.34196496, 13.34759498, 13.35322499, 13.35885501, 13.36448503,
   13.37011504, 13.37574506, 13.38137507, 13.38700509, 13.39263511,
   13.39826512, 13.40389514, 13.40952516, 13.41515517, 13.42078519])

y = -np.asarray([ -1.72620157,  -3.71811187,  -6.01282809,  -6.98874144,
    -8.36645166, -14.31787771, -23.3688849 , -26.14679496,
   -18.85970983, -10.73888697,  -7.10763373,  -5.81176637,
    -5.44146953,  -5.37165105,  -4.6454408 ,  -2.90307138,
    -1.66250349,  -1.66096343,  -1.8188269 ,  -1.41980552])

mean = sum(x * y) / sum(y)
sigma = np.sqrt(sum(y * (x - mean)**2) / sum(y))

def Gauss(x, a, x0, sigma):
    return a * np.exp(-(x - x0)**2 / (2 * sigma**2))

popt,pcov = curve_fit(Gauss, x, y, p0=[max(y), mean, sigma])

plt.plot(x, y, 'b+:', label='data')
plt.plot(x, Gauss(x, *popt), 'r-', label='fit')
plt.legend()

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