近似の方が正確？

質問

確率変数bが平均m、分散σ^2 の正規分布に従うとして、1/bの分布を考えます。mはσより十分大きいとして、bのとる値を
　b = (1 + ε)m
とするとε^２までとった時
　1/b ≒ (1 - ε + ε^2)/m
1/b^2 ≒ (1 - 2ε + 3ε^2)/m^2
すると1/bの平均と分散は
　E[1/b] = (1 + E[ε^2])/m
　E[1/b^2] = (1 + 3E[ε^2])/m^2
　σ^2[1/b] = E[1/b^2] -(E[1/b])^2
bの変動係数をCV=σ/m とすると E[ε^2] = CV^2 なので
　E[1/b] = (1 + CV^2)/m
　σ^2[1/b] = ( 1 + 3 CV^2 - (1 + CV^2)^2 )/m^2
例えばm=100, σ=10として乱数と上の式で計算したものを比較すると

平均
　乱数：0.01010367　上の式：0.0101
標準偏差　
　乱数：0.001043964　上の式：0.0009949874
とほぼ一致する結果が得られました。ところが
　E[1/b]
= (1/√(2πσ^2))∫db exp[-(b-m)^2/2σ^2]/b
として計算しようとすると積分がb～0のところで発散してしまって計算できません。そこでやっと質問です。mが10σより大きければbが０以下になることなど実際上ありません。そのような重要でない領域が計算上は大変な影響を与えてしまう。近似をするとちゃんと計算できるのに厳密にやろうとするとできないというのはどのように考えたらよいのでしょうか。またεの全てのオーダーまでとったとき1/bの平均と標準偏差をmとσで表わすとどうなるのでしょうか。

回答No.7

siegmund · ベストアンサー

siegmund です．

grothendieck さん：

> bは1/bの特異点の近傍でもランダムな値をとります。
> とすれば極限をとるのにコーシーの主値のように強い制限をつけるのではなく、
> 通常の広義積分の定義による方がランダムネスを忠実に表わしている様に思います。

私は No.3 でコーシーの主値のことを書いたのですが，
ランダムなサンプリングは grothendieck さんの言われるように
通常の広義積分の定義に対応する様に思われます．

> (2)の左辺は存在しないのに右辺は存在すると言う不思議なことになっています。

grothendieck さんがご自身で書かれていますように，
E[ε^n] = (n-1)!! (σ/m)^n　　(nが偶数の時)
ですので，(2)の右辺は数学的には存在しないというのが正しいでしょう．
E[ε^n] と E[ε^(n+2)] の比を取ってみればすぐにわかりますように，
今は σ/m=1/10 なので，n=100 で E[ε^n] が最小になります．
実際，
　　n=98 　E[ε^n] = 2.75292×10^(-22)
　　n=100 　E[ε^n] = 2.72539×10^(-22)
　　n=102 　E[ε^n] = 2.75265×10^(-22)
です．
まともな(？)値まで回復するにはかなり大きな n が必要で
　　n=200 　E[ε^n] = 6.66631×10^(-14)
　　n=250 　E[ε^n] = 0.00004
　　n=300 　E[ε^n] = 3.75327×10^6
です．
つまり，n を増やしてゆくとき，かなり長い間 E[1/b] は
安定であるように見えます．

grothendieck さんの最大サンプリングは 50000000 回で，
累積分布が 1/50000000 になるのは b=-5.49σ のところです．
したがって，
中心から 5.5σ以上離れた部分は事実上サンプリングされておらず，
10σ離れた付近から来る E[1/b] の異常性は観測できません．
サンプリング回数に関しても「安定領域」が存在します．

grothendieck さんの求められた
E[1/b] = 0.01010325 or 0.01010315
はこれらの「安定領域」での E[1/b] の(見かけの)値であると考えられます．
εによる展開であまり高次は取らない，
は中心からあまりはずれたところは考慮しない，
ということに相当するでしょう．
また，数値計算の方も，
中心から 5.5σ 以上離れた部分は考慮しない，
ということになっていますので，
両者で E[1/b] が一致するのはなるほどと言う気がします．

大きな n の E[ε^n] に効く b の値の領域も調べてみました．
x^n exp(-x^2/2σ^2) の極大がどこで起きるか調べますと(x=b-m)，
x=10√n とわかりますので
n=10 ・・・ x=3.16σ
n=30 ・・・ x= 5.47σ
n=50 ・・・ x= 7.07σ
n=100 ・・・ x= 10σ
n=200 ・・・ x= 14.14σ
n=250 ・・・ x= 15.81σ
n=300 ・・・ x= 17.32σ
積分値には極大付近が効くでしょうから，
中央から5.5σ程度まで取ると，
n≧30 の E[ε^n] は正しく求められないことになります．
上に書きましたように，解析的値は
　　n=200 　E[ε^n] = 6.66631×10^(-14)
　　n=250 　E[ε^n] = 0.00004
　　n=300 　E[ε^n] = 3.75327×10^6
でしたから n=300 で調べてみます
(n=200 程度じゃどうせ効きませんので)．

