numpy.random.normal другое распределение: выбор значений из распределения

У меня есть степенное распределение энергий, и я хочу выбрать n случайных энергий на основе распределения. Я пытался сделать это вручную, используя случайные числа, но это слишком неэффективно для того, что я хочу сделать. Мне интересно, есть ли метод в numpy (или другом), который работает как numpy.random.normal, за исключением того, что вместо использования нормального распределения может быть указано распределение. Итак, на мой взгляд, пример может выглядеть так (похоже на numpy.random.normal):

import numpy as np

# Energies from within which I want values drawn
eMin = 50.
eMax = 2500.

# Amount of energies to be drawn
n = 10000

photons = []

for i in range(n):

    # Method that I just made up which would work like random.normal,
    # i.e. return an energy on the distribution based on its probability,
    # but take a distribution other than a normal distribution
    photons.append(np.random.distro(eMin, eMax, lambda e: e**(-1.)))

print(photons)

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

Я не уверен, существует ли такой метод, но кажется, что он должен. Надеюсь, понятно, что я хочу сделать.

РЕДАКТИРОВАТЬ:

Я видел numpy.random.power, но мой показатель степени равен -1, поэтому я не думаю, что это сработает.


python numpy random distribution
person davly    schedule 07.07.2014    source источник
comment
какой именно pdf вам нужен? распределение мощности — это особый случай бета-версии, можете ли вы использовать его вместо docs.scipy.org/doc/numpy/reference/generated/?   -  person wim    schedule 07.07.2014
comment
@wim, я считаю, что мне нужна кусочная функция, которая f (x) = 0 за пределами моего энергетического диапазона и f (x) = x ** a (где a может быть значением от -5 до 5) внутри. Я не вижу, как здесь будет работать бета.   -  person davly    schedule 08.07.2014
comment
@davly обновил мой ответ фрагментом кода на случай, если это будет полезно   -  person John Greenall    schedule 08.07.2014

Ответы (3)


arrow_upward
1
arrow_downward

Выборка из произвольных PDF-файлов на самом деле довольно сложна. Существуют большие и объемные книги о том, как эффективно и точно делать выборки из стандартных семейств дистрибутивов.

Похоже, вы, вероятно, могли бы обойтись с помощью специального метода инверсии для примера, который вы привели.

person Robert Kern    schedule 07.07.2014
comment
Как мне реализовать собственный метод инверсии? Я не вижу такого, как здесь - person davly; 08.07.2014
comment
Выведите обратную функцию CDF. Используйте random_sample(), чтобы получить значения, равномерно распределенные между 0 и 1. Пропустите их через обратный CDF, чтобы получить значения, соответствующие желаемому распределению. В вашем случае обратный CDF равен lambda u: eMin*(eMax/eMin)**u. - person Robert Kern; 09.07.2014

arrow_upward
1
arrow_downward

Если вы хотите сделать выборку из произвольного распределения, вам нужна обратная функция кумулятивной плотности (а не pdf).

Затем вы равномерно выбираете вероятность из диапазона [0,1] и вводите ее в обратную функцию cdf, чтобы получить соответствующее значение.

Часто бывает невозможно получить cdf из pdf аналитически. Однако, если вы хотите аппроксимировать распределение, вы можете сделать это, вычислив f (x) через равные промежутки времени по его области, а затем выполнив общую сумму по этому вектору, чтобы получить приближение cdf и из этого аппроксимировать обратное.

Грубый фрагмент кода:

import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate

def f(x):
   """
   substitute this function with your arbitrary distribution
   must be positive over domain
   """
   return 1/float(x)


#you should vary inputVals to cover the domain of f (for better accurracy you can
#be clever about spacing of values as well). Here i space them logarithmically
#up to 1 then at regular intervals but you could definitely do better
inputVals = np.hstack([1.**np.arange(-1000000,0,100),range(1,10000)])

#everything else should just work
funcVals = np.array([f(x) for x in inputVals])
cdf = np.zeros(len(funcVals))
diff = np.diff(funcVals)
for i in xrange(1,len(funcVals)):
   cdf[i] = cdf[i-1]+funcVals[i-1]*diff[i-1]
cdf /= cdf[-1]

#you could also improve the approximation by choosing appropriate interpolator
inverseCdf = scipy.interpolate.interp1d(cdf,inputVals)

#grab 10k samples from distribution
samples = [inverseCdf(x) for x in np.random.uniform(0,1,size = 100000)]

plt.hist(samples,bins=500)
plt.show()
person John Greenall    schedule 07.07.2014

arrow_upward
-1
arrow_downward

Почему бы вам не использовать eval и не поместить дистрибутив в строку?

>>> cmd = "numpy.random.normal(500)"
>>> eval(cmd)

вы можете манипулировать строкой по своему усмотрению, чтобы установить распределение.

person BigBrownBear00    schedule 07.07.2014
comment
извините, я неправильно понял ваш вопрос. - person BigBrownBear00; 07.07.2014