いろいろな分布の尖度

統計学用語辞典

いろいろな分布の尖度

• 左右対称な分布を考える。
• 平均値を０，標準偏差が１となるように調整する。
• 正規分布より尖っている分布としてロジスティック分布がある。
• 正規分布より平らな分布として一様分布がある。
• 三角分布も正規分布より平らである。
• これらの分布の確率密度関数のグラフを描くには，R が便利である。
• 三角分布，一様分布の確率密度関数は R には用意されていないが，簡単にかける。
• これらの分布の尖度を計算するのは面倒なので，シミュレーションの結果を援用する。
  三角分布の尖度は -0.6，一様分布の尖度は -1.2 というのは簡単な積分で計算できるが，ロジスティック分布の尖度の計算が面倒だった。
  ロジスティック分布の尖度は 1.2 である。

# 尖度を求める関数を定義しておく。 method="ordinary" を使う。

kurt <- function(x, method = c("Excel", "ordinary"))
{
	method <- match.arg(method)	# 省略可能なパラメータの処理。標準は Excel（SPSS） と同じ計算法
	n <- length(x)			# データ数
	if (method == "Excel") {		# Excel（SPSS）と同じ計算法
		n*(n+1)*sum(scale(x)^4)/(n-1)/(n-2)/(n-3)-3*(n-1)^2/(n-2)/(n-3)	# scale は元のデータを標準化する関数
	}
	else {
		sum(((x-mean(x))/sqrt((n-1)*var(x)/n))^4)/n-3	# 標準化は分散の平方根を取った標準偏差による
	}
}

# シミュレーション本体。
# 各分布に従う乱数を n 個（デフォールトで 3000 個）発生させ，尖度を求める。

sim <- function(n=3000, plt= FALSE)
{
	x <- seq(-3.5, 3.5, 0.1)
	
	# ロジスティック分布
	y <- dlogis(x, scale=sqrt(3)/pi)	# 確率密度
	d <- rlogis(n, scale=sqrt(3)/pi)	# 乱数
	k <- kurt(d, method="ordinary")		# 尖度
	
	# 標準正規分布
	y2 <- dnorm(x) 				# 確率密度
	d2 <- rnorm(n) 				# 乱数
	k2 <- kurt(d2, method="ordinary")	# 尖度
	
	# 三角分布
	d3 <- (runif(n)-runif(n))*sqrt(6)	# 乱数
	k3 <- kurt(d3, method="ordinary")	# 尖度
	
	# 一様分布
	d4 <- runif(n, min=-sqrt(3), sqrt(3))	# 乱数
	k4 <- kurt(d4, method="ordinary")	# 尖度

	if(plt) { # plt が TRUE ならグラフを描く
		plot(x, y, type="l")
		lines(x, y2, type="l", col="red")
		lines(c(-sqrt(6), 0, sqrt(6)), c(0, 1/sqrt(6), 0), type="l", col="blue")
		lines(c(-sqrt(3), -sqrt(3), sqrt(3), sqrt(3)), c(0, 1/sqrt(12), 1/sqrt(12), 0), type="l", col="magenta")
		abline(h=0, v=0)
		op <- par(cex=1.4)
		text(1, 0.45, paste("logistic:", k), pos=4, col="black")
		text(1, 0.40, paste("normal:", k2), pos=4, col="red")
		text(1, 0.35, paste("triangular:", k3), pos=4, col="blue")
		text(1, 0.30, paste("uniform:", k4), pos=4, col="magenta")
		par(op)
	}
	c(k, k2, k3, k4) # 尖度を要素として持つベクトルを結果として返す
}

# シミュレーションの実行

n <- 500 # 500 セット行う
result <- matrix(0, ncol=4, nrow=n)	# 結果の保管場所を確保
for (i in 1:n) result[i,] = sim()	# 各行に毎回の結果を保管
apply(result, 2, mean)			# n 個の尖度の平均値を計算する