積分範囲　　数値積分結果
±10σ　　　E[ε^300] = 7.61932×10^(-24)
±15σ　　　E[ε^300] = 1374.25
±20σ　　　E[ε^300] = 3.75285×10^6
±25σ　　　E[ε^300] = 3.75327×10^6

という具合になりました．
したがって，中央から5.5σ程度までしか取らないことにすると，
E[ε^300] も効きません．

まとめてみます．

(a) grothendieck さんの数値計算は
中心値から±5.5σ程度までをサンプリングしていることになります．
したがって，b=0=m-10σ 付近から来る E[1/b] の異常性は観測されません．

(b) grothendieck さんの ε 展開は発散しますが，
n=250 程度までは一定値に収束するように見えます．
それ以後は発散してゆくのが見えますが，
大きな n に対する E[ε^n] には中心値から離れた場所からの寄与が効きますから，
例えば中心値から±5.5σ 程度しか取らないとすると，
一定値に収束したままです．

(c) (a)(b)ですべてコンシステントに説明できたと考えています．

(d) ポアソン分布で代用した結果と少し違うのは，
代用の結果でしょうか．
ポアソン分布と正規分布が実質的に同じになるのは
m>1000 のときなどと言われています．

(e) 統計力学で相転移関係の数値計算でも似たような事情が生ずることがあり，
サンプリングの回数，あるいは系の大きさ，などと
異常の生ずる相転移点に(例えば)温度がどれくらい近いが，
関連していたりします．

お礼

ご回答ありがとうございます。(σ/m)が１より大きい時1/bの平均は試行からは読み取ることができず、(σ/m)が１より小さいとき1/bの平均は始めは収束するように見えますが、あるところから先は収束しません。その意味では「1/bの平均は定義されない」という厳密な理論は正しいわけです。しかしm=100,σ=10 のとき誰でも1/bの平均としておよそ0.010103 という値を読み取るでしょう。1/bを計算したいとき、bがある誤差の範囲内でしか分かっていないということは現実には多いと思われます。
http://oshiete1.goo.ne.jp/kotaeru.php3?q=2382839

このとき「1/bの平均は定義されません。1/bの誤差は評価できません」という答えしか与えてくれない「厳密な」理論よりはεで展開して始めの数項をとった方が良いように感じられます。
　E[1/b]～(1 + E[ε^2] + E[ε^4] +…)/m　(2)
の右辺は収束級数としては存在しないが、漸近級数として存在すると考えています。量子電気力学でも「摂動展開は漸近級数だから理論値は存在しない」とはされていないと思います。

回答No.11

stomachman

No.10、書き間違い訂正です。

> 時刻tに於けるyを推定するのにこの式を使うのは全く妥当で、その推定誤差はNに従うと期待できます。しかし、yからtを推定するのに

は

時刻tに於けるyを推定するのにこの式を使うのは全く妥当で、その推定誤差は正規分布に従うと期待できます。しかし、yからtを推定するのに

に訂正です。

お礼

ご回答ありがとうございます。
[2] 正規分布に従う物理量bについて、bを観測したデータからx=1/bを計算したとき、xの分布や期待値を考える

ことは少なくないと思います。例えば光速cを測定した時、cの測定誤差が微細構造乗数α=e^2/hc にどの程度の誤差を与えるか知りたいというようなことは当然の要求と思います。「厳密な」理論によれば「cに誤差がある時1/cの値はさっぱり分かりません。微細構造乗数の誤差は評価できません」ということで、おかしな話だと思います。

回答No.10

stomachman

stomachmanです。No.9へのコメントについてです。

> 「xの分布をξ(n,x)とすると、ξはxについて偶関数になるからξ(n,x)の平均は0」とは必ずしも言えないと思います。

　仰る通り、∀x(x≠0 → f(x)=f(-x)) であっても、fが「超関数の意味で偶関数」でなければ（つまりδ関数の奇数次の微分の項が入っていれば） ∫ {x=-∞～∞} xf(x)dx は0になりません。でもξにはそんな成分はありませんので、超関数の意味でも偶関数であり、積分は0になります。

> 非負の関数を積分をして負になると言うのは理解できません。

「非負の関数を定積分したら常に非負である」のなら、　a＜0＜bのとき、

∫ {x=a～b} (1/x^2) dx = 1/a-1/b ＜0

は誤りですか？
　左辺を素直に（x=0を避けて）数値積分しても、この積分はどうしても正にしかなりません。「ただしく」数値積分するには

∫ {x=-∞～∞} (1/x^2) dx = 0

（これもお説によれば理解できない式）を使って

∫ {x=a～b} (1/x^2) dx = -∫ {x=-∞～a} (1/x^2) dx -∫ {x=b～∞} (1/x^2) dx

