【一様分布・二項分布・正規分布・ポアソン分布】統計モデリングで使おう

SHARE

こんにちは、たくろー(@takuro_data )です！

前回は、変数の相関について学習し、データの特徴をとらえる平均、バラツキ、相関とひと通り終えました。

【相関係数・順位相関係数】2変数の関係を把握しよう

今回は、データの背後にある母集団の特徴を標本からとらえる推測統計の基礎となる確率分布について学習します。

要点
  1. 一様分布は、各事象の起きる確率（生起確率）が等しい分布
  2. ２項分布は、成功・失敗といった事象についての分布
  3. 正規分布は、平均値を中心としたつり鐘型の分布
  4. ポアソン分布は、まれだが一定の時間内にはある程度の頻度で起こる事象の数の分布

一様分布：サイコロ

各事象の起きる確率（生起確率）が等しい分布です。変数xが区間[a,b]の間にあるときの、平均は\(\mu=\frac{a+b}{2}\)、分散は\(\sigma^2=\frac{(b-a)^2}{12}\)です。

#乱数シミュレーション
n <- c(10,100,1000,10000)
par(mfrow=c(2,2))
for(i in n){
  dice_sim <- floor(runif(i,1,7))
  hist(dice_sim,breaks=c(0,1,2,3,4,5,6),probability=T,col="gray", main = paste("n = ", i))
}
mean(dice_sim)
var(dice_sim)
sd(dice_sim)

#期待値と分散の理論値
a <- 1
b <- 6
mu <- (a+b)/2
var <- ((b-a)^2/12)
sd <- sqrt(var)
mu
var
sd
> mean(dice_sim)
[1] 3.5329
> var(dice_sim)
[1] 2.877205
> sd(dice_sim)
[1] 1.696233
> mu
[1] 3.5
> var
[1] 2.083333
> sd
[1] 1.443376

２項分布：コイン投げ

２項分布は、成功・失敗といった事象についての分布です。成功・失敗のように結果が２種類しかない試行を、ベルヌーイ試行といいます。\(n\)回の試行で\(m\)回成功する確率は、次式で計算します。平均は\(np\)で、分散は\(np(1-p)\)です。

$$_{n}C_{m}･p^m･(1-p)^{n-m}$$

乱数シミュレーション

# 乱数シミュレーション
size <- 10
p <- 0.5
n <- c(10,100,1000,10000)
par(mfrow=c(2,2))
for(i in n){
  coin_sim <- rbinom(i,size,p)
  hist(coin_sim,probability=T,col="gray",main = i)
}
mean(coin_sim)
var(coin_sim)
sd(coin_sim)

# 期待値と分散の理論値
mu <- size * p
var <- size * p * (1-p)
sd <- sqrt(var)
mu
var
sd
> mean(coin_sim)
[1] 4.9842
> var(coin_sim)
[1] 2.518002
> sd(coin_sim)
[1] 1.586821
> mu
[1] 5
> var
[1] 2.5
> sd
[1] 1.581139

繰り返しの数を増やすと、次で説明する正規分布に近づきます。

確率密度・累積分布・確率点

# 確率密度関数: size回のうちx回成功する確率
size <- 10
p <- 0.5
x <- 3
dbinom(x,size,p)
#計算確認
choose(size, x) * p^x * p^(size-x)

# 累積確率
pbinom(x,size,p)
#計算確認
ruiseki = 0
for(i in 0:x){
  ruiseki <- ruiseki + dbinom(i,size,p)
}
ruiseki

# 確率点
qbinom(ruiseki,size,p)

# 関数から描画
size <- 50
p <- 0.5
x<-0:40
par(mfrow=c(1,2))
plot(x,dbinom(x,size,p),type="h",lwd=3,col="gray")
plot(x,pbinom(x,size,p),type="h",lwd=3,col="gray")
> dbinom(x,size,p)
[1] 0.1171875
> choose(size, x) * p^x * p^(size-x)
[1] 0.1171875
> pbinom(x,size,p)
[1] 0.171875
> ruiseki
[1] 0.171875
> qbinom(ruiseki,size,p)
[1] 3

確率密度はdxxx、累積分布はpxxx、確率点はqxxxです。詳細は下記記事をご覧ください。

参考 60. 確率分布と乱数R-Source

正規分布：つり鐘型の分布

正規分布は、平均値を中心としたつり鐘型の分布です。検定などで正規分布が前提とされることが多く、統計学で最も重要な分布です。ガウス分布とも呼ばれます。平均\(\mu\)、標準偏差\(\sigma\)の確率密度関数は次式のとおりです。

$$f(x)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$$

乱数シミュレーション

#乱数シミュレーション
mu<-10
sigma<-20
par(mfrow=c(1,1))
normal_sim <- rnorm(1000, mean = mu,sd = sigma)
hist(normal_sim)

sample_mean <- mean(normal_sim)
sample_sd <- sd(normal_sim)
sample_mean
sample_sd
> sample_mean
[1] 8.616459
> sample_sd
[1] 19.79836

確率密度

