R言語
一般線形回帰（OLS）でコスト分析をしたら係数がマイナスの項目が出てきた、工数をかけるほど費用が減っていく･･･と言う嘆きを見かけたのですが、説明変数の相関が高いと推定される係数の標準誤差が大きくなって推定量がおかしくなる、多重共線性による問題…
ループでスカラ演算を繰り返すよりもベクトル演算をした方が速いRですが、標準構成ではNetlib BLASと言う線形代数ライブラリBLASのリファレンス実装を使っていて、行列演算、ベクトル演算も速いとは言えないものとなっています。チューニングされたOpenBLAS…
プロセッサのマルチコア化が進んだ現代なので、時間のかかる計算では並列処理をするコードを書くのが望ましいです。スレッドやセマフォの制御を直接プログラミングすると骨が折れるのですが、シングルコア向け逐次処理コードを並列処理に変換してくれるOpenM…
Berndt-Hall-Hall-Hausmanアルゴリズムは、化石的な最適化アルゴリズムで、実証分析を行っている経済学徒でお世話になったことがある人は、そこそこいると思います。TSPと呼ばれる統計解析パッケージ*1では標準的に使えたのですが、RではmaxLikパッケージ*2…
ある化石的な最適化アルゴリズムを書くための準備*1なのですが、黄金分割探索をRで実装してみましょう。計算するだけならばoptimize関数を関数を用いればよいのですが、建築や美術など見た目以外でも黄金比が実用的に使えるのを実感できます。黄金分割探索は…
感染症の罹患経験など、無作為抽出のyes/noの観測値は二項分布に従うことが知られていますが、その場合も正規分布を用いて区間推定を行うことが多いです*1。中心極限定理で正規分布に漸近することを利用しているのですが、稀にしか観測されない事象で、全体…
glmコマンドに頼り切って二項分布（と言うかロジスティック回帰）やポアソン回帰をしている日々ですが、速度狂は最初から何でもC++で･･･と言う話がある*1ので一般化線形回帰（GLM: Generalized Linear Model）の手順を少し詳しく見てみましょう。最尤法を使…
Twitterで、ぬ氏が「二つの日本語の単語列（10-20文字ぐらい）があって、なるべく共通の単語が多いとか意味合い的に近いものを自動でマッチしたいときって何かいい方法ありますか？」と言うお題を出していたので、お気楽な方法を考えてみました。わかち書き…
GitHubを使えって絶対に言うな！！！えーっと、BCGワクチン接種が新型コロナウイルス感染症（COVID-19）に有効説検証論文に付随した話で、同一世界にワクチン接種ありとなしの集団がいたらどうなんの？ — と言う御題があって、（ドラクエのマップのような縦…
データを男女などで二つ以上の群にわけて、比較したいときは色々とあります。平均や分散の違いをt検定やF検定などで比較するのが教科書的かつ堅実な方法ですが、統計学に必ずしも詳しくない人々に説明する場合は、ビジュアライゼーションも考えていかないと…
他人への質問を眺めていても学ぶことはあります。昨日は「クロス集計の独立性検定を示した表で、セルごとにアスタリスクがあるのは何故か？」と言う質問がされていましたが、私も理由が思いつきませんでした。以下のような頻度表、時折、見てきた気がするの…
計量の教科書には載っているものの非線形モデルには使えず、内生性を制御できないので顧みられることの少ない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)# …
