平均と標準偏差,範囲を指定した乱数発生方法

論文に出てる記述統計からシュミレーションをしてみようと思った。いろいろググってみたのだけれど,なかなかうまい方法がなかったのでメモ。

Rでは,様々な分布に従う乱数を発生させることができる。一番有名なのはrnormというやつ。正規分布に従う乱数を発生させる。平均と標準偏差を指定して,こんなかんじで。

>dat <-rnorm(100, mean=50, sd=5) #平均値が50で標準偏差が5の乱数を100個

> dat
[1] 47.49772 48.21094 56.25741 49.96995 44.78259 48.71722 50.17432 55.17608
[9] 44.92424 46.83734 48.06748 53.36903 46.53421 50.15734 40.72731 48.05790
[17] 55.66313 53.06009 50.12512 46.79404 54.13993 45.70641 45.05168 51.28307
[25] 50.83705 48.71708 48.01274 49.00765 56.62082 49.58216 52.19713 63.85808
[33] 48.90386 56.36748 48.44979 53.96921 43.93527 50.40558 46.51421 48.81157
[41] 39.57205 51.43287 50.73080 52.36210 44.21936 59.07660 52.05478 54.53502
[49] 52.26897 54.84681 45.80431 53.86563 53.97093 55.28294 55.68036 57.21981
[57] 48.21196 42.46615 44.65935 46.11213 52.64849 47.62270 46.36596 52.57540
[65] 47.27981 63.84164 52.44592 43.90385 50.75210 51.00874 57.80445 50.17152
[73] 47.95725 50.36983 46.72289 47.57410 44.48939 53.98110 59.12215 47.21243
[81] 42.14203 40.78201 50.64352 48.90662 53.63576 55.65404 46.23628 48.17049
[89] 51.76097 52.69611 49.41447 53.25395 52.08408 41.06224 53.22397 41.33806
[97] 58.23008 44.21629 53.06959 53.20374

ただしこのやり方では,SDが大きくなると負の数値が出てしまったりする。また,例えばテストの点数とかの乱数を作りたい場合,100より上もいらないのに100を超える数値も出る。

> dat1 <- rnorm(100, mean=60, sd=30)
> dat1
[1] 44.494760 89.170093 27.292283 66.450161 62.407602 66.882294
[7] 61.546427 67.186071 49.042632 63.244362 84.543721 60.576172
[13] 51.848288 83.655269 41.251925 116.468853 64.299882 88.146247
[19] 18.110042 47.436176 56.263604 96.849951 89.739825 23.088249
[25] 34.150477 2.516926 4.948293 56.340174 25.665198 95.580383
[31] 83.448675 113.616310 -3.453262 108.500719 56.467913 100.086726
[37] 66.637288 43.201111 59.247994 45.427133 39.968471 37.975361
[43] 16.056423 22.321241 29.502743 39.415555 68.695397 94.841580
[49] 49.784221 82.842063 53.112870 7.554488 34.891446 38.734127
[55] 105.762073 39.075296 89.386547 56.279559 28.455354 84.684350
[61] 84.719927 86.950250 72.399872 48.875458 73.908410 41.047765
[67] 38.112460 38.006715 25.781026 64.453294 63.488881 62.025456
[73] 88.706781 53.570347 71.571211 37.240862 32.510384 20.557917
[79] 48.065483 94.325533 94.798381 64.708053 108.530045 26.543721
[85] 58.603449 52.491463 85.727697 106.653768 47.484039 102.625380
[91] 49.518054 36.609159 27.157635 59.725756 39.526326 82.433037
[97] 90.007551 97.304846 29.513355 80.113184

これと同じようなことは,Excelを使うと次のような式でできる。

=NORMINV(RAND(),60, 30)

ただしこれでも先ほどのように負の値と100位上の値ができてしまう。うーむ。数学的なことがわからないので困った。Excelを使う場合,ある程度の数値ならばここに書いてあるような方法を使ってできるみたい。

http://www016.upp.so-net.ne.jp/sige-lab/before2002/1995make_rand.pdf

ただしこのやり方もSDが広すぎると対応できない。と,ぐぐっていたらこんなのをみつけた。

R言語について。正規分布に従う疑似乱数を発生させる「rnorm()」関数で、正の値だ… http://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q1177484151?fr=pc_tw_share_q #

