餡子付゛録゛

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

RでBEEP音を鳴らす方法

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

ggplot2やplotlyを使わないbubble plot

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

DiagrammeR/Mermaid.jsで描いた図をreveal.jsに何とか埋め込む方法

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

Rでバイナリファイルを置換

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

Rによる解析結果をShinyでインタラクティブにプレゼンテーションしよう

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

Rでnames(x) <- "a name"的な代入を定義するには

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

Rの関数の引数が参照渡しから値渡しに変わるとき

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

Rでエッジワース・ボックスを描こう

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

Rのリスト処理のざっとした説明

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

前世紀風のRの日付処理

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

Rで線形混合モデル×多重代入法

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

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)に有効説検証論文に付随した話で、同一世界にワクチン接種ありとなしの集団がいたらどうなんの? — と言う御題があって、(ドラクエのマップのような縦…