Блог

scipy.stats.rv_continuous.fit, позволяет вам фиксировать параметры при подгонке распределения, но это зависит от выбора scipy параметризации. Для гамма-распределения используется параметризация k, theta (shape, scale), поэтому его было бы легко подогнать, например, сохраняя тета-константу. Я хочу соответствовать набору данных, в котором я знаю среднее значение, но наблюдаемое среднее значение может отличаться из-за ошибки выборки. Это было бы легко, если бы scipy использовал параметризацию, которая использует mu = k * theta вместо theta. Есть ли способ заставить scipy сделать это? А если нет, есть ли другая библиотека, которая может?

Вот некоторый пример кода с набором данных, который имеет наблюдаемое среднее значение 9,952, но я знаю, что фактическое среднее значение базового распределения равно 11:

 from scipy.stats import gamma
observations = [17.6, 24.9, 3.9, 17.6, 11.8, 10.4, 4.1, 11.7, 5.7, 1.6, 
                8.6, 12.9, 5.7, 8.0, 7.4, 1.2, 11.3, 10.4, 1.0, 1.9, 
                6.0, 9.3, 13.3, 5.4, 9.1, 4.0, 12.8, 11.1, 23.1, 4.2, 
                7.9, 11.1, 10.0, 3.4, 27.8, 7.2, 14.9, 2.9, 5.5, 7.0, 
                3.9, 12.3, 10.6, 22.1, 5.0, 4.1, 21.3, 15.9, 34.5, 8.1, 
                19.6, 10.8, 13.4, 22.8, 27.6, 6.8, 5.9, 9.0, 7.1, 21.2, 
                1.0, 14.6, 16.9, 1.0, 6.5, 2.9, 7.1, 14.1, 15.2, 7.8, 
                9.0, 4.9, 2.1, 9.5, 5.6, 11.1, 7.7, 18.3, 3.8, 11.0, 
                4.2, 12.5, 8.4, 3.2, 4.0, 3.8, 2.0, 24.7, 24.6, 3.4, 
                4.3, 3.2, 7.6, 8.3, 14.5, 8.3, 8.4, 14.0, 1.0, 9.0]
shape, _, scale = gamma.fit(observations, floc = 0)
print(shape*scale)
  

и это дает

 9.952
  

но я хотел бы подгонку так, чтобы shape*scale = 11.0

fit Метод распределений SciPy обеспечивает оценку максимального правдоподобия параметров. Вы правы, что оно обеспечивает только подгонку формы, местоположения и масштаба. (На самом деле, вы сказали форму и масштаб, но SciPy также включает параметр местоположения. Иногда это называется гамма-распределением с тремя параметрами.)

Для большинства распределений в SciPy fit метод использует числовой оптимизатор для минимизации отрицательного логарифмического правдоподобия, как определено в nnlf методе. Вместо того, чтобы использовать fit метод, вы могли бы сделать это самостоятельно, используя всего пару строк кода. Это позволяет вам создать целевую функцию только с одним параметром, скажем, с формой k , и внутри этой функции установить theta = mean/k , где mean желаемое среднее значение, и вызвать gamma.nnlf для оценки отрицательного логарифмического правдоподобия. Вот один из способов, которым вы могли бы это сделать:

 import numpy as np
from scipy.stats import gamma
from scipy.optimize import fmin


def nll(k, mean, x):
    return gamma.nnlf(np.array([k[0], 0, mean/k[0]]), x)


observations = [17.6, 24.9, 3.9, 17.6, 11.8, 10.4, 4.1, 11.7, 5.7, 1.6,
                8.6, 12.9, 5.7, 8.0, 7.4, 1.2, 11.3, 10.4, 1.0, 1.9,
                6.0, 9.3, 13.3, 5.4, 9.1, 4.0, 12.8, 11.1, 23.1, 4.2,
                7.9, 11.1, 10.0, 3.4, 27.8, 7.2, 14.9, 2.9, 5.5, 7.0,
                3.9, 12.3, 10.6, 22.1, 5.0, 4.1, 21.3, 15.9, 34.5, 8.1,
                19.6, 10.8, 13.4, 22.8, 27.6, 6.8, 5.9, 9.0, 7.1, 21.2,
                1.0, 14.6, 16.9, 1.0, 6.5, 2.9, 7.1, 14.1, 15.2, 7.8,
                9.0, 4.9, 2.1, 9.5, 5.6, 11.1, 7.7, 18.3, 3.8, 11.0,
                4.2, 12.5, 8.4, 3.2, 4.0, 3.8, 2.0, 24.7, 24.6, 3.4,
                4.3, 3.2, 7.6, 8.3, 14.5, 8.3, 8.4, 14.0, 1.0, 9.0]


# This is the desired mean of the distribution.
mean = 11

# Initial guess for the shape parameter.
k0 = 3.0
opt = fmin(nll, k0, args=(mean, np.array(observations)),
           xtol=1e-11, disp=False)

k_opt = opt[0]
theta_opt = mean / k_opt
print(f"k_opt:     {k_opt:9.7f}")
print(f"theta_opt: {theta_opt:9.7f}")
  

Этот скрипт печатает

 k_opt:     1.9712604
theta_opt: 5.5801861
  

В качестве альтернативы, можно изменить условия первого порядка для экстремума логарифмического правдоподобия, показанного в википедии, так, чтобы был только один параметр, k . Тогда условие для экстремального значения может быть реализовано как скалярное уравнение, корень которого можно найти, скажем, scipy.optimize.fsolve с помощью . Ниже приведен вариант вышеупомянутого скрипта, который использует эту технику.

 import numpy as np
from scipy.special import digamma
from scipy.optimize import fsolve


def first_order_eq(k, mean, x):
    mean_logx = np.mean(np.log(x))
    return (np.log(k) - digamma(k)   mean_logx - np.mean(x)/mean
            - np.log(mean)   1)


observations = [17.6, 24.9, 3.9, 17.6, 11.8, 10.4, 4.1, 11.7, 5.7, 1.6,
                8.6, 12.9, 5.7, 8.0, 7.4, 1.2, 11.3, 10.4, 1.0, 1.9,
                6.0, 9.3, 13.3, 5.4, 9.1, 4.0, 12.8, 11.1, 23.1, 4.2,
                7.9, 11.1, 10.0, 3.4, 27.8, 7.2, 14.9, 2.9, 5.5, 7.0,
                3.9, 12.3, 10.6, 22.1, 5.0, 4.1, 21.3, 15.9, 34.5, 8.1,
                19.6, 10.8, 13.4, 22.8, 27.6, 6.8, 5.9, 9.0, 7.1, 21.2,
                1.0, 14.6, 16.9, 1.0, 6.5, 2.9, 7.1, 14.1, 15.2, 7.8,
                9.0, 4.9, 2.1, 9.5, 5.6, 11.1, 7.7, 18.3, 3.8, 11.0,
                4.2, 12.5, 8.4, 3.2, 4.0, 3.8, 2.0, 24.7, 24.6, 3.4,
                4.3, 3.2, 7.6, 8.3, 14.5, 8.3, 8.4, 14.0, 1.0, 9.0]


# This is the desired mean of the distribution.
mean = 11

# Initial guess for the shape parameter.
k0 = 3
sol = fsolve(first_order_eq, k0, args=(mean, observations),
             xtol=1e-11)

k_opt = sol[0]
theta_opt = mean / k_opt
print(f"k_opt:     {k_opt:9.7f}")
print(f"theta_opt: {theta_opt:9.7f}")
  

Вывод:

 k_opt:     1.9712604
theta_opt: 5.5801861