に置き換える必要があります。(No.2で参考linkを示しました。）
　奇関数の積分の話はもうちょっとややこしくなりますが、同じようなことです。
　これは単なる机上のお遊びではありません。まさに上記の積分が実用されている例として、医療に使われるCTの画像再構成があります。データf(nΔx)（nは整数）に対して

g(x) = f(x) ＊　(1/x^2)

（＊はconvolution）をg(nΔx)について計算するんですが、これを実際どうやってやるかというと、

f(x) ＊　((1/x^2)＊b(x))
b(x)= |x|<(Δx/2)のとき　1/Δx, それ以外は0

に置き換えて、括弧内の＊は解析的に計算して　(1/((nΔx)^2)＊b(nΔx)を算出し（これはデータによらないので数表にしておけます）、外側の＊は数値積分します。(b(x)は多様な「平滑化関数」で置き換えることが出来ますが、本質的には同じ事です。）n=0における計算に注目です。無限大になっちゃう極限には（事実上生じ得ないほどの近く以上には）近づかないで、その周りで粗視化して丸めてしまう、という処方は「くりこみ」の最も素朴な説明とそっくりです。

> 積分を通常の広義積分とする理論は観測される値を与えてくれないと言うことを書きましたが、積分を正則化する理論は逆に観測されない値を特定してしまいます。

　母平均・母分散と標本平均・標本分散が乖離してしまう話(No.9)は、仰るのとまさに同じ事を言っています。
　しかし、むしろ、
[1] 物理量xについて、1/xを観測したデータbからx=1/bと推定し、しかもbは正規分布に従う、と仮定する
[2] 正規分布に従う物理量bについて、bを観測したデータからx=1/bを計算したとき、xの分布や期待値を考える
のどちらも、やってることがどこかおかしい（整合性のある物理的意味づけができない）のだと思います。
　たとえ話になりますけれど、
　y = A t + B
という理論に従う変量yをいろんなtについて測定し、その測定誤差が正規分布に従うとします。測定値から最小自乗法でA,Bを計算できたとしましょう。時刻tに於けるyを推定するのにこの式を使うのは全く妥当で、その推定誤差はNに従うと期待できます。しかし、yからtを推定するのに
　t = (y-B)/A
という手は使えません。A,Bを出すのに使った同じ測定値を用いて、最小自乗法で
　t = C y + D
のC,Dを計算すると、C≠1/A, D≠-B/Aになります。
　誤差を含んだデータと、誤差のないモデルとを単純に同一視できない、という事情は、こんな簡単な例に既に現れています。ですから、上記[1]の場合、「b=1/xだからx=1/b」というロジックは危ないのです。

> 漸近展開は(σ/m)が１より大きい時と１より小さいときの区別をします。

　そうでしょうか。打ち切る次数を上げれば広い範囲で収束すると思います。実際、m=0の正規分布を多項式展開するのは容易ですよ？

> こうしてみると近似であるはずの漸近展開が一番良いような気がしませんか。

　確かに、物理の常套手段ですね。ですが、超関数はこのあたりのもやもやを乗り越えるために発明されたとは、言えませんでしょうか。

お礼

ご回答ありがとうございます。Cauchy分布の密度関数
　f(x) = c/π(c^2 + (x-m)^2)
でパラメータmを通常は母平均とは呼ばないと思います。超関数は「一点を除いて０だが積分すると１になる関数」などを合理的に定義したいということが動機の一つだと思いますが、Cauchy分布などではむしろ平均は定義されないとする方が自然だと思います。確かに超関数論が威力を発揮するような問題もあるのでしょうが、少なくとも統計とは関係が少なそうに思います。

回答No.9

stomachman

No.6へのコメントについてです。

> 「断面積が負になるはずがないから」という理由で負の測定値を無視することはしないと思います（そんなことをしたら負の側のノイズは捨てて、正の側のノイズだけ残していることになるので断面積を過大評価することになります）。

仰る通りだと思います。No.6はこのお考えに沿って書いています。

　m=0の場合には、解析的には1/bの期待値は0になります。しかしRを使った数値実験が変な結果になる。これは、b=0に近い正と負のサンプルが均等には出現しないからでしょう。サンプル数を増やせば増やすほど、|b|が0にごく近いサンプルがたまたま含まれて来て、それだけで平均値を大きく動かしてしまう。
b<- rnorm(100,mean=0,sd=1);mean(1/b);max(1/b);min(1/b)
[1] -1.770347
[1] 16.01663
[1] -147.8735
てな具合に。

　ですから、ご質問は「有限個のサンプルによる平均値が、母集団の平均値(期待値）から大きくずれてしまう」ということを問題になさっているのだろうと思います。これは「有限個のサンプルに基づいて計算した1/bの平均の分布」がどうなるか、という話(No.7でsiegmund教授が論じていらっしゃるのもこのことだろうと思います）であり、母集団の期待値とは区別して議論する必要がありそうです。
　「平均0の正規分布に従う確率変数bのn個のサンプルから計算した1/bの平均xの分布」をξ(n,x)とすると、ξはxについて偶関数になるからξ(n,x)の平均は0。しかし分散を考えると、n=1において、ξの分散は（定数倍を無視して）
σ^2 = ∫{-∞,∞}exp(-(b^2))/(b^2) db
であり、この積分（が存在する事は既にNo.4のフーリエ変換の議論で示した通りですが、）は負の値になっちゃう。当然ξ(n,x)の分散(σ^2)/nも負であり、つまりξには普通の意味では分散が存在しません。数値実験を相当の回数繰り返せば分布の形が見えて来るでしょうけれども、それと共に、幾らでも絶対値の大きい値が現れて来る。
　だから数値実験で(1/b)の期待値を出そうとすると、幾らやっても収束しない、という事情になるんでしょう。

　となると、(1/b)の期待値を計算することは、それ単独でホントに物理的に意味あるの？という方から攻める手かな、と想像します。例えば、物理的に意味のある確率変数x(＞0)があって、xの逆数b=1/xを観測する実験において、ノイズ成分を出来るだけ除いたらb≦0になることもある。このとき、x=1/bと推定するのが妥当とは限りません。これは、測定値bにはまだ誤差が含まれているということを考慮してxの最尤推定をどう行うか、という話でしょう。もちろん、その誤差の性質によって、またxに関する先験分布をどの程度知っているかによっても、確率モデルを作る処方は違ってくるでしょう。

お礼

ご回答ありがとうございます。
「xの分布をξ(n,x)とすると、ξはxについて偶関数になるからξ(n,x)の平均は0」
とは必ずしも言えないと思います。奇関数の積分をすべて０とすれば平均を持たないとされるCauchy分布でも平均を持つことになります。
「σ^2 = ∫{-∞,∞}exp(-(b^2))/(b^2) db
は負の値になっちゃう。」
非負の関数を積分をして負になると言うのは理解できません。
下に積分を通常の広義積分とする理論は観測される値を与えてくれないと言うことを書きましたが、積分を正則化する理論は逆に観測されない値を特定してしまいます。また通常の広義積分、積分正則化のどちらも(σ/m)が１より大きい時と１より小さいときの区別をしませんが、この二つは質的に異なるように思われます。漸近展開は(σ/m)が１より大きい時と１より小さいときの区別をします。そして(σ/m)が１より小さいときは観測されるおよその値を与えます（およそで良いのです。それしか観測されないのですから）。
こうしてみると近似であるはずの漸近展開が一番良いような気がしませんか。

回答No.8

siegmund

siegmund です．
ちょっと考察に不足がありました．

∫{-∞～∞} (1/b) f(b) db
f(b) ≡ exp[-(b-m)^2/2σ^2]
を考えるときに，積分に効く b の幅をあまりちゃんと考えていませんでした．

(1/b) f(b) の極大は b=98.9898 で起こり，
このときの (1/b) f(b) の値は 4.009×10^(-4) です．
これは b = m (= 100)のときの
(1/m) f(m) = 3.989×10^(-4)
とほとんど変わりません．
つまり，割合ブロードで積分に効きやすい．

一方，(1/b) f(b) は b=0 で発散しますが，
積分に効く領域は非常に狭いことになっています．
f(0) = 7.694×10^(-24)
ですから，(1/m) f(m) と同程度の値になる b は
b = 10^(-20) 程度になります．
乱数による数値計算で b=0 付近からの異常性を観測するには，
そもそも確率分布が f(0) = 7.694×10^(-24) と小さい付近で
-10^(-20) ～ 10^(-20) をピンポイントでヒットしないといけないわけで，
その実現確率は
7.694×10^(-24) × 2 × 10^(-20) ～ 10^(-43)
程度です．
１秒間に１テラ回(10^12)回のサンプリングを宇宙開闢以来続けても，
宇宙の寿命が 137億年～4.3×10^17 秒ですから，
サンプリング回数は 10^29 程度で，
問題の領域には全然ヒットしません．

お礼

ご回答ありがとうございます。私はこれまで漸近展開して始めの数項を取れば良いと楽観論的だったのですが、悲観論的になってきました。まれに大きく外れた値が観測されれば異常値として捨てられる場合が多いと思いますが、実はそれは異常値ではなくて、ほぼ収束しているように見えていたものが本当は収束していなかったということに過ぎないのですね。実験で真の値は（近似的にすら）分からない、あなたの捨てたその観測値は異常値ではありませんということで、実験家にはショッキングなことではないでしょうか。でも私は漸近展開が捨てきれない気がします。

回答No.6

stomachman

どうも違和感があります。話を変な方向へ先走らせすぎているのではないかという…ご質問の真意を再度確認すべきかと思います。

> E[1/b]= (1/√(2πσ^2))∫db exp[-(b-m)^2/2σ^2]/b
として計算しようとすると積分がb～0のところで発散してしまって計算できません。

と仰っている。
　正規分布の確率密度関数φ(x) （平均m、標準偏差σ)は十分滑らかなので、x=0の周りで級数展開して、少なくとも|x|≪1のとき

φ(x)=Σ{k=0～∞}a[k](x^k)

とできます。-cからcの範囲(0＜c≪1)での定積分

I = ∫{-c,c}(φ(x)/x)dx

を考えたとき、高校のレベルでフツーに計算すると

I = a[0]∫{-c,c} (1/x)dx +Σ{k=1～∞}a[k]∫{-c,c} (x^(k-1)) dx
= a[0]((ln |c| - ln |-c|)+Σ{k=1～∞}a[k]((c^k)-((-c)^k))/k
= a[0]((ln |c| - ln |-c|)+Σ{k=1～∞}a[k]((c^k)-((-c)^k))/k
= 2Σ{k=1～∞}a[2k-1](c^(2k-1))/(2k-1)

となって、特異点らしきものはどこにもない。

　でも、こんなの計算できないと仰る。そればかりか、（この展開を使うとI=0になるのが自明であるところの）m=0の場合こそが一番オカシイと仰る。

　うーむむ、だとすると、これはきっと、上の計算では「公式通り」で通過したa[0]の項

∫{-c,c} (1/x) dx = (ln |c| - ln |-c|) = 0

というのが本当に辻褄が合っているのか、という所にご質問の焦点があるんだろう、と考えていたんです。

　しかしです。物理の難しい計算をガンガンやってらっしゃるgrothendieckさんが「今頃になって」

∫{-c,c}奇関数dx = 0

に引っかかりを感じるというのはなんだか変な気もします。こりゃ、stomachmanはポイントを外しているのではないか。

お礼

ご回答ありがとうございます。通常の標本平均はコーシーの主値に対応はしていないと思います。m=0の場合にRで1/bの標本平均を計算してみましたが、やはり下のように0に収束する気配は感じられません。

b <- rnorm(10000,mean=0,sd=10)
> mean(1/b)
[1] 0.04400062
> rm(b)
> b <- rnorm(100000,mean=0,sd=10)
> mean(1/b)
[1] -0.03495033
> rm(b)
> b <- rnorm(1000000,mean=0,sd=10)
> mean(1/b)
[1] -4.213309
> rm(b)
> b <- rnorm(10000000,mean=0,sd=10)
> mean(1/b)
[1] -9.555608

したがって、m=0では定義されず、m=100,σ=10 のときには 1/b の平均として 0.010103…を与えるような方法が望ましいのではないかと思います。
ニュートリノか何かの稀な反応の断面積を測定するとします。このときいくつかのカウントが観測され、それからバックグラウンドを差し引いたら断面積が負になってしまったということは十分にあり得ることです。このとき「断面積が負になるはずがないから」という理由で負の測定値を無視することはしないと思います（そんなことをしたら負の側のノイズは捨てて、正の側のノイズだけ残していることになるので断面積を過大評価することになります）。そう考えると、積分をb=0の手前で切断することにも問題がありそうです。やはり元は近似であったはずの漸近展開を適当なところまでとるのが一番良さそうに思えます。

回答No.5

stomachman

No.2へのコメントを拝見しました。

　No.4にも書きましたけれど、1/xを超関数として扱うなら全然問題ない、というのが一つのアプローチです。物理をお考えなら、超関数を持ち込むのが駄目ということもないんじゃありませんでしょうか。

　超関数の定義の仕方はいろいろありますが、ここでは
1. ガウス関数を、x∈{-∞,∞}の定積分を1に保ったまま幅を狭めた極限がδ関数
2. 超関数gに収束する関数列{g[n]}について、{g'[n]}が超関数に収束するとき、それがgの導超関数g'
を元に、一群の具体的な超関数を作り出すというやりかたで十分です。こうすれば、作られた超関数は全部フーリエ変換できる。その上、{-∞,∞}の定積分に関して部分積分ができるという性質を持っています。任意の「良い」関数T（全ての自然数Nについて、|x|→∞においてTの全ての導関数とTがO(1/|x|^N)である関数）に対して
∫{-∞,∞}g'(x)T(x)dx=-∫{-∞,∞}g(x)T'(x)dx
と言える、ということです。（ガウシャンはまさに「良い」関数の代表選手です。）

　超関数としての1/xは、概略以下のように導入できます。関数方程式
x f(x) = 1
を考えると、一般解として任意定数Cを含む
f(x) = (d/dx)(ln |x|) + C δ(x)
が出ます。そこでさらにf(x)が奇の超関数（すなわち、任意の「良い」偶関数G(x)について、f(x)G(x)をx∈{-∞,∞}で定積分すると0になるもの）であることを要求すると、解は一意に決まって
f(x) = (d/dx) ln |x|
x≠0のとき f(x) = 1/x
∫{-∞,∞}f(x)dx = 0
などが言えます。
　つまり、単にf(x)=1/xと言っただけでは、超関数の観点からはC δ(x)の分だけの不定性を含んでいる。普通の関数だと思って見ているうちはx=0を特異点として除外しているのですが、その特異点にδ関数が入り込んでいる。超関数としての1/xはこのδ関数を取り除いたものであり、いわば「pureな1/xを取り出した」に過ぎません。No.2で「格別オカシなことをしているわけじゃありません」と申し上げたのはそういう意味です。
　また、f(x)=1/xをフーリエ変換したものFで考えれば、Fは普通の関数なのだけれど任意定数項Cを含んでいる。そこでさらにFが奇関数であることを要求することによってC=0にしたことになります。
　なお、ガウス分布の平均がm=0の場合の(1/b)の平均は、（No.4のフーリエ変換の議論をするまでもなく、(1/x)を奇の超関数とした時点で直ちに、）0と決まります。「奇関数と偶関数の積は奇関数だから積分したら0」と素朴に言っちゃって構わないわけです。

　No.2で括弧付きで「有限部分」と書いちゃったのは、正確にはご指摘の通り-1乗の特異性をもつ積分に関するCauchy主値のことです。しかし超関数としての1/xならば、x=0の周りでの極限操作は不要です。C δ(x)が予め取り除いてあるためにx=0を避ける必要がない、と解釈することができます。

回答No.4

stomachman

siegmund先生、ご無沙汰してました。ここに来るとついつい遊び過ぎてしまうので、ちょっと自粛していたんです。で、ご回答No.3に関して、

> 何度測定しても長さや質量が負になることはありえません．

　ガウスの誤差論から出て来る正規分布は、無数の誤差要因が合成されて生じる、という前提がありますから、体重計で測った質量が負にならないのは当たり前、というルールの成り立つ世界とは相容れません。その意味でNo.2で「正規分布は極限として存在するのであって、そのままムキダシで現れて来るものじゃない」なんて書きました。これはご指摘の通りです。
　しかし、量子の世界の話となると、非局所性ということがある。無限遠まで広がった確率密度が瞬時に現れるコトがホントにありそうな気もしますんで悩ましい。（逆に、シュレーディンガー方程式に何か非線形項が抜けてるんじゃないか、という議論もありますが。）

> E[1/|b|] や E[1/b^2] だとそういうこともできません．

　いや、1/b^2の場合は同じようなもんです。1/b, 1/b^2 は、それぞれ超関数としてフーリエ変換で -πi sgn(ω), -(π^2)|ω| に写ります(sgnは符号関数）。正規分布の方は、ヨコにずらして平均を0に持って来てからフーリエ変換するとガウス関数のまま。ですから、どちらの場合も積を取れば（たとえば-(π^2)|ω|×ガウス関数(ω)）、|ω|が大きいところで速やかに0に落ちてくれます。ということは、この積はおとなしい関数のフーリエ変換になっているのであり、逆変換したものは特異点などない普通の関数g(x)になる。そして、g(x)は元々ガウス関数と1/bあるいは1/b^2との畳み込み積だから、正規分布の平均値がxであるとき（分散は一定のまま）の（1/bなり1/(b^2)なりの）期待値を表す関数になってる筈です。
　でも1/|b|の場合はこの手も駄目。フーリエ変換すると -2(ln(|ω|)+C) なので、ω≒0のところでやっぱり吹っ飛んでいますから、ガウス関数と積を取って逆変換したら、出てくるのは特異点を持つ超関数になっちゃいます。

> 具体的にどうされたかは不明ですが，正規分布ということを前面に出すなら，

　stomachmanも、おそらくgrothendieckさんの計算はsiegmund先生の仰るやり方でなさったんだろうと思います。数値計算でご質問の積分を調べるには、1/bの大きさごとにb軸の方を区間に分けて、それぞれで積分を考える、いわばルベーグ積分方式が良いかと思います。b≒0近辺での打ち消し合いがはっきり観察できるでしょう。

> Poisson 分布でλが大きいときは正規分布で近似できることが知られています

　b≦0では確率密度が断然0であるポアソン分布が、決して0にならない正規分布で「近似できる」とはどういうことか、というところがちょっと問題アリで、その問題アリはb≦0の部分の累積確率密度の程度のアリだろうと思います。すなわちNo.2の最後に例示した
f(b) = max(0, exp[-(b-m)^2/2σ^2] -c)
という確率密度関数と丁度同程度の問題アリ、ではないでしょうか。
　もちろん、この近似に問題アリなのか、それともそもそも正規分布だと仮定したことに問題アリなのか、というところこそがポイントになる訳です。言い換えると、もし1/bの代わりに1/|b|でもいい、という話であるのなら、b＜0での確率密度は0とすべきじゃないですかというコトです。

回答No.3

siegmund

grothendieck さん，stomachman さん，お久しぶりです．
特に stomachman さん，復活うれしい限りです．

純粋に数学的に見る限り E[1/b] が計算できないのは grothendieck さんご自身，
および totoro7683 さん，stomachman さんの言われるとおりです．
E[1/b] の場合は b=0 付近からの log 発散は
totoro7683 さん，stomachman さんのように
主値積分と解釈することによって一応回避できますが，
E[1/|b|] や E[1/b^2] だとそういうこともできません．

stomachman さんの二番煎じになりますが，
実際の問題ではある量が正規分布に従うとしていても，
それが使えるのは分布の中心付近(数倍のσ程度)でしょう．
測定論では正規分布が基本ですが
> 例えばm=100, σ=10として
が長さや質量の測定だとすると，
何度測定しても長さや質量が負になることはありえません．

さて，grothendieck さんは乱数で E[1/b] や σ[1/b] を計算されました．
具体的にどうされたかは不明ですが，正規分布ということを前面に出すなら，
(1)　　f(b) = exp[-(b-m)^2/2σ^2]
に従う確率で b のサンプリングを多数回おこない，
そのときの 1/b あるいは 1/b^2 の値を求めて加え合わせて平均，
ということになると思います．
こういうサンプリングをすると，
発散の原因になる b=0 付近は f(b) の値が非常に小さいため，
事実上サンプリングされないことになります．
実際，m=100, σ=10 としますと
(2)　　f(0)/f(m) = exp(-50) ≒ 2×10^(-22)
です．
したがって，grothendieck さんの乱数による数値計算は正規分布の中央付近だけ
見ていることになっているのでしょう．
こう考えると，grothendieck さんの近似計算が乱数計算結果を非常によく再現するのも
納得できます．

今は b=0 のところが(2)のように実現確率が極めて低くなっていますので
期待値の発散の様子が見えません．
m，σを適当に選んで実用的なサンプリング回数で b=0 付近がヒットするようにすれば，
サンプリング回数と共に期待値が発散してゆくのが見えるのではないでしょうか．
E[1/b] ですと正負のキャンセルがありますから E[1/b] よりは E[1/|b|] や E[1/b^2] の方が
素直に見えそうです．

なお，
(3)　　E[1/|b|] = ∫db f(b)/|b|
などを計算する際に，b が一様分布として f(b)/|b| の期待値を求めるという考えもあります
(積分区間は適当にカットしないといけませんが)．
これは普通の数値積分に近い手法です．
これだとたちどころに発散が引っかかるでしょう．

正規分布的アプローチ，一様分布的アプローチ，
どちらも数式的には同じことですが，
有限回の試行という点からは大きな違いがあります．

grothendieck さんの例は m=100, σ=10で， m=σ^2 の関係があります．
この関係からすぐ思い出すのは Poisson 分布
(4)　　P(n;λ) = exp(-λ)λ^n / n!
です(n=0,1,2,...)．
Poisson 分布は平均も分散もλですから，λ=100 と思えばよいわけです．
Poisson 分布でλが大きいときは正規分布で近似できることが知られていますから，
逆に今の問題を Poisson 分布で代用してみます．
正規分布は連続分布で Poisson 分布は離散分布ですから，
よく使われる処方箋にしたがって前者の 1/b 平均は後者の 1/(n+1/2) 平均に
対応するものとします．
そうしますと
(5)　　Σ{n=0 →∞} P(n;100) {1/(n+1/2)} = 0.0100508
(6)　　Σ{n=0 →∞} P(n;100) {1/(n+1/2)^2} = 1.0418458255157529×10^(-6)
から(mathmatica による)，
平均値が 0.0100508，標準偏差が 0.00102071 となり
grothendieck さんの結果
> 平均
> 　乱数：0.01010367　上の式：0.0101
> 標準偏差　
> 　乱数：0.001043964　上の式：0.0009949874
をよく再現します．

お礼

ご回答ありがとうございます。乱数の計算はお書きになっている方法でしました。試行回数をN、平均100標準偏差を10として統計ソフトRで書くと

　b <- rnorm(N,100,10)
　mean(1/b)
となります。するとNを大きくした時、ある値に近づくように見えます。

　　N　　1/bの標本平均
10000　0.01008869
100000　0.01009949
1000000　0.01010202
10000000　0.01010350
50000000　0.01010325

他方、
　1/b = 1/(1+ε)m
= (1 - ε + ε^2 - ε^3…)/m　　(1)
と展開するとεの次数を上げた時、乱数計算の値に近づくように見えます。

εの次数　(1)の平均
　0　　0.01
　2　　0.0101
　4　　0.010103
　6　　0.01010315

もっとも1/bの正確な平均は存在しないのでNをいくら大きくしてもある所から先は収束しなくなるのでしょう。εの方も
　E[ε^n] = (n-1)!! (σ/m)^n　　(nが偶数の時)
は(σ/m)が１より小さくても n→∞ で発散するので(1)は収束級数ではありません。
1/bの正確な平均は存在しないが、およその平均あるいは「事実上の平均」は存在し、その値は0.010103…になるのではないかと考えています（先生のされたPoissonn分布の結果と一致しないのは気になりますが）。そのまま乱数にすると精度が上がりにくいので計算するためにはご教示下さった一様分布でのf(b)/b期待値の計算が良い方法だと思います。ただ私の質問の眼目は

　E[1/b]～(1 + E[ε^2] + E[ε^4] +…)/m　(2)

のように漸近展開される実体は何かと言うことで、(2)の左辺は存在しないのに右辺は存在すると言う不思議なことになっています。

回答No.2

stomachman

[1] まず、数式の上での話としては、totoro7683さんと同意見です。
正規分布の確率密度関数は至る所0でないので、いくらmが大きくたって、b=0で0にはならない。だから1/bを掛けて積分すると、ややこしいことが起こる。しかしこれは正規分布の側の問題じゃなく、単に1/bの積分の問題です。
　単に1/bを積分したって、b=0付近で変なことが起こる。でも、1/bの積分を公式通り計算するとそんなの問題にならない。普段問題にならないから気にしてなかったものが、今回はアカラサマに見えて来た、という事情なのだと思われます。（こちらもご参考に→http://oshiete1.goo.ne.jp/kotaeru.php3?q=32716）。

f(b) = exp[-(b-m)^2/2σ^2]
∫[p,q]db f(b)/b 　(∫[p,q]db は、b=pからb=qまでの定積分）
を計算するにあたって、
∫[p,q]db f(b) = ∫[p,-ε]db f(b)/b + ∫[-ε,ε]db f(b)/b + ∫[ε,q]db f(b)/b
と分割します。
一方、f(b)は滑らかな関数ですから、ε→+0のとき
∫[-ε,ε]db f(b)/b → ∫[-ε,ε]db f(0)/b = 0
従って、ε→+0のとき
∫[p,q]f(b) db = ∫[p,-ε]db f(b)/b + ∫[ε,q]db f(b)/b
となる。積分すると、右辺の二つの項にはそれぞれεが残るけれども、打ち消し合ってくれます。残ったのが「有限部分」。怪しげな名前ですが、格別オカシなことをしているわけじゃありません。

[2]計算のテクニックで問題が処理できたって、そもそもそういう変な状況が生じるってこと自体がおかしいじゃないか、とgrothendieckさんならお考えだろうと思います。
　1/bは僅かな確率で、めちゃくちゃ巨大な値や負の値を取るんですよね。そして

> bが０以下になることなど実際上ありません。そのような重要でない領域

と仰っている。となると、はて、何をモデル化したらこんなものが出て来たんでしょうか。そこんとこを、ちょっと考えてみてはどうだろう、と思います。
　逆に言うと、「bが正規分布に従う」というのがそもそも仮定でしょう。裾野が幾らでも広がっている正規分布というものが実現するような確率現象なんてありますかね。正規分布は極限として存在するのであって、そのままムキダシで現れて来るものじゃない。そもそもが「近似用」なんだと考えることもできるわけです。
　1/bにゲンジツ問題における意味解釈を与えられるのは、例えばbがポアソン分布だとかベータ分布だとかに従っているようなときじゃないでしょうか。それを正規分布で近似して計算したら、裾野のところが邪魔になった。さて、何のために近似したかというと計算をラクにするためである。なのに計算で困ってしまったのなら、それは近似のしかたがまずかった、ということに他なりません。　このスタンスでは、

> εの全てのオーダーまでとったとき

には興味がない。
　ご質問にあるような多項式で近似するのも良いですし、裾野が邪魔なら切り取ったもので近似したって良いわけです。例えば
f(b) = max(0, exp[-(b-m)^2/2σ^2] -c)
でも良いでしょう。

お礼

ご回答ありがとうございます。確かに正の値しかとらない確率変数を負の値もとりうる正規分布で近似しているような場合もあるでしょう。しかし厳密に正規分布に従う場合もあると思います。例えば自由場の汎関数積分の測度はガウシアンになります。積分を正則化すれば良いのかも分かりません。bは1/bの特異点の近傍でもランダムな値をとります。とすれば極限をとるのにコーシーの主値のように強い制限をつけるのではなく、通常の広義積分の定義による方がランダムネスを忠実に表わしている様に思います。もっと困ったこともあります。平均mが1/bの特異点に一致し、m=0であるとします。このとき1/bの標本平均は収束しないでしょう。また通常の定義での1/bの平均も存在しません。ところが積分を正則化するとこの場合も平均が存在してしまいます。

回答No.1

totoro7683

統計のことは知りませんが、つまり発散積分から何らかの意味のある有限確定値を取り出すことが出来るかということですね。
これについてはアダマールの発散積分の有限部分という考えがあります。正の部分の積分を取り出すと、ε→0∫（ε→∞)db exp[-(b-m)^2/2σ^2]/bはbが０付近で発散しますが発散項を引いて有限部分を取り出します。bが負の部分も同じ。関数A(ｂ)=exp[-(b-m)^2/2σ^2]とおくとε→0(∫（ε→∞)dbA(b)/b+logεA(ε))
は有限確定値を定める。
これについてはアダマールの偏微分方程式という本に書いてあります。
E[1/b]に意味を持つのはこの有限部分で発散項は無視してよいのではないでしょうか。おそらくこれが求める真の値ではないでしょうか。なぜ無限大を無視して有限部分をとらなければいけないのかは分かりません。

お礼

ご回答ありがとうございます。後で書くように、例えばbが平均100、標準偏差10の正規分布に従うとき、1/bの平均は
0.010103...
になるというのが私の予想です。アダマールの発散積分の有限部分をとると、この値がでるのでしょうか。またコーシーの主値とは違うのでしょうか。教えて頂ければ幸いです。