Как представить 1D-вектор как сумму кривых Гаусса с помощью scipy/numpy?

UPD: Спасибо, все работает.

У меня есть 1D-вектор, который представляет собой гистограмму. Это выглядит как сумма нескольких гауссовых функций: введите здесь описание изображения

Я нашел curve_fit образец кода на SO, но не знаю, как его изменить, чтобы получить больше гауссовских кортежей (mu, sigma). Я слышал, что «curve_fit» оптимизирует только одну функцию (в данном случае одну гауссову кривую).

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

def estimate_sigma(hist):
    bin_edges = np.arange(len(hist))
    bin_centres = bin_edges + 0.5

    # Define model function to be used to fit to the data above:
    def gauss(x, *p):
        A, mu, sigma = p
        return A*numpy.exp(-(x-mu)**2/(2.*sigma**2))

    # p0 is the initial guess for the fitting coefficients (A, mu and sigma above)
    p0 = [1., 0., 1.]

    coeff, var_matrix = curve_fit(gauss, bin_centres, hist, p0=p0)

    # Get the fitted curve
    hist_fit = gauss(bin_centres, *coeff)

    plt.plot(bin_centres, hist, label='Test data')
    plt.plot(bin_centres, hist_fit, label='Fitted data')

    print 'Fitted mean = ', coeff[1]
    coeff2 =coeff[2]
    print 'Fitted standard deviation = ', coeff2

    plt.show()

Эта функция находит одну гауссову кривую, тогда как визуально их 3 или 4: введите здесь описание изображения

Пожалуйста, не могли бы вы посоветовать некоторые функции numpy/scipy для достижения gmm-представления 1D vector в форме ([m1, sigma1],[m2, sigma2],..,[mN,sigmaN])?


python numpy scipy gmm
person Olha    schedule 14.12.2016    source источник
comment
Вы можете добавить вторую функцию, которая правильно складывает несколько gauss(x, *p), например. gauss(x, *p[0:2]) + gauss(x, *p[2:4]) + ... и подходит к этому. (С дополнительным параметром амплитуды, возможно, для каждого термина.)   -  person tBuLi    schedule 14.12.2016
comment
@tBuLi Спасибо, работает!!!   -  person Olha    schedule 14.12.2016
comment
@lanery Да, код почти такой же: stackoverflow.com/a/26954881/2567725. Поэтому я отмечаю его как дубликат.   -  person Olha    schedule 15.12.2016

Ответы (1)


arrow_upward
0
arrow_downward

В соответствии с рекомендациями tBuLi я передал дополнительные коэффициенты гауссовых кривых в gauss, а также в curve_fit. Теперь подобранная кривая выглядит так: введите здесь описание изображения

Обновленный код:

def estimate_sigma(hist):
    bin_edges = np.arange(len(hist))
    bin_centres = bin_edges + 0.5

    # Define model function to be used to fit to the data above:
    def gauss(x, *gparams):
        g_count = len(gparams)/3
        def gauss_impl(x, A, mu, sigma):
            return A*numpy.exp(-(x-mu)**2/(2.*sigma**2))
        res = np.zeros(len(x))
        for gi in range(g_count):
            res += gauss_impl(x, gparams[gi*3], gparams[gi*3+1], gparams[gi*3+2])
        return res

    # p0 is the initial guess for the fitting coefficients (A, mu and sigma above)
    curves_count = 4
    p0 = np.tile([1., 0., 1.], curves_count)

    coeff, var_matrix = curve_fit(gauss, bin_centres, hist, p0=p0)

    # Get the fitted curve
    hist_fit = gauss(bin_centres, *coeff)

    plt.plot(bin_centres, hist, label='Test data')
    plt.plot(bin_centres, hist_fit, label='Fitted data')

    # Finally, lets get the fitting parameters, i.e. the mean and standard deviation:
    print coeff

    plt.show()
person Olha    schedule 14.12.2016
comment
Выглядит неплохо! Лично я бы поместил определение gauss_impl вне определения gauss, но кроме этого прекрасного решения. Рад, что смог помочь :). Вам разрешено отвечать на свои вопросы, так что вы можете принять это как ответ. - person tBuLi; 15.12.2016