Хорошо, вы, кажется, довольно запутались в нескольких вещах. Давайте начнем с самого начала: вы упомянули «многомерную функцию», а затем перейдете к обсуждению обычной гауссовой кривой с одной переменной. Это не многомерная функция: когда вы интегрируете ее, вы интегрируете только одну переменную (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), вы можете легко вычислить интеграл любого гауссиана. Этот интеграл появляется так часто, что имеет специальное имя: error error 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().
Надеюсь, это демистифицировало некоторую путаницу и объяснило, почему функция ошибок является правильным ответом на вашу проблему, а также как это сделать самостоятельно, если вам интересно. Если какое-либо из моих объяснений не было ясным, я предлагаю сначала взглянуть на Википедию; если у вас все еще есть вопросы, не стесняйтесь спрашивать.