Подгонка методом наименьших квадратов кумулятивной функции плотности к некоторым квантилям

Question

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

# a list of (p,x) tuples, where P(X<x)=p
percentiles = [(0.2,8),(0.4,12),(0.5,16),(0.9,30)]

Пользователь указывает семейство распределения (например, Нормальное). Когда имеется более 2 процентилей, система уравнений переопределена, поэтому мы захотим найти параметры, которые лучше всего подходят для ввода, на основе наименьших квадратов.

У меня проблемы с реализацией этого. В этом минимальном примере ниже curve_fit просто возвращает значение по умолчанию 1 для обоих параметров. Что я делаю не так?

from scipy import stats
from scipy import optimize

# a list of (p,x) tuples, where P(X<x)=p
percentiles = [(0.2,8),(0.4,12),(0.5,16),(0.9,30)]

fit = optimize.curve_fit(
            lambda x,mu,sigma: stats.norm(mu,sigma).cdf(x),
            xdata=[x[1] for x in percentiles],
            ydata=[x[0] for x in percentiles])
print(fit[0])

Edit : Johan C указал, что мне важны первоначальные предположения. Вот (грубый) метод, который я использовал, чтобы получить начальное предположение о параметрах, введенных пользователем:

def percentiles_to_list(percentiles):
    out =[]
    i = 1
    c = 1
    for p,q in percentiles:
        if c == len(percentiles):
            number_to_append = int(100 - i)
        else:
            number_to_append = int(p*100-i)
        out += [q]*number_to_append
        i = p*100
        c += 1
    return out

def initial_guess(percentiles):
    lis = percentiles_to_list(percentiles)
    mean = np.mean(lis)
    stdev = np.std(lis)
    return mean,stdev

Answer 1

Для получения точного соответствия важны первоначальные предположения для curve_fit. Здесь данный 50% -ный процентиль является идеальным начальным предположением для mu. Для sigma мы можем принять любое разумное первоначальное предположение. Все предположения по умолчанию - 1, что в данном случае кажется слишком далеким. Помимо (или вместо ) первоначального предположения, также могут быть установлены границы параметров для curve_fit.

Обратите внимание, что из-за функции strict monotoni c мы могли переключить роль из x и y, чтобы проверить, даст ли это более удовлетворительное соответствие.

Следующий код пробует оба соответствия и показывает результат графически. Оба совпадают довольно сильно.

from matplotlib import pyplot as plt
import numpy as np
from scipy import stats
from scipy import optimize

# a list of (p,x) tuples, where P(X<x)=p
percentiles = [(0.2, 8), (0.4, 12), (0.5, 16), (0.9, 30)]

fit_params_ppf, _fit_covariances = optimize.curve_fit(lambda x, mu, sigma: stats.norm(mu, sigma).ppf(x),
                                                      xdata=[percent for percent, percentile in percentiles],
                                                      ydata=[percentile for percent, percentile in percentiles],
                                                      p0=[16, 8])
plt.scatter([percentile for percent, percentile in percentiles], [percent for percent, percentile in percentiles],
            label='given percentiles')
xs = np.linspace(stats.norm(*fit_params).ppf(0.001), stats.norm(*fit_params).ppf(0.999), 500)
plt.plot(xs, stats.norm(*fit_params).cdf(xs), 'b--',
         label=f'fit ppf: $\\mu={fit_params_ppf[0]:.2f}, \\sigma={fit_params_ppf[1]:.2f}$')

fit_params_cdf, _fit_covariances = optimize.curve_fit(lambda x, mu, sigma: stats.norm(mu, sigma).cdf(x),
                                                      xdata=[percentile for percent, percentile in percentiles],
                                                      ydata=[percent for percent, percentile in percentiles],
                                                      p0=[16, 8])
plt.plot(xs, stats.norm(*fit_params).cdf(xs), 'r:',
         label=f'fit cdf: $\\mu={fit_params_cdf[0]:.2f}, \\sigma={fit_params_cdf[1]:.2f}$')

plt.legend()
plt.show()

итоговый сюжет