餡子付゛録゛

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

R言語

観察された相対リスクが全て交絡バイアスによるものだとしても、未測定交絡因子の相対リスクはそんなに高くなくても良いかも知れない

「受動喫煙防止法について論点整理①:受動喫煙による健康リスク・死亡者数の推定はどのくらい信用できるか?」と言うブログのエントリーで、受動喫煙の相対リスク1.22倍が未測定交絡因子によって引き起こされたバイアスだとして、未測定交絡因子の相対リスク…

Rの関数をラッピング

Rは対話型インターフェイスなので、長々と関数に引数を指定するのは手間隙なので、グローバル・オプションを設定したくなります。しかし、発展途上のライブラリが多いため、グローバル・オプションをつけると上手く動作しなくなったり、そもそもグローバル・…

マクロ経済学の動的計画法説明用コードのMatlabからRへの移植

一橋大学の阿部教授が院生向けのレクチャー・ノートをアップロードしていたのを見つけ、拝読して勉強させて頂いていたのですが、説明用のコード(P.16--20)がMatlabだったのでRに移植してみました。プロプライエタリなMatlab嫌いな人に役立つかも知れないの…

Rで機械学習(SVM)

人間の直観的な認識力を手軽に模倣するフレームワークである機械学習がもてはやされるようになってもう10年以上経つと思いますが、世間では第3次人工知能ブームでまだまだ注目されているようです。扱いやすいパッケージが拡充されているのはもちろん『情報処…

結婚市場からの離脱を考慮した夫婦の社会的階層格差

@kabutoyama_taro氏に、「ネット界隈が想起する"上昇婚"のイメージ」は・・・数学的にありえないと言われた*1のですが、ネット界隈の仮定に基づけばあり得るので、シミュレーションをして例示したいと思います。数式を書いて検証すると煩雑になりますし、「…

Rの並列処理(doParallel/foreach)で繰り返し処理を中断する方法

確認してみたのですが、試行錯誤してみた結果、結論は「ありません」でした*1。仕方が無いので、OpenMP 4.0以前の代替措置*2と同じ事をしてみたいと思います。速度的には改善になるはずです。 library("doParallel")cl <- makeCluster(detectCores()) regist…

Rでテキストファイルを読み込んで分解

RでCSVファイルを上手く読み込めないときは、read.tableにあわせてテキスト・ファイルを書き換えた方が手軽な気がしますが、それが都合で出来ないときは、テキスト・ファイルとしてCSVファイルを読み込んでから、カンマで分離することもできます。 con <- fi…

お気軽にAR(1)の構造転換点を区間推定してみる

学術的な意義は無いのですが、時系列データの構造転換点を求めたいときがあります。記述的に転換点を仮定して推定することが多いわけですが、主観的になりやすく言い合いになることがあるからです。統計学者に殺されそうな荒業ですが、尤度を使って構造転換…

特定月の分散が大きいときの大数の法則

別のブログで使った大数の法則を示すシミュレーションのソースコードです。 # 乱数から分析データを作る 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:…

Rで各集団における大きい方から2番目の値を調べる

ツイッターで見かけた御題なのですが、(学校の課題などで)ありそうなので、置いておきます。 # お試しデータを作る 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月の方が行楽客が少ないとしても、行楽客が減少し出したとは言えないでしょう。そもそも月ごとに日数も異な…

「年/月」の形式を月末日でDate型に変更

月次データで「年/月」と言う形式(e.g. 2013/11、2014/2)は良く見かけると思いますが、Rのread.table関数でデータフレームに読み込むと文字列型*1になってしまいます。日付型でないと、subset関数などで絞り込むときに不便ですね*2。しかし、as.Date関数は…

パネルデータで全員が一斉に説明変数の値を変化させたときのタイムダミーの係数

立命館大学の柴田悠氏への私信なのですが、パネルデータ分析の固定効果モデルに年次などタイムダミーが入っているときに、ある瞬間に個々が同じような行動をとって、説明変数の値を一斉に変化させたときに、それがタイムダミーの係数として観測されない事を…

