餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

OpenBLASをリンクしたWindows版R 4.2 PR

標準のNetlib BLASよりも高速な線形代数ライブラリOpenBLASをリンクしたR 4.1.xを使ってきたのですが、今年から4.2系が標準になるようです。古いRを使い続けていると、更新されたパッケージを使うときに警告で出て目障りなので、バージョンアップをしてみま…

反復測定分散分析を主観ベイズ推定に置き換えよう

TJO氏にふられた質問で、反復測定分散分析(repeated measures anova)をベイズ推定に置き換えられないかと言う話があったので、質問文から想定される状況を整理し、欠損値補定をFully Bayesian Imputation Method *1で、モデル比較をBayes Factorで行なうコ…

ベイズ推定の事前分布で多重共線性に対処しよう

線形回帰などでベイズ統計学を用いる利点は何なのかと言う御題が流れていました。ざっと思いつくところで、 主観的事前分布で、データ以外の情報を加味できる (同じことだが)逐次ベイズ推定にできる 信頼区間よりも信用区間の方が解釈がしやすい 有意水準…

完全情報最尤推定法(FIML)で操作変数のある連立方程式モデルを推定

計量経済学では同時方程式を考えることが多く、そのために一般最小二乗法(OLS)の発展である操作変数法(IVM)とそのバリエーションである二段階最小二乗法(2SLS)*1や三段階最小二乗法(3SLS)、モーメント法の発展である一般化モーメント法(GMM)がよく…

MCMCpackのMCMCmetrop1RでBayes Factorを雑に計算

