--------(--)

スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

2012-08-04(Sat)

正規分布に従う乱数発生の方法

ほぼ以下を元ネタにしているのでこちらの方が詳細に書いてあると思う。
60. 確率分布と乱数

3 変量正規分布に従う乱数を生成する関数 r3norm() を定義する
> r3norm <- function(mu, A, n) {
> U <- svd(A)$u; V <- svd(A)$v; D <- diag(sqrt(svd(A)$d))
> B <- U %*% D %*% t(V) # 行列 A の平方根
> w <- c()
> for (i in 1:n) w <- append(w, list(mu + B%*%cbind(rnorm(3))))
> return(w)
> }
> mu3 <- cbind(c(1,1,1)) # 平均ベクトル(縦ベクトル)
> A3 <- array(c(2,1,1,1,2,1,1,1,2), dim=c(3,3))
> w3 <- r3norm(mu3, A3, 2000)
> head(w3)
[[1]]
[,1]
[1,] 3.760969
[2,] 1.506681
[3,] 2.271915
---- 以下省略


では、4変量正規分布に従う乱数を生成する関数 r4norm() を定義すると
> r4norm <- function(mu, A, n) {
> U <- svd(A)$u; V <- svd(A)$v; D <- diag(sqrt(svd(A)$d))
> B <- U %*% D %*% t(V) # 行列 A の平方根
> w <- c()
> for (i in 1:n) w <- append(w, list(mu + B%*%cbind(rnorm(4))))
> return(w)
> }
> mu4 <- cbind(c(1,1,1,1)) # 平均ベクトル(縦ベクトル)
> A4 <- array(c(2,1,1,1,1,2,1,1,1,1,2,1,1,1,1,2), dim=c(4,4))
> w4 <- r4norm(mu4, A4, 2000)
> head(w4)
[[1]]
[,1]
[1,] 0.9252319
[2,] -0.4436436
[3,] -0.3861779
[4,] -0.1728883
---- 省略


この4乱数が本当に正規分布を描いているか確認してみる。

# 一つ目のデータだけ取り出す
> y<-c()
> x<-c(1:2000)
> for (i in 1:2000) y<-append(y, w4[[i]][1])
> head(y)
[1] 0.9252319 -0.6451168 3.6471124 -1.1930477 0.4600912 1.0727863

# Summaryデータから平均を取得
> summary(y)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-3.48500 0.03779 0.99330 0.98500 1.96500 5.68700
> quantile(y) # クォンタイル点
0% 25% 50% 75% 100%
-3.48457453 0.03778559 0.99330414 1.96540774 5.68672752
# 標準偏差を取得
> variance <- function(x) var(x)*(length(x)-1)/length(x)
> sqrt(variance(y))
[1] 1.389023

> # グラフ描画
> plot(y,dnorm(y, mean=0.9850483, sd=1.38937), type="n")
> curve(dnorm(x, mean=0.9850483, sd=1.38937), type="l",add=T)

正規分布1
確かにちゃんと正規分布が描かれています。納得。

関連記事
スポンサーサイト

コメントの投稿

管理者にだけ表示を許可する

コメント

プロフィール

kumagonjp2

Author:kumagonjp2
Python,Django,R,Mongo,MySQL,Struts,Spring,データマイニングなどサーバー関係のメモを残していきます。

最新記事
最新コメント
最新トラックバック
月別アーカイブ
カテゴリ
雪が3Dで降るブログパーツ ver2

マウスで見る方向変えられます

検索フォーム
RSSリンクの表示
リンク
ブロとも申請フォーム

この人とブロともになる

QRコード
QR
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。