統計モデルに観測値と観測値の割り算値を入れたときのバイアス

「データ解析のための統計モデリング入門」の著者の久保拓弥氏へのお返事なのですが、別所の「統計モデルに観測値と観測値の割り算値を入れても問題ない」に対して以下のようなコメントをもらったので、具体的な例を作ってみました*1。 …まあ,みなさん割算…

Rの拡張でライフゲームを作ってみる

わざわざ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. モデル ポアソン分布はパラメーターで決定されるわけですが、を説明変数…

Rのグラフに回帰分析の結果である数式を書き込む

データをプロットしたグラフに、推定された予測値や、推定されたパラメーターを数式ごと書き込みたいときはありませんか? そう難しい話では無いのですが、幾つかコツがあるので例を作ってみました。 1. データセット グラフを描くにはデータセットが要るの…

Rのグラフ中文字列の上下左右の寄せ

text関数でプロット領域に文字を書き込むことができるわけですが、指定座標と文字位置の関係を示すadjパラメーターの意味を整理してみます。 グラフの中心に表示座標をとり、adjだけ色々と変えてみたグラフが以下です。 c(1, 0)だと、文字列の右下に表示座標…

Rのグラフで余白の調整

Rのプロットはデフォルトでは余白が多いです。しかし、TeXの方でタイトルや注釈を入れる場合は、プロット領域を大きく取りたくなるので、この余白を削りたくなります。par関数のmar、mgpオプションで余白を調整できるので覚えておきましょう。 1. marオプシ…

Rで多体問題を数値解析

古典物理学と言えばニュートンの運動方程式で、何かの軌道をプロットするのが定番です。三体以上の問題は特殊なケースしか閉じた解が出てこないので、数値解析が行われています。銀河にある星の数をシミュレーションすると、スーパーコンピューターな数値演…

RでKdV方程式を数値解析でプロットしてみる

浅水波を表すKdV方程式は、今では可積分系として閉じた解が知られている非線形偏微分方程式ですが、これは意外に最近知られるようになった事で、二十世紀の中ごろにザブスキーとクルスカルと言う物理学者がシミュレーションをしてソリトンと言う孤立波を見つ…

Rでグラフの軸ラベルを回転させる

計量分析やシミュレーションに関係の無いデータで、Rでグラフを描いている人は少数派だと思います。 データ入力を表計算ソフトなどで行い、そこからテキスト・ファイルで保存するのに手間がかかる*1と言うのもあると思いますが、慣れていないと上手くレイア…

動学マクロ経済学と言う名の非線形連立方程式を解いてみる

動学マクロ経済学に関する蓮見氏のレクチャーノートを拝読していて、もう少しRのソースコード部分に補足説明があったほうが分かりやすい気がしたので、ラムゼー・モデル(Ramsey–Cass–Koopmans model)に関して紹介してみたいと思います。このモデルは連続関…

Rのグラフにスプライン曲線を描きこむ

TTFなどのフォントもそうですし、CADなどグラフィックス分野の方が使う事が多いと思いますが、間欠的なデータから連続した線を描画したいときは実用上は多々あります。ラグランジュ補間、ニュートン補間などバリエーションは幾つもありますが、Rではスプライ…

Rで偏微分方程式

速度的に実用性が低い感じもするのですが、Rで偏微分方程式の数値解を計算することもできます。微分方程式一般に言えることですが、せっせと計算して閉じた解を求めるのは苦労するので、数値解は実用的には重要です。 deSolveパッケージを使えば三元連立式ぐ…

Rで常微分方程式

自発的には揺れないボロビルの力学的振動は、単純化すると以下のように表されます。 が揺れ幅(変位)、が固有振動数で、が抵抗です。は時間を表します。これをプロットしてみましょう。 1. 閉じた解を求める 定係数二階同次方程式なので、数学的に解くこと…

Rでパスワードを作る

最近はウェブサービスが乱立している事もあり、パスワードが漏洩する事件も良く聞くようになりました。この時に複数のサイトでパスワードが共通にしている人は、被害が他のサイトにも拡大していく可能性があります。 それを避けるために複数サイトで異なるパ…

