無作為な標本(データ)を得ようとするために,過剰な(不要な)操作をしすぎ。
pas <- runif(10,0.2,0.8)
pbs <- runif(10,0.2,0.8)
s <- sample(1:10,1)
pa <- pas[s]
pb <- pbs[s]
というのは結局は
pa <- runif(1, 0.2, 0.8)
pb <- runif(1, 0.2, 0.8)
と同じ。
つまり,
無限母集団である 0.2~0.8 の範囲の一様乱数を 10 個得て,
その中の無作為な 1 個 s <- sample(1:10,1) を取り出す pa <- pas[s]
というのは,
無限母集団である 0.2~0.8 の範囲の一様乱数を 1 個得るのと同じであるということ。
シミュレーションは,現実を「シミュレート」するものではあるが,馬鹿丁寧に[シミュレート]する必要はない。
まあ,実行時間は問題になるようなレベルの話ではないが,プログラムを読む人の負担を考えれば,どうかな~と思う。
まあ,最適化されたプログラムの方が分かりにくいではないかといわれれば,しかたない。対話できませんね。
「東北大学プレスリリースについての疑問と再分析(pdf)」にて,トレンドを含む変数同士の回帰分析では偽の相関関係が観察されているのではないかということについて,書かれている。
難しいことをやっているように見えるが,偏相関係数を考えて見れば物事はもう少し簡単になるかもしれない。
使用する変数は,年,若年世代1人あたりの新規国債発行額,若年世代の投票率
年 若年世代投票率 若年一人あたり新規国債発行額
年 0.937 -0.769 0.897
若年世代投票率 -0.610 0.784 -0.622
若年一人あたり新規国債発行額 0.837 0.241 0.904
上三角行列は(普通の)単相関係数
下三角行列は偏回帰係数
対角要素は重相関係数
「若年世代投票率」と「若年一人あたり新規国債発行額」の相関係数は -0.622 でかなり強い負の相関
しかし,「年」の影響を取り除いた「若年世代投票率」と「若年一人あたり新規国債発行額」の偏相関係数は 0.241 と,弱いが正の相関関係があるということになっている。
つまり,この結果からいうと,「若者が投票すると新規国債発行額が増える」,棄権しろということか(笑)
ちなみに,同様にして,高齢世代投票率を見てみると。
年 高齢世代投票率 若年一人あたり新規国債発行額
年 0.905 -0.373 0.897
高齢世代投票率 -0.263 0.385 -0.291
若年一人あたり新規国債発行額 0.889 0.105 0.899
高齢世代投票率が下がるとやはり新規公債発行額は増える(相関係数は -0.291 なので,若年世代投票率ほどの相関ではないが)。
偏相関係数は 0.105 で,やはり正の相関。つまり,高齢世代の投票率が上がっても,新規公債発行額は増える。
いずれにしても,少なくとも観察された限り,なにがどうであれ,新規公債発行額はどんどんどんどん増えている(行く)のだ。
世代別投票率と新規国債発行額
図3では,曲線を当てはめているが,その曲線は漸近指数曲線かな?
いずれの世代の投票率も新規国債発行額とは負の相関。
投票率はまあ,曲線にあてはまっているといっても良いが,国債発行額は,単純な曲線にあてはまっているとはいいがたい。
上の図を描くための R プログラム。
このブログの特性への対策のため,以下のプログラムでは,付値は = に変えた。
par(mar=c(3, 3, 1, 4), mgp=c(1.8, 0.8, 0))
plot(若年世代投票率 ~ 年, data=d, pch=15, col="red", ylim=c(35, 85), ylab="投票率")
points(高齢世代投票率 ~ 年, data=d, pch=17, col="darkgreen")
年 = seq(1967, 2012, length=500)
(ans.red = nls(若年世代投票率 ~ a*b^(年-1967)+c, data=d, start=list(a=57, b=0.98, c=19)))
p = ans.red$m$getPars()
lines(predict(ans.red, newdata=list(年=年)) ~ 年, col="red")
text(1990, 58, sprintf("y = %.3f * %.3f^(年-1967) %s %.3f", p[1], p[2], ifelse(p[3] > 0, "+", "-"), abs(p[3])), pos=2, col="red")
(ans.green = nls(高齢世代投票率 ~ a*b^(年-1967)+c, data=d, start=list(a=-7, b=1.01, c=85)))
p = ans.green$m$getPars()
lines(predict(ans.green, newdata=list(年=年)) ~ 年, col="darkgreen")
text(1994, 77, sprintf("y = %.3f * %.3f^(年-1967) %s %.3f", p[1], p[2], ifelse(p[3] > 0, "+", "-"), abs(p[3])), pos=1, col="darkgreen")
par(new=TRUE)
plot(新規国債発行額 ~ 年, data=d, pch=16, col="blue", axes=FALSE, ylab="")
mtext("新規国債発行額", 4, line=1.8)
axis(4, at=0:5*10, lab=0:5*10)
(ans.blue = nls(新規国債発行額 ~ a*b^(年-1967)+c, data=d, start=list(a=100, b=1.1, c=0)))
p = ans.blue$m$getPars()
lines(predict(ans.blue, newdata=list(年=年)) ~ 年, col="blue")
text(1975, 2, sprintf("y = %.3f * %.3f^(年-1967) %s %.3f", p[1], p[2], ifelse(p[3] > 0, "+", "-"), abs(p[3])), pos=4, col="blue")