マルコフ連鎖モンテカルロ法(MCMC)を使いたいベイジアンの統計ユーザーはStanを使っている人が多く、R使いもフロントエンドRStanを通して利用していることがほとんどだと思うのですが、MCMCpackと言うのもあります。このMCMCpackでベイズ因子(Bayes Facto…

Rの主成分回帰で多重共線性を制御

一般線形回帰(OLS)でコスト分析をしたら係数がマイナスの項目が出てきた、工数をかけるほど費用が減っていく・・・と言う嘆きを見かけたのですが、説明変数の相関が高いと推定される係数の標準誤差が大きくなって推定量がおかしくなる、多重共線性による問題…

VSCodeで最小限のLaTeX環境を整備する

Visual Studio Code(以下、VSCode)1.5でLaTeX環境を整備したのですが、参考にしたページの設定が心理的な参入障壁になりそうなほど頑張っていたので、備忘録をかねてWindows 10でplatex/pbibtexを使うだらっとした最小限の設定を紹介します。 1. LaTeXのイ…

OpenBLASをリンクしたWindows版R 4.1 PR

ループでスカラ演算を繰り返すよりもベクトル演算をした方が速いRですが、標準構成ではNetlib BLASと言う線形代数ライブラリBLASのリファレンス実装を使っていて、行列演算、ベクトル演算も速いとは言えないものとなっています。チューニングされたOpenBLAS…

Rの拡張でOpenMPを使ってみる

プロセッサのマルチコア化が進んだ現代なので、時間のかかる計算では並列処理をするコードを書くのが望ましいです。スレッドやセマフォの制御を直接プログラミングすると骨が折れるのですが、シングルコア向け逐次処理コードを並列処理に変換してくれるOpenM…

疎な巨大対称行列の一般化固有値問題のC++/Eigenでの効率的な解き方

「水素原子に対するSchrödinger方程式を有限要素法で数値的に解いてみる(C++のソースコード付き)」で紹介されているC++/Eigenのコードを、(ページ著者が)Juliaに移植してみたら12倍くらい速くなったと言うツイートが流れてきて、気になってC++の方のソー…

FASTGLを使ってC++で積分

数値解析における数値積分法の代表としてガウス求積と言う技法があり、それに必要なガウス点と重みの計算でJuliaがC++より桁違いに速いと言う話が流れていたのですが、Juliaで使われているFastGaussQuadrature.jlが使っているアルゴリズムのC++版を使ったと…

モンテカルロ法で円周率を計算するベンチマーク

モンテカルロ法で円周率を計算するベンチマークで、Juliaと比較してgccが悪い言われているのですが、本当なのか検証してみましょう。CとJuliaでこの分岐が単純なマイクロベンチマークを実行すると、乱数生成の速度が結果を大きく左右しえます。後方互換性の…

SIMD拡張命令セットでgcc用のsum関数を書いてみる

とあるJuliaを使った教材に、JuliaやNumPyのsumはSIMD拡張命令セットを使うから速いけれども、C言語でsumを書くとそうではないから遅いような記述があって、CからSIMDを叩けば良いだろうと思ったので、私にしては珍しくC99のコードで、単浮動小数点のsumを書…

RでべたべたとBHHH法

Berndt-Hall-Hall-Hausmanアルゴリズムは、化石的な最適化アルゴリズムで、実証分析を行っている経済学徒でお世話になったことがある人は、そこそこいると思います。TSPと呼ばれる統計解析パッケージ*1では標準的に使えたのですが、RではmaxLikパッケージ*2…

Rで黄金分割探索

ある化石的な最適化アルゴリズムを書くための準備*1なのですが、黄金分割探索をRで実装してみましょう。計算するだけならばoptimize関数を関数を用いればよいのですが、建築や美術など見た目以外でも黄金比が実用的に使えるのを実感できます。黄金分割探索は…

信頼区間がマイナスに入らない二項分布の区間推定

感染症の罹患経験など、無作為抽出のyes/noの観測値は二項分布に従うことが知られていますが、その場合も正規分布を用いて区間推定を行うことが多いです*1。中心極限定理で正規分布に漸近することを利用しているのですが、稀にしか観測されない事象で、全体…

Rと手作業で覚えるGLM

glmコマンドに頼り切って二項分布(と言うかロジスティック回帰)やポアソン回帰をしている日々ですが、速度狂は最初から何でもC++で・・・と言う話がある*1ので一般化線形回帰(GLM: Generalized Linear Model)の手順を少し詳しく見てみましょう。最尤法を使…

使っている名詞で2つの文の相関係数を取ってみる

Twitterで、ぬ氏が「二つの日本語の単語列(10-20文字ぐらい)があって、なるべく共通の単語が多いとか意味合い的に近いものを自動でマッチしたいときって何かいい方法ありますか?」と言うお題を出していたので、お気楽な方法を考えてみました。わかち書き…

トーラス世界の感染進行モデル

GitHubを使えって絶対に言うな!!!えーっと、BCGワクチン接種が新型コロナウイルス感染症(COVID-19)に有効説検証論文に付随した話で、同一世界にワクチン接種ありとなしの集団がいたらどうなんの? — と言う御題があって、(ドラクエのマップのような縦…

Rでガブリエル比較区間を計算する方法

データを男女などで二つ以上の群にわけて、比較したいときは色々とあります。平均や分散の違いをt検定やF検定などで比較するのが教科書的かつ堅実な方法ですが、統計学に必ずしも詳しくない人々に説明する場合は、ビジュアライゼーションも考えていかないと…

クロス集計表のセルごとの有意確率を出してみよう

他人への質問を眺めていても学ぶことはあります。昨日は「クロス集計の独立性検定を示した表で、セルごとにアスタリスクがあるのは何故か?」と言う質問がされていましたが、私も理由が思いつきませんでした。以下のような頻度表、時折、見てきた気がするの…

行列を使った計算で、SURの感じを掴んでみよう

計量の教科書には載っているものの非線形モデルには使えず、内生性を制御できないので顧みられることの少ないSUR *1 をふと思い出したので、手順を確認してみました。 1. 2つの方程式からなるSURモデル 理屈はGreeneのEconometric Analysisに詳しく載ってい…

相関行列から一定以上の大きさの変数の組み合わせを抽出

100や1000もある説明変数同士の相関係数を見るはだるいと言う話があって、確かに100×100や1000×1000どころか、10×10の行列でも目視で確認すると見落としが出そうです。相関係数0.5以上をリストするコードを書いてみましょう。 # テスト用のデータフレームを…

よくある操作変数法(IV)のための練習データとOLSによる推定バイアス

内生性について上手くツイートでは問題点を説明できないので、念頭に置いているコードを出します。 # # 教科書例の需要と供給のデータを作る # 真のモデル: # S = 2 + 3*p + 4*z + ν # D = 1 - 1*p + μ # S = D # p:価格, S:供給, D:需要, z:気候か何か, νと…

FE-IV練習用データ生成から、one-wayクラスター頑強標準誤差の計算まで

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…

拡張Dickey-Fuller検定を横断面方向に拡張する

同一のデータ生成プロセスを踏んでいる複数の系列 y1, y2, y3, … を、それぞれ個別に単位根検定にかけると、定常であるのに標本数が不足して単位根の棄却に失敗しがちになる一方、c(y1, y2, y3, …) と単純に接続して単位根検定した場合、接続面に断続が生じ…

打ち切りリセット付単位根過程のAR(1)での推定

ちょっと特異なデータを、そこそこ正しく推定するための検討用です。結論は、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…

オッズを使ってロジスティック回帰してみる

ロジット・モデルと言うと個票に応用する二項選択モデルを最尤法などで推定することを思い浮かべる人が多いと思いますが、都道府県ごとの率など集計データに一般線形回帰をかける事で推定する事もできます。観測数不明の集計データしか取れない場合も多いの…