BEEP音を鳴らしたかったのですが、検索するとbeeprパッケージを使え、alarm関数を使えと言うコレじゃないソリューションや、system関数でOSのコマンドを叩けと言うベタな手が引っかかります。どうやらRでBEEP音を鳴らしたい人々はほとんどいないようです。C…

Rではggplot2やplotlyを使うと簡単にbubble plotを描画できますが、標準のplotを使ってもそんなに手間暇でもないです。より楽ができて、見栄えのするggplot2かplotlyの利用をお勧めしますが。 1. データセット 以下のデータセットをプロットすることを考えま…

Rで簡単にフローチャートやシーケンス図やER図を描いてくれるDiagrammeRパッケージで生成した図を、標準的な方法でR Markdown Format for reveal.js Presentationsで使おうとしたら、微妙に挙動不審なhtmlが出来上がったので、問題の回避策を記しておきたい…

なるべく避けた方が良い気がするのですが、Rでもバイナリのデータを扱えます。0から255の範囲を整数のpack形式と、0か1のunpack形式のraw型として処理されるのですが、C言語あたりから入った人だとややこしく感じるかも知れません*1。 1. pack/unpack形式のr…

数理モデルの数値解析や統計解析では、パラメーターを色々といじったり、変数や分析期間を変えたりして計算をやり直したり、結果の一部分を切り取って参照したいときがあります。プログラムをぽちぽちと変えていけば実現できるわけですが、ぱっと見ではその…

今まで一度たりとも作る必要を感じたことがなく、言語定義を見ても特段何も書いていないようなので、関数呼び出しへの代入（？）はユーザー定義はでき無いのかなと思っていた*1のですが、ヘルプを見ていたら*2 The functions 'dim' and 'dim などと書いてあ…

実用上はRの関数の引数は値渡しと理解しておけばよいのですが、引数を他の変数に代入したり、関数の中の関数呼び出しで使ったりする前は、参照渡しになっています。実際、 fn1 <- function(x) x^2 fn2 <- function(x, fn){ rm("fn1", envir = parent.env(env…

Edgeworth Box以前に描いたコードを整理しなおしたと言うか、手計算で微分して均衡や契約曲線を求めていたのでどちらかと言うと作り直したものですが、applyやwhichやexpressionなどに慣れるのに良さそうだったので。数値解析で均衡点を求めているので、初期…

説明用です。他のプログラミング言語の連想配列/マップ/ハッシュと同様に使えるわけですが、[ ]と[[ ]]の使い分けがややこしいですね。 # 空リストを作成 lst <- list() # リストを作成 lst <- list("a"=123, "b"=456, "c"=789) # ベクターをリストに変えて…

以前書いたエントリーを見直していて、操作を手早く紹介する例をつくっておきたくなりました。文字列型や数値型の日付データを、タイムゾーンUTCでのUNIX秒を保持するクラスPOSIXctか、日付構造体になっているリストPOSIXlt型に変換して操作して、文字列型に…

反復測定分散分析をlme4パッケージとmiceパッケージを使った線形混合モデル*1×多重代入法で行なう作業を確認したので公開します*2。以前のエントリー「反復測定分散分析を主観ベイズ推定に置き換えよう」で、こっちの方が楽と書いたのですが、欠損値補定がや…

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

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

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

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

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

ループでスカラ演算を繰り返すよりもベクトル演算をした方が速い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 …