ここ10年間ぐらいで多重代入法を用いた欠損データ処理が一般化しましたが、多重代入法を使わなくても多重代入法と同様の結果が得ることもできます。Fully Bayesian Imputation Methodと呼ばれている手法ですが*1、最尤法による単一代入（imputation model）…

マルコフ連鎖モンテカルロ法（MCMC）には様々なアルゴリズムがあり、メトロポリス・ヘイスティングス法は尤度関数があれば推定ができる使い勝手のよい方法だと知られています。しかし、高次元で効率が悪いことが知られており、また多峰の確率分布には上手く…

仕様上はできそうなので試してみたらできたと言うだけの話なのですが、Rで呼び出したFortranのプロシージャからRの関数を呼ぶ方法を記しておきます。技術的にどうこうと言うのは無いのですが、流れはややこしくて以下のようになります。 .CallでRのコールバ…

Rもだいたい線形代数の範囲で片付けている限りは、そこそこ速いことがわかります。Fortranとの速度差が出ないというか、R vs R/FortranだとRの方が微妙に速く*1。遅くなるのになぜFortranでRの対数尤度関数とその一階微分を用意したかと言うと、Fortranでハ…

（多くのプログラミング言語の標準ライブラリのものよりは）新しく高速な乱数生成アルゴリズムdSFMTのFortranでの使い方を確認してみました。Fortran標準の乱数生成プロシージャrandom_numberは遅くて偏りのある線形合同法なので、シミュレーションで使うと…

野球の出塁率と長打率を足し合わせた値OPSの評価をするために、打点と得点の合計値を他の指標で回帰するコードです。非負制約をつけるために、ベイズ推定をしています。なお、結果は以下で、OLSで推定してもほとんど同様の値になりました。線形回帰モデルな…

「Modern fortranでルンゲクッタ法」のFortranのコードをRから呼び出すようにしてみました。 RからFortranを呼び出す方法は前にまとめたものがあるのと、ソースコードはMercurialのリポジトリで公開しているので詳細は省きますが、 Rから呼び出すFortranのサ…

「ルンゲクッタ法で様々な連立微分方程式を解く数値計算例(C言語)」のラングフォードの方程式のプロット結果をgifアニメにしてみたので、コードを残しておきます。 Langford Reactor Cで常微分方程式を計算 Rでプロット Cで常微分方程式を計算 gnuplotあたり…

数値解析の速度は申し分がないFortranですが、プロットなどのI/O周りの使い勝手は悪いです。Fortranでサブルーチンを書いて、Rから呼び出すのが良いかなとは思うのですが、Fortranの計算結果をファイルに書き出して、Rに読ますのもひとつの方法です。テキス…

計算機でmapと言うと、リストの要素を無名関数で一括処理することや（Python）、ハッシュテーブルのような構造（Java）を思い浮かべる人が多いかも知れませんが、mmapはファイルに仮想メモリーを割り当てるOSの提供するAPIです。巨大データを処理したり、高…

MCMCpackで順序ロジットのベイズ推定をするかなと思ったのですが、用意されていたのは順序プロビットのMCMCoprobitだったので、順序プロビットをします。順序ロジットの尤度関数はそんなに複雑ではないのですが、微積分をしない限りは大差ないので。 限界効…

前のエントリーを書いているときに、Rのordinalパッケージに限界効果を計算する関数が用意されていないことに気づきました。MASSパッケージとererパッケージを使えばよいわけですが*1、ererパッケージの更新が微妙な感じがしたので、ordinalパッケージの結果…

最近は理工系はもちろん、社会科学系の大学教員も研究で統計解析をしていることが多いので、Rで成績評価する方法を解説しても需要はあまり無さそうな気もするのですが、人文系でニーズがあるようなので。 名簿と試験の成績とレポートの評価の3つのCSVファイ…

ここ5年ぐらい代表的ブラウザーのGoogle Chromeが非SSL/TLSであるhttp://のページに警告を出すようになっていて、Bloggerで独自ドメイン www.anlyznews.com で書いているブログのほうもLet's Encryptを使ってSSL/TLS化したのでメモ書きを。ブログごとの設定…

ノイズが多い時系列データに対しては、観測値から観測誤差を取り除いて、分析者の興味関心である観察不能な状態を取り出すカルマンフィルターがよく使われています。しかし、カルマンフィルタは線形近似になるため、より一般的に用いることができるアルゴリ…

だいぶ前にパッケージを使わず固定効果モデルで推定する例は書いていたのですが、変量効果モデルで推定する例も書いてみました。計量経済学で言う変量効果モデルは、観測値ごとの誤差に加えて、個体ごとの誤差がある線形回帰モデルです。は被説明変数、は説…

「欠測処理でRのmiceパッケージのmice()を使うとして、処理時間の問題皆さんどうしてますかね。」と言う質問に、「多重代入法はM個の代入を独立して実行するのが望ましいので、マルチコアを使える演算環境ならば、並列処理することができます…Flexible Imput…

Kloke and McKean (2012) "Rfit: Rank-based Estimation for Linear Models," The R Journal, Vol.4(2)にアルゴリズムの概要と利用方法の詳細が書いてあるので説明する必要はないのですが、U検定やBrunner-Munzel検定を使ってきた外れ値がある観測値がコーシ…

教科書の分散比のF検定をつかった等分散性の検証は、F検定が正規分布に依存しているので利用できる場合が限られるわけですが、それは分布の裾野の厚さ、尖度が問題であって、尖度が分かればそれで自由度を補正してより広い分布に使えることが知られています…

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。以前のエントリー「反復測定分散分析を主観ベイズ推定に置き換えよう」で、こっちの方が楽と書いたのですが、欠損値補定がや…