2011年10月31日

n次元確率ベクトルの1次関数の平均と分散

いやはや、数理統計学でここ一週間ずうっと悩んでいたものがようやく解決した。
やはり、独立というものをきちんと理解しておかないといけない。

a,Xがそれぞれベクトルと行列だとして、cov(Xi, Xj)=0のとき(X1,...,Xnが互いに独立のとき)、
var(aX) = a1^2*var(X1) + a2^2*var(X2) +, ..., + an^2*var(Xn)
となります。これは、a*var(X)'*aと一致するはずなのですが、一致しません。なんでだろう、と思って色々やってみていたのだけれど、要は平均すれば一致する、ということのようです。シミュレーションしてみました。

> set.seed(1)
> n <- 100
> S <- 40
> a <- seq(1, n)
> X <- X2 <- X3 <- NULL
> temp <- temp2 <- numeric(s)
> for (s in 1 : S) {
+ X <- matrix(runif(n^2), nrow=n, ncol=n)
+ X2 <- matrix(NA, nrow=n, ncol=n)
+ X3 <- matrix(0, nrow=n, ncol=n)
+
+ ###期待値ベクトル###
+ ###ほぼ0だが、微妙に差が出てしまうようだ###
+ ###aが横ベクトルであることに注意###
+ round(apply(a*X, 1, mean)-
+ a*apply(X, 1, mean))
+
+ ###共分散行列###
+ for(j in 1 : ncol(X)){
+ for (i in 1 : nrow(X)) {
+ X2[j, i] <- cov(X[j, ], X[i, ])
+ }}
+ ###対称行列となっている###
+ #X2==t(X2)
+ #var(X)==t(var(X))
+ #X2==var(t(X))
+
+ ###独立とは言っても一回一回はばらつく###
+ temp[s] <- a%*%X2%*%a-
+ (a^2)%*%apply(X, 1, var)
+
+ ###対角だけ取ってきて無理矢理独立にすれば(?)いつでも差は0###
+ diag(X3) <- diag(X2)
+ temp2[s] <- a%*%X3%*%a-
+ (a^2)%*%diag(X3)
+ }
> jpeg("variance_independent.jpeg")
> hist(temp)
> dev.off()
windows
2
> temp;temp2
[1] -173.4482 1113.3236 -7706.5664 502.3490 -7383.8686 -993.2250
[7] 1724.8721 2393.0135 5236.6691 3939.5105 827.2226 -1489.3155
[13] -2922.6024 447.2575 -1754.0834 -1780.8433 1820.3705 -10123.5030
[19] -2082.8554 2079.5110 -3439.7196 -2314.1255 -5338.1769 6276.1852
[25] -2316.7343 1050.0900 -2247.3146 -3779.2307 8175.3117 1395.5665
[31] -3659.9512 -2296.1468 -7303.6068 116.1688 -1407.6821 -4242.2000
[37] -1149.0380 6796.3341 952.2544 797.3517
[1] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 -3.637979e-12
[6] -3.637979e-12 0.000000e+00 0.000000e+00 -3.637979e-12 0.000000e+00
[11] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[16] 0.000000e+00 3.637979e-12 0.000000e+00 0.000000e+00 0.000000e+00
[21] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 -3.637979e-12
[26] 0.000000e+00 -3.637979e-12 0.000000e+00 3.637979e-12 -3.637979e-12
[31] 0.000000e+00 3.637979e-12 0.000000e+00 7.275958e-12 0.000000e+00
[36] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00

出来た図。
variance_independent.jpeg
確かに、何度も繰り返してやれば一応イコールで結ばれるようです。なるほどなぁ。独立って難しい。
posted by しばきん at 00:17| 神奈川 ☔| Comment(0) | TrackBack(0) | 統計 | このブログの読者になる | 更新情報をチェックする

2011年10月30日

神様のパズル

神様のパズル



昨日、初めてお台場とやらに行って、日本科学未来館というところに行ってきた。基本的に物理の面白いお話を多く紹介していて、要所要所にボランティアの解説員がいる。一人のおじさんとニュートリノのコーナーで熱く議論をしてしまい、人の話が終わるまで質問するな、と怒られた(笑)。だってエネルギーが質量となるなんて(E=mc^2)よく分からなかったんだもの。ってか80歳だそうだが、あそこまで数式を巧みに操れるものなのか、と思った。うーん、僕も頑張らねば。

