ここでは、R言語で確率分布をプロットしたり、確率分布からの乱数を生成したりする方法を解説していきます。
確率密度や確率を算出するための関数さらに乱数を生成する関数を紹介し、関数の使用例をまとめています。
今回紹介するプログラミングコードも以下のzipファイル中のRスクリプトに保存しています。
ポチップ
確率分布のための関数
確率密度・質量、分布関数、確率点、乱数の算出
確率分布の各点の確率密度や分布関数、確率点、乱数を算出する関数を以下にまとめました。
表の左が確率分布で右側が各分布に対応する関数になります。
なお、関数の引数についてですが全て説明すると膨大になるので省略します。関数の引数は各分布のパラメータを指定するものがあり、様々なパラメータに関する確率分布の密度・質量、分布関数、確率点、乱数を算出することができます。
| 分布 | 確率密度・確率質量 | 分布関数 | 確率点 | 乱数 |
| ベータ分布 | dbeta() | pbeta() | qbeta() | rbeta() |
| 二項分布 | dbinom() | pbinom() | qbinom() | rbinom() |
| コーシー分布 | dcauchy() | pcauchy() | qcauchy() | rcauchy() |
| カイ二乗分布 | dchisq() | pchisq() | qchisq() | rchisq() |
| 指数分布 | dexp() | pexp() | qexp() | rexp() |
| F分布 | df() | pf() | qf() | rf() |
| ガンマ分布 | dgamma() | pgamma() | qgamma() | rgamma() |
| 幾何分布 | dgeom() | pgeom() | qgeom() | rgeom() |
| 超幾何分布 | dhyper() | phyper() | qhyper() | rhyper() |
| ロジスティック分布 | dlogis() | plogis() | qlogis() | rlogis() |
| 対数正規分布 | dlnorm() | plnorm() | qlnorm() | rlnorm() |
| 負の二項分布 | dnbinom() | pnbinom() | qnbinom() | rnbinom() |
| 正規分布 | dnorm() | pnorm() | qnorm() | rnorm() |
| ポアソン分布 | dpois() | ppois() | qpois() | rpois() |
| t分布 | dt() | pt() | qt() | rt() |
| スチューデント化された分布 | dtukey() | ptukey() | qtukey() | rtukey() |
| 一様分布 | dunif() | punif() | qunif() | runif() |
| ワイブル分布 | dweibull() | pweibull() | qweibull() | rweibull() |
| ウィルコクソンの順位和統計量の分布 | dwilcox() | pwilcox() | qwilcox() | rwilcox() |
| ウィルコクソンの符号和統計量の分布 | dsignran() | psignran() | qsignran() | rsignran() |
関数の使用例
上に記載した関数の使いかたを紹介していきます。
ここでは、確率密度及び確率質量を算出する関数（dnormやdunif）と乱数を生成する関数（rnormやrunif）の実行例について触れていきます。
正規分布
関数dnormを使って正規分布の確率密度関数を描画する方法を以下にまとめました。
1 2 3 4 5 6 7 | #正規分布の確率密度関数 curve(dnorm(x, mean = 0, sd = 1), xlim = c(-10, 10) , xlab = "", ylab = "pdf", col = "firebrick") par(new = TRUE) curve(dnorm(x, mean = 0, sd = 2), xlim = c(-10, 10) , xlab = "", ylab = "pdf", xaxt = "n", yaxt = "n", col = "royalblue") par(new = TRUE) curve(dnorm(x, mean = 0, sd = 3), xlim = c(-10, 10) , xlab = "", ylab = "pdf", xaxt = "n", yaxt = "n", col = "green") legend("topleft", legend = c("μ=0, σ=1", "μ=0, σ=2", "μ=0, σ=3"), col = c("firebrick", "royalblue", "green"), lty = c(1, 1, 1)) |
上から平均0、分散1の正規分布、平均0、分散4の正規分布、平均0、分散9の正規分布を描画しています。
次の例は、関数rnormから生成される正規乱数が実際に正規分布に従っているかをプロットしています。
1 2 3 4 5 6 7 | #正規乱数 par(mfrow = c(1, 3)) sampleSizes <- c(10, 100, 1000) for (n in sampleSizes) { hist(rnorm(n, 0, 1), xlim = c(-4, 4), xlab = "") } par(mfrow = c(1, 1)) |
サンプルサイズを大きくすることでヒストグラムがより正規分布の形に近づいていることが見て取れます。
一様分布
次に一様分布に関する例をみていきます。
次のようにして一様分布の確率密度関数をプロットすることが可能です。
1 2 3 4 5 6 7 | #一様分布の確率密度関数 curve(dunif(x, 0, 1), xlim = c(-1, 4), ylim = c(0, 1.5), xlab = "", ylab = "pdf", col = "firebrick") par(new = TRUE) curve(dunif(x, 0, 2), xlim = c(-1, 4), ylim = c(0, 1.5), xlab = "", ylab = "pdf", xaxt = "n", yaxt = "n", col = "royalblue") par(new = TRUE) curve(dunif(x, 0, 3), xlim = c(-1, 4), ylim = c(0, 1.5), xlab = "", ylab = "pdf", xaxt = "n", yaxt = "n", col = "green") legend("topleft", legend = c("a=0, b=1", "a=0, b=2", "a=0, b=3"), col = c("firebrick", "royalblue", "green"), lty = c(1, 1, 1)) |
各区間に対応する一様分布を描いていることが分かります。
次に、一様乱数を生成する関数runifの使用例についてみていきます。
次のように一様乱数を用いることで円周率を近似的に求めることができます。
1 2 3 4 5 6 7 8 9 10 11 12 13 | #モンテカルロ法・円周率 par(mfrow = c(1, 3)) sampleSizes <- c(100, 1000, 10000) for (n in sampleSizes) { x <- runif(n, 0, 1) y <- runif(n, 0, 1) distances <- sqrt(x^2 + y^2) surface <- 4 * length(distances[distances < 1]) / n plot(x[distances < 1], y[distances < 1], xlab = "x", ylab = "y", col = "firebrick", sub = paste0("π=", surface)) par(new = TRUE) plot(x[distances >= 1], y[distances >= 1], xlab = "", ylab = "", xaxt = "n", yaxt = "n", col = "royalblue") } par(mfrow = c(1, 1)) |
一様乱数を用いて、面積が1の正方形中の座標をランダムに抽出します。半径1以内の円内（ここでは扇内）の座標の数と正方形の面積との比を用いることで円周率を求めています。サンプル数を大きくすると3.14に近づいていることが分かります。
モンテカルロ法を用いた円周率のシミュレーションゲームは次の記事で紹介しています。 円周率を近似するシミュレーションゲームを一日で作ってみました。 普段はR言語を用いた統計解析を解説していますが、学校の課題でアプリの作成をすることになったので一日で作りました。 一日で作成したのでクオ ... 続きを見る
円周率シミュレーターを一日で作ってみた【Unity】
次の例は、標本平均の分散の収束をプロットしています。
1 2 3 4 5 6 7 | #標本平均の収束 xbars <- NULL for (n in seq_len(500)) { x <- runif(n, 0, 1) xbars <- append(xbars, mean(x)) } plot(xbars, xlab = "index", ylab = "sample mean", type = "l", col = "royalblue") |
サンプルサイズが大きくなることで標本平均の分散が小さくなることが分かります。
ポチップ二項分布
次に二項分布に関する関数の使いかたを紹介していきます。
次のようにして各パラメータの二項分布の確率関数を描画することが可能です。
1 2 3 4 5 6 7 | #二項分布の確率関数 plot(0 : 50, dbinom(0 : 50, size = 50, prob = 0.5), xlim = c(0, 50) , xlab = "", ylab = "pmf", col = "firebrick") par(new = TRUE) plot(0 : 50, dbinom(0 : 50, size = 50, prob = 0.7), xlim = c(0, 50) , xlab = "", ylab = "pmf", xaxt = "n", yaxt = "n", col = "royalblue") par(new = TRUE) plot(0 : 50, dbinom(0 : 50, size = 50, prob = 0.9), xlim = c(0, 50) , xlab = "", ylab = "pmf", xaxt = "n", yaxt = "n", col = "green") legend("topleft", legend = c("p=0.5, n=50", "p=0.7, n=50", "p=0.9, n=50"), col = c("firebrick", "royalblue", "green"), lty = c(1, 1, 1)) |
次に関数rbinomの使用例を紹介します。パラメータp=0.5の二項分布のnを大きくすると、正規分布に近づくことをプロットで確かめています。
1 2 3 4 5 6 7 | #二項乱数 par(mfrow = c(1, 3)) sampleSizes <- c(250, 500, 750) for (n in sampleSizes) { hist(rbinom(n, size = n, prob = 0.5), ylab = paste0("histogram (n=", n, ")")) } par(mfrow = c(1, 1)) |
サンプルサイズnが750を超えるとヒストグラムはほぼ正規分布になっています。
ポアソン分布
最後にポアソン分布に関する関数を紹介します。
次のようにしてポアソン分布の確率関数を描くことができます。
1 2 3 4 5 6 7 | #ポアソン分布の確率関数 plot(0 : 50, dpois(0 : 50, lambda = 1), xlim = c(0, 50), ylim = c(0, 0.5), xlab = "", ylab = "pmf", col = "firebrick") par(new = TRUE) plot(0 : 50, dpois(0 : 50, lambda = 5), xlim = c(0, 50), ylim = c(0, 0.5), xlab = "", ylab = "pmf", xaxt = "n", yaxt = "n", col = "royalblue") par(new = TRUE) plot(0 : 50, dpois(0 : 50, lambda = 9), xlim = c(0, 50), ylim = c(0, 0.5), xlab = "", ylab = "pmf", xaxt = "n", yaxt = "n", col = "green") legend("topright", legend = c("λ=1", "λ=5", "λ=9"), col = c("firebrick", "royalblue", "green"), lty = c(1, 1, 1)) |
関数rpoisの使いかたではないですが、rbinomの引数pを小さくすることで二項分布がポアソン分布に近づくことを確かめています。
1 2 3 4 5 6 7 | #二項乱数・ポアソン分布への収束 plot(0 : 100, dbinom(0 : 100, size = 100, prob = 0.2), xlim = c(0, 50), ylim = c(0, 0.5), xlab = "", ylab = "pmf", col = "firebrick") par(new = TRUE) plot(0 : 100, dbinom(0 : 100, size = 100, prob = 0.1), xlim = c(0, 50), ylim = c(0, 0.5), xlab = "", ylab = "pmf", xaxt = "n", yaxt = "n", col = "royalblue") par(new = TRUE) plot(0 : 100, dbinom(0 : 100, size = 100, prob = 0.01), xlim = c(0, 50), ylim = c(0, 0.5), xlab = "", ylab = "pmf", xaxt = "n", yaxt = "n", col = "green") legend("topright", legend = c("p=0.2", "p=0.1", "p=0.01"), col = c("firebrick", "royalblue", "green"), lty = c(1, 1, 1)) |
パラメータpを小さくすることで、二項分布の確率関数が上で紹介したポアソン分布の確率関数に近づいていることが分かります。
まとめ
R言語で確率密度を得たり乱数を発生させたりする方法を紹介しました。
この記事で見てきたように、簡単に確率分布をプロットしたり確率点を算出することが可能です。
また、任意の分布の乱数を簡単に発生させることが可能であり、シミュレーションで非常に役立ちます。
