6

scipy のクワッド メソッドを使用してガウスを積分しようとすると (ガウスという名前のガウス メソッドがあるとしましょう)、必要なパラメーターをガウスに渡し、クワッドを残して正しい変数の積分を行うのに問題がありました。多次元関数で quad を使用する方法の良い例はありますか?

しかし、これにより、一般的にガウス分布を統合するための最良の方法について、より壮大な疑問が生じました。ガウス分布が scipy に統合されているのを見つけられませんでした (驚いたことに)。私の計画は、単純なガウス関数を作成し、それをクワッド (または、現在は固定幅の積分器) に渡すことでした。あなたならどうしますか?

編集:固定幅は、固定dxを使用して曲線の下の面積を計算するtrapzのようなものを意味します。

これまでのところ、quad に入ることができるラムダ関数を返すメソッド 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、mu、およびシグマで呼び出す必要があります)を渡して、クワッドのような値をいくつか入力しようとしたとき

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

パラメータ 10、2、および 0 は、必ずしも N=10、sigma=2、mu=0 と一致するとは限らなかったため、より拡張された定義が求められました。

scipy.special の erf(z) では、最初に t が何であるかを正確に定義する必要がありますが、そこにあることを知っておくと便利です。

4

5 に答える 5

34

わかりました、あなたはいくつかのことについてかなり混乱しているようです。最初から始めましょう。「多次元関数」について言及しましたが、その後、通常の 1 変数ガウス曲線について説明します。これは多次元関数ではありません。積分すると、1 つの変数 (x) だけが積分されます。真の多次元関数である「多変量ガウス分布」と呼ばれるモンスターがあり、統合する場合、2つ以上の変数を統合する必要があるため、区別することが重要です(前述の高価なモンテカルロ法を使用します)。しかし、あなたは通常の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は であることに注意してください: 他の 3 つのパラメーターは、特定のガウス分布に対して固定されています。

ここで、数学的な事実として、すべてのガウス曲線が同じ形状を持っていることは確かに真実であり、わずかにずれているだけです。したがって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)の正規分布の積分です。誤差関数の定義により、次のようになります。musigmanz

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

この事実の詳細または証明が必要な場合は、通常の CDF に関するウィキペディアの記事を参照してください。

わかりました、それで十分な背景説明になるはずです。(編集済みの) 投稿に戻ります。「scipy.special の erf(z) では、最初に 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 CPU サイクルがかかり、除算には約36、累乗は約 200 です。これは数年前のことです。それでも、それはそれらの相対的な複雑さを示しています。同様にexp(x)、力ずくで計算するのは非常に悪い考えです。の優れた実装を作成するときに、一般的なスタイルのべき乗exp(x)よりも大幅に高速かつ正確にすることができるトリックがあります。a**b

定数 pi と e の numpy バ​​ージョンを使用したことはありません。私はいつも昔ながらの数学モジュールのバージョンに固執してきました。どちらかを好む理由がわかりません。

quad()電話で何をしようとしているのかわかりません。quad(gen_gauss, -inf, inf, (10,2,0))再正規化されたガウスをマイナスの無限大からプラスの無限大まで積分する必要があり、ガウスは実線上で 1 に積分されるため、常に 10 (正規化係数) を出力する必要があります。10 からかけ離れた答え (正確には10とは思わないでしょうquad()) は、どこかで何かがおかしくなっていることを意味します...実際の戻り値と、おそらく の内部の仕組みを知らずに、何がおかしくなったかを言うのは難しいですquad()

うまくいけば、混乱の一部が解明され、エラー関数が問題に対する正しい答えである理由と、興味があればそれをすべて自分で行う方法が説明されました. 私の説明が明確でない場合は、まずウィキペディアをざっと見てみることをお勧めします。まだ質問がある場合は、遠慮なく質問してください。

于 2009-02-04T12:32:51.073 に答える
15

scipyには、「エラー関数」、別名ガウス積分が付属しています。

import scipy.special
help(scipy.special.erf)
于 2009-02-04T04:11:03.257 に答える
4

ガウス分布は正規分布とも呼ばれます。scipy norm モジュールの 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

于 2011-01-18T13:49:59.387 に答える
3

多変量ガウス分布を扱っていると思います。もしそうなら、SciPyはすでにあなたが探している関数を持っています:それはMVNDIST(「多変量正規分布」)と呼ばれます。SciPyのドキュメントは相変わらずひどいので、関数がどこに埋め込まれているかさえわかりませんが、ドキュメントはSciPyの最悪の部分であり、過去に終わりのない私を苛立たせてきました。

単一変数ガウス分布は、古き良き誤差関数を使用するだけであり、その多くの実装が利用可能です。

一般的な問題への攻撃に関しては、そうです、James Thompsonが述べているように、独自のガウス分布関数を記述して、それをquad()にフィードしたいだけです。ただし、一般化された積分を回避できる場合は、そうすることをお勧めします。特定の関数に特化した積分手法(MVNDISTが使用するような)は、標準のモンテカルロ多次元積分よりもはるかに高速になり、非常に遅くなる可能性があります。高精度のために。

于 2009-02-04T04:10:36.423 に答える
2

-infinity から +infinity への積分を常に行って、常に答えがわかるようにしないのはなぜですか? (冗談です!)

私の推測では、SciPy に定型化されたガウス関数がまだない唯一の理由は、それが簡単な関数であるということです。独自の関数を作成し、それを quad に渡して統合することについてのあなたの提案は、素晴らしい音に聞こえます。これを行うために、承認された SciPy ツールを使用します。最小限のコード作業で済み、他の人が SciPy を見たことがない場合でも非常に読みやすくなっています。

固定幅積分器とは正確には何を意味しますか? QUADPACK が使用しているものとは異なるアルゴリズムを使用するということですか?

編集:完全を期すために、平均が 0 で標準偏差が 0 から +infinity までの 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 )

Gaussian 関数は少し長いため、これは少し醜いですが、それでも書くのはかなり簡単です。

于 2009-02-04T03:59:32.280 に答える