Rでlmコマンドの結果からStd. Errorを取得する

たまに回帰分析の係数の区間推定を求めたいときがあります。 lmコマンドの推定結果から標準誤差(Std. Error)を取得するときは、以下のようにすると簡単です。 r <- lm(expression) # expressionは任意の推定式です coef(summary(r))[, 2] 例えば3番目の係…

Rで解くFizzBuzz問題

前のエントリーでFizzBuzz問題がもっと速くなるのでは無いかと指摘したら、“エレガントではないけどと速いコード”を書いて頂けたので、他のコードとあわせて紹介したいと思います。遅い順に合計7種類。 1. @kos59125 氏 Type I @kos59125 氏が説明用に示した…

Rのベクター処理にバグ? - ではありませんでした

FizzBuzz問題を解いていて気付いたのですが、Rにバグがあるかも知れません。と思ったら、無かったです。 ちょっと紛らわしいコードを書きます。ベクターに添字をすることで、ベクターの一部分の文字列を取得しているのですが、1から11まで一度に計算したとき…

RでFizzBuzz問題を高速化

「裏 RjpWiki」で「R らしく FizzBuzz」と言うエントリーがあるのですが、ワンライナーですっきりしていますが、速度的にちょっと彼らしくない書き方になっています。 ifelse()を減らせば、もっと速度が改善しますよ。手元の計測では裏 RjpWiki版が4.96秒、i…

Rで解いていないFizzBuzz問題

「Rで解くFizzBuzz問題」を見ていて理解できる事は、Rで高速なコードを書くのはトリッキーだと言う事ですね。素直にCで実装した方が速いかも知れないと言うことで、試してみました。 速度的には、Rで最も速かった「裏 RjpWiki 変態チック版」の約1.5倍という…

「Rと手作業で覚える最尤法」のコードに関して補足

どうもコメント欄で粗雑なコードで御迷惑をおかけしましたと書いたら、『裏 RjpWiki』と言うブログで「Rと手作業で覚える最尤法」に関して、『ある意味「えれがんと」だなあ』と皮肉を言われています(ダメ出し:簡明直截なプログラムを書こう)。ダメ出しを…

M-Hアルゴリズムによるポアソン分布推定コードのチューニング

「Rでマルコフ連鎖モンテカルロ法を試す」で書いたMCMCのコードが、NR法による最尤推定とのつながりで混乱していたせいか、それをチューニングしようと言うエントリーが出てしまったので、ちょっと補足をします。 以前のエントリーでは対数尤度関数を使って…

-1*f()より-f()

