Common Lispで統計関数を車輪の再発明(再訪)

だいぶ前に書いたCommon Lispのプログラムを発掘しました。

marui.hatenablog.com

まだapplymapcarを知らなかったので、再帰だけを使って計算していました。たとえばベクトル(というか数値の入ったリスト)の和を計算するのにも、以下のように再帰を使っていました。

(defun sum (x)
  (if (null x)
      0
      (+ (first x) (sum (rest x)))))

(sum '(1 2 3 4))                        ;=> 10

applyを使えば一発で書けます。

(defun sum (x)
  (apply '+ x))

(sum '(1 2 3 4))                        ;=> 10

同様に幾何平均を計算するところも、以下のようにまどろっこしいことをしていました。(とはいえ(expt (* 1 2 3 4) (/ 1 4))みたいに直接乗算を使っておらず、logexpを使ってオーバーフローが生じにくいようにしてはいますね)

(defun geomean-sub (x)
  (if (null x)
      0
      (+ (log (first x)) (geomean-sub (rest x)))))

(defun geomean (x)
  (exp (/ (geomean-sub x) (length x))))

(geomean '(1 2 3 4))                    ;=> 2.213364

mapcarを使うと次のように書き換えられます。ただ、算術平均を計算するmean関数も欲しかったので、そちらを作って下請けとして使っています。

(defun mean (x)
  (/ (sum x) (length x)))

(defun geomean (x)
  (exp (mean (mapcar 'log x))))

(geomean '(1 2 3 4))                    ;=> 2.213364

長いデータ列の処理をすることを考えると、再帰でスタックを消費するのはよくありませんから(データが多いと計算できなくなってしまう)、applymapcarなどを使うほうが安全そうです。こんなかんじで、プログラミング言語の勉強のために、ゆっくり車輪の再発明をしています。

Matlablinspaceに対応するものもCommon Lispで作ってみました。(linspace x1 x2)とすると、x1からx2の範囲を等間隔に100個作ってくれます。個数の指定をしたいときにはオプションで個数を書いてあげるとOK。

;; linspace
;; (linspace -2.0 +1.0)                 ;=> (-2.0 -1.969697 -1.939394 ... 0.969697 1.0)
;; (linspace 1 5 5)                     ;=> (1 2 3 4 5)
;;
(defun linspace (x1 x2 &optional (n 100))
  (let ((x (loop for x from 0 to (1- n)
                 collect (/ x (1- n)))))
    (mapcar #'(lambda (v) (+ v x1))
            (mapcar #'(lambda (v) (* v (- x2 x1))) x))))

追記 データはリストではなくベクトルにしたほうがいいのかも。あるいはRのようなデータフレーム型を作ったほうがいいか。……規模が大きくなってきたら考えよう。(というより、遊びで作っているだけなので規模が大きくならないようにしよう。統計処理をするツールはLisp-Statとかがあるし。)

本日の給油(CRF1100Lアフリカツイン/2BL-SD10)

給油日 オドメーター (km) 給油量 (L) 単価 (円/L) 燃費 (km/L) 距離単価 (円/km)
2021-10-24 10 9.69 155.01 - -
2021-10-24 48 7.62 154.99 - -
2021-10-30 267 13.46 162.04 16.27 9.96
2021-11-06 507 15.29 157.03 15.70 10.00
2021-11-14 737 13.01 157.03 17.68 8.88
2021-11-27 924 13.07 157.00 14.31 10.97
2021-12-25 1188 15.90 155.97 16.60 9.39
2022-01-10 1398 13.71 153.03 15.32 9.99
2022-02-05 1740 18.22 158.01 18.77 8.42
2022-03-07 2007 16.46 162.03 16.22 9.99
2022-04-10 2283 14.23 162.97 19.40 8.40
2022-06-26 2628 18.97 166.00 18.19 9.13

途中に海外出張が入ったり雨の週末が多かったりして乗る機会が減ってしまい、2〜3ヶ月ぶりの給油になりました。とはいえ距離は乗っているかも。

今日は取手市国道294号線沿いにある家系ラーメン王道家の本店に行ってきました。まだ6月なのに気温は37度になっており、走行中はアスファルトからの照り返しもあるのでバイク内蔵の温度計は41度を指していました。フル装備のプロテクターだと熱中症になりそうで、途中で何回かコンピにに寄って給水をしました。

今回はタンクがだいぶ軽くなるまで走ってからの給油だったのと、このところの燃料高で、これまでで一番ガソリン代が高かったかな。

本日の給油(スーパーカブ110プロ/JA07)

給油日 オドメーター (km) 給油量 (L) 単価 (円/L) 燃費 (km/L) 距離単価 (円/km)
2021-03-12 12042.2 3.08 133.77 45.88 2.92
2021-05-03 12177.4 2.46 139.84 54.96 2.54
2021-06-06 12290.5 1.99 142.21 56.83 2.50
2021-07-31 12417.3 2.62 148.09 48.40 3.06
2021-09-25 12563.7 2.81 148.04 52.10 2.84
2021-11-04 12711.0 2.44 156.97 60.37 2.60
2022-01-03 12814.3 2.18 153.21 47.39 3.23
2022-02-23 12910.1 2.55 158.04 37.57 4.21
2022-03-06 13058.0 2.83 160.07 52.26 3.06
2022-05-26 13182.1 2.46 158.10 50.45 3.09
2022-06-22 13278.4 1.99 165.83 48.39 3.43

燃費の良いカブをときどき走らせているくらいだから困りはしませんが、さらにガソリン単価が上がりました。この単価になったのは2014年以来かと思います。仕事で使っている人や燃費の悪い車に乗っている人は大変だろうなぁ。

線形混合モデルを(Juliaで)理解したい

以前「線形混合モデルを理解したい」という記事を読んで、Rで線形混合モデルを試したことがあるのですが、今回は、Juliaで同じことをやってみたときにハマったところを書き残しておきます。

やりたかったことは、前述の記事からもってきた以下のようなコードです。palmerpenguinsというパッケージに含まれているpenguinsというデータから、NAが含まれるデータを削除して、lme4でモデルを作る、というものです。もちろんRCallを使えば以下のように何も考えずに実行でき、「Julia最高!」となります。

using RCall                            # JuliaからRを呼び出す
R"library(palmerpenguins)"
R"summary(penguins)"
R"library(tidyverse)"                  # Rでパイプ演算子(%>%)を使うために必要
R"usedata <- penguins %>% na.omit()"   # NAデータを削除
R"glmm_model <- lme4::lmer(bill_length_mm ~ flipper_length_mm
    + (1 + flipper_length_mm|species),
    data = usedata)"                   # モデルの作成と分析
R"summary(glmm_model)"
R"lme4::ranef(glmm_model)"             # 回帰係数を表示する

これを、Rを最低限しか使わずにJuliaで実行します。まず、palmerpenguinsパッケージのデータを読み込みます。

using RCall, CSV, DataFrames

R"library(palmerpenguins)"
filepath = R"""palmerpenguins::path_to_file("penguins.csv")"""
penguins = CSV.read(convert(String, filepath), DataFrame; missingstring="NA")
penguins = dropmissing(penguins)

RDatasetsパッケージにはpalmerpenguinsは含まれていません。RDatasetsを使えば同じマシンにインストールされているRのパッケージを参照できるものだと思い込んでいたので、最初にここでハマりました。というわけで、R側にインストールされたパッケージからCSVファイルを読み込む必要があります。Rにインストールされているパッケージを参照して、データファイルへのパスを取ってきてCSV.read()します。path_to_file()の中でダブルクォートを使っているので、ダブルクォート三つで式全体を囲んでいます。式の中のダブルクォートをバックスラッシュでエスケープしてもいいでしょう。

さて、ハマりポイント二つ目が、NAの扱いでした。Juliaでは欠損値をmissingとして扱います。なので、ファイルから読み込むときに欠損値がどのような表記になっているかをmissingstringで指定してあげます。複数の表し方があるときにはmissingstring=["NA", "-999"]といった書き方ができます。

線形混合モデルの分析にはMixedAnova.jllme()が使えそうです。このパッケージはGLM.jlMixedModels.jlFixedEffectModels.jlなどを統合してANOVAをくっつけたような便利パッケージです。以下のように使いましたが、ANOVAを使わないのであればMixedModels.jlだけでも同じことができるかもしれません。

using MixedAnova
mm_init()

lmm1 = lme(@formula(bill_length_mm ~ flipper_length_mm
    + (1 + flipper_length_mm|species)), penguins)

出力されたのは次のようなモデルです。Rのlme4を使ったものと値がやや違いますが、許容範囲なのかな。Rで実行したときはパラメータが収束していないというwarningが出るので、そのせいかもしれません。とはいえ「Juliaは誤りが多くて使えない(Why I no longer recommend Julia)」なんて意見もあるので気を付けないといけませんね。

Linear mixed model fit by maximum likelihood
 bill_length_mm ~ 1 + flipper_length_mm + (1 + flipper_length_mm | species)
   logLik   -2 logLik     AIC       AICc        BIC    
  -795.0174  1590.0348  1602.0348  1602.2924  1624.8836

Variance components:
               Column        Variance   Std.Dev.    Corr.
species  (Intercept)        143.4020628 11.9750600
         flipper_length_mm    0.0038364  0.0619387 -0.95
Residual                      6.5214560  2.5537141
 Number of obs: 333; levels of grouping factors: 3

  Fixed-effects parameters:
─────────────────────────────────────────────────────────
                       Coef.  Std. Error      z  Pr(>|z|)
─────────────────────────────────────────────────────────
(Intercept)        -0.059232   8.14618    -0.01    0.9942
flipper_length_mm   0.221575   0.0416848   5.32    <1e-06
─────────────────────────────────────────────────────────

係数を得るにはlme4と同様にranef()が使えますが、raneftables()のほうが情報が多くて理解しやすいです。

julia> raneftables(lmm1)
(species = (species = String15["Adelie", "Chinstrap", "Gentoo"],
(Intercept) = [9.780118308163642, 4.4103374755445595, -14.1904557836386],
flipper_length_mm = [-0.06844124963200851, 0.0054061201200947765, 0.06303512951155885]),)

以上、環境はJulia 1.7.3 / Windows 10でした。

julia> versioninfo()
Julia Version 1.7.3
Commit 742b9abb4d (2022-05-06 12:58 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i5-10400 CPU @ 2.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS =

本日の給油(スーパーカブ110プロ/JA07)

給油日 オドメーター (km) 給油量 (L) 単価 (円/L) 燃費 (km/L) 距離単価 (円/km)
2021-03-12 12042.2 3.08 133.77 45.88 2.92
2021-05-03 12177.4 2.46 139.84 54.96 2.54
2021-06-06 12290.5 1.99 142.21 56.83 2.50
2021-07-31 12417.3 2.62 148.09 48.40 3.06
2021-09-25 12563.7 2.81 148.04 52.10 2.84
2021-11-04 12711.0 2.44 156.97 60.37 2.60
2022-01-03 12814.3 2.18 153.21 47.39 3.23
2022-02-23 12910.1 2.55 158.04 37.57 4.21
2022-03-06 13058.0 2.83 160.07 52.26 3.06
2022-05-26 13182.1 2.46 158.10 50.45 3.09

いつもだいたい同じセルフスタンドで給油していますが、ここ1年間の間にガソリンの値上げがすごいことになっていますね。20円も上がってしまった。