#確率密度関数
norm_dist <- function(x,mu,sigma){
  1/(sqrt(2*pi)*sigma)*exp(-(x-mu)^2/(2*sigma^2))
}
norm_dist(0,0,1)
dnorm(0)
> norm_dist(0,0,1)
[1] 0.3989423
> dnorm(0)
[1] 0.3989423
#いろいろな正規分布描画
x<-seq(-4,4,0.1)
curve(dnorm(x,0,0.5),from=-4,to=4,type="l",col="blue")
curve(dnorm(x,0,1),add=T,col="pink")
curve(dnorm(x,0,1.5),add=T,col="green")
curve(dnorm(x,-2,1),add=T,col="orange")

標準正規分布

標準化とは、データの平均値を0に、標準偏差を1に変換することです。変換後のデータは標準化変量と呼ばれます。標準化した正規分布が標準正規分布（z分布）です。

#標準化:x-mu / sigma
n_sample <- (normal_sim - sample_mean) / sample_sd
hist(n_sample)
mean(n_sample)
sd(n_sample)
> mean(n_sample)
[1] 1.316568e-17
> sd(n_sample)
[1] 1
#標準正規分布よく使うタイル点
#2.5%,5%,95%,97.5%タイル点
qnorm(0.025)
qnorm(0.05)
qnorm(0.95)
qnorm(0.975)
> qnorm(0.025)
[1] -1.959964
> qnorm(0.05)
[1] -1.644854
> qnorm(0.95)
[1] 1.644854
> qnorm(0.975)
[1] 1.959964
#シグマ区間
pnorm(1) - pnorm(-1)
pnorm(2) - pnorm(-2)
pnorm(3) - pnorm(-3)
> pnorm(1) - pnorm(-1)
[1] 0.6826895
> pnorm(2) - pnorm(-2)
[1] 0.9544997
> pnorm(3) - pnorm(-3)
[1] 0.9973002

データが正規分布に従うとすると、１シグマ（標準偏差）内におさまる確率が68.3%、２シグマ内が95.4%、３シグマ内が99.7%です。

歪度と尖度

歪度と尖度は、標本の分布の形が正規分布からどの程度離れているかを測るための指標です。

歪度は、分布が左右対称か、右裾が長いか、左裾が長いか、分布のゆがみを表す指標で、次式で計算します。正の場合、右裾が長いです。

$$S_{w}=\frac{1}{n}\sum_{i=1}^n (\frac{x_{i}-\bar{x}}{s})^3$$

尖度は、分布のとんがり度合いを表す指標で、次式で計算します。正規分布は３で、３より大きいか小さいかで判断します。

$$S_{k}=\frac{1}{n}\sum_{i=1}^n (\frac{x_{i}-\bar{x}}{s})^4$$

株価の収益率がマイナスのショックが大きいため、左裾が長く、尖っている傾向があるので、RのサンプルデータEuStockMarketsで用いて実行してみます。EuStockMarketsは、欧州の主要株価指数（ドイツDAX、スイスSMI、フランスCAC、英国FTSE）の1991年〜1998年の日次データです。

## 歪度と尖度
#歪度と尖度を含むパッケージを読み込み
library(moments)

data("EuStockMarkets")
str(EuStockMarkets)
dax <- EuStockMarkets[,1]
r_dax <- diff(log(dax)) *100
hist(r_dax)
mean(r_dax)
sd(r_dax)
skewness(r_dax)
kurtosis(r_dax)
> mean(r_dax)
[1] 0.06520417
> sd(r_dax)
[1] 1.030084
> skewness(r_dax)
[1] -0.5540533
> kurtosis(r_dax)
[1] 9.279689

ポアソン分布：稀な事象の分布

試行回数がとても多く、事象発生確率（生起確率）がとても小さいときの２項分布です。不良品の数、交通事故数などに用いられます。平均\(\lambda\)（試行回数\(n\)×確率\(p\)）の事象の確率が起こる回数\(x\)の確率は、次式で計算します。分散は平均と同じ値です。

$$f(x)=\frac{e^{-\lambda}\lambda^x}{x!}$$

乱数シミュレーション

#乱数シミュレーション
n <- c(100,1000,5000,10000)
p <- 0.002
#lambda <- n*p
par(mfrow=c(2,2))
for (i in n) {
  pois_sim <- rpois(i,i*p)
  hist(pois_sim,probability=T,col="gray", main = paste("n=",i))
}
mean(pois_sim)
sd(pois_sim)

mu <- n[4]*p
var <- mu
sd <- sqrt(var)
mu
sd
> mean(pois_sim)
[1] 20.0253
> sd(pois_sim)
[1] 4.472568
> mu
[1] 20
> sd
[1] 4.472136

\(n\)が増える（\(p\)が大きくなる）と、\(\lambda\)が大きくなり、平均・分散も大きくなり、正規分布に近づきます。

確率密度

#確率密度関数
pois_dist <- function(x,lambda){
  exp(-lambda)*lambda^x/factorial(x)
}
lambda <- 20
x <- 15
pois_dist(x,lambda)
dpois(x,lambda = lambda)
> pois_dist(x,lambda)
[1] 0.05164885
> dpois(x,lambda = lambda)
[1] 0.05164885

さいごに

今回は、よく使う確率分布として一様分布、２項分布、正規分布、ポアソン分布を扱いました。

次回は、検定でよく使う確率分布として、カイ二乗分布、F分布、t分布を扱います。

【カイ二乗分布・F分布・t分布】仮説検定で使おう

Rで統計学入門

前の記事

【相関係数・順位相関係数】2変数の関係を把握しよう

次の記事

【カイ二乗分布・F分布・t分布】仮説検定で使おう