Лучший способ написать функцию Python, которая интегрирует гауссову функцию?

Question

При попытке использовать метод quad scipy для интегрирования гауссовой величины (скажем, существует гауссов метод с именем gauss), у меня возникли проблемы с передачей необходимых параметров в gauss и отказом от Quad выполнять интегрирование по правильной переменной.Есть ли у кого-нибудь хороший пример использования Quad с многомерной функцией?

Но это привело меня к более важному вопросу о том, как вообще лучше всего интегрировать гауссову функцию.Я не нашел гауссовского интегрирования в scipy (к моему удивлению).Мой план состоял в том, чтобы написать простую функцию Гаусса и передать ее в Quad (или, может быть, теперь в интегратор фиксированной ширины).Что бы вы сделали?

Редактировать:Фиксированная ширина означает что-то вроде ловушки, которая использует фиксированное значение dx для расчета площадей под кривой.

На данный момент я пришел к методу make___gauss, который возвращает лямбда-функцию, которая затем может перейти в Quad.Таким образом, я могу создать нормальную функцию со средним значением и дисперсией, которые мне нужны перед интегрированием.

def make_gauss(N, sigma, mu):
    return (lambda x: N/(sigma * (2*numpy.pi)**.5) *
            numpy.e ** (-(x-mu)**2/(2 * sigma**2)))

quad(make_gauss(N=10, sigma=2, mu=0), -inf, inf)

Когда я попытался передать общую функцию Гаусса (которую нужно вызывать с помощью x, N, мю и сигма) и заполнить некоторые значения, используя квадрат, например

quad(gen_gauss, -inf, inf, (10,2,0))

параметры 10, 2 и 0 НЕ обязательно соответствовали N=10, сигма=2, mu=0, что потребовало более расширенного определения.

erf(z) в scipy.special потребовал бы от меня точно определить, что такое t изначально, но приятно знать, что оно есть.

Answer 1

Хорошо, вы, кажется, очень озадачены некоторыми вещами.Начнем с самого начала:вы упомянули «многомерную функцию», но затем переходите к обсуждению обычной гауссовой кривой с одной переменной.Это нет многомерная функция:когда вы интегрируете его, вы интегрируете только одну переменную (x).Это различие важно проводить, поскольку существует является монстр, называемый «многомерным распределением Гаусса», который представляет собой настоящую многомерную функцию и, если ее интегрировать, требует интегрирования по двум или более переменным (при этом используется дорогостоящий метод Монте-Карло, о котором я упоминал ранее).Но вы, похоже, просто говорите об обычном гауссове с одной переменной, с которым гораздо проще работать, интегрировать и все такое.

Распределение Гаусса с одной переменной имеет два параметра: sigma и mu, и является функцией одной переменной, которую мы обозначим x.Похоже, вы также носите с собой параметр нормализации. n (что полезно в паре приложений).Параметры нормализации обычно нет включаются в вычисления, так как их можно просто прикрепить обратно в конце (помните, что интегрирование — это линейный оператор: int(n*f(x), x) = n*int(f(x), x) ).Но мы можем носить его с собой, если хочешь;мне нравится обозначение нормального распределения: тогда

N(x | mu, sigma, n) := (n/(sigma*sqrt(2*pi))) * exp((-(x-mu)^2)/(2*sigma^2))

(читай это как «нормальное распределение x данный sigma, mu, и n дается...") Пока все хорошо;это соответствует вашей функции.Обратите внимание, что единственное истинная переменная вот x:остальные три параметра зафиксированный для любого конкретного гауссиана.

Теперь математический факт:доказуемо верно, что все кривые Гаусса имеют одинаковую форму, они просто немного сдвинуты.Таким образом, мы можем работать с N(x|0,1,1), называемое «стандартным нормальным распределением», и просто переведем наши результаты обратно в общую кривую Гаусса.Итак, если у вас есть интеграл N(x|0,1,1), вы можете тривиально вычислить интеграл любой гауссианы.Этот интеграл встречается настолько часто, что имеет специальное название:тот функция ошибки erf.Из-за некоторых старых соглашений это не так. точно erf;также существует пара аддитивных и мультипликативных факторов.

Если Phi(z) = integral(N(x|0,1,1), -inf, z);то есть, Phi(z) — интеграл стандартного нормального распределения от минус бесконечности до z, то по определению функции ошибок верно, что

Phi(z) = 0.5 + 0.5 * erf(z / sqrt(2)).

Аналогично, если Phi(z | mu, sigma, n) = integral( N(x|sigma, mu, n), -inf, z);то есть, Phi(z | mu, sigma, n) - интеграл нормального распределения при заданных параметрах mu, sigma, и n от минус бесконечности до z, то по определению функции ошибок верно, что

Phi(z | mu, sigma, n) = (n/2) * (1 + erf((x - mu) / (sigma * sqrt(2)))).

Взгляни на статья в Википедии об обычном CDF если вы хотите более подробную информацию или доказательство этого факта.

Ладно, этого должно быть достаточно.Вернемся к вашему (отредактированному) сообщению.Вы говорите: «erf(z) в scipy.special потребует от меня точно определить, что такое t изначально».Я понятия не имею, что вы имеете в виду;где же t (время?) вообще в это ввязываться?Надеемся, что приведенное выше объяснение немного прояснило функцию ошибок, и теперь стало яснее, почему функция ошибки является подходящей функцией для этой работы.

Ваш код Python в порядке, но я бы предпочел замыкание, а не лямбду:

def make_gauss(N, sigma, mu):
    k = N / (sigma * math.sqrt(2*math.pi))
    s = -1.0 / (2 * sigma * sigma)
    def f(x):
        return k * math.exp(s * (x - mu)*(x - mu))
    return f

Использование замыкания позволяет предварительно вычислить константы. k и s, поэтому возвращаемой функции придется выполнять меньше работы при каждом вызове (что может быть важно, если вы ее интегрируете, а это означает, что она будет вызываться много раз).Кроме того, я избегал любого использования оператора возведения в степень **, что медленнее, чем простое возведение в квадрат, и вытащило деление из внутреннего цикла и заменило его умножением.Я вообще не рассматривал их реализацию на Python, но с тех пор, как я в последний раз настраивал внутренний цикл на чистую скорость с использованием сырой сборки x87, я, кажется, помню, что сложение, вычитание или умножение занимают около 4 циклов ЦП каждое, деление примерно 36, а возведение в степень около 200.Это было пару лет назад, так что относитесь к этим цифрам с недоверием;тем не менее, это иллюстрирует их относительную сложность.Так же, расчет exp(x) метод грубой силы — очень плохая идея;есть приемы, которые можно использовать при написании хорошей реализации exp(x) что делает его значительно быстрее и точнее, чем обычный a**b экспоненциализация стиля.

Я никогда не использовал numpy-версию констант pi и e;Я всегда придерживался старых версий математического модуля.Я не знаю, почему вы можете предпочесть любой из них.

Я не уверен, что ты собираешься делать с quad() вызов. quad(gen_gauss, -inf, inf, (10,2,0)) должен интегрировать перенормированную гауссиану от минус бесконечности до плюс бесконечности и всегда должен выдавать 10 (ваш коэффициент нормализации), поскольку гауссиан интегрируется до 1 по реальной линии.Любой ответ, далекий от 10 (я бы не ожидал точно 10 с тех пор quad() это всего лишь приближение, в конце концов) значит что-то где-то напортачено...трудно сказать, что случилось, не зная фактического возвращаемого значения и, возможно, внутренней работы quad().

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

Answer 2

scipy поставляется с «функцией ошибки», также известной как интеграл Гаусса:

import scipy.special
help(scipy.special.erf)

Answer 3

Я предполагаю, что вы имеете дело с многомерными гауссианами;если да, то в SciPy уже есть функция, которую вы ищете:оно называется MVNDIST («Многомерное нормальное РАСПРЕДЕЛЕНИЕ»).Документация SciPy как всегда ужасна, поэтому я даже не могу найти, где запрятана функция, но это где-то там.Документация, безусловно, худшая часть SciPy, и в прошлом она меня бесконечно расстраивала.

Гауссианы с одной переменной просто используют старую добрую функцию ошибок, для которой доступно множество реализаций.

Что касается решения проблемы в целом, да, как упоминает Джеймс Томпсон, вы просто хотите написать свою собственную функцию гауссова распределения и передать ее в quad().Однако если вы можете избежать обобщенной интеграции, это хорошая идея — специализированные методы интеграции для конкретной функции (например, использование MVNDIST) будут намного быстрее, чем стандартная многомерная интеграция Монте-Карло, которая может быть чрезвычайно медленной. для высокой точности.

Answer 4

Распределение Гаусса также называют нормальным распределением.Функция cdf в модуле нормы scipy делает то, что вы хотите.

from scipy.stats import norm
print norm.cdf(0.0)
>>>0.5

http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html#scipy.stats.norm

Answer 5

Почему бы просто не выполнять интегрирование всегда от -бесконечности до +бесконечности, чтобы всегда знать ответ?(шучу!)

Я предполагаю, что единственная причина, по которой в SciPy еще нет стандартной функции Гаусса, заключается в том, что ее написать тривиально.Ваше предложение написать собственную функцию и передать ее в Quad для интеграции звучит превосходно.Для этого он использует общепринятый инструмент SciPy, для вас требуется минимум усилий для написания кода, и он очень удобен для чтения для других людей, даже если они никогда не видели SciPy.

Что именно вы подразумеваете под интегратором фиксированной ширины?Вы имеете в виду использование другого алгоритма, чем тот, который использует QUADPACK?

Редактировать:Для полноты картины вот что-то вроде того, что я бы попробовал для гауссовой функции со средним значением 0 и стандартным отклонением 1 от 0 до + бесконечности:

from scipy.integrate import quad
from math import pi, exp
mean = 0
sd   = 1
quad(lambda x: 1 / ( sd * ( 2 * pi ) ** 0.5 ) * exp( x ** 2 / (-2 * sd ** 2) ), 0, inf )

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