別のブログで使った大数の法則を示すシミュレーションのソースコードです。
# 乱数から分析データを作る
set.seed(20150209)
x.sd <- 16.96534 # 想定標準偏差
n <- 24 # 分析年数
x <- as.numeric(n*12) # 分析データ
for(i in seq(1, (n-1)*12, 12)){
x[i:(i+11)] <- c(rnorm(3, sd=x.sd), rnorm(1, sd=x.sd*2), rnorm(8, sd=x.sd))
}
# 分析データの状態を確認
plot(x, type="l", main="集計前データ", xlab="", ylab="")
# n年分を合計してプロットする
x.sum <- numeric(12)
for(i in 1:(n-1)){
x.sum <- x.sum + abs(x[(i*12-11):(i*12)])
}
max <- (as.integer(max(x.sum)/n/5)+1)*5
B <- barplot(x.sum/n, names.arg=sprintf("%d月", 1:12), ylim = c(0, max), main=sprintf("シミュレーション（%d年分）", n))
# 偏差値を計算してみる
y.mean <- mean(x.sum/n)
y.sd <- sd(x.sum/n)
sprintf("4月の偏差値: %.2f", ((x.sum/n)[4]-y.mean)/y.sd*10+50)