ここで回答者さんが平均値を指定して,負の値を出さないようにする乱数発生のRスクリプトを書いてくださっている。引用します。

>zz<-rnorm(10000,mean=3.34)

>while(any(zf<-zz<0)){#ひたすら粘る(お勧めできない)#
>zz<-rnorm(10000,mean=3.34)}

>abs(zz)#原点で折り返す#

>zz[zz>=0]#負値は省く#

>(zz>0)*zz#0に圧縮#

>ifelse(zz<0,3.34,zz)*zz#平均を強制#

>while(any(zf<-zz<0)){#正値で上書き#
zz[zf]<-rnorm(sum(zf),mean=3.34)}

これのrnormの引数でSDも指定してやればいいのではないか,と。

>zz<-rnorm(100,mean=60. sd=30)

>while(any(zf<-zz<0)){#ひたすら粘る(お勧めできない)#
zz<-rnorm(100,mean=60, sd=30)}

>abs(zz)#原点で折り返す#

>zz[zz>=0]#負値は省く#

>(zz>0)*zz#0に圧縮#

>ifelse(zz<0,60,zz)*zz#平均を強制#

>while(any(zf<-zz<0)){#正値で上書き#
zz[zf]<-rnorm(sum(zf),mean=60, sd=30)}

ただしこれだけでは100以上の値が含まれるので,

>ifelse(zz>100, 60, zz)

として0以下の値の処理と同じように平均に置き換えてあげる。シュミレーションした平均値と標準偏差を出すと,こんな感じ。

> mean(zz)
[1] 65.13294
> sd(zz)
[1] 29.33622

平均が5点ほど上にずれたけどsdはまあまあ。このやり方は,範囲の外を取る値を平均値で置き換えるやり方だけれど,それらの値を境界の0と100で置き換えるということもありなんじゃないだろうか。やってみる。

> zz<-rnorm(100,mean=60,sd=30)

while(any(zf<-zz<0)){#ひたすら粘る(お勧めできない)#
zz<-rnorm(100,mean=60, sd=30)}

abs(zz)#原点で折り返す#

zz[zz>=0]#負値は省く#

(zz>0)*zz#0に圧縮#

ifelse(zz<0,0,zz)*zz#負の値を0に置き換え#

while(any(zf<-zz<0)){#正値で上書き#
zz[zf]<-rnorm(sum(zf),mean=60,sd=30)}

ifelse(zz>100, 100, zz) #100以上の値を100に置き換え#

> mean(zz)
[1] 60.1755
> sd(zz)
[1] 27.46259

標準偏差は少し小さくなったが平均はさきほどよりも近くなった。もともとこういうやり方では最初に指定したものとぴったり一致することはないにせよ,このやり方だとまあまあの精度でシュミレーションができそう。ただまあサンプルサイズがこれより少なくなると推定の精度は下がる。逆に増やせば上がる。

最後に。ここまでやってきたやり方は勝手に標本分布を正規分布と仮定してやっているのでその点は注意が必要。t検定なんかやってるんだったらまあ正規性を仮定してのことなので良いけれど,ポアソン分布とかカイ二乗分布とかだったらこのやりかたはできない。そういう分布では平均と標準偏差の意味するものが違ってくる。そういう場合も分布のパラメータがわかればシュミレーションはできるのかな。またなにか機会があったらそっちでも遊んでみよう。

そんな感じで昨日の基礎研ではシュミレーションしたデータを使って,論文のRQに答える検定を自分で考えてやってみるということをやろうとしたのでした。バタバタと準備したのでなんかつまらなくなってしまったけれど。もうちょっとちゃんと準備してやればもっと勉強にもなったよなと思うので,また基礎研で僕の番が回ってきたら同じようなことにチャレンジしてみたい。今度は,テキストファイルのデータをExcelにまとめてそれをRでもlangtestでもSPSSでなんでもいいから使って結果を解釈するところまで。

なんかあの準備不足で色々あたふたする感じが僕の普段の授業がどうなっているかを表している気がしないでもなかったり…(遠い目

というわけで,そうちゃんと戯れてこようと思います(じゅるり

なにをゆう たむらゆう

おしまい

 

広告

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中