ベイズ統計 第1回演習

  • ここでは、慶應大学SFCでの開講科目「ベイズ統計」での演習で用いている分析例やRコード・JAGSコードなどを、授業履修者のために示しています。
  • 必ずしも全てのバージョンのRやOSで動作確認を行っていません。
  • Rパッケージやデータの更新などにより、本書執筆時とはデータやRコードを予告無しに変更する場合があります。 本書記載のR分析例以外は、R最新版等で動作確認できていないRコードも存在しますが、参考のために記載しておきます(旧バージョンで動作することがあります)。

第1回R演習の内容

  • Rを使ったことがない人は、まずRをダウンロードしインストールしてください。
  • 今回はまず、確率論の基本的な試行の例として用いられる、コイントスとサイコロのシミュレーションを行います。それにより、理論的確率と経験的確率の違いを改めて理解します。
  • 次に、授業で出てきた様々な確率分布を表示し、その性質を理解します。ここでは、一様分布、二項分布、ベータ分布、正規分布、ガンマ分布を扱います。

Rとパッケージのインストール

  • RはCRANからダウンロードし、インストール可能です。
  • インストールと基本的な操作は、こちらも参考にできます。
  • ベイズ統計関連パッケージは、CRAN Task View: Bayesian Inferenceを参考にしてください。

コイントス

  • 100回のコイントス
    # コイントス回数
    N <- 100
    # 確率0.5で表と裏がでる確率をN回ランダムに行う
    cointoss1 <- sample(x=c(0,1), prob=c(0.5,0.5), size=N, replace=TRUE)
    # 全てのデータを見てみる
    cointoss1
    # 最初の20回分をみてみる
    head(cointoss1, 20)
    # 最後の20回分をみてみる
    tail(cointoss1, 20)
    # 図をプロットする際の日本語フォントを指定
    par(family="HiraMaruProN-W4")
    # plot(x, y, ...):図をプロット
    # cumsum():累積和
    plot((1:N), cumsum(cointoss1)/(1:N) , 
    # 点をoで描画
    type="o",
    # x軸とy軸の範囲
    xlim=c(1,N) , ylim=c(0.0,1.0),
    # 図タイトルおよびx軸とy軸の名前
    main="コイントスで表が出る確率 (N=100)" , xlab="コイントス回数" , ylab="表が出る確率",
    # xy軸、xy軸名前、図タイトルの大きさの指定
    cex.axis=1.2, cex.lab=1.2, cex.main=1.5) 
    # 確率0.5の線の描画
    # col="":色を指定
    lines(c(1,N),c(0.5,0.5), col="red", cex=1.2)
    head(cumsum(cointoss1)/(1:N), 20)
    tail(cumsum(cointoss1)/(1:N), 20)
  • 10000回のコイントス
    N <- 10000
    cointoss2 <- sample(x=c(0,1), prob=c(0.5,0.5), size=N, replace=TRUE)
    par(family="HiraMaruProN-W4")
    plot((1:N), cumsum(cointoss2)/(1:N) ,  type="o", log="x", xlim=c(1,N) , ylim=c(0.0,1.0),
    main="コイントスで表が出る確率 (N=10000)" , xlab="コイントス回数" , ylab="表が出る確率",
    cex.axis=1.2, cex.lab=1.2, cex.main=1.5) 
    lines(c(1,N),c(0.5,0.5), col="red", cex=1.2)
    head(cumsum(cointoss2)/(1:N), 20)
    tail(cumsum(cointoss2)/(1:N), 20)
  • 10000回のコイントスを5回行った結果を重ねてみる
    N <- 10000
    cointoss2 <- sample(x=c(0,1), prob=c(0.5,0.5), size=N, replace=TRUE)
    par(family="HiraMaruProN-W4")
    plot((1:N), cumsum(cointoss2)/(1:N) ,  type="l", xlim=c(1,N) , ylim=c(0.0,1.0),
    main="コイントスで表が出る確率 (N=10000)" , xlab="コイントス回数" , ylab="表が出る確率",
    cex.axis=1.2, cex.lab=1.2, cex.main=1.5, col=1, lwd=5) 
    cointoss2 <- sample(x=c(0,1), prob=c(0.5,0.5), size=N, replace=TRUE)
    lines((1:N), (cumsum(cointoss2)/(1:N)), col=2, lwd=3)
    cointoss2 <- sample(x=c(0,1), prob=c(0.5,0.5), size=N, replace=TRUE)
    lines((1:N), (cumsum(cointoss2)/(1:N)), col=3, lwd=1.5)
    cointoss2 <- sample(x=c(0,1), prob=c(0.5,0.5), size=N, replace=TRUE)
    lines((1:N), (cumsum(cointoss2)/(1:N)), col=4, lwd=1)
    cointoss2 <- sample(x=c(0,1), prob=c(0.5,0.5), size=N, replace=TRUE)
    lines((1:N), (cumsum(cointoss2)/(1:N)), col=5, lwd=1)

