Подгонка данных с использованием scipy truncnorm

Question

У меня есть данные, которые следуют за распределением Гаусса. Однако данные действительно гауссовы только для диапазона значений [xa, xb], поэтому я хочу подогнать усеченное нормальное распределение, используя scipy.stats.truncnorm , используя тот факт, что я знаю диапазон [xa , хь]. Моя цель - найти местоположение и масштаб.

Я не понимаю, как исправить XA и XB в посадке. Параметры формы - это «a» и «b», но они зависят от местоположения и масштаба, которые мне неизвестны. Более того, кажется невозможным сделать начальное предположение о «a» и «b» (они могут быть заморожены только с помощью fa и fb?). Когда я делаю:

par = truncnorm.fit(r, a=a_guess, b=b_guess, scale= scale_guess, loc = loc_guess)

Я получаю

Неизвестные аргументы: {'a': 0.0, 'b': 2.4444444444444446}.

Кроме того, припадки, которые я получаю, очень нестабильны. Вот пример:

from scipy.stats import truncnorm
import matplotlib.pyplot as plt

xa, xb = 30,250 
loc, loc_guess = 50, 30
scale, scale_guess = 75, 90
a,b = (xa-loc)/scale, (xb-loc)/scale

fig, ax = plt.subplots(1, 1)
x = np.linspace(xa,xb,10000)    
ax.plot(x, truncnorm.pdf(x, a, b, loc=loc, scale=scale),
        'r-', lw=5, alpha=0.6, label='truncnorm pdf')

r = truncnorm.rvs(a, b, loc=loc, scale=scale, size=10000)
par = truncnorm.fit(r, scale= scale_guess, loc = loc_guess)
ax.plot(x, truncnorm.pdf(x, *par),
        'b-', lw=1, alpha=0.6, label='truncnorm fit')
ax.hist(r, density=True, histtype='stepfilled', alpha=0.3)
plt.legend()
plt.show()

1-й пример 2-й пример

У меня также часто есть это предупреждение:

/ home / elie / anaconda2 / envs / py36 / lib / python3.6 / site-packages / scipy / stats / _continuous_distns.py: 5823: RuntimeWarning: деление на ноль, встречающееся в журнале self._logdelta = np.log (self._delta)

Answer 1

Как вы обнаружили, проблема в том, что параметры, которые вы хотите сохранить фиксированными, xa и xb, не являются собственными параметрами truncnorm. truncnorm имеет параметры формы a и b, которые определяют форму путем установки интервала x для нормального распределения standard . Затем эта форма смещается и масштабируется параметрами loc и scale. Отношение

xa = a*scale + loc
xb = b*scale + loc

Чтобы исправить xa и xb, вы можете использовать один из минимизаторов SciPy, который принимает ограничения равенства. Здесь я буду использовать scipy.optimize.fmin_slsqp. (Вместо этого вы можете использовать функцию «omnibus» scipy.optmize.minimize, которая включает решатель SLSQP в качестве одной из своих опций.)

Вот скрипт, который демонстрирует, как использовать fmin_slsqp для этой проблемы. Функция func является целевой функцией, которую нужно минимизировать. Это просто оболочка для truncnorm.nnlf, функция отрицательного логарифмического правдоподобия. Функция constraint возвращает массив, содержащий два значения. Эти значения равны 0, когда ограничение выполнено.

import numpy as np
from scipy.stats import truncnorm
from scipy.optimize import fmin_slsqp

import matplotlib.pyplot as plt


def func(p, r, xa, xb):
    return truncnorm.nnlf(p, r)


def constraint(p, r, xa, xb):
    a, b, loc, scale = p
    return np.array([a*scale + loc - xa, b*scale + loc - xb])


xa, xb = 30, 250 
loc = 50
scale = 75

a = (xa - loc)/scale
b = (xb - loc)/scale

# Generate some data to work with.
r = truncnorm.rvs(a, b, loc=loc, scale=scale, size=10000)

loc_guess = 30
scale_guess = 90
a_guess = (xa - loc_guess)/scale_guess
b_guess = (xb - loc_guess)/scale_guess
p0 = [a_guess, b_guess, loc_guess, scale_guess]

par = fmin_slsqp(func, p0, f_eqcons=constraint, args=(r, xa, xb),
                 iprint=False, iter=1000)

xmin = 0
xmax = 300
x = np.linspace(xmin, xmax, 1000)

fig, ax = plt.subplots(1, 1)
ax.plot(x, truncnorm.pdf(x, a, b, loc=loc, scale=scale),
        'r-', lw=3, alpha=0.4, label='truncnorm pdf')
ax.plot(x, truncnorm.pdf(x, *par),
        'k--', lw=1, alpha=1.0, label='truncnorm fit')
ax.hist(r, bins=15, density=True, histtype='stepfilled', alpha=0.3)
ax.legend(shadow=True)
plt.xlim(xmin, xmax)
plt.grid(True)

plt.show()

Вот сюжет, который он генерирует. Данные выборки случайны, поэтому график будет отличаться при каждом запуске.

Примечание: иногда генерируется случайный набор данных, для которого fmin_slsqp завершается с ошибкой с "обнаруженным недопустимым значением" во время вычисления. Я не исследовал это дальше, но вы можете столкнуться с этим с вашими данными.