R言語
計量の教科書には載っているものの非線形モデルには使えず、内生性を制御できないので顧みられることの少ないSUR *1 をふと思い出したので、手順を確認してみました。 1. 2つの方程式からなるSURモデル 理屈はGreeneのEconometric Analysisに詳しく載ってい…
100や1000もある説明変数同士の相関係数を見るはだるいと言う話があって、確かに100×100や1000×1000どころか、10×10の行列でも目視で確認すると見落としが出そうです。相関係数0.5以上をリストするコードを書いてみましょう。 # テスト用のデータフレームを…
内生性について上手くツイートでは問題点を説明できないので、念頭に置いているコードを出します。 # # 教科書例の需要と供給のデータを作る # 真のモデル: # S = 2 + 3*p + 4*z + ν # D = 1 - 1*p + μ # S = D # p:価格, S:供給, D:需要, z:気候か何か, νと…
dとsの式でpが内生変数、zが操作変数、iが個体を表す番号ですが、こんなんで。 noi <- 20 # 個体数 t <- 5 # 観測期間 obs <- t*noi i <- rep(1:(noi), each=t) fe <- runif(noi, min=0, max=100) a0 <- rep(fe, each=t) a1 <- -1 b0 <- 2 b1 <- 3 b2 <- 4 z …
ゲーム理論で考える、じゃんけんの拡張: ニュースの社会科学的な裏側で使ったコードです。通常のゲーム理論の数値演算で扱われる計算手順に習っていないと言うか、フィーリングで描いたので、何か勘違いがあるかも知れないです。 ##########################…
観測数100ぐらいでも、誤差項の分散など次第では多重共線性が推定結果に影響を与える例を作ってみました。 多重共線性が無いケース set.seed(2103) n <- 100 x1 <- 0.3 * runif(n, min=0, max=10) + rnorm(n, sd=1) x2 <- 0.7 * runif(n, min=0, max=10) + r…
同一のデータ生成プロセスを踏んでいる複数の系列 y1, y2, y3, … を、それぞれ個別に単位根検定にかけると、定常であるのに標本数が不足して単位根の棄却に失敗しがちになる一方、c(y1, y2, y3, …) と単純に接続して単位根検定した場合、接続面に断続が生じ…
ちょっと特異なデータを、そこそこ正しく推定するための検討用です。結論は、AR(1)にしておけば無問題でした。 ### データ生成 ### genSample <- function(n=250){ ex <- rnorm(n) x <- rep(NA, length = n) x[1] <- ex[1] * 0.2 for (j in 2:n) { # 係数0.9…
ロジット・モデルと言うと個票に応用する二項選択モデルを最尤法などで推定することを思い浮かべる人が多いと思いますが、都道府県ごとの率など集計データに一般線形回帰をかける事で推定する事もできます。観測数不明の集計データしか取れない場合も多いの…
「受動喫煙防止法について論点整理①：受動喫煙による健康リスク・死亡者数の推定はどのくらい信用できるか？」と言うブログのエントリーで、受動喫煙の相対リスク1.22倍が未測定交絡因子によって引き起こされたバイアスだとして、未測定交絡因子の相対リスク…
Rは対話型インターフェイスなので、長々と関数に引数を指定するのは手間隙です。グローバル・オプションを設定したくなります。しかし、発展途上のライブラリが多いため、グローバル・オプションをつけると上手く動作しなくなったり、そもそもグローバル・オ…
一橋大学の阿部教授が院生向けのレクチャー・ノートをアップロードしていたのを見つけ、拝読して勉強させて頂いていたのですが、説明用のコード（pp.16—20）がMatlabだったのでRに移植してみました。プロプライエタリなMatlab嫌いな人に役立つかも知れないの…
人間の直観的な認識力を手軽に模倣するフレームワークである機械学習がもてはやされるようになってもう10年以上経つと思いますが、世間では第3次人工知能ブームでまだまだ注目されているようです。扱いやすいパッケージが拡充されているのはもちろん『情報処…
@kabutoyama_taro氏に、「ネット界隈が想起する"上昇婚"のイメージ」は・・・数学的にありえないと言われた*1のですが、ネット界隈の仮定に基づけばあり得るので、シミュレーションをして例示したいと思います。数式を書いて検証すると煩雑になりますし、「…
確認してみたのですが、試行錯誤してみた結果、結論は「ありません」でした*1。仕方が無いので、OpenMP 4.0以前の代替措置*2と同じ事をしてみたいと思います。速度的には改善になるはずです。 library("doParallel")cl <- makeCluster(detectCores()) regist…
RでCSVファイルを上手く読み込めないときは、read.tableにあわせてテキスト・ファイルを書き換えた方が手軽な気がしますが、それが都合で出来ないときは、テキスト・ファイルとしてCSVファイルを読み込んでから、カンマで分離することもできます。 con <- fi…
学術的な意義は無いのですが、時系列データの構造転換点を求めたいときがあります。記述的に転換点を仮定して推定することが多いわけですが、主観的になりやすく言い合いになることがあるからです。統計学者に殺されそうな荒業ですが、尤度を使って構造転換…
別のブログで使った大数の法則を示すシミュレーションのソースコードです。 # 乱数から分析データを作る 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:…
ツイッターで見かけた御題なのですが、（学校の課題などで）ありそうなので、置いておきます。 # お試しデータを作る set.seed(20150208) n <- 30 df <- data.frame(type=c("A", "B", "C")[round(runif(n, min=0.5, max=3.5))], value=60 - ((1:n)-10)^2)# …
時系列データには、季節バイアスが入りがちです。年中行事はもちろんのこと、天候の変化も周期的に発生します。ゴールデンウィークのある5月よりも6月の方が行楽客が少ないとしても、行楽客が減少し出したとは言えないでしょう。そもそも月ごとに日数も異な…
月次データで「年/月」と言う形式（e.g. 2013/11、2014/2）は良く見かけると思いますが、Rのread.table関数でデータフレームに読み込むと文字列型*1になってしまいます。日付型でないと、subset関数などで絞り込むときに不便ですね*2。しかし、as.Date関数は…
立命館大学の柴田悠氏への私信なのですが、パネルデータ分析の固定効果モデルに年次などタイムダミーが入っているときに、ある瞬間に個々が同じような行動をとって、説明変数の値を一斉に変化させたときに、それがタイムダミーの係数として観測されない事を…
「データ解析のための統計モデリング入門」の著者の久保拓弥氏へのお返事なのですが、別所の「統計モデルに観測値と観測値の割り算値を入れても問題ない」に対して以下のようなコメントをもらったので、具体的な例を作ってみました*1。 …まあ，みなさん割算…
わざわざC言語で書くほどのものではなく、単に行列の受渡しを確認しただけですが、Rの中で使えるライフゲームを作ってみたので公開します。 以下のソースコードを LifeGame.c と名づけて保存します。 #include <R.h> #include <Rinternals.h>SEXP LifeGame(SEXP m) { SEXP ans, </rinternals.h></r.h>…
某所の(1)ポアソン回帰モデルの説明が、(2)対数変換OLSと同じになっている気がします。違うものだと思うのですが、シミュレーションをして(1)と(2)の推定をして確認してみました。 1. モデル ポアソン分布はパラメーターで決定されるわけですが、を説明変数…
データをプロットしたグラフに、推定された予測値や、推定されたパラメーターを数式ごと書き込みたいときはありませんか？ そう難しい話では無いのですが、幾つかコツがあるので例を作ってみました。 1. データセット グラフを描くにはデータセットが要るの…
text関数でプロット領域に文字を書き込むことができるわけですが、指定座標と文字位置の関係を示すadjパラメーターの意味を整理してみます。 グラフの中心に表示座標をとり、adjだけ色々と変えてみたグラフが以下です。 c(1, 0)だと、文字列の右下に表示座標…
Rのプロットはデフォルトでは余白が多いです。しかし、TeXの方でタイトルや注釈を入れる場合は、プロット領域を大きく取りたくなるので、この余白を削りたくなります。par関数のmar、mgpオプションで余白を調整できるので覚えておきましょう。 1. marオプシ…
古典物理学と言えばニュートンの運動方程式で、何かの軌道をプロットするのが定番です。三体以上の問題は特殊なケースしか閉じた解が出てこないので、数値解析が行われています。銀河にある星の数をシミュレーションすると、スーパーコンピューターな数値演…
浅水波を表すKdV方程式は、今では可積分系として閉じた解が知られている非線形偏微分方程式ですが、これは意外に最近知られるようになった事で、二十世紀の中ごろにザブスキーとクルスカルと言う物理学者がシミュレーションをしてソリトンと言う孤立波を見つ…