サイコロ投げ

  • 100回のサイコロ投げ
    N=100
    # seq():連続した整数
    # rep(x,n): xをn回繰り返す
    diceroll <- sample(x=seq(1,6),prob=rep(1/6,6),size=N,replace=TRUE)
    head(diceroll,20) # 最初の20個を表示する
    diceroll1 <- diceroll[diceroll==1] # 1の目が出たものを取り出す
    diceroll1 # すべて表示してみる
    length(diceroll1) # 個数を数える
    table(diceroll) # クロス集計
    prop.table(table(diceroll)) # 割合を計算
    # ヒストグラムで度数を表示
    par(family="HiraMaruProN-W4")
    hist(diceroll, breaks=seq(0,6), freq=TRUE, ylim=c(0,50), 
    main="サイコロ投げ 100回", xlab="サイコロの目", ylab="度数")
    lines(c(0,6),c(1/6*N,1/6*N), col="red", cex=1.2)
    text(5.5,1/6*N+1, "理論確率", col="red")
    # ヒストグラムで確率を表示
    par(family="HiraMaruProN-W4")
    hist(diceroll, breaks=seq(0,6), freq=FALSE, ylim=c(0,0.4), 
    main="サイコロ投げ 100回", xlab="サイコロの目", ylab="確率")
    lines(c(0,6),c(1/6,1/6), col="red", cex=1.2)
    text(5.5,0.18, "理論確率", col="red")
  • 10000回のサイコロ投げ
    N=10000
    diceroll <- sample(x=seq(1,6),prob=rep(1/6,6),size=N,replace=TRUE)
    # head(diceroll,20) # 最初の20個を表示する
    table(diceroll) # クロス集計
    prop.table(table(diceroll)) # 割合を計算
    # ヒストグラムで度数を表示
    par(family="HiraMaruProN-W4")
    hist(diceroll, breaks=seq(0,6), freq=TRUE, ylim=c(0,(1/6*N+1000)), 
    main="サイコロ投げ 10000回", xlab="サイコロの目", ylab="度数")
    lines(c(0,6),c(1/6*N,1/6*N), col="red", cex=1.2)
    # ヒストグラムで確率を表示
    par(family="HiraMaruProN-W4")
    hist(diceroll, breaks=seq(0,6), freq=FALSE, ylim=c(0,0.3), 
    main="サイコロ投げ 10000回", xlab="サイコロの目", ylab="確率")
    lines(c(0,6),c(1/6,1/6), col="red", cex=1.2)
    text(5.5,0.18, "理論確率", col="red")
  • ガリレオのサイコロ問題をシミュレーションで確かめてみる
    N=10000
    dice.G <- matrix(0, ncol=4, nrow=N)
    colnames(dice.G)=c("dice1", "dice2", "dice3", "Total")
    dice.G[,1] <- sample(x=seq(1,6),prob=rep(1/6,6),size=N,replace=TRUE) #1個目のサイコロ
    dice.G[,2] <- sample(x=seq(1,6),prob=rep(1/6,6),size=N,replace=TRUE) #2個目のサイコロ
    dice.G[,3] <- sample(x=seq(1,6),prob=rep(1/6,6),size=N,replace=TRUE) #3個目のサイコロ
    dice.G[,4] <- rowSums(dice.G[,1:3]) # 行和
    table(dice.G[,4]) # クロス集計
    prop.table(table(dice.G[,4])) # 割合を計算

一様分布

  • 区間[0, 1]で1となる一様分布
    x <- seq(-1,2,length.out=1000)
    par(family="HiraMaruProN-W4")
    plot(x, dunif(x, min=0,max=1), xlim=c(-1,2), ylim=c(0,2), type="l", main="一様分布")
  • 区間[0, 0.5]で2となる一様分布
    x <- seq(-1,2,length.out=1000)
    par(family="HiraMaruProN-W4")
    plot(x, dunif(x, min=0,max=0.5), xlim=c(-1,2), ylim=c(0,2), type="l", main="一様分布")

