■
ずっと気になっている（しかしマイナな）問題が解決しないので，そのメモを記録しておきます．
所得分布の生成モデルは，投資行動の反復から対数正規分布を導出する数理モデルである．仮定は次のとおりである（Hamada 2003; 2004; 2016）.
- 人々が最初に初期資産
y_0を持っている． - 資産のうち割合
bを投資して，確率pで成功すると，資産が増える．失敗すると投資した資産を失う． - 資産の一定割合を所得として得る．人々が投資行動を
n回繰り返すと，資産と所得の分布は対数正規分布にしたがう．
このモデルの結論の一つは，導出された対数正規分布のパラメータがモデルの外生的パラメータであるy_0, b, p, nによって内生的に決まる，というものである．
具体的には，資産の分布は
となる．証明はHamada (2003; 2004; 2016）をみよ．
これは2項分布が正規分布
で近似できる，という性質を利用した命題である（確率変数を基準化した場合にはド・モアブル=ラプラスの定理となる）．
さて
このモデルから得た命題をシミュレーションによって再現できない，というのが問題だった．
どういうことか説明しよう．
まず次の関数によって，確率pで資産y_0 bを増やし，1-pで資産y_0 bを失う仮定を表現する．コードは例によって一番楽に書けるMathematica．
gmi[y0_, b_, p_, n_] := Module[{y = y0}, Do[(* repeat n time *) If[Random[] < p, y = y + (y*b), y = y - (y*b)], {n}]; y];
一方で，理論的に導出される定まる確率分布のパラメータを使って，対数正規分布に従う疑似乱数をコンピュータに作らせる．
(* generate samples from a pdf with parameters derived from Hamada model *) lognsamples[y0_, b_, p_, n_, pop_] := Module[{m, s, a, B }, a = Log[(1 + b)/(1 - b)]; B = Log[y0] + n Log[1 - b]; m = B + a n p; s = Sqrt[n p (1 - p) a^2]; RandomVariate[LogNormalDistribution[m, s], {pop}]];
モデルが正しく対数正規分布を導出できているのであれば，ミクロレベルでの資産生成プロセスを集約した結果と，分布からサンプリングした結果は，ほぼ等しくなるはずである．
modelfit[y0_, b_, \[Gamma]_, n_, pop_] := Module[{dat0, dat1}, (* generate samples from a model *) dat0 = Table[gmi[y0(*ini.c*), b(*interest*), p(*success*), n(*repeat*)], {pop(* population *)}]; dat1 = lognsamples[y0(*ini.c*), b(*interest*), p(*success*), n(*repeat*), pop]; Histogram[{dat0, dat1}, Automatic, "Probability", ChartStyle -> {Red, Blue}]];
ところが,この関数を実行すると，「理論値から得たサンプルデータ」と「疑似乱数で作ったサンプルデータ」がズレる．
ズレの解消の仕方として，疑似乱数でサンプルデータを作るのではなく，確率密度関数を理論値で得た離散的所得階級値間のレンジで積分して，確率密度を少し「まとめて」やれば，うまくいくことは分かっている．
だが，いちいち積分させるのもコードが煩わしい．
先日，ふと気づいた．疑似乱数から生成したデータベクトルの中には同じ値がdouble精度で1つもないのに，理論に基づいて生成した方は，double精度で同じ値がいくつも含まれているのだ．ようするに，理論値のベクトルの方が疑似乱数のベクトルよりも分散が小さいのである．
これは一体，どういうことなのか？
理論値は理論値で間違っていない．疑似乱数ベクトルも指定したパラメータを持つ分布からランダムにサンプルを抽出している．
分布を導出しているはずの数理モデルから生成したサンプルが，なぜPCで発生させた確率分布疑似乱数と分布がずれてしまうのか？
そこで理論値は厳密すぎて，疑似乱数のランダムさを失っていると考えてみる．
確かに関数gmi[y0_, b_, p_, n_]の中で一様乱数を使っているので投資結果はランダムに変化するのだが，獲得した資本額がある決まった値しかとらない．だから理論値にノイズを足してやればいいのでは？ というわけだ．
別の言い方をすれば疑似乱数の方は連続確率変数の実現値を表現しているのに，理論モデルは離散的確率変数の実現値しか表現できていない．だから確率密度関数を区間毎に積分する操作が，実現値の離散化に対応していたのである．よって修正はただノイズ
RandomVariate[NormalDistribution[0, 0.1]]
を加えるだけでよいのでは，と考えた.
修正したコードはこうなる
lognsamples[y0_, b_, p_, n_, pop_] := Module[{m, s, a, B}, a = Log[(1 + b)/(1 - b)]; B = Log[y0] + n Log[1 - b]; m = B + a n p; s = Sqrt[n p (1 - p) a^2]; RandomVariate[LogNormalDistribution[m, s], {pop}]]; gmi[y0_, b_, p_, n_] := Module[{y = y0}, Do[(* repeat n time *) If[Random[] < p, y = y + (y*b), y = y - (y*b)], {n}]; y]; modelfit[y0_, b_, p_, n_, pop_] := Module[{dat0, dat1}, dat0 = Table[gmi[y0(*ini.c*), b(*interest*), p(*success*), n(*repeat*)] + RandomVariate[NormalDistribution[0, 1/Sqrt[n]]], {pop(* population *)}]; (* add randome error since the model is too rigorous *) dat1 = lognsamples[y0(*ini.c*), b(*interest*), p(*success*), n(*repeat*), pop]; Histogram[{dat0, dat1}, ChartStyle -> {Red, Blue}]];
これで，だいぶマシにはなる．だがパラメータの値によってノイズの分散に微妙に調整が必要なので 完璧ではない．
多分，記録しておかないと数日したら忘れちゃうので，書き留めておきます．
順番を決めるコード
新学期がはじまるので，演習報告担当者の順番決め用のあみだくじコード（Mathematica）を再掲しておきます．
(* 動的表現 サイズ固定*) Manipulate[ Dynamic[Refresh[RandomSample[Range[n]], UpdateInterval -> x]], {{x, Infinity, "シャッフル"}, {0 -> "Go", Infinity -> "Stop"}}, {{n, 20, "人数"}, 2, 100, 1}, ContentSize -> {600, 150}]
ダイナミックを使ってるので，動的に更新します． この短さでGUIまで自動的に生成するとは驚きである． やっぱりMathematicaで書くとコードがシンプルでええなあ．
探索的分析の自動化実験
内容
データフレームを突っ込むと，データに含まれる全変数の中から1つを順番に応答変数として選び，残りの変数に回するような分析を自動的にやってくれるコードを書いてみた．
「本格的な分析の前に，変数間の関係を探索的・網羅的に知りたい」場合に使うことを想定している． 「こういう頭を使わない機械的な分析はやってはいけない」という例を示すための実験でもある．
目的
データがn個の変数を含む場合に，
が応答変数で，それ以外の全ての変数が説明変数
が応答変数で，それ以外の全ての変数が説明変数
･･･
が応答変数で，それ以外の全ての変数が説明変数
という分析を自動化したい．
コードの作成
基本的なアイデアは次の2つのステップに区分できる．
一つは，テキストで式を自動生成してから，fomulaオブジェクトに変換してlmに突っ込むという方法である．
これは，Rの型キャストを利用した方法であり，オブジェクトの使い方がやや邪道な気もするが，型を意識した
コードの理解には適している．
xnam <- paste("x", 1:25, sep="") fmla <- as.formula(paste("y ~ ", paste(xnam, collapse= "+")))
もう一つは，応答変数を先頭列にもってきて，データフレームごと，lmに突っ込むというワイルド（？）なやり方である．
d2 <- data.frame(y,x1,x2,x3) lm(d2) Call: lm(formula = d2)
この場合，データフレームの先頭列をrotateするコードが別途必要になる．ちなみに上の二つのコードはこちらからの引用である．
上記のコードを元に，自動データ分析関数のプロトタイプを作ってみた．
analyzer <- function(dataframe){ n <- length(dataframe)#data.frame内の変数の個数を取得 no.col <- 1:n #ローテーション用ベクトル n1 <- n-1 #ループエンド用 for(i in 0:n1){ x <- no.col[1] # print(x) no.col[1] <- no.col[i+1] no.col[i+1] <- x #print(no.col) dataframe0 <- dataframe[,no.col] out = lm(formula = dataframe0) print(summary(out)) }}
計算例
なにかテキトーなデータ行列を使って試してみる．
dat0 = data.frame(swiss)#data行列の定義 head(dat0) #中身確認 analyzer(dat0)
アウトプットはこんな感じになる．
Call: lm(formula = dataframe0) Min 1Q Median 3Q Max -15.2743 -5.2617 0.5032 4.1198 15.3213 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 66.91518 10.70604 6.250 1.91e-07 *** Agriculture -0.17211 0.07030 -2.448 0.01873 * Examination -0.25801 0.25388 -1.016 0.31546 Education -0.87094 0.18303 -4.758 2.43e-05 *** Catholic 0.10412 0.03526 2.953 0.00519 ** Infant.Mortality 1.07705 0.38172 2.822 0.00734 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 7.165 on 41 degrees of freedom Multiple R-squared: 0.7067, Adjusted R-squared: 0.671 F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10 Call: lm(formula = dataframe0) Min 1Q Median 3Q Max -29.958 -9.107 1.287 10.805 24.504 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 135.1861 22.7578 5.940 5.26e-07 *** Fertility -0.7410 0.3027 -2.448 0.018727 * Examination -0.9627 0.5117 -1.881 0.067050 . Education -1.4908 0.4119 -3.619 0.000804 *** Catholic 0.1674 0.0762 2.197 0.033745 * Infant.Mortality -0.3609 0.8637 -0.418 0.678268 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 14.87 on 41 degrees of freedom Multiple R-squared: 0.618, Adjusted R-squared: 0.5715 F-statistic: 13.27 on 5 and 41 DF, p-value: 1.045e-07 Call: lm(formula = dataframe0) Residuals: Min 1Q Median 3Q Max -9.0429 -2.8491 -0.5308 3.0459 6.8874 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 24.57200 8.23969 2.982 0.00480 ** Fertility -0.09523 0.09371 -1.016 0.31546 Agriculture -0.08254 0.04387 -1.881 0.06705 . Education 0.33037 0.12858 2.569 0.01392 * Catholic -0.06778 0.02108 -3.215 0.00254 ** Infant.Mortality 0.09730 0.25297 0.385 0.70250 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 4.353 on 41 degrees of freedom Multiple R-squared: 0.7346, Adjusted R-squared: 0.7022 F-statistic: 22.7 on 5 and 41 DF, p-value: 7.624e-11
中略．
このように，関数を実行するとlmの結果がずらずら出てくる．
さて，以上の結果をさらにstep関数に渡してモデル選択まで自動でやろうとすると，エラーが出てしまう．
どうやら
formula = dataframe0
の部分が，step関数のオブジェクトとして不適切らしい．
しかたないので，一旦データフレームの列をrotateしてから，順番の入れ替わった変数名リストを取り出し，
テキスト型でy ~ x1 + x2 + ... + xnの形にしてからformulaにキャストする，という荒技(?)
を使ってみる．ついでにlmの実行の度に応答変数が変わるので，名前が分かるようにprintしておく．
#step関数を使い, AIC基準でモデル選択まで実行 analyzer101 <- function(dataframe){ n <- length(dataframe)#data.frame内の変数の個数を取得 no.col <- 1:n#ローテーション用ベクトル n1 <- n-1#ループエンド用 for ( i in 0:n1 ){ x <- no.col[1] # print(x) no.col[1] <- no.col[i+1] no.col[i+1] <- x #ベクトルのrotate #print(no.col) dataframe0 <- dataframe[,no.col]#データフレーム列のrotate variable = names(dataframe0) xnam = variable[2:n] #説明変数のラベル ynam = variable[1] #oucomeのラベル model = as.formula(paste(paste(ynam,"~") , paste(xnam, collapse= "+"))) out1 = lm(model,data=dataframe0) print(paste("Outcome is", names(dataframe0)[1])) out2 <- step(out1,trace=1)#AIC基準でモデル選択 print(out2) }}
結果を確認してみよう. データフレームはさっきと同じ．
> analyzer101(dat0) [1] "Outcome is Fertility" Start: AIC=190.69 Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality Df Sum of Sq RSS AIC - Examination 1 53.03 2158.1 189.86 <none> 2105.0 190.69 - Agriculture 1 307.72 2412.8 195.10 - Infant.Mortality 1 408.75 2513.8 197.03 - Catholic 1 447.71 2552.8 197.75 - Education 1 1162.56 3267.6 209.36 Step: AIC=189.86 Fertility ~ Agriculture + Education + Catholic + Infant.Mortality Df Sum of Sq RSS AIC <none> 2158.1 189.86 - Agriculture 1 264.18 2422.2 193.29 - Infant.Mortality 1 409.81 2567.9 196.03 - Catholic 1 956.57 3114.6 205.10 - Education 1 2249.97 4408.0 221.43 Call: lm(formula = Fertility ~ Agriculture + Education + Catholic + Infant.Mortality, data = dataframe0) Coefficients: (Intercept) Agriculture Education Catholic Infant.Mortality 62.1013 -0.1546 -0.9803 0.1247 1.0784 [1] "Outcome is Agriculture" Start: AIC=259.3 Agriculture ~ Fertility + Examination + Education + Catholic + Infant.Mortality Df Sum of Sq RSS AIC - Infant.Mortality 1 38.59 9101.4 257.50 <none> 9062.9 259.30 - Examination 1 782.31 9845.2 261.20 - Catholic 1 1066.74 10129.6 262.53 - Fertility 1 1324.81 10387.7 263.72 - Education 1 2895.46 11958.3 270.33 Step: AIC=257.5 Agriculture ~ Fertility + Examination + Education + Catholic Df Sum of Sq RSS AIC <none> 9101.4 257.50 - Examination 1 816.46 9917.9 259.54 - Catholic 1 1076.03 10177.5 260.76 - Fertility 1 1947.58 11049.0 264.62 - Education 1 3091.51 12192.9 269.25 Call: lm(formula = Agriculture ~ Fertility + Examination + Education + Catholic, data = dataframe0) Coefficients: (Intercept) Fertility Examination Education Catholic 132.6210 -0.7989 -0.9802 -1.5193 0.1681
以下，略．
いまのところnumericなデータにしか使えないが，気が向いたらfactorが応答変数になった場合は自動的に 多項ロジットを実行する機能などを追加してみる．
繰り返すが，こんなAICの使い方は意味が無いし，「有意性を見つけてから理屈を考える」というアドホックな分析は，研究としてやってはいけない． 上記の関数からは，理論的に無意味なアウトプットもたくさん出てくるので，使用はリコードの確認や 探索的分析にとどめよう．
mathjax実験
はてぶろがマークダウンに対応しているのを知ったとき，へえそう，くらいにしか思ってなかった.
が...
いやいや，とっくにMathjax使えるようになっとるやん．
とか
とか．
まあ，Rmarkdownやそのものより使いづらいけど，少なくとも見た目はRpubsより綺麗やん．
これからは全ての資料はで書いて，WEBにアップするときだけマークダウンする方式でいこうかな・・・．
というか，そういう変換コード，誰か既に作ってるよな．．．
p値の分布を計算する
人工的に作成したデータを分析すると，統計の理解が捗ります．「p値ちっちゃい，やったー」みたいな意味のない分析するより，こっちのほうが統計初学者には勉強になるのではないか，と思います．
次のような例を考えてみましょう．
仮説検定を繰り返してp値の分布がどうなっているかを調べる
最初に正規分布に従う乱数を作ります
rnorm(10, 0, 1)
たとえばこのコードは平均0標準偏差1の正規分布に従う乱数を10個作ります．
次に，正規分布(平均0，標準偏差sd)に従う乱数列を二つ作り，その平均値をt検定する関数を定義します関数の名前はcoin1です．
#coin1 coin1<-function(n,sd){ x <- rnorm(n, 0, sd) y <- rnorm(n, 0, sd) z<-t.test(x,y,var.equal=T) z$p.value}
関数testは関数coin1をk回繰り返します．
test<-function(k,n,sd){ x <- c() # 空の（要素数ゼロの）リストを作る for (i in 1:k) { y <- coin1(n,sd) #コイン1000回投げt検定の結果p値のみ返す x <- c(x, y) } x }
最後にp値の分布を確認します．
hist(test(k=1000,n=10000,sd=5), breaks = 25)
仮説検定は通常，データに対して一回しか用いません．だから計算されたp値は1個です．
しかし計算された，たった一つのｐ値をみても，正しい仮説が何かは分かりません．帰無仮説が正しいときp値は一様分布に従います．
