RとRubyで統計学実験

目次へ

データを作る

一様乱数の重ねあわせ

$\left[0, 1 \right]$ の区間の一様乱数(Uniform-Distributed Random Number, 区間内のどこにも等しい確率で落ちる乱数)を $n$ 個発生させ，その平均値（標本平均）がどのように分布するかを知りたい．多数のデータが足しあわされると，ベル型の正規分布が生じることを調べたいのである．

      #!/usr/local/bin/ruby
      # add_unif.rb
      r = Random.new
      n_sampling = 4800
      size = 3 #標本の大きさ
        n_sampling.times do
        sum = 0.0
        size.times do 
          sum += r.rand()
        end
        puts sum / size
      end

上のプログラムで，標本の大きさsizeを 1 から 20 程度まで変えて実行し，その結果をグラフ化してみよう(最初のグラフのようすがちがったら，こちらを参照)．

注意: こういう作業を進めるときに困ってしまうことのひとつは，「いくつもあるファイルの名前を決めるのが面倒だ！」ということだ．この文書に貼ってある画像の名前は HTML のソースを見るとすぐにわかるので，参考にするといい（賢い作業重要！）．

これらの図を見ると，一見してサイズが小さい時には正規分布とは程遠いことが分かる．大きいとなんとなく正規分布っぽい．それでは，本当に正規分布になるとしたら，どのへんを境にして正規分布と言えるようになるのだろうか？それを調べていこう．

R でプロットする

R にはデータをプロットする機能が豊富に用意されている．ここでは単純な一次元データの分布をヒストグラムとしてプロットする機能を，手順を追って試してみよう．今，各行にデータが1個だけ書き込まれているファイルを sample.dat として進めてみる．

ターミナルで R として，コマンドラインモードで R を起動する．あるいは（より合理的に） hoge.R のようなファイルを作って，Rscript で実行する．

データを読み込んで数値型のデータ集合に変換する

 
      DT <- read.table("sample.dat",header=FALSE) # 先頭行からデータが始まるので FALSE 指定
      str(DT)  # DT の構造を確認して， V1 という名前の数値変数であることをチェック	    
      x <- DT$V1  
      summary(x)  ## データ x の最小値，第一四分位数，中央値，平均，第三四分位数，最大値が表示される

plot関数で表示してみる．
```
      plot(x)
    
```
これはわけのわからない図を作り出す．

ヒストグラムでデータを表示

   
       library("ggplot2")  ## 次の qplot を使うために ggplot2 ライブラリを読み込む
       qplot(x,geom="histogram",binwidth=0.05) ## 柱（ビン）の幅を0.05としてヒストグラムを表示

データを検定する

データ集合に関する仮説検定

R は１次元データの集合がどのような確率分布に従っているのかをテスト（検定）する機能がそなわっている．下にそれらを箇条書きにして示してある．ここで使いたいのは，正規分布のテストに使われるシャピロ・ウィルク検定だ．

シャピロ・ウィルク検定 shapiro.test(x): データが正規分布に従っているかどうかを検定する．p値が 0.1 程度よりも小さければ，正規分布から有意に外れている．
分散の違いのF検定 F.test(a,b): ２組のデータ a, b が同じ母集団から抽出されたとみなせる（等しい）分散を持つかどうかを検定する．t値が小さければ有意に異なる．
t検定 t.test(x,mu=$v$): データが t分布に従っているかどうかを検定する．$v$ は仮説で置いている平均．p値が 0.1 程度よりも小さければ合格．
2標本t検定 t.test(a,b,var.equal=TRUE): ２組のデータ a, b が同じ母集団から抽出されたとみなせる（等しい）平均値を持つかどうかを検定する．t値が小さければ有意に異なる．

手持ちのデータは正規分布に従っているか？

ヒストグラムで可視化して目視によって検定する

GNUPLOT やR を使ってヒストグラムを描くだけで，分布が正規分布に近いかどうかの判断がある程度できる．

Ｒを使って統計的検定を行なう．

