# 尖度を求める関数を定義しておく。 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 個の尖度の平均値を計算する