2012-03-22 25 views
2

分布のn番目の瞬間を計算したい。all.moments function奇妙な結果

library(moments) 
r<-rnorm(10000) 
rr<-all.moments(r,order.max=4) 
rr 

[1] 1.000000000 0.002403360 0.962201478 -0.022694670 2.852696159 

これは私が、3-ことを知っている原因、それは、真実ではないように私には思える:私は、私は、このようにall.momentsをテストしているR.にライブラリ「瞬間」のall.moments機能を使用するようにしようとしています第4および第4モーメントは、正規分布において0でなければならない。 私のミスはどこですか?

答えて

6

第3の瞬間は、の歪度です。あなたは正しいです:正規分布の場合、これはゼロです。あなたは正規分布からサンプリングするだけなので、結果はほぼゼロになります。

4次の瞬間は、のkurtosisです。正規分布の場合、これは3σ^ 4です。この場合、σは1なので、結果は3になるはずです。


見積もりの​​精度を向上させるには、サンプルサイズを改善します。 1E7観測のサンプルについて:

> library(moments) 
> r <- rnorm(1e7) 
> all.moments(r,order.max=4) 
[1] 1.0000000000 0.0004028138 0.9995373115 0.0007276404 2.9976881271 
1

ことが期待で唯一の真のだからではなく、正確に、かつ高いモーメントが大きなばらつきを持っているので?

(第四モーメント(以下V5)がゼロに近くさえない理由のために、同様Andrieの答え@を参照してください。)

> library(moments) 
> R <- t(replicate(50,all.moments(rnorm(1e4),order.max=4))) 
> summary(R) 
     V1   V2     V3    V4    
Min. :1 Min. :-0.024921 Min. :0.9714 Min. :-0.0987174 
1st Qu.:1 1st Qu.:-0.009527 1st Qu.:0.9911 1st Qu.:-0.0341950 
Median :1 Median : 0.001021 Median :0.9994 Median : 0.0067138 
Mean :1 Mean :-0.001047 Mean :1.0006 Mean :-0.0002613 
3rd Qu.:1 3rd Qu.: 0.004711 3rd Qu.:1.0147 3rd Qu.: 0.0299731 
Max. :1 Max. : 0.023356 Max. :1.0398 Max. : 0.1283456 
     V5  
Min. :2.775 
1st Qu.:2.921 
Median :3.005 
Mean :3.007 
3rd Qu.:3.092 
Max. :3.325 
0

私はpythonで働いて、同じ問題を抱えて、私はと思いますコンピュータ上の浮動小数点演算の精度は限定されています。多数の力を持っていても助けにはならない。 Pythonコードと結果を切り取って貼り付けます。このコードは、標準法線の最初の20分を計算しようとします。要するに、分布の高次モーメントを数値的に計算することは容易ではないと思います。ここでは「高次」とは10を超える意味です。別の実験では、より多くのサンプルを描画することで18日目に得られる分散を減らそうとしましたが、それは私の「普通の」コンピュータでは実用的ではありませんでした。

N = 1000000 
w = np.random.normal(size=N).astype("float128") 

for i in range(20): 
    print i, mean(w**i) # simply computing the mean of the data to the power of i 

はあなたに与える:

0 1.0 
1 0.000342014729693 
2 1.00124397377 
3 0.000140133725668 
4 3.00334304532 
5 0.00506625342217 
6 15.0227401941 
7 0.0238395446636 
8 105.310071549 
9 -0.803915389936 
10 948.126995798 
11 -34.8374820713 
12 10370.6527554 
13 -1013.23231638 
14 132288.117911 
15 -26403.9090218 
16 1905267.02257 
17 -658590.680295 
18 30190439.4783 
19 -16101299.7354 

しかし、正しい瞬間があります:1、0、1、0、3、0、15、0、105、0、945、0、10395、0 、135135、0、2027025、0、34459425、0、654729075、

+0

数値の不正確さの問題ではありません。高次モーメントの分散が爆発するからです。 –

+0

しかし、ここでは、標準法線から数えてN個の乱数を与え、例えば、その18乗を計算してから、これの平均を計算します。この手順では、Nを大きくすると、18日目の真の値に近づくべきではないでしょうか? – Frank

+0

はい、ただし、*ロット*を増やす必要があります。 –