Блог

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

В конечном счете, я хочу имитировать шум, но я пытаюсь понять, почему, если у меня есть два гауссиана со средними значениями, которые действительно далеки друг от друга, мой curve_fit должен возвращать их среднее значение. Он этого не делает.

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

N_tot = 1000
# Draw from the major gaussian. Note the number N. It is
# the main parameter in obtaining your estimators.
mean = 0; sigma = 1; var = sigma**2; N = 100
A = 1/np.sqrt((2*np.pi*var))
points = gauss.draw_1dGauss(mean,var,N)

# Now draw from a minor gaussian. Note Np
meanp = 10; sigmap = 1; varp = sigmap**2; Np = N_tot-N
pointsp = gauss.draw_1dGauss(meanp,varp,Np)
Ap = 1/np.sqrt((2*np.pi*varp))      

# Now implement the sum of the draws by concatenating the two arrays.
points_tot = np.array(points.tolist() pointsp.tolist())
bins_tot = len(points_tot)/5
hist_tot, bin_edges_tot = np.histogram(points_tot,bins_tot,density=True)
bin_centres_tot = (bin_edges_tot[:-1]   bin_edges_tot[1:])/2.0

# Initial guess
p0 = [A, mean, sigma]

# Result of the fit
coeff, var_matrix = curve_fit(gauss.gaussFun, bin_centres_tot, hist_tot, p0=p0)

# Get the fitted curve
hist_fit = gauss.gaussFun(bin_centres, *coeff)
plt.figure(5); plt.title('Gaussian Estimate')
plt.suptitle('Gaussian Parameters: Mu = '  str(coeff[1])  ' , Sigma = '   str(coeff[2])   ', Amplitude = '   str(coeff[0]))
plt.plot(bin_centres,hist_fit)
plt.draw()        

# Error on the estimates
error_parameters = np.sqrt(np.array([var_matrix[0][0],var_matrix[1][1],var_matrix[2][2]]))
  

Возвращаемые параметры по-прежнему сосредоточены около 0, и я не уверен, почему. Он должен быть сосредоточен вокруг 10.

Редактировать: изменены части целочисленного деления, но по-прежнему не возвращается подходящее значение. Я должен получить среднее значение около ~ 10, поскольку большинство моих точек берутся из этого дистрибутива (т. Е. Второстепенного дистрибутива)

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

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

Когда распределение является таким однобоким (90% выборок относятся к большему из двух пиков), условия ошибки на главном пике уничтожают локальные минимумы на меньшем пике и минимуме между пиками.

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

С помощью следующих фрагментов я могу запустить ваш код:

 bin_centres = bin_centres_tot

def draw_1dGauss(mean,var,N):
    from scipy.stats import norm
    from numpy import sqrt
    return scipy.stats.norm.rvs(loc = mean, scale = sqrt(var), size=N)

def gaussFun(bin_centres, *coeff):
    from numpy import sqrt, exp, pi
    A, mean, sigma = coeff[0], coeff[1], coeff[2]
    return exp(-(bin_centres-mean)**2 / 2. / sigma**2 ) / sigma / sqrt(2*pi)

plt.hist(points_tot, normed=True, bins=40)