Rで「-1*f()が気持ち悪い」と言う指摘を受けて、-f()だとSyntax errorになりそうで逆に気持ち悪いんじゃないかと思っていたら、年寄りの戯言でした。 どういう事かと言うと、この二つがRではS式が異なるわけです。 f <- function(n){ n + 1 } as.list(quote(…

加算後の除算≿除算後の加算

かなり大量に演算しないと目立ちませんが、500万サンプルで比較すると、0.10秒と0.16秒で差が計測できました。 obs <- 5000000 x <- runif(obs) y <- runif(obs) ### まとめて割り算 ### gc();gc() system.time({ n <- x + y m <- n / 2 print(sum(m)) }) ##…

Rは行列演算≿apply関数≿ループ演算

Rのプロビット回帰を用いたマイクロ・ベンチマークを書いたので、ソースコードを公開しておきます。 1. セットアップ 推定結果が常に同じになるように、セットアップします。 ### セットアップ ### set.seed(20120209) obs crt x1 x2 y frml # 説明変数を行…

ifelse()と論理値の掛け算

Rでは論理値TRUEとFALSEは、それぞれ1と0として演算する事ができます。ゆえにyとzがxの値に関わらず計算可能な値になるのであれば、ifelse()の代わりに論理値で演算を行う事ができます。 ifelse()はベクター演算にならないので、論理値を使う方が速いです。…

配列の拡張とパフォーマンス

動的に配列サイズが変更されるほど、パフォーマンスは劣化します。 num = 50000### 配列を拡張していく ### gc(); gc() system.time({ x <- numeric(0) i <- 1 while(i <= num){ x[i] <- i i <- i + 1 } })### 配列を一度に初期化 ### gc(); gc() system.tim…

Rでトービット・モデルによる打ち切りデータの推定

金額で測った不動産の保有量を考えましょう。大半の貧乏人はゼロでありますが、一定以上の所得がある人から保有量が増加してきます。集合住宅、持ち家、そして別荘や投資物件など増加していくわけです。 こういう説明変数の値に対して、一定まではゼロ(もし…

地獄のミサワ的Rユーザーのソースコード

別のブログでアイナちゃんの課題をミサワ君とサエコ先輩が助ける話を書いたのですが、検証用にその物語中のソースコードを公開しておきます。 1. データセットの作成 データセットは課題で出されたモノになりますが、ファイルを用意するのも面倒なので、今回…

Rでマルコフ連鎖モンテカルロ法を試す

地味にここ5年間ぐらい、マルコフ連鎖モンテカルロ法(MCMC)が流行っているようです。汎用的な分布でベイズ推定を行う時に有用な数値解析アルゴリズムの総称で、Metropolis-Hastings algorithm(M-Hアルゴリズム)などが主要なメソッドとして使われています…

Rによる数値微分

Rはderiv()関数で記号/微分ができるものの、和や積、つまりsum()やprod()があると微分できません。以下のようなコードを実行すると、『関数 'sum' は導関数の表中にありません』のようにエラーメッセージが出ます。 # ベクトルxを代入 x # 変数bを2とする b …

Rと手作業で覚える最尤法

OLSより進んだ統計手法で最初に覚えるのは最尤法だと思います。大半の人はツールとして知っていて、あまり中身を意識していない気がするのですが、「尤度」の説明無しで『尤度が最大になるパラメーターを求める方法』と言う説明が横行しているのは、問題があ…

Rでニュートン・ラフソン法

文系の学生などで特に目的無くプログラミングを覚えたい人には、Rが最適です。統計解析に使わなくても、ちょっと計算した気になるプログラミングの練習をしてみましょう。 1. よくある近似解を求める 関数f(x)=0があるとして、テーラー展開をしてみます。f(x…

Rで年次ダミーを作成

実際の計量分析ではダミー変数は無くてはならないテクニックで、技術進歩の計測でもしない限りは、複数年にまたがるデータには年次ダミーが含まれている事が通例です。Rでは簡単にダミー変数を作る事ができます。 例えば以下のようなデータフレームdfがある…

Rでデータフレームに、t-1期のデータを追加する(ミクロデータ向き)

被説明変数が説明変数に影響を与えている事を同時性と言います。アイドルのCMやドラマ露出時間と、アイドルの人気度の関係などが例にあげられるでしょう。単回帰では問題ないのですが、重回帰分析では分析に深刻な影響を及ぼす可能性があります。 これを解決…

Rでバランスド・データの作り方

複数年に渡るミクロデータでは、個体が入ったり抜けたりします。例えば地域住人のデータでは、転居や死亡などで個体が減りますし、逆に転入や誕生で個体が増えたりします。サンプル数が膨大な場合は誤差になりますが、限定されたサンプルでは推定結果に影響…

Rでミクロデータにマクロデータを紐付ける

ミクロデータへのマクロデータへの紐付けはデータ分析では良くある作業です。例えば企業の財務データを、その年のインフレ率で割り引く作業があります。Rの演算は各年各企業のデータにインフレ率が無いと演算しにくいので、Excelのvlookup関数で紐付けを済ま…

今日は第何週?

R言語は、他の人気スクリプト言語に比較すると、ベクトル演算以外は記述が煩雑な気がします。特に日付関係は特殊なオブジェクトが乱立しており、少なくとも他の言語とは扱い方がちょっと違います。 例えば、今日が第何週かを調べるには以下のようにすれば良…