その後、ぶーらぶらするがお土産コーナーでこの本を発見。文庫版だったが、面白くて面白くて、昨日一気に読み終わってしまった。宇宙は作れるのか、というお話から始まっているから物理一色の話かと思いきや、途中から哲学的要素も含みだす。ディベートで、宇宙は人が作れる派と作れない派で分かれるのだが、作れない派が負けそうになったときに、「宇宙は人が作っていないなら、誰が作ったのか」という旨の問いをするくだりがある。そこで今まで作れないとしていた派が分裂して、神を容認する立場、科学のみを信じる立場、科学で説明できないことがあることは認めるが、神も容認できない立場、に分かれる。その後も、
「彼」という表現で、何度も宇宙を創った人、が登場する。このあたり、非常に僕は好きだ。途中の専門用語はちんぷんかんぷんだったが、文庫でここまで数学や物理のことを組み込んだ本は初めて読んだ。数式は一切出てこないし、恋愛色も薄いが、夢中になれてしまった。次はどうなるどうなる?と非常にわくわくした。大変オススメの一冊。ちなみに本のなかでも、ヒロインの天才少女はE=mc^2が実感としてよくわからんと言っていた。こういうの見るとほっとするあたり、いかんな。
posted by しばきん at 10:27| 神奈川 ☁| Comment(0) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする

2011年10月26日

明日のゼミ

明日はゼミでドクターの進捗状況の発表。
とにかく分かりやすくを心がけたが、少し崩しすぎたかもしれない。

それにしてもこの時期にまだ論文が出ていないのがまずい。焦っても仕方ないが、今後は最適性の原理に従った行動をしなければならないだろう。
「与えられた初期値が何であれ残りの決定は最初の決定から生じた状態に対して最適でなければならない」
そのために坊主にしたのだし。頑張らねば。
posted by しばきん at 17:50| 神奈川 ☀| Comment(0) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする

2011年10月25日

ワードでコメントが多すぎた

最近自分の論文に時間を割かれて、ほとんど他人の論文が読めていない。きちんと時間管理をしなければこうなるのだな、と反省。また少しずつ前のペースで論文が読めるよう調整していきたい。ところで先ほど、ワードでの変更履歴の記録が多すぎて、ワードがクラッシュした。おまけに、今日半日分の修正箇所が全て復元されなかった。今それを直しているが、今度からあまり変更履歴をため込まないよう細心の注意を払いたい。
posted by しばきん at 19:29| 神奈川 ☀| Comment(0) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする

2011年10月14日

対数正規分布と正規分布

たまに両者の関係を忘れることがあるので復習。

> ###対数正規分布と正規分布の関係###
> set.seed(1)
>
> ###LN(2, 0.5^2)に従う乱数を発生###
> y <- rlnorm(1000, mean=2, sd=0.5)
>
> ###平均値と標準誤差###
> mean(y);sd(y)
[1] 8.394254
[1] 4.644934
> hist(y)
>
> ###対数化して平均値と標準誤差を算出###
> y_log <- log(y)
> hist(y_log)
>
> ###ちゃんと平均値が2, 標準誤差が0.5にほぼなっている###
> mean(y_log);sd(y_log)
[1] 1.994176
[1] 0.5174579
>
> ###対数化した値から、対数化する前の平均値と標準誤差を推定###
> ###mean(y);sd(y)とほぼ同じ値を返す###
> exp(mean(y_log) + sd(y_log)^2/2)
[1] 8.398535
> sqrt(exp(2*mean(y_log) + sd(y_log)^2)*(exp(sd(y_log)^2) - 1))
[1] 4.653705
ラベル:統計学 R
posted by しばきん at 20:24| 神奈川 ☔| Comment(0) | TrackBack(0) | 研究生活 | このブログの読者になる | 更新情報をチェックする
×

この広告は1年以上新しい記事の投稿がないブログに表示されております。