最初のシミュレーションに戻って，連続一様分布に従うデータの平均が正規分布とみなせるかどうかを調べたい．上述のようにシャピロ・ウィルク検定が使える． R による上の作業で x を得たあと，次のようにすればよい．

  shapiro.test(x)

この結果，p値としてたとえば 0.1 より小さい値になったら，x は危険率 0.1 で正規分布と有意な違いはないという仮説は棄却される（この言明の仕方に注意して，勘違いしないようににすること）．

やりがちなミス

size = 1 の図が下のようになる人が多いはず．「あれ？一様分布のはずなのにぜんぜんバラバラじゃない！！」と感じるようなら，君はえらい．そしたら GNUPLOT の設定を工夫しよう．よくやらかしてしまうことなので，グラフを見るときには縦軸と横軸の数値や単位に気をつけること．

Q-Q プロットで正規分布かどうかを知る

Q-Q プロットというのは，与えられたデータ集合の各分位点の位置と，そのデータと同じ平均，分散を持つ生正規分布の分位点とを縦軸横軸にプロットしたものである．もしもデータが正規分布していれば，プロットは直線になるし，そうでないと歪むことになる．

  #!/usr/bin/env Rscript
  library(MASS) # MASS ライブラリを使う
  P <-  read.table("sample.dat",header=F)
  str(P) # P がどんな名前と型の変数を持っているか調べる(V1 で numericのはず) 
  x <- as.numeric(as.character(P$V1))
  png("NdistCheck.png")
  qqnorm(y1,main = "Main Title",ylab="hogehoge",ylim=c(750,1150))
  qqline(y1)

課題

上で説明した連続一様分布乱数の平均が作る分布を， $n = 1, 2, 3, 5, \ldots, 60$ までのケースについて作ってみて，どのあたりで正規分布と言えるようになるかを調べなさい．すべてについてやる必要はないので，適当に「間引き」して調べなさい．
ポアンカレはフランスの天才数学者の1人であり，「ポアンカレ予想」という難問を後世に贈って，現代に至るまで数学者たちをきりきりまいさせたことで名前を知られている．彼にはこんな逸話がのこされている．ただし後世の作り話らしい．
ポアンカレは，地方に単身赴任しているとき，毎日買ってきたパンの重量を量って記録をとっていた．着任から1年後，彼はそのデータを持って警察に出向き，パン屋を告発した．ポアンカレによると，1000g という表示のパンの重さの平均値は 950 g しかなかったというのである．その結果，パン屋は警察に呼ばれて絞られるはめになった．

その1件以来，ポアンカレに渡されるパンはすべて 1000 g以上のものばかりになった．しかし，ふたたび1年後，ポアンカレはデータを手に警察に出向き，パン屋を告発した．確かにパンは 1000 g 以上だったのだが，その分布をみるとパン屋が作っているパンは前と同じものであるにも関わらず，その中から 1000 g 以上のものだけをポアンカレに売っていたことが，重量の分布から明らかだったからである．
1. 最初に1年間にポアンカレが買った365個のパンの重量の分布を，正規分布乱数を使ったシミュレーションによって再現しなさい．パンの重量の平均は 950 g，標準偏差は 50 g であったとする．得られた結果のヒストグラムを描き，Q-Qプロットとウィルク・シャピロ検定を行なって正規分布性を確認しなさい．正規分布乱数の発生の仕方については，資料を参考にすること．
2. 翌年にポアンカレが買った365個のパンの重量の分布を，シミュレーションによって再現しなさい．そのあと，上と同様の検査を行いなさい.
結果の例
この課題で得られる結果の例は，次のようになるはずだ．もちろん乱数を使ったシミュレーションだから，結果は毎回異なる．
シャピロ・ウィルク検定の結果|. p 値が重要
```
	Shapiro-Wilk normality test
data:  y1
W = 0.9974, p-value = 0.8303

	Shapiro-Wilk normality test
data:  y2
W = 0.8833, p-value = 4.959e-16
```
Q-Qプロットの結果

資料

Ruby と Java の乱数発生