餡子付゛録゛

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

Bloggerで独自ドメインを利用している場合のSSL/TLS化

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

粒子フィルタ(particle filter)で状態とパラメーターを同時推定

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

Rで計量経済学の変量効果モデルを推定

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

miceの並列処理で時間短縮

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

コーシー分布も扱えるランク基準線形モデル推定(Rank-based Estimation for Linear Models)

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

尖度で自由度を補正したF検定による等分散性の検証

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

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を書…