二項分布

  • n=10回の試行でs=5 (r=0.5)となるベルヌーイ試行をシミュレートする⇒r=0.5となるベルヌーイ試行を10回乱数生成する
    rbinom(10,1,0.5)
  • n=10回の試行でs=1 (r=0.1)となるベルヌーイ試行
    par(family="HiraMaruProN-W4")
    plot(dbinom(seq(0,10,1), 10, 0.1), type="l", ylab="確率")
    text(3, 0.37, labels="B(10,0.1)")
  • n=10回の試行でs=3 (r=0.3)となるベルヌーイ試行
    par(family="HiraMaruProN-W4")
    plot(dbinom(seq(0,10,1), 10, 0.3), type="l", ylab="確率")
    # lines(dbinom(seq(0,10,1), 10, 0.3))
    text(6, 0.2, labels="B(10,0.3)")
  • n=10回の試行でs=5 (r=0.5)となるベルヌーイ試行
    par(family="HiraMaruProN-W4")
    plot(dbinom(seq(0,10,1), 10, 0.5), type="l", ylab="確率")
    # lines(dbinom(seq(0,10,1), 10, 0.5))
    text(8, 0.2, labels="B(10,0.5)")
  • n=10回の試行でs=7 (r=0.7)となるベルヌーイ試行
    par(family="HiraMaruProN-W4")
    plot(dbinom(seq(0,10,1), 10, 0.7), type="l", ylab="確率")
    # lines(dbinom(seq(0,10,1), 10, 0.7)) 
    text(6, 0.2, labels="B(10,0.7)")

ベータ分布

  • ベータ分布Β(0.5, 0.5)
    x <- seq(0,1,length=100)
    par(family="HiraMaruProN-W4")
    curve(dbeta(x, 0.5,0.5), type="l", xlim=c(0,1), ylim=c(0, 8), ylab="確率", lwd=2)
  • Β(0.5, 1) ⇒ 単調減少分布となる
    curve(dbeta(x, 0.5, 1), type="l", xlim=c(0,1), col="black", add=TRUE, lty=2, lwd=2)
  • Β(1, 2) ⇒ 単調増加分布となる
    curve(dbeta(x, 1, 2), type="l", xlim=c(0,1), col="grey", add=TRUE, lwd=2)
  • Β(1, 1) ⇒ 一様分布となる
    curve(dbeta(x, 1, 1), type="l", xlim=c(0,1), col="red", add=TRUE, lwd=2)
  • Β(5, 5) ⇒ 0.5をピークとする釣鐘型の分布となる
    curve(dbeta(x, 5, 5), type="l", xlim=c(0,1), col="blue", add=TRUE, lwd=2)
  • Β(10, 2)
    curve(dbeta(x, 10, 2), type="l", xlim=c(0,1), col="green", add=TRUE, lwd=2)
  • Β(10, 12)
    curve(dbeta(x, 10, 12), type="l", xlim=c(0,1), col="yellow",  add=TRUE, lwd=2)
  • 凡例
    legend(.7,8,c("Β(0.5, 0.5)","Β(0.5, 1)","Β(1, 2)","Β(1, 1)","Β(5, 5)","Β(10, 2)","Β(10, 12)"), 
    col=c("black","black","grey","red","blue","green","yellow"),lwd=c(2,2,2,2,2,2,2),  lty=c(1,2,1,1,1,1,1))

ポアソン分布

  • ポアソン分布Po(1)
    par(family="HiraMaruProN-W4")
    plot(dpois(0:50, lambda=1), type="b", xlim=c(0,20), lwd=2, ylab="確率", xlab="x")
    text(3.5, 0.35, "Po(1)")
  • Po(3)
    lines(dpois(0:50, lambda=3), type="b", lwd=2, col=2)
    text(5, 0.24, "Po(3)", col=2)
  • Po(5)
    lines(dpois(0:50, lambda=5), type="b", lwd=2, col=3)
    text(7, 0.19, "Po(5)", col=3)
  • Po(10)
    lines(dpois(0:50, lambda=10), type="b", lwd=2, col=4)
    text(11, 0.14, "Po(10)", col=4)
  • 2013年度サッカーJリーグ1部の全306試合の得点分布をポアソン分布で近似した場合、Po(1.44)となる
    # 306試合×2チーム=612チーム分のポアソン分布
    J2013.p <- dpois(0:10, lambda=1.44)*612
    # 実際の得点分布
    J2013.s <- c(132,227,154,66,23,6,4,0,0,0,0)
    # 図のプロット
    par(family="HiraMaruProN-W4")
    plot(J2013.p, type="h", xlim=c(0,11), ylim=c(0,250), lwd=3, xlab="得点", ylab="チーム数", xaxt="n")
    # 自分でx軸のラベルを付ける⇒plot()でxaxt="n"としておく
    axis(1, at=1:11, labels=0:10)
    # 別のplotを重ねあわせる
    par(new=TRUE)
    # 実際の得点分布をプロット
    plot(J2013.s, type="b", xlim=c(0,11), ylim=c(0,250), col="blue", xlab="", ylab="", xaxt="n", yaxt="n")

正規分布

  • 区間[-100, 100]で平均0, 標準偏差20の正規分布
    x <- seq(-100,100,length=2000)
    plot(x, dnorm(x, 0,20), type="l", ylim=c(0,0.02), xlim=c(-100,100), lwd=2)
  • 平均10, 標準偏差20の正規分布
    curve(dnorm(x, 10,20), type="l", ylim=c(0,0.02), xlim=c(-100,100), lwd=2, col="red", add=TRUE)
  • 平均0, 標準偏差40の正規分布
    curve(dnorm(x, 0,40), type="l", ylim=c(0,0.02), xlim=c(-100,100), lwd=2, col="blue", add=TRUE)
  • 凡例
    legend(30,0.02,c("N(0,20^2)","N(10,20^2)","N(0,40^2)"), col=c("black","red","blue"),lwd=c(2,2,2))
  • NFL2011-2013シーズンのホームチームとアウェーチムの得失点差が平均5.42、標準偏差13.54の正規分布に従うとすると
    x <- seq(-40, 60, length=2000)
    # 確率分布
    plot(x, dnorm(x, 5.42, 13.54), type="l", ylim=c(0,0.04), xlim=c(-40, 60))
  • NFL2011-2013シーズンのホームチームとアウェーチムの得失点差のデータをここからダウンロードして下さい。ユーザ名とパスワードは授業中にアナウンスします。
    # データの読み込み
    NFLscore <- read.table("NFLscore.csv", sep=",", header=TRUE)
    NFL <- c(NFLscore[,1],NFLscore[,2],NFLscore[,3])
    # 度数によるヒストグラム
    par(family="HiraMaruProN-W4")
    hist(NFL,breaks=seq(-40,65,5), main="NFL2011-2013シーズン得失点差分布(Home-Away)",
    xlab="得失点差", ylab="試合数", xlim=c(-40,60), ylim=c(0,150), col="grey")
    hist(rnorm(768,5.42,13.44), add=TRUE, border="red", breaks=seq(-40,65,5))
    # 相対度数によるヒストグラム
    par(family="HiraMaruProN-W4")
    hist(NFL,breaks=seq(-40,65,5), main="NFL2011-2013シーズン得失点差分布(Home-Away)",
    xlab="得失点差", ylab="試合数", xlim=c(-40,60), ylim=c(0,0.05), col="grey", probability=TRUE)
    lines(x,dnorm(x,5.42,13.54), col="red") 
  • 標準正規分布N(0,1)
    x <- seq(-3,3,by=0.01)
    par(family="HiraMaruProN-W4")
    plot(x, dnorm(x, 0, 1), type="l", main="標準正規分布", ylab="確率")
    # y=0 ライン
    abline(h=0)
    # 区間[-1,1]を表示
    reg.x <- x[x>=-1&x<=1]
    reg.y <- dnorm(reg.x, 0, 1)
    pol.x <- c(reg.x[1],reg.x,tail(reg.x,1))
    pol.y <- c(0, reg.y,0)
    polygon(pol.x, pol.y, density=8, col="red")

ガンマ分布

  • ガンマ分布Γ(1,1)
    # Γ(α, λ) ⇒ dgamma(x, shape, rate, scale=1/rate)
    # α⇒shape, λ=rate
    x <- seq(0,10,by=0.01)
    # α=1, λ=1⇒指数関数と同じになる
    par(family="HiraMaruProN-W4")
    plot(x, dgamma(x, 1, 1), type="l", lwd=2, ylab="確率") 
  • Γ(1,3)
    curve(dgamma(x, 1, 3), type="l", lwd=2) 
  • Γ(5, 1)
    x <- seq(0,20,by=0.01)
    plot(x, dgamma(x, 5, 1), type="l", lwd=2)
  • Γ(40,1)
    plot(x, dgamma(x, 40,1), type="l",  lwd=2)
  • Γ(2, 0.5)
    x <- seq(0,20,by=0.01)
    plot(x, dgamma(x, 2, 0.5), type="l", lwd=2)

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2014-06-23 (月) 10:13:54 (1947d)