Подбор распределения по гистограмме с использованием scipy

Я могу что-то неправильно понять, но я считаю, что подгонка к гистограмме — это именно то, что вам следует делать: вы пытаетесь аппроксимировать плотность вероятности. И гистограмма максимально приближена к базовой плотности вероятности. Вам просто нужно нормализовать его, чтобы интеграл был равен 1, или разрешить вашей подобранной модели содержать произвольный префактор.

import numpy as np
import scipy.stats as stats
import scipy.optimize as opt
import matplotlib.pyplot as plt

orig_hist = np.array([10, 5, 3, 2, 1])
norm_hist = orig_hist/float(sum(orig_hist))

popt,pcov = opt.curve_fit(lambda x,c: stats.weibull_min.pdf(x,c), np.arange(len(norm_hist)),norm_hist)

plt.figure()
plt.plot(norm_hist,'o-',label='norm_hist')
plt.plot(stats.weibull_min.pdf(np.arange(len(norm_hist)),popt),'s-',label='Weibull_min fit')
plt.legend()

Конечно, для вашего ввода подгонка Вейбулла будет далеко не удовлетворительной:

Обновлять

Как я упоминал выше, Weibull_min плохо подходит для вашего примера входных данных. Большая проблема заключается в том, что это также плохо соответствует вашим фактическим данным:

orig_hist = np.array([ 23., 14., 13., 12., 12., 12., 11., 11., 11., 11., 10., 10., 10., 9., 9., 8., 8., 8., 8., 8., 8., 8., 8., 8., 8., 8., 7., 7., 7., 7., 7., 7., 7., 7., 7., 7., 7., 7., 7., 6., 6., 6., 6., 6., 6., 6., 6., 6., 6., 6.], dtype=np.float32)

У этой гистограммы есть две основные проблемы. Первый, как я уже сказал, заключается в том, что он вряд ли соответствует распределению Weibull_min: оно максимально вблизи нуля и имеет длинный хвост, поэтому ему нужна нетривиальная комбинация параметров Вейбулла. Кроме того, ваша гистограмма явно содержит только часть распределения. Это означает, что мое предложение по нормализации, приведенное выше, гарантированно потерпит неудачу. Вы не можете избежать использования произвольного параметра масштаба в вашей подгонке.

Я вручную определил масштабированную функцию подбора Вейбулла в соответствии с формулой из Википедии:

my_weibull = lambda x,l,c,A: A*float(c)/l*(x/float(l))**(c-1)*np.exp(-(x/float(l))**c)

В этой функции x — независимая переменная, l — lambda (параметр масштаба), c — k (параметр формы), а A — префактор масштабирования. Слабым преимуществом введения A является то, что вам не нужно нормализовать гистограмму.

Теперь, когда я поместил эту функцию в scipy.optimize.curve_fit, я обнаружил то же, что и вы: она на самом деле не выполняет подгонку, а придерживается исходных параметров подгонки, какие бы вы ни установили (используя параметр p0; все догадки по умолчанию равны 1 для каждого параметр). И curve_fit, похоже, считает, что примерка сходится.

После более чем часового биения головой о стену я понял, что проблема в том, что сингулярное поведение на x=0 отбрасывает нелинейный алгоритм наименьших квадратов. Исключив самую первую точку данных, вы получите фактическое соответствие вашим данным. Я подозреваю, что если мы установим c=1 и не позволим этому подходить, то эта проблема может исчезнуть, но, вероятно, будет более информативно разрешить это подгонку (поэтому я не проверял).

Вот соответствующий код:

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

orig_hist = np.array([ 23., 14., 13., 12., 12., 12., 11., 11., 11., 11., 10., 10., 10., 9., 9., 8., 8., 8., 8., 8., 8., 8., 8., 8., 8., 8., 7., 7., 7., 7., 7., 7., 7., 7., 7., 7., 7., 7., 7., 6., 6., 6., 6., 6., 6., 6., 6., 6., 6., 6.], dtype=np.float32)

my_weibull = lambda x,l,c,A: A*float(c)/l*(x/float(l))**(c-1)*np.exp(-(x/float(l))**c)

popt,pcov = opt.curve_fit(my_weibull,np.arange(len(orig_hist))[1:],orig_hist[1:]) #throw away x=0!

plt.figure()
plt.plot(np.arange(len(orig_hist)),orig_hist,'o-',label='orig_hist')
plt.plot(np.arange(len(orig_hist)),my_weibull(np.arange(len(orig_hist)),*popt),'s-',label='Scaled Weibull fit')
plt.legend()

Результат:

In [631]: popt
Out[631]: array([  1.10511850e+02,   8.82327822e-01,   1.05206207e+03])

окончательные подобранные параметры находятся в порядке (l,c,A) с параметром формы около 0.88. Это соответствует расходящейся плотности вероятности, что объясняет, почему появляется несколько ошибок, говорящих

RuntimeWarning: в мощности обнаружено недопустимое значение

и почему нет точки данных из фитинга для x=0. Но, судя по визуальному совпадению данных и подгонки, можно оценить, приемлем результат или нет.

Если вы хотите переборщить, вы, вероятно, можете попробовать сгенерировать точки, используя np.random.weibull с этими параметрами, а затем сравнить полученные гистограммы со своими.