Хорошо, вы, кажется, довольно запутались в нескольких вещах. Давайте начнем с самого начала: вы упомянули «многомерную функцию», но затем перейдем к обсуждению обычной кривой Гаусса с одной переменной. Это не многомерная функция: при ее интегрировании вы интегрируете только одну переменную (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.
Я никогда не использовал пустую версию констант pi и e; Я всегда придерживался старых версий математических модулей. Я не знаю, почему вы можете предпочесть любой из них.
Я не уверен, что вы собираетесь делать со звонком quad(). quad(gen_gauss, -inf, inf, (10,2,0)) должен интегрировать ренормализованную гауссиану от минус бесконечности до плюс бесконечности и всегда должен выдавать 10 (ваш коэффициент нормализации), поскольку гауссиана интегрируется до 1 по реальной линии. Любой ответ, далекий от 10 (я бы не стал ожидать точно 10, поскольку quad() в конце концов — это только приблизительное значение) означает, что где-то что-то накосячило... трудно сказать, что накосячило, не зная фактического результата ценность и, возможно, внутреннюю работу quad().
Надеюсь, это развеяло некоторую путаницу и объяснило, почему функция ошибок является правильным решением вашей проблемы, а также как сделать все это самостоятельно, если вам любопытно. Если какое-либо из моих объяснений было неясным, я предлагаю сначала взглянуть на Википедию; если у вас еще остались вопросы, не стесняйтесь спрашивать.
person
kquinn
schedule
04.02.2009
