キャンプ道具を揃えた2023（更新版）

遊び

キャンプ道具を揃えたシリーズの2023年更新版です。これまで購入したものは以下。

そして、2023年に購入したものたち。

amzn.toで始まるリンクはアフィリエイト、それ以外のリンクはアフィではありません。避けたい人は気を付けてください。

寝袋

今年に入るまでは寝袋は家族のものを借りていましたが、キャンプに行くようになると自分用が欲しくなります。モンベルの「ホローバッグ#3（https://amzn.to/43GmlGG）」を購入しました。ファミリーバッグ#3と迷ったのですが、マミー型っぽくも使えるということでホローバッグのほうを選びました。使用可能限界温度が摂氏2度ということで、暑がりの僕なら温度がマイナスに行かない範囲の冬キャンプにも使えそうかもと思いましたが、真冬に使うには寒いと思いますよ……。

もうホローバッグはキッズ用しかカタログに載っていないので、ちょうど終売のタイミングで購入できたようです（アウトレット品はまだ入手可能かもしれません）。

webshop.montbell.jp

バッグ

雨の中のバイク移動に備えて、キャンプ道具を揃えた2022（バッグと小物類）でも書いていた防水バッグを買いました。トロリーシートバッグ DH-746（https://amzn.to/3RCGvME）の上に60 Lの防水シートバッグを載せてもアフリカツインなら楽々です。モールシステムなどは何もついていないシンプルなものですが、容量は大きくて良いです。できるだけ荷物を減らしたいのですが、だんだんと増えてしまいますね。

www.daytona.co.jp

このバッグを購入したこともあって、アフリカツインのパニアケースにも食指が動いたのでした。

marui.hatenablog.com

飯盒のための色々

僕はとくにキャンプ地で何か特別な活動をしたい人ではありません。屋外炊爨してぼーっと過ごすということがしたいので、調理器具についつい手が伸びてしまいます。できるだけコンパクトに収納できる高機能なものを探していたのですが、戦闘飯盒2型が良さそうに思えました。すでにフェニックスライズの戦闘鉄板チタン FB-2 Ti（https://amzn.to/3L4u4H1）は持っていましたし、飯盒にあわせた大きさのFUZZチタン Ti-02（https://amzn.to/3DdlRyh）もCRUNCH VOX（https://amzn.to/3Ozf9YR）も揃えていましたので、実のところソラマメ型の飯盒を導入するのはほぼ確定していたのでした。

camphack.nap-camp.com

ソラマメ型の2合炊き飯盒は、ロスコ、エバニュー、ムースルームワークスなど、各社から似たようなものがいくつか販売されています。でも、買うからには自衛隊の仕様書（防衛省仕様書改正票「飯ごう，2形」DSP S 2036C(1)）に完全準拠したものがいいです。ロスコの「戦闘飯盒2型」とエバニューの「山岳飯盒弐型」は取っ手の端の形状が仕様と異なりますし、ムースルームワークスの「戦闘飯盒2型」は耳金の塗装がされていません。

いろいろと探していたら、東海理研株式会社の製品情報に陸上自衛隊仕様の携帯シャベルと飯ごうの記載があるのを見つけました。一般販売はされていなさそうなのですが、岐阜県関市へふるさと納税すると返礼品として飯盒2型がもらえるという情報を得ました。ふたに「GIFU SEKI」と書かれていてダサい郷土愛を感じるんですが、品物じたいは実際に自衛隊に納入されているものと同等に見えます。すでに関市の返礼品としての扱いは終わっていますが、群馬県太田市が同等品を扱っています。

www.furusato-tax.jp

実際に、この飯盒を入手してからは他の調理なべを使わなくなってしまいました。鍋もフライパンも湯沸かしもこれひとつですべてまかなえます。

アルコールストーブとライター

Amazonで半値以下の販売価格になっていたので、トランギアのアルコールバーナー（https://amzn.to/3RjjmBm）を衝動買いしてしまいました。真鍮製なので、アルミと接触させないように気を付けないといけません。真鍮は銅と亜鉛の合金ですが、銅とアルミを接触させるとアルミに腐食が起こるらしいです。それは真鍮とアルミでも同じなのだとか。チタン製のアルストを狙っていたのですが、衝動買いは事故のようなものなので、しょうがないですね。

www.iwatani-primus.co.jp

フェニックスライズのCRUNCH VOXに突っ込んで飯盒炊爨に使ったり、ミニ五徳をのせて湯を沸かしたりしています。

火をつける道具はイワタニのアウトドアトーチバーナーとマッチしか持っていませんでした。イワタニのコンパクトジュニアバーナーには着火装置がついているので、火をつける専用品は不要だったのですね。一方、アルコールバーナーには着火装置がついていないので、別途SOTOのスライドガストーチ ST-487（https://amzn.to/3PdRflJ）を購入しました。ライター用のガスを注入して使えます。（てっきりCB管から直接ガスを注入できるのかと思っていましたが、差し込み穴の径が異なるようでガスが入りませんでした）

ランタン

小さくて光量があると評判のゴールゼロの価格が落ち着いてきたので導入しました。いざというときに乾電池も使えるレッドレンザー ML4（https://amzn.to/43w1K7W）とも迷ったのですが、最終的には充電池の容量が決め手になりました。レッドレンザーの750 mAhに対してゴールゼロは2600 mAhと3倍以上の差があります*1。容量が大きいほど長時間使えるわけなので、ゴールゼロを選びました。一泊二日のキャンプで電池切れになったことはありません。

スマホにもフラッシュ兼用のライトがついていますが、光量の面では専用品のほうが断然良いです。小さくてどこにでも持ち運べるので、キャンプ以外の場面でも重宝しています。

arinomi.co.jp

ペグとハンマーとグローブ

房総と猪苗代でのキャンプで困ったのが、テント付属のアルミペグが地面に刺さらないという問題でした。どちらも友人と一緒に行っていたので鍛造ペグを数本借りて急場をしのぎましたが、ソロのときにペグが刺さらないと危険です。テントが飛ばされないように、最低限の本数でいいので強いペグが欲しくなりました。また、小型のハンマーではうまくいかないことが明らかになったので、ハンマーも専用設計のものを買うことにしました。

チタン好きなので、ペグについてもチタン一択です。ハンマーは鍛造ペグで有名なエリッゼステークのものにしました。バイクで行くのでショート版。ステンレスヘッドのやつが外装に傷があるとかで安くなっていたので、貯まっていたAmazonポイントで購入。

両方そろうと使い勝手の差は明らかでした。ランタンもそうですが、とくにハンマーはなぜ最初からキャンプ専用品にしなかったのか……。

ペグうちなどの作業をするときにはモンベルの軍手を使っていましたが、いつのまにか失くしてしまったので、手袋専門メーカーが出しているイカツイやつにしてみました。熱いものも触れます。最初のうちは革のにおいが気になるので、陰干ししておくといいです。

*1:iPhone SE 3のバッテリー容量は2018 mAhなので、小型スマホよりも容量が大きいことになります。

ピンクノイズを生成するLISP関数

LISP 音

以前PythonやJSFX用のものも作ったのですが、ピンクノイズはいろいろなところで活躍するので、Common Lisp用にも作っておきます。

Voss-McCartney法

DSP generation of Pink (1/f) Noiseには、(1)ホワイトノイズにフィルタをかけてピンクにする方法、(2)Voss-McCartney法、の二つがあると紹介されています。Voss-McCartney法はホワイトノイズを低域に重みをつけて重ね合わせていく方法です。

具体的には、いくつかの乱数生成器を用意して、それらで1, 2, 4, 8, 16, ...サンプルごとに新しい乱数を作らせ、すべてを加算して出力するというものです。たとえば4サンプルごとの生成をする乱数生成器は、4サンプルに一度だけ乱数を更新しますが、それ以外の3サンプルではその前に生成した乱数を返します。より低い頻度（ここでは1/4の標本化周波数）での乱数生成をしていることになりますね。つまり、1, 2, 4, 8, 16, ...サンプルごと生成して加算するということは、オクターブ低いホワイトノイズを次々と作って足し合わせてピンクノイズを作っていくことになります。2nが大きくなればなるほど低い周波数まで伸びるピンクノイズになります。

Common Lispでの実装

とりあえず標本化周波数を決めておきます。

(defparameter +samprate+ 48000)

Voss-McCartney法の本体です。指定した秒数ぶんのピンクノイズを生成します。返り値ですが、多値を使ってピンクノイズ振幅値の列と各サンプルの時刻の列の二つを返しています。どちらも1次元配列（simple-vector）です。

(defun synth-pinknoise (dur &optional (fs +samprate+))
  "Generates pinknoise (wideband noise having -3dB/oct slope) of specified duration (seconds). Uses Voss-McCartney algorithm. http://www.firstpr.com.au/dsp/pink-noise/

Example:
(wavwrite \"pinknoise.wav\" (synth-pinknoise 2.0))"
  (let* ((lowfreq 20)
         (levels (ceiling (log (/ fs lowfreq) 2)))
         (tt (sample-time dur fs))
         (out (make-array (length tt)))
         (x (make-randn-array levels)))
    (dotimes (n (length tt))
      (dotimes (m levels)
        (if (= (1- (expt 2 m)) (mod n (expt 2 (1+ m))))
            (setf (aref x m) (randn)))
        (setf (aref out n) (+ (randn) (vecsum x)))))
    (values (normalize out) tt)))

sample-timeは、指定した秒数の各サンプルの時刻を計算する関数です。横軸に時間・縦軸に振幅をとってグラフを書く時などに使います。標本化周波数が48000 Hzのときに1秒間の時刻列は、#(0 1/48000 2/48000 ... 47999/48000)というようになります。以下のように定義しています。rangeはPythonの、linspaceはMatlabの関数をそれぞれ真似しています。

(defun sample-time (dur &optional (fs +samprate+))
  "Make a vector of time values of each sample in seconds for specified duration."
  (coerce (mapcar #'(lambda (x) (/ x fs))
                  (range (round (* dur fs)))) 'vector))

(defun range (n)
  "Generates a list of integers from 0 to N-1."
  (linspace 0 (1- n) n))

(defun linspace (x1 x2 &optional (n 100))
  "Generates a list of N (100 by default) equally spaced numbers between X1 and X2."
  (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))))

randnとmake-randn-arrayは標準正規分布に従う乱数を作る関数で、以下のエントリに書いたものです。

marui.hatenablog.com

vecsumはベクトル内の数の総和を計算します。

(defun vecsum (vec)
  (loop as x across vec
        sum x))

ピンクノイズは音波形として使用したいので、出力時には±1の範囲に収めたいです。これをノーマライズといい、以下のnormalizeで実行しています。±1ビチビチに入れてしまうと（0 dBFSにすると）、オーディオ編集ソフトによっては「レベルオーバーしてるよ」という表示になることがあるので-1 dBFSをデフォルトにしていますが、-3 dBFSとかのほうが安全かもしれません。

(defun normalize (x &optional (dB -1.0))
  "Makes a sampled signal normalized to specified level (in dB)."
  (let ((res (make-array (array-dimensions x)
                         :initial-element 0.0))
        (gain (db->amp dB))
        (tmpmax (loop for n from 0 to (1- (length x))
                      maximize (abs (aref x n)))))
    (dotimes (n (length x))
      (setf (aref res n) (* gain (/ (aref x n) tmpmax))))
    res))

(defun db->amp (dB &optional (p0 1.0))
  (* (expt 10 (/ dB 20)) p0))

できあがったピンクノイズをオーディオ編集ソフトに取り込んで30分にまとめたのが以下の動画です。YouTubeにアップロードするときに聴覚符号化で圧縮されてしまうので、本来の音とは少し違いますが……。

youtu.be

同じスペクトルを持つのにクレストファクタが違う波形（Julia Advent Calendar 2023）

JuliaLang Advent Calendar

Julia Advent Calendar 2023の22日目の記事です。毎年Juliaで音に関係したことを書いていて、これまでにこのブログで書いたアドカレ記事は以下の通りです。

今年は、以下の論文を参考にして「同じスペクトルを持つのにクレストファクタが違う波形」を作ってみます。

  • Schroeder, M. R.. “Synthesis of Low-Peak-Factor Signals and Binary Sequences with Low Autocorrelation.“ IEEE Transactions on Information Theory, pp.85–89. January 1970. DOI: 10.1109/TIT.1970.1054411.

クレストファクタ

クレストファクタ（波高率）は「振幅と実効値の比」で計算されます。振幅の幅はゼロから波形ピーク値で測られます（文献によっては負側のピークから正側のピークまでとしているものもあります）。一方で実効値は波形の時間平均に近いもので、以下の式 (root mean square) で計算できます。音の波形のような交流波において振幅値の算術平均を計算すると正負の値が打ち消しあってしまうので、それを防ぐために二乗してから平均をとり平方根でもとの単位に戻すという流れで計算しています。

\displaystyle x_\mathrm{rms} = \sqrt{\frac{1}{N}\sum_{i=1}^{N}\left(x_i\right)^2}

もし波形が±1の範囲に正規化された正弦波（サイン波）だったとすると、振幅値のピークは1、実効値は 1/√2 になります。クレストファクタは振幅と実効値の比なので 1 ÷ 1/√2 = √2 です。正規化波形に対しての実効値は必ず1以下になるので、クレストファクタは1以上の値になります。最もクレストファクタが小さくなる（1になる）のは矩形波のように-1と+1の信号しか含まない信号です。

一般に、定常的な持続音（オルガンやクラリネットなど）のほうがクレストファクタは小さくなりがちで、逆に、大きな振幅値が時々しか現れない音（太鼓のリズムとか）はクレストファクタが大きくなりがちです。正規化波形で考えるので実際に音を聞いたときに感じる音の大きさとは異なりますが、クレストファクタは音ファイルから音量感を計算する指標の一つとして使うこともできます（たとえば以下の論文参照）。

パルス列の生成

例として余弦波の重ね合わせでパルス列を作っていくことを考えます。

\displaystyle r(t) = \sum_{k=1}^{N}\left(\frac{p_{k}}{2}\right)^{1/2} \cos\left(\frac{2\pi k t}{T} + \theta_k\right)

このとき、各余弦波の位相をあわせて（つまり 𝜃ₖ をすべて同じ値にして）重ね合わせると、高いクレストファクタ（論文中ではpeak-factorと書いてある）になります。ここで、𝑝ₖ は 𝑘 番目の倍音のパワーで、∑ 𝑝ₖ = 1を満たす必要があります。以下のコードでは論文中の数式に従って、第16倍音まで計算しています。

# number of harmonics
N = 16

# relative power of each harmonics
p = zeros(N)
for n in 1:N
    p[n] = 1/8 * (sin(π * (2n - 1) / 32))^2
end
p ./= sum(p)

# phase angles
ϕ = zeros(N)
for n in 1:N
    ϕ[n] = -π/2
end

# signal generation (Eq.1)
fs = 1000
T = 1
t = range(-π/2, stop=π/2, step=1/fs)
r = zeros(size(t))
for k in 1:N
    r .+= (p[k] / 2)^(1/2) * cos.(2π * k * t / T .+ ϕ[k])
end

plot(t, r, xlabel="Time", ylabel="Amplitude", legend=false, ylim=(-2.5, 2.5))

位相を調整する

測定信号や通信用の信号を生成したいときには、一般にクレストファクタは低ければ低いほどよいです。許された振幅を十分に使えるということは、信号・雑音比（S/N）を低くすることができるし、たとえばFM変調で伝送するときを考えると周波数帯域幅をまんべんなく使えるということですからね。

できるだけ全体にわたってピーク値が均等になるように、各倍音の位相を倍音の強さに応じて調整します。

\displaystyle \phi_n = \phi_1 - 2\pi\sum_{l=1}^{n-1}\left(n-l\right)p_l

それをプログラムにすると以下のようになります。

# phase angles
ϕ = zeros(N)
ϕ[1] = -π/2
for n in 2:N
    tmp = 0
    for l in 1:(n-1)
        tmp += (n - l) * p[l]
    end
    ϕ[n] = ϕ[1] - 2π * tmp
end

# signal generation (Eq.1)
fs = 1000
T = 1
t = range(-π/2, stop=π/2, step=1/fs)
r = zeros(size(t))
for k in 1:N
    r .+= (p[k] / 2)^(1/2) * cos.(2π * k * t / T .+ ϕ[k])
end

plot(t, r, xlabel="Time", ylabel="Amplitude", legend=false, ylim=(-2.5, 2.5))

パワースペクトルを調べると、各倍音の強さは両方とも同じになっています。位相がずれているので異なる信号なのですが、定常音についてヒトは位相に鈍感なので、実際にスピーカから再生すると同じ音に聞こえます。波形はまったく違うのに、同じ音に聞こえるって不思議ですね。

アフリカツインにパニアケース購入

CRF1100Lアフリカツイン バイク

少しずつ快適なキャンプになるようにと考えていると、なぜかキャンプ道具が増えていきます。容量45リットルのシートバッグの上に、容量60リットルのソフトな防水シートバッグを載せて使っていましたが、車体重心が高くなってしまうことがネックでした。下の写真右側のように、さらにその上に30リットルのバッグを付けたときは、さすがにヤバかったので短距離しか走行しませんでした（リアキャリアから左右15 cmは出っ張ってないし、地上2 m以内だし、合計60 kgも超えてないので、法令違反にはなってないはず）。

シートバッグの上にシートバッグを重ねると重心が高くなって危険

やはり安全かつしっかりと積載したいので、ぼっち女campさんの積載動画を参考にしつつパニアケースを導入することにしました。ぼっち女campさんはYamaha YZF-R3とKTM 250 DUKEなので、ちょっと真似するだけでアフリカツインのほうが多く積載できるはずです。

購入したパニアケースは以下の、バイク本体と鍵が共有できるのが特徴のホンダ純正品です。なかなかお高いうえに、錠の設定をしてもらう必要があるので工賃もかかります（工賃はそれほどでもないですが）。

hondago-bikegear.jp

このケースはGIVIによるアフリカツイン専用設計とのことです。もともとアフリカツインは、バイク本体のリアキャリア付近に差し込み穴が開いており、ケースが取り付けられるようになっています。そこにパニアケースから出ているアームをひっかけてロックできる仕組みです。そのためケースを取り付けるためのステーが不要なので、パニアケースを使わないときには純正の形のままで乗ることができます。

容量は、左側が40リットル、右側が30リットルです。右側が小さいのは、マフラーに干渉しないような形状になっているからです（上の写真でも右側パニアが若干えぐれているのがわかります）。マフラーからの熱が心配なので、熱に耐えられるものしか入れないようにしたほうがいいですね。また、ケースのフタは横向き（箱の下部にヒンジがあって、箱の上部が割れてフタが車体から左右に離れるように）に開くので、車体につけたまま開くときには中身が落ちないように気を付けないといけません。

このプラスチック製ケースとは別に、アルミケースも純正オプション設定されているんですが、あちらはフタが上向きに開くので、パニアケースにかぶさるくらい大型のシートバッグを付けた状態だとフタが開かないのではないかと思いました。また、アルミ特有の白錆が出てくるという話を聞いたのと、取り付けのためにステーが必要なのも残念ポイントでした。見た目はアドベンチャーバイクっぽさが向上するアルミのほうがいいんですけどね。

使ってみた感想としては、まず「大きくてかっこいい」「荷物の重心が低くなるので安定感がある」です。また、パニアケースは積載装置と考えてよいそうなので、パニアケースの上に乗るくらいのシートバッグも積載できそうです。手で持つためのハンドルがついているのも便利なポイントです。

逆に、気になったことのは「開錠してるときには鍵が抜けない」というのがあります。バイクを動かすときにはケースは施錠されていないといけないので、安全面を考えると良いことではあるのですが、左右両方のケースを同時に開けたいときにどうすればいいのか……。（→フタを開けたまま施錠状態にすれば鍵が抜けるので、左右ケースともに開けっ放しにして対応しました）

重回帰分析 (Lisp Advent Calendar 2023)

LISP Advent Calendar

Lisp Advent Calendar 2023の第1日目の記事です。Lispについては全くのヘナチョコなので1日目は恐縮なのですが、カレンダーの先頭が空になっているのも幸先が悪いので……。

さて、車輪の再発明を続けて、ようやく逆行列の計算までたどりつきました。今回は逆行列を使った応用をしたいと思います。

重回帰分析

逆行列の応用例として、重回帰分析を取り上げます。

重回帰分析は、ひとつの量的な応答変数を、複数の量的な説明変数で予測・説明するための線形モデルです。例えば、アパートの家賃（応答変数、単位：円）を、駅からの徒歩時間（説明変数、単位：分）と延べ床面積（説明変数、単位：m2）から推測するようなときに用いることができます。数式は 𝑦 = 𝛽₀ + 𝛽₁ 𝑥₁ + 𝛽₂ 𝑥₂ + ... + 𝛽ₙ 𝑥ₙ + 𝜖 のように書かれます。家賃の例の場合、家賃 = 係数0 + 係数1×徒歩時間 + 係数2×面積 + 誤差のようになるでしょうか。係数0は賃貸物件としてのベースラインの金額、係数1は徒歩時間が1分増えるごとに増える（あるいは減る）金額、係数2は延べ床面積が1 m2増えるごとに上がる金額を表します。また、誤差には時間と面積だけでは説明しきれない要素（築年数、階数、インターネット回線の品質など）や物件ごとの個体差が含まれます。多くのデータをもとに係数の値を求めるのが重回帰分析です。係数0と係数1だけを求めればいい場合（説明変数がひとつしかないとき）を単回帰分析と区別することもありますが、やっていることに違いはありません。

係数が決まったとして、そこから計算できる予測値は実際の値からは多少ずれます。そのずれの総量を最小化すれば最適な回帰式が得られるはずです。その発想で係数を求めるのが最小二乗法です。イメージとしては、下のグラフのように各データ点から近似直線までの距離（つまり誤差）の二乗和を最小にするように近似直線を定めるかんじです。

誤差の二乗和を最小にするためには微分を使って計算を進めることになります。やたらとめんどくさそうですが、実は簡単な行列計算に帰着します。英語版WikipediaのLinear Regressionの項に詳しく書かれていますが、まず説明変数𝑥₁ ... 𝑥ₙを以下のような行列にまとめます。ここで𝑛は説明変数の数、𝑚はデータの数となっています。行列の1列目に並んでいる1は𝛽₀と掛け算されるために追加されています。このようにすると 𝑦 = 𝛽₀1 + 𝛽₁ 𝑥₁ + 𝛽₂ 𝑥₂ + ... + 𝛽ₙ 𝑥ₙ + 𝜖 を 𝑦 = 𝑋𝛽 + 𝜖 と書け、𝑋𝛽 の部分は行列とベクトルの積で一発計算ができるようになるわけです。

\displaystyle X = \begin{bmatrix} \mathbf{x}^\mathsf{T}_1 \\ \mathbf{x}^\mathsf{T}_2 \\ \vdots \\ \mathbf{x}^\mathsf{T}_m \end{bmatrix} = \begin{bmatrix} 1 & x_{11} & \cdots & x_{1n} \\ 1 & x_{21} & \cdots & x_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{m1} & \cdots & x_{mn} \end{bmatrix}

ここで予測誤差である 𝜖 = 𝑦 - 𝑋𝛽 の総和を最小にしたいです。𝜖 の値は正にも負にもなるので、単純に総和を計算するとゼロになってしまいます。そこで、𝜖 の二乗和を使います。つまり、𝑋𝛽 で予測した値と、データで与えられている応答変数 𝑦 との差の二乗和ということです。以下の式を微分して最小になる 𝛽 を探ります。

\displaystyle \|X\vec{\beta} - Y\|^2

ここでは導出は省略しますが、実は以下のような行列の転置、逆行列、行列の積という操作だけで 𝛽 の最適値が求まってしまいます。導出については前述のWikipedia記事を参照してください。

\displaystyle \vec{\hat{\beta}} = \left(X^\textsf{T} X\right)^{-1} X^\textsf{T} Y

Lisp実装

この式を単純にCommon Lispコードにしたのが以下です。

(defun mra (y X)
  (let ((XX (make-array `(,(array-dimension X 0) ,(1+ (array-dimension X 1))) :initial-element 1)))
    (loop for i from 0 below (array-dimension X 0)
          do (loop for j from 0 below (array-dimension X 1)
                   do (setf (aref XX i (1+ j)) (aref X i j))))
    (matrix-multiply (matrix-multiply (inv (matrix-multiply (matrix-transpose XX) XX)) (matrix-transpose XX)) y)))

なんの工夫もない、単純なコードですね。まず、第1列に1が並び、それ以降にXの内容が入っているXX行列を作り、それに上記式を適用した結果を返すだけです。matrix-multiplyとmatrix-transposeはこのブログでは紹介したことはありませんが、二次元配列で与えられた行列についての積と転置の単純な実装になっています。

実行例

南風原朝和『心理統計学の基礎 https://amzn.to/48K5fKY』（有斐閣アルマ，2002年）の第8章に掲載されている例を使って実行してみましょう。この本はt検定、分散分析、回帰分析、因子分析などの背景について、数式を使いつつも直感的に理解できるように書かれている良書です。

p.226に掲載されている「母親の協調性価値得点 (𝑥₁)」「通年年数 (𝑥₂)」「子どもの協調性得点 (𝑦)」を抜き出しました。まずは、Rで実行してみます。重回帰分析はモデル式をlmで評価させるだけです。

kachi <- c(12, 12, 7, 17, 14, 9, 10, 13, 15, 12, 12, 15, 11, 14, 17, 17, 16, 15, 15, 10, 12, 9, 12, 12, 19, 11, 14, 15, 15, 15, 16, 15, 12, 10, 11, 12, 15, 13, 15, 12, 12, 12, 13, 17, 13, 11, 14, 16, 12, 12)
nensu <- c(2, 2, 2, 3, 2, 2, 3, 3, 3, 1, 3, 3, 2, 2, 4, 2, 4, 3, 4, 2, 2, 1, 2, 2, 4, 2, 3, 2, 3, 3, 2, 3, 2, 2, 3, 1, 2, 3, 2, 2, 2, 3, 3, 3, 2, 3, 2, 4, 2, 2)
kyouchou <- c(6, 11, 11, 13, 13, 10, 10, 15, 11, 11, 16, 14, 10, 13, 12, 15, 16, 14, 14, 8, 13, 12, 12, 11, 16, 9, 12, 13, 13, 14, 12, 15, 8, 12, 11, 6, 12, 15, 9, 13, 9, 11, 14, 12, 13, 9, 11, 14, 16, 8)
res <- lm(kyouchou ~ kachi + nensu)
summary(res)

以下のように、子どもの協調性得点 = 5.1716 + 0.3027×母親の協調性価値得点 + 1.1259×通園年数という結果になりました。

Call:
lm(formula = kyouchou ~ kachi + nensu)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.0563 -1.4780 -0.0563  1.4751  4.9437 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   5.1716     1.6647   3.107   0.0032 **
kachi         0.3027     0.1467   2.064   0.0446 * 
nensu         1.1259     0.4713   2.389   0.0210 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.12 on 47 degrees of freedom
Multiple R-squared:  0.3137,    Adjusted R-squared:  0.2845 
F-statistic: 10.74 on 2 and 47 DF,  p-value: 0.000144

今回のLisp実装でこの3つの係数が得られれば成功ということにします。

(defparameter *kachi* #(12 12 7 17 14 9 10 13 15 12 12 15 11 14 17 17 16 15 15 10 12 9 12 12 19 11 14 15 15 15 16 15 12 10 11 12 15 13 15 12 12 12 13 17 13 11 14 16 12 12))
(defparameter *nensu* #(2 2 2 3 2 2 3 3 3 1 3 3 2 2 4 2 4 3 4 2 2 1 2 2 4 2 3 2 3 3 2 3 2 2 3 1 2 3 2 2 2 3 3 3 2 3 2 4 2 2))
(defparameter *kyouchou* #(6 11 11 13 13 10 10 15 11 11 16 14 10 13 12 15 16 14 14 8 13 12 12 11 16 9 12 13 13 14 12 15 8 12 11 6 12 15 9 13 9 11 14 12 13 9 11 14 16 8))

(mra (matrix-transpose *kyouchou*) (cbind (matrix-transpose *kachi*) (matrix-transpose *nensu*)))

matrix-transposeは行列の転置ですが、入力がsimple vectorの場合はn×1のarrayにします。cbindは二つの2次元arrayを横方向につなげる関数です。この結果、#2A((64133/12401) (11263/37203) (13962/12401))が得られました。64133÷12401 = 5.171599、11263÷37203 = 0.3027444、13962÷12401 = 1.1258769なので、Rでの計算結果と同じになりました。

ここでは決定係数（R2）の計算はしていませんが、𝑋𝛽の分散を𝑦の分散で割ると得られます。今回のデータでは0.31367722です。

逆行列の計算

LISP

Common Lispで逆行列の計算をするための準備をしてきました。

今日はついに逆行列の計算です。とは言っても、効率の良い計算をするわけではなく、単純な𝑂(𝑛³)での実装です。しかもsolve関数を呼ぶたびにLUP分解が実行されるのでかなり効率が悪いです。でも、まずは効率のことは無視して進めましょう。

前回までの記事で、LUP分解を使って連立方程式を解けるようになりました。連立方程式は、行列𝐴とベクトル𝒃に対して𝐴𝒙=𝒃が成立するベクトル𝒙を求めるということになります。ここで、𝐴の逆行列𝑋があったとすると、𝐴𝑋=𝐸ₙが成立します（𝐸ₙは𝑛×𝑛の単位行列）。𝐸ₙの一列目（[1; 0; 0; ...; 0]）は、𝐴と𝑋の1列目の積になります。つまり𝐴𝒙₁=𝒆₁ということなんですが、これはすべての列にあてはまるので𝐴𝒙ₖ=𝒆ₖとなります。この連立方程式を解けば、逆行列𝑋の𝑘列目である𝒙ₖが求まるというわけです。これを𝑛列すべてに適用すれば逆行列全体が求まります。

(defun inv (A)
  "Compute an inverse matrix of A.
Reference: Cormen, Leiserson, Rivest, and Stein. Introduction to Algorithms, 4ed. MIT Press. 2022."
  (assert (/= (det A) 0) nil "Input matrix is singular.")
  (let* ((n (array-dimension A 0))
         (x (make-array n))
         (Ainv (make-array `(,n ,n) :initial-element 0)))
    (loop for k from 0 below n
          do (let ((ek (make-array n :initial-element 0)))
               (setf (aref ek k) 1)
               (setf x (solve A ek)))
             (loop for j from 0 below n
                   do (setf (aref Ainv j k) (aref x j))))
    Ainv))

さっそく使ってみましょう。

CL-USER> (inv #2A((2 5) (1 3)))
#2A((3 -5) (-1 2))

Matlabだと以下のようになります。同じ結果ですね。

>> inv([2 5 ; 1 3])
ans =
     3    -5
    -1     2

さて、試してみて驚いたことに、このプログラムに手を加えずに複素数が入った行列の計算もできるようなのです。

CL-USER> (inv #2A((2.0d0 5.0d0) (1.0d0 #C(3.0d0 4.0d0))))
#2A((#C(0.5384615384615384d0 -0.3076923076923077d0)
     #C(-0.07692307692307693d0 0.6153846153846154d0))
    (#C(-0.015384615384615385d0 0.12307692307692308d0)
     #C(0.03076923076923077d0 -0.24615384615384617d0)))

>> inv([2 5 ; 1 3+4i])
ans =
  0.538461538461538 - 0.307692307692308i -0.076923076923077 + 0.615384615384615i
 -0.015384615384615 + 0.123076923076923i  0.030769230769231 - 0.246153846153846i

Matlabと同じくらいの精度。こんなことってある？　ビックリしました。

前述したように、そもそも計算オーダーは𝑂(𝑛³)ですし、1回しか呼び出さなくていいLUP分解を𝑛回も呼び出しているので無駄が多いです。でも、今回はここまでにします。

本日の給油（スーパーカブ110プロ／JA07）

スーパーカブ110プロ
給油日 オドメーター (km) 給油量 (L) 単価 (円/L) 燃費 (km/L) 距離単価 (円/km)
2022-06-22 13278.4 1.99 165.83 48.39 3.43
2022-07-27 13405.0 2.21 157.92 57.29 2.76
2022-08-31 13546.7 2.92 156.16 48.53 3.22
2022-11-05 13699.4 3.21 157.01 47.57 3.30
2023-02-06 13848.1 3.26 157.06 45.61 3.44
2023-03-02 13973.9 2.60 156.92 48.38 3.24
2023-04-01 14095.9 2.99 156.86 40.80 3.84
2023-07-13 14218.0 2.73 163.00 44.73 3.64
2023-09-25 14380.5 3.60 176.94 45.14 3.92
2023-11-25 14529.5 2.93 162.12 50.85 3.19

ガソリン代は落ち着いてきたけどまだ高いですね。長距離やパワーが必要なときにはアフリカツインを使えばいいという気持ちなので、のんびりトコトコ走るのにスーパーカブを使っています。2台持ちだとそういった使い分けができていいですね。