ガウスを統合するPythonの関数を記述するための最良の方法?
質問
(ガウスという名前のガウスの方法があると言うことができます)ガウスを統合するためにscipyのダウンロードのクワッドメソッドを使用しようと、私は問題ガウスするために必要なパラメータを渡して、正しい変数の上に統合を行うためにクワッドを残しを持っていました。誰もが多次元関数/ wのクワッドを使用する方法の良い例がありますか?
しかし、これは、一般的にガウスを統合するための最良の方法についての詳細壮大な質問に私を導きました。私は(私の驚きに)scipyのダウンロードに統合ガウスを見つけることができませんでした。私の計画は、単純なガウス関数を記述し、クワッド(または多分今固定幅インテグレータ)に渡すことでした。あなたならどうしますか?
編集:曲線下の面積を計算するために固定されたDXを使用していますtrapzのようなものを意味固定幅
。私はこれまでのところに来ていることは、その後のクワッドに行くことができますラムダ関数を返すメソッドmake___gaussです。この方法で私は統合前に必要な平均と分散を持つ正規の関数を作ることができます。
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、μ= 0、。
scipy.specialでERF(z)はtは最初は正確に定義するために私を必要とするが、それは素敵なことがあると知ってます。
解決
さて、あなたはいくつかのことについてかなり混乱しているようです。最初に起動してみましょう:あなたは、「多次元機能」を述べたが、その後通常の1変数のガウス曲線を議論するために行きます。あなたがそれを統合するとき、あなたは唯一の変数(X)を統合:これはのないの多次元関数です。 が存在しているため区別は、作ることが重要である。のモンスターは、真の多次元関数である「多変量ガウス分布」と呼ばれ、統合された場合、高価なモンテを使用する(二つ以上の変数の上に統合する必要が私は前に述べたモンテカルロ法)。しかし、あなたはただ、で動作する統合する方がはるかに簡単です通常の1変数のガウス、について話しているように見える、そしてすべてのこと。
1変数のガウス分布は、2つのパラメータ、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 の上のWikipediaの記事をご覧くださいあるいは、この事実の証明ます。
さて、それは十分な背景説明をする必要があります。戻るあなた(編集済み)ポストへ。あなたは「scipy.specialでERF(z)はtは最初は正確に定義するために私を必要とする」と言います。私はあなたがこのことで何を意味するのか見当がつかない。どこt
は(時間?)で、このすべてに入るのでしょうか?うまくいけば、上記の説明では、誤差関数を少し詳解しており、それは誤差関数は、仕事のための適切な機能である理由を、今より明確です。
あなたのPythonコードはOKですが、私はラムダ以上の閉鎖を希望ます:
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
の事前計算を可能にするので、返された関数ウィルlはあまり仕事に(あなたはそれが何度も呼び出されますを意味し、それを統合している場合、重要であり得る)が呼び出されるたびに行う必要があります。また、私はちょうど二乗を書くよりも遅い指数演算子の**
、のいずれかの使用を避け、そして内側のループの外に除算を掲揚し、乗算とそれを交換しました。私は、Pythonでの実装では全く見ていないが、生のx87アセンブリを使用して、純粋なスピードのための内部ループをチューニング私の最後の時から、私はそれが追加されます覚えているようだ、減算、または乗算は約4 CPUサイクルごとを取る、およそ分割します36、そして数年前だったので、塩の粒とこれらの数字を取る約200累乗。まだ、それは彼らの相対的な複雑さを示しています。同様に、exp(x)
ブルートフォースの道を計算することは非常に悪い考えです。一般exp(x)
スタイルの累乗よりも、それは非常に速く、より正確にa**b
の良い実装を書くときに、あなたが取ることができるトリックがあります。
私は、定数PIとEのnumpyのバージョンを使用したことがありません。私はいつも、昔ながらの数学モジュールのバージョンで立ち往生してきました。あなたはどちらを好むかもしれない理由を私は知らない。
私はあなたがquad()
呼び出しとするために行っているかわからないんだけど。 quad(gen_gauss, -inf, inf, (10,2,0))
はプラス無限大にマイナス無限大から再正規化ガウスを統合するべき、とガウスは、実際のラインの上に1に統合するので、常に、(あなたの正規化因子)10を吐き出す必要があります。これまで10から任意の答えは(quad()
は、すべての後、唯一の近似値であるので、私はの正確の10期待していない)、実際のリターンを知らずに台無しに何を言うのは難しい...何かがどこかにめちゃくちゃにされることを意味します値とquad()
の可能性内部動作
うまくいけば、それは混乱の一部を詳解し、誤差関数は、あなたが興味があれば、すべてを自分で行う方法と同様に、あなたの問題に正解である理由を説明しています。私の説明のいずれかが明確ではなかった場合は、私が最初にウィキペディアをざっと見てみる提案します。あなたはまだ質問がある場合は、お気軽にお問い合わせください。
他のヒント
"誤差関数" とscipyのダウンロード船、通称ガウスの積分ます:
import scipy.special
help(scipy.special.erf)
私はあなたが、多変量ガウス分布を処理していると仮定します。もしそうなら、scipyのダウンロードはすでにあなたが探している機能を持っています。それはMVNDIST(「多変量正規分布)と呼ばれていますscipyのダウンロードのドキュメントでは、これまで通り、ひどいですので、私も機能が埋葬された場所を見つけることが、それはどこかにそこにありますA>。ドキュメントを簡単にscipyのダウンロードの最悪の一部であり、過去には終わりに私をイライラしています。
シングル変数ガウス分布は、ちょうど多くの実装が利用可能であるかの、古き良き誤差関数を使用します。
ジェームス・トンプソンが言及として一般的な問題を攻撃に関しては、はい、あなただけの()独自のガウス分布関数を書き、クワッドにそれを送りたいです。あなたは一般的な統合を避けることができる場合は、しかし、そうすることをお勧めします - 特定の機能に特化した統合技術は(MVNDISTが使用するように)極端に遅くなることが標準モンテカルロ多次元統合、よりもはるかに高速であることを行っています高精度のために。
ガウス分布も正規分布と呼ばれます。 scipyのダウンロード規範モジュール内のCDF関数が何をしたいんます。
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する
なぜ、常に、-Infinityから+無限大にあなたの統合をしませんか? (冗談!)
私の推測では、缶詰のガウス関数は、scipyのダウンロードには、既に存在しないことを唯一の理由は、それは書くための些細な機能だということであるということです。独自の関数を記述し、優れたサウンドを統合するクワッドに渡すについてあなたの提案。それは、これを行うための受け入れscipyのダウンロードツールを使用して、それはあなたのための最小限のコードの努力だし、それは彼らがscipyのダウンロードを見たことがない場合でも、他の人のために非常に読みやすいです。
正確にあなたが固定幅の積分器とはどういう意味ですか?あなたが使用しているものは何でもQUADPACKとは異なるアルゴリズムを使用して意味するのですか?
編集:完全を期すために、ここで私は0の平均値と0から+無限大への1の標準偏差を有するガウスのために試してみた何のようなものです。
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 )
これはまだ書くのは非常に簡単ガウス関数は少し長いので少し醜いですが、。