Накопительная функция распределения с гамма-функцией в Python

Question

Я имею дело с функцией яркости Шехтера, которая выглядит следующим образом:

phi(L)dL = norm. Factor * (L/Lstar)^(a) * exp (L/Lstar) d(L/Lstar)

Скажем, L/Lstar это л.

Аналитическое решение его кумулятивной функции распределения дается гамма-функцией: N = коэффициент нормы * Гамма (a+1, l).

Это неполная гамма-функция, поскольку пределы интегрирования составляют от L до бесконечности.

Теперь я пытаюсь построить cdf в Python. Я использовал:

import scipy.special as ss
si= [ss.gammainc(a+1, l[i]) for i in a]  #cdf  

(где l[I] - массив, который я сделал из случайных чисел)

Итоговый график равен 1 и выглядит как cdf. Но теперь я хочу рандомизировать это. Итак, вместо cdf = 1, я устанавливаю cdf = случайное число (сгенерированное равномерно Python.) Теперь, если я хочу построить гистограмму с числом отсчетов против L со случайной выборкой, мне нужно инвертировать гамма-функцию.

Мой вопрос: как мне инвертировать гамма-функцию в Python?

Вот что у меня сейчас:

u= [random.uniform(0,1) for i in a]

l= [ss.gammaincinv(a+1, u[i]) for i in a]

plt.plot(l, u, '.')

plt.show()

plt.hist(l, bins=50,rwidth= 1.5,histtype='step', lw= 0.7, normed= True, range=(-0.5, 1))

plt.show()

Компилятор не жалуется, но гистограмма неправильной формы. Я думаю, что случайно выбранная гистограмма cdf должна восстановить форму PDF.

Что я делаю неправильно? По-видимому, версия неполной гамма-функции Сципи "регуляризована", что означает, что она делится на полную гамма-функцию. Так что, если я умножу gammainc(a+1, u[I])* gamma(a+1), он все равно не будет работать.

Оси имеют логарифмическое масштабирование.

Какие-либо предложения?

Итог: мне нужно сделать гистограмму cdf функции светимости Шехтера путем случайной выборки.

Answer 1

Первая попытка:

Функция - это отображение домена в диапазон. Таким образом, вы можете написать это так:

def function(x):
    # ...

Domain = list(range(0, 1000)) # [0,1000)
mapping = {}
inverse_mapping = {}
for x in Domain:
    y = function(x)
    mapping[x] = y
    inverse_mapping[y] = x

def inverse_function(y):
    return inverse_mapping[y] # not a continuous function. needs improvement

Если что-то вроде этого у вас на уме, дайте мне знать. Мы можем улучшить его для монотонных функций, таких как cdf.