Maximaを活用した数学学習

in [数ナビの部屋]

「Maxima」を活用した数学学習を取りまとめました.

「Maxima」を活用して数学の学習ロードを駆け抜けよう!

  • MathJaxを使用しているので、 スマホでは表示に時間がかかることがあります。
  • モバイル利用(Android)でのメニュー選択は、 SiteMapを利用するか、 または「長押し」から「新しいタブを開く」を選択してください。
■ 数式処理ソフト「Maxima」を活用した数学学習


[御案内] 「Maxima」(マキシマ)は,フリーの数式処理ソフトです. 有料の Mathematica や Maple に劣らないレベルの数式処理が可能であり, Linux,Windows,MacOSのみならず,Android版もあります. ここでは,数学学習での Maxima の活用法について解説します.

[お知らせ] スマホ(Android)版Maximaの解説本を出版しました. 計算問題やグラフの確認をするときに非常に重宝します. フリーソフトなので一度試してみてください. PC版のコマンドレファランスとしても利用できます。
いつでも・どこでも・スマホで数学!
試し読み 森北出版 amazon 楽天 honto 7net 紀伊國屋電子版読書メーター

■数学学習での活用
  • 以下では,「TeXmacs」+「Maxima」の画面で基本的な使い方を解説します.
    詳しい解説は,リンク集に登録したサイトを参照してください.
  • 世の中の多くのことは、 「正規分布」ではなく「ベキ分布」に従っているようです。ベキ分布:リンク集」も参照してください。
確率・統計
 統計に特化したソフトウェアとしては,SPSSSASRなどが有名ですが, Maximaにも統計向けのパッケージが備わっています.

■回帰分析

 回帰分析では, 互いに関連性があると思われる2種類のデータ \(\small X, Y\) があるとき, それらの関係を最もよく近似すると思われる関数 \(\small Y=f(X)\) を 求めようとします. その関数を1次式で求めるのが単回帰分析です. 説明変数が2変数以上のときは重回帰分析と呼ばれます. いずれも1次式による線形回帰です。パッケージ「stats」を読み込むと, これらの計算を行うコマンド「linear_regression」を利用することができます. 指数関数などの1次式ではない関数で近似する場合は, 別のコマンド「lsquares_estimate」を利用します.

単回帰分析

 単回帰分析では, 目的変数を \(\small Y\), 説明変数を \(\small X\) とするとき, 与えられたデータを「最も良く近似する1次式」を求めます. 「最も良く近似する」ということをどのように考えるかが問題ですが, 設定した式で推定される値と実際の値との差を求めて, その平方和が最小になるように定めます. この考え方は「最小二乗法」と呼ばれ,いろいろな場面で利用されます.

最小二乗法
 コマンド(linear_regression)に データ(行列)を入力すると分析結果が得られますが, ブラックボックスとして利用するのではなく, どのようにして結果が得られているのかも 理解した上で利用すべきでしょう. ここでは,基本となる「最小二乗法」の考え方について解説します.

 今,\(\small (X, Y)\) の具体的なデータを \(\small (x_1, y_1), (x_2, y_2), \ldots, (x_n, y_n)\) とすると,\(\small Y=b+mX\) による \(\small x=x_i\) のときの \(\small y_i\) の推測値は \(\small \hat{y}_i=b+mx_i\) です. 最小二乗法は,実際の値との差の平方 \(\small (y_i-\hat{y}_i)^2\) を合計して

\[\small \sum_{i=1}^{n}\left\{y_i-(b+mx_i)\right\}^2\] が最小になるように \(\small b, m\) の値を定めます. \(\small (x_i, y_i)~(1\leq i\leq n)\) の値は与えられるので, この式は \(\small b, m\) の関数です. それを \(\small F(b, m)\) とおくと, 結局は2変数関数 \(\small F(b, m)\) の最小値を求める問題になり, そのような点では極値になっていると考えられます. 1変数関数 \(\small f(x)\) の場合は, \(\small x=a\) で極値を取れば \(\small f'(a)=0\) が成り立ちます. 偏微分法(多変数関数の微分法)で学ぶように, 2変数関数でも同様のことが成り立ちます. つまり,\(\small F(b, m)\) が極値をとれば, その点で \(\small b, m\) で微分するといずれも偏導関数の値は\(0\)になり \[\small \frac{\partial}{\partial b}F(b, m)=0,\quad \frac{\partial}{\partial m}F(b, m)=0\] が成り立ちます. \(\small \partial\) は,2変数以上の場合の微分を表す記号です.
 ここで,実際に \(\small F(b,m)\) を微分すると,合成関数の微分法により \[\begin{align*} \small \frac{\partial F}{\partial b} &\small =2\sum_{i=1}^{n}\left\{y_i-(b+mx_i)\right\}\cdot (-1)~ (1)\\ &\small =-2\left(\sum_{i=1}^{n}y_i-bn-m\sum_{i=1}^{n}x_i\right)\\ &\small =0\\ \small \frac{\partial F}{\partial m} &\small =2\sum_{i=1}^{n}\left\{y_i-(b+mx_i)\right\}\cdot (-x_i)\\ &\small =-2\left(\sum_{i=1}^{n}x_iy_i-b\sum_{i=1}^{n}x_i\right.\\ &\small    \left. -m\sum_{i=1}^{n}x_i^2\right)=0\quad (2) \end{align*} \] となります.したがって,\(\small b, m\) は次の連立1次方程式を満たします. \[ \begin{cases}\small\displaystyle m\sum_{i=1}^{n}x_i^2+b\sum_{i=1}^{n}x_i=\sum_{i=1}^{n}x_iy_i ~ (3)\\ \small \displaystyle m\sum_{i=1}^{n}x_i+bn=\sum_{i=1}^{n}y_i\qquad (4) \end{cases}\] これを \(\small m\) について解くと次のようになります. \[\small m=\frac{\displaystyle n\sum_{i=1}^{n}x_iy_i-\sum_{i=1}^{n}x_i\sum_{i=1}^{n}y_i} {\displaystyle n\sum_{i=1}^{n}x_i^2-\left(\sum_{i=1}^{n}x_i\right)^2} \]  ここで, 共分散 \(\small s_{xy}=\displaystyle \frac1{n}\sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y})\) は, \[\small s_{xy}=\frac1{n}\sum_{i=1}^{n}x_iy_i-\bar{x}\bar{y} \] と変形することができ,分散 \(\small\displaystyle \sigma_x^2= \frac1{n}\sum_{i=1}^n(x_i-\bar{x})^2\) は \[\small \sigma_x^2=\frac1{n}\sum_{i=1}^{n}x_i^2-(\bar{x})^2\] と表すことができることに注意すると, \(\small m\) の分子・分母を \(\small n^2\) で割ることにより \[\small m=\frac{s_{xy}}{\sigma_x^2}\] と簡潔な式で表されます.
 一方, \(\small b\) の値は,(4)より \[\begin{align*}\small b&\small =\frac1{n}\sum_{i=1}^{n}y_i-m\cdot \frac1{n}\sum_{i=1}^{n}x_i\\ &\small =\bar{y}-m\cdot\bar{x}\end{align*}\] となります. つまり,2変量 \(\small X, Y\) について, \(\small \displaystyle \sum_{i=1}^nx_i, \sum_{i=1}^ny_i, \sum_{i=1}^nx_i^2, \sum_{i=1}^nx_iy_i\) の値が分かれば, \(\small b, m\) の値を計算できることになります.
 ここで, 実際の値と推測値との差の平方の合計を \(\small n\) で割った平均の式 \[\begin{align*} &\small \frac1{n}\sum_{i=1}^{n}(y_i-\hat{y}_i)^2\qquad (5)\\ &\small =\frac1{n}\sum_{i=1}^{n}\left\{y_i-(b+mx_i)\right\}^2 \end{align*}\] に,今求めた \(\small b, m\) の値を代入して \[\small \frac1{n}\sum_{i=1}^{n}\left\{(y_i-\bar{y}) -\frac{s_{xy}}{\sigma_x^2}(x_i-\bar{x})\right\}^2\] を計算すると, \[\begin{align*} &\small \sigma_x^2=\frac1{n}\sum_{i=1}^{n}(x_i-\bar{x})^2\\ &\small \sigma_y^2=\frac1{n}\sum_{i=1}^{n}(y_i-\bar{y})^2\\ &\small s_{xy}=\frac1{n}\sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y}) \end{align*} \] であることなどから,(5)は \[\begin{align*} \small \sigma_y^2-\frac{s_{xy}^2}{\sigma_x^2} =\sigma_y^2\left\{1-\left(\frac{s_{xy}}{\sigma_x\sigma_y}\right)^2\right\} \end{align*}\] と変形することができます. ここに現れた \(\small \frac{s_{xy}}{\sigma_x\sigma_y}\) が 相関係数 \(\small r\) です.\(\small r\) の値が 1 に近いほど 推測値と実際の値との差が小さいことなり, データが一つの直線 \(\small y=b+mx\) の近くに分布していることを 示しています.この直線が回帰直線であり, その係数 \(\small b, m\) が回帰係数です.
 以上の詳細は, 統計の教科書やWebサイトを 参照してください.

▲戻る(トップメニュー)

回帰係数の計算
 具体的なデータで,この回帰直線 \(\small y=b+mx\) を求めてみましょう. 手計算での検算がしやすいように \(\small X, Y\) の値を次の数値で考えます.

x 1 2 3 4 5
y 2 1 4 2 4

 ここでは, Maximaを紙と鉛筆の代わりに利用して \(\small b, m\) を求めてみましょう.

 最初に, 「linear_regression」が登録されているパッケージ「stats」を読み込み, サンプルデータを(%i4)(%i4)によりリストで入力します. データ数が多いときは, あらかじめ csv ファイルで作成して「load」で 読み込みます. そして,Maximaの入力書式に合わせて,(%i5)で転置した行列を「dat」とします. 「linear_regression」に入力するデータは, 最終列を目的変数とする行列である必要があります.
 次に,このデータから基本統計量を一つずつ計算していきましょう. 「いろいろな統計量」で解説したように, 平均は「mean」,分散は「var」により求めることができます.

 (%i6)では,データの大きさを求めています.それは,要するに行の数です. (%i7)(%i8)では,\(\small \Sigma{x_i}, \Sigma{y_i}\) を求めています. (%i9)(%i10)では,\(\small \Sigma{x_i^2}, \Sigma{y_i^2}\) を求めています. 「xdat^2」により各要素が2乗されるので,その平均を5倍すれば2乗和を 求めたことになります. (%i11)は積和 \(\small \sum{x_iy_i}\) を求めています. リストとリストの積では,対応する要素の積が計算されます.

 (%i12)は,xdatの分散を求めています. (%i13)では,同じ分散を \(\small \frac1{n}\sum{x_i^2}-\bar{x}^2\) により 求めて一致することを確認しています.「std」は標準偏差で分散の平方根を とったものです.

 (%i17)では,共分散をその定義である \(\small \frac1{n}\sum(x_i-\bar{x})(y_i-\bar{y})\) により求めています. (%i19)では,\(\small \frac1{n}\sum{x_iy_i}-\bar{x}\bar{y}\) により求めて 一致することを確認しています.
 ここではコマンドの使い方の練習として 各列の平均や分散を個別に求めましたが, 元データである行列「dat」を入力すると, (%i20)(%i21)のように,各列の平均や分散が一気に求められます.

 (%i24)(%i25)では,前述で求めた回帰係数の式 \(\small m=s_{xy}/\sigma_x^2, b=\bar{y}-m \bar{x}\) から \(\small m, b\) を求めています. したがって,求める回帰直線は \(\small y=\frac{11}{10}+\frac12x\) です. 下記により描画されます.

 統計ソフトを利用すると,この回帰直線はコマンド一つで求められます. それをどのようにして求めているのかを理解する上では, 以上のような個々の変量の基本統計量などを自分で計算しながら求めることも 必要と思います.その計算過程は,Maximaを, まさに「紙と鉛筆」の代わりに利用していることになります.

▲戻る(トップメニュー)

linear_regression
 次に,「stats」パッケージに登録されている 線形回帰を求めるコマンド 「linear_regression」を利用して求めてみましょう.

 「Stats」パッケージは,すでに(%i1)で読み込み済みです. 数式処理ソフトであっても,このコマンドに入力する行列は 具体的な数値である必要があります. いろいろな結果が小数で表示され,デフォルトでは16桁で表示されます. そのままだと桁数が多すぎて結果が読みにくいので, ここでは桁数を制限しています.表示桁数は, (%i28)にある「fpprintprec」により指定することができます.
 線形回帰を行うには, (%i29)のように(%i2)で定めたデータ行列「dat」を入力するだけです. このコマンドで3列以上の行列を入力すると重回帰分析を行います. いずれも,行列の最後の列を目的変数とする必要があります. 結果は(%o29)のような多岐の内容で表示されますが, これは重回帰分析に対応した出力結果です. マニュアルには記載されていませんが, 単回帰分析 \(\small y=a+bx\) の場合は, 後述する 「simple_linear_regression」を利用した方が 出力結果が分かりやすいようです. また,(%o29)は,「Texmacs+Maxima」での出力です. 最新版のwxMaximaでは,このうち 1,4,5,6,7,8,9,11 の値だけが 表示されるように修正されています.
 なお,(%o29)の内容について, マニュアルには次のように記載されています.
  1. b_estimation: 回帰係数推定
  2. b_covariances: 回帰係数推定の共分散行列
  3. b_conf_int:回帰係数の信頼区間
  4. b_statistics:係数テストの統計
  5. b_p_value:係数テストのp値
  6. b_distribution:係数テストの確率分布
  7. v_estimation:不偏分散推定量
  8. v_conf_int:分散信頼区間
  9. v_distribution:分散テストの確率分布
  10. residuals:残差
  11. adc:調整済み決定係数
  12. aic:赤池情報量基準
  13. bic:Bayes情報量基準
 以上のような項目名しか書かれていないので, 統計初心者には何のことなのか全く意味不明と思われます. マニュアルはコマンドの使い方を解説しているだけであり, そこでの数学の内容を解説しているわけではありません. 回帰係数以外の項目は,係数や残差の区間推定や依拠した確率分布に 関する情報なので,推定・検定に関する知識が必要になります.
 さしあたっては,回帰係数 \(\small b, m\) の値が求められればよいかと 思います.それが,最初に書かれている「1. b_estimation」の箇所にあります. リストで [b, m] の結果が書かれています. 回帰直線は求められても, それが実際の値をどの程度近似しているのかが気になります. それは,「10. residuals」(残差)の箇所に書き込まれています. 予測値を \(\small \hat{y}\) とするとき, \(\small y-\hat{y}\) の値がリストで出力されています.たとえば, \(\small y_1=2\), \(\small \hat{y}_1=11/10+1/2\cdot 1=8/5\) なので, \(\small y_1-\hat{y}_1=2-8/5=2/5=0.4\) です.
 個別の誤差だけ示されても,全体としての近似の良さはどの程度なのでしょうか. それを表す指標として決定係数と呼ばれるものがあります. それは,予測値を \(\small \hat{y}\) とするとき \[\small R^2= 1-\frac{\sum_{i=1}^{n}(y_i-\hat{y}_i)^2}{\sum_{i=1}^{n}(y_i-\bar{y})^2}\] により定義される値です. 分母は実測値の変動を表し,分子は予測値との「ずれ」の変動を表すので, \(\small R^2\) の値が 1 に近いほど分子(つまり「ずれ」)が小さくて 良い予測式といえます.分母は全変動,分子は残差変動と呼ばれています. この値は, \[\small R^2= \frac{\sum_{i=1}^{n}(\hat{y}_i-\bar{y})^2}{\sum_{i=1}^{n}(y_i-\bar{y})^2}\] と表すこともできます.分子は予測値の変動を表します.
 この2つの式が一致することは,次のようにして分かります. まず,回帰直線を求めるための前提とした式(1)(2)から \[\begin{align*} &\small \sum_{i=1}^{n}(y_i-\hat{y}_i)=0\\ &\small \sum_{i=1}^{n}(y_i-\hat{y}_i)x_i=0 \end{align*}\] であることに注意すると,\(\small \hat{y}_i=b+mx_i\) に対して \[\small \sum_{i=1}^{n}(y_i-\hat{y}_i)\hat{y}_i=0\] が得られます.さらに, \[\begin{align*} &\small (y_i-\bar{y})^2\\ &\small =\left\{(y_i-\hat{y}_i)+(\hat{y}_i-\bar{y})\right\}^2\\ &\small =(y_i-\hat{y}_i)^2+2(y_i-\hat{y}_i)(\hat{y}_i-\bar{y})\\ &\small \hspace{10mm}+(\hat{y}_i-\bar{y})^2 \end{align*}\] と変形できて,右辺の第2項を \[\small 2(y_i-\hat{y}_i)\hat{y}_i-2(y_i-\hat{y}_i)\bar{y}\] と分けて総和を取ると 0 になることから, 次の式が成り立ちます. \[\begin{align*} \small \sum_{i=1}^{n}(y_i-\bar{y})^2 &=\sum_{i=1}^{n}(y_i-\hat{y}_i)^2\\ &\small \qquad+\sum_{i=1}^{n}(\hat{y}_i-\bar{y})^2 \end{align*}\]  右辺の第1項を左辺に移項して \(\small \sum{(y_i-\bar{y})^2}\) で 割れば,\(\small R^2\) の2つの式が一致することが分かります. 特に, 単回帰(\(\small y=b+mx\)) の場合は相関係数の2乗と一致します (参照1参照2). ただし,説明変数を増やしていけば「ずれ」は小さくなるので, 定数項を除く説明変数の数を \(\small p\) として, \[\small R^2= 1-\frac{\frac1{n-p-1}\sum_{i=1}^{n}(y_i-\hat{y}_i)^2} {\frac1{n-1}\sum_{i=1}^{n}(y_i-\bar{y})^2}\quad (6)\] で考える場合もあります.これを,自由度調整済みの決定係数といいます. 「11. adc」は,この値を表します(参照). 12〜13の情報量基準も,この回帰式の良さの程度を表す数値です. 他の 2〜8 は,求めた回帰係数の区間推定や「ずれ」の分散に関する情報です. これらの数値を理解するには,推定・検定に関する知識が必要なので, さしあたっては気にしないことです.
 この表示結果から値を取り出すコマンドもあります. それを行うコマンドは「take_inference」です. 「結果のフォーマット」の箇所を参照してください.

▲戻る(トップメニュー)

「stats.mac」のコード
 Maximaの内部では,これらの値をどのようにして計算しているのでしょうか. その計算過程は「stats」パッケージに書かれています. 「load(stats)$」ではなく「load(stats);」とすると, 読み込むパッケージの保存場所が表示されます.最後に書かれているファイル名 「satats.mac」が推定・検定の統計処理を記述しているファイルです. Lisp言語ではなくMaximaのプログラムとして記述されているので, プログラムの知識がある方が見れば, 何をやっているのかを何となく把握することができます. マニュアルでは,線形回帰に関して「linear_regression」しか記載されて いませんが,sats.macをみると, 単回帰を行う「simple_linear_regression」というコマンドもあります. (%o29)の出力は,重回帰に対応した出力です.
 単回帰を行う「simple_linear_regression」は後半の1300行以降に記述されており, その中の「estimation」の箇所に単回帰の場合の計算過程が記述されています. 「dat」を入力したデータ行列,回帰直線を \(\small y=a+bx\) として, 次の手順で計算が進められています.
/* estimations */
means: mean(dat),
covar: cov(dat),
corr: covar[1,2] / sqrt(covar[1,1] * covar[2,2]),
b: covar[1,2] / covar[1,1],
a: means[2] - b * means[1],

/* computing predictions, residuals, residual variance and adc */
pred: transpose(a + b * col(dat,1))[1],
res: transpose(col(dat,2))[1] - pred,
sig2: mean(res^2), /* ml estimator for sigma^2 */
resvar: n * sig2 / (n-2), /* residual variance */
adc: 1 - (1-1/n) * resvar / covar[2,2],
下段は,私の理解による「simple_linear_regression」の コードに関するコメントです.
  1. means: mean(dat)
    /* 各列の平均からなるリストを作成して「means」とする.*/
  2. covar: cov(dat)
    /* データ行列から「cov」により共分散行列を作成して「covar」とする. 2変量の場合は, \(\small \begin{pmatrix} \sigma_x^2& s_{xy}\\ s_{yx}&\sigma_y^2\end{pmatrix}\) が生成され, この行列から2変量の分散と共分散が求められる. */
  3. corr: covar[1,2]/sqrt(covar[1,1] * covar[2,2])
    /* 相関係数を計算して「corr」とする. 相関係数は,\(\small \frac{s_{xy}}{\sigma_x\sigma_y}\) で定義される値です.*/
  4. b:covar[1,2]/covar[1,1]
    /* 回帰直線 \(\small y=a+bx\) の係数 \(\small b\) の値を計算する. それは,\(\small b=\frac{s_{xy}}{\sigma_x^2}\) により計算され, \(\small s_{xy}\) はcovar[1,2]により参照できる.*/
  5. a:means[2]-b*means[1]
    /* 回帰直線 \(\small y=a+bx\) の係数 \(\small a\) の値を計算する. それは,\(\small a=\bar{y}-a\bar{x}\) により計算され, \(\small \bar{x}, \bar{y}\) はそれぞれmeans[1], means[2]により参照できる. */
  6. pred:transpose(a + b * col(dat,1))[1]
    /* 求めた回帰式による予測値からなるリストをpredとする. 予測値からなる1列の行列「a+b\(\small *\)col(dat,1)」を 転置して(transpose)1行からなる行列にして, その1行目をリストとして取り出しています. 行列Mに対してM[k]とすると,k行目がリストで取り出されます.*/
  7. res:transpose(col(dat,2))[1]-pred
    /* 実際の値と予測値との誤差(残差)からなるリストをresとする.*/
  8. sig2:mean(res^2)
    /* 残差の平方(res^2)の平均をsig2とする.*/
  9. resvar:n*sig2/(n-2)
    /* (残差)^2 の合計(残差平方和)を \(\small n-2\) で割った値をresvarとする.
    この値は,線形回帰式 \(\small y=a+bx\) において, 残差の分散の不偏推定量であることが知られています (参照). */
  10. adc:1-(1-1/n)*resvar/covar[2,2]
    /* この定義式を計算すると, \[\begin{align*} &\small {\rm adc}\\ &\small =1-\left(1-\frac1{n}\right)\frac{\rm resvar}{{\rm covar}[2,2]}\\ &\small =1-\frac{n-1}{n}\frac{\frac1{n-2}\sum_{i=1}^{n}(y_i-\hat{y}_i)^2} {\frac1{n}\sum_{i=1}^{n}(y_i-\bar{y})^2}\\ &\small =1-\frac{\frac1{n-2}\sum_{i=1}^{n}(y_i-\hat{y}_i)^2} {\frac1{n-1}\sum_{i=1}^{n}(y_i-\bar{y})^2} \end{align*}\] となり,adcは説明変数が1つだけの場合の自由度調整済み決定係数(6)を 表しています. */
下図は,ここでの1〜10の順番に(%i2)の「dat」を計算したものです.

 上図では,平均,共分散行列,相関係数,そして回帰係数の順に計算されています.

 上図では,予測値,残差,残差平方の平均,残差の不偏推定量, そして自由度調整済み決定係数が計算されています. stats.macでは,この後は回帰係数の信頼区間や「ずれ」の不偏推定量などが 求められています.

▲戻る(トップメニュー)

単回帰による表示結果
 下図は,「simple_linear_regression」による データ行列「dat」の分析結果です. 前述の計算手順により計算されたものです.

 (%o43)をみると,この分析結果は(%o29)よりも分かりやすいといえるでしょう. 解説するまでもないとは思いますが, 回帰直線を \(\small y=a+bx\) とするとき, 以下の内容が示されています.スマホ対応の都合上, 画面の右半分を省略しました。上図では 後半の3つは左半分だけが示されています.
 なお,最新版のwxMaximaでは, (%o43)の半分以下の内容しか表示されません. おそらく,初等的な項目や専門性の高い項目は非表示にしたものと思われます. 表示はされなくても,計算はされています. どのような内容が計算されているのかは 「items_inference」で知ることができます. 「結果のフォーマット」を参照してください.
  1. model: モデル式 \(\small 0.5x+1.1\).
  2. means: 説明変数と目的変数の平均.
  3. variances:説明変数と目的変数の分散.
  4. correlation:相関係数の値.
  5. adc:自由度調整済み決定係数.
  6. a_estimation:回帰係数 \(\small a\) の値.
  7. a_con_int:回帰係数 \(\small a\) の信頼区間.
  8. b_estimation:回帰係数 \(\small b\) の値.
  9. b_con_int:回帰係数 \(\small b\) の信頼区間.
  10. hypotheses:回帰係数 \(\small b\) に関する仮説検定.
    帰無仮説は「H\(\small _0: b=0\)」, 対立仮説は「H\(\small _1: b\ne 0\)」の両側検定.
  11. statistics:検定統計量の値.
  12. distribution:確率分布は自由度3の \(\small t\) 分布.
  13. p_value:棄却域の \(\small p\) 値.
    「statistic」の値が「p_value」を超えているので, 帰無仮説は棄却されることが分かる.
  14. v_estimation:残差の不偏分散.
  15. v_conf_int:残差分散の信頼区間.
  16. cond_mean_conf_int:回帰直線の信頼区間.
  17. new_pred_conf_int:回帰直線による予測区間.
  18. residualas:個々の値に対して, \(\small [(y_i-\hat{y}_i)^2, y_i-\hat{y}_i]\) からなるリスト.
 7〜17の箇所は,推定・検定に関する知識が必要です. 18では,個々の値に対して,残差のみならずその平方の値も示されています. 16〜17は信頼区間と予測区間に関するものだと思われます. 上図では区間の全体が示されていませんが, 16を \(\small [m_1, m_2]\),17を \(\small [p_1, p_2]\) とすると,次の式が示されています. \[\begin{align*} &\small m_1: \\ &\small 0.5x-3.983(0.1(3.0-x)^2+0.2)^{0.5}\\ &\small \hspace{30mm}+1.1\\ &\small m_2:\\ &\small 0.5x+3.983(0.1(3.0-x)^2+0.2)^{0.5}\\ &\small \hspace{30mm}+1.1\\ &\small p_1:\\ &\small 0.5x-3.983(0.1(3.0-x)^2+1.2)^{0.5}\\ &\small \hspace{30mm}+1.1\\ &\small p_2:\\ &\small 0.5x+3.983(0.1(3.0-x)^2+1.2)^{0.5}\\ &\small \hspace{30mm}+1.1 \end{align*}\]  下図は,これらの曲線をデータプロット図に書き入れたものです.

 16,17の計算式はstats.macに記載されているのですが, その式は通常の統計の教科書には記載されていません. これらは,(回帰係数ではなく) 回帰直線自体の信頼区間や予測値の範囲を示しているようです. 1回だけの調査ではなく, 何度も調査してデータを得ることができたとすると, そのたびに異なる回帰直線が得られます. その回帰直線は m1 と m2 の間にあること, そして予測される値は p1 と p2 の間にあることを示しているようです.
 以上のことを解析するには, 誤差の分布の仕方について幾つかの仮定をする必要があります. 詳細は下記を参照してください.

▲戻る(トップメニュー)

結果のフォーマット
 線形回帰の表示結果を利用してさらなる計算を進めるとき, 結果の数値を取り出すには「take_inference」を利用します. たとえば,(%o29)の回帰係数「b_estimation」の値を[b,m]に割り当てるには, 取り出す項目の先頭に「'」をつけて下図の(%i56)のようにします. 同様にして,出力結果を自分の使用する文字に割り当てることができます. また,表示はされなくても,どのような結果が保持されているかをみるには, 結果を(%oN)とすると「items_inference(%oN);」とします(図は略).

 (%o29)の結果の出力様式自体を自分で指定することもできます. それには「inference_result」を利用します. たとえば,結果を次のような形で出力させたいとします.
"Linear Regression 'y=mx+b'"
yosoku_siki=\(\small \cdots\)
zansa_hendou=\(\small \cdots\)
soukan_keisu=\(\small \cdots\)
kettei_keisu=\(\small \cdots\)
 この場合は,次のように指定することになります.

 「inference_result」は3つの引数を持ちます. 第1引数でタイトルを,第2引数で表示したり事後に取り出し可能な値を, そして第3引数で表示させる内容を指定します. 表示させる内容は「'word」の形で指定し, その値はMaximaのコマンドで指定します. そこで指定する値は,stats.mac の「linear_regression」や 「simple_linear_regression」の箇所で利用されている変数です. したがって,このような指定を行うには, stats.mac の内容をある程度読みこなせる必要があります. 第3引数は, 第2引数で指定した内容をどういう順番で表示するかを指定しています. ここでは指定した順番になっていますが, たとえば [1,3] とすると式と相関係数だけが表示されます。 事後に利用する必要がある値は第2引数で多めに指定すれば, 第3引数で結果を表示させなくても「take_inference」で取り出すことができます.
 (%o61)での「yosoku_siki」の係数が同じになってます. satats.mac では,予測式を \(\small y=a+bx\) として計算がなされています. (%i56)は前述の解説で予測式を \(\small y=b+mx\) として \(\small [b,m]\) に \(\small [a,b]\) を割り当てたので, 同じ値が割り当てられてしまいました. 改めてやり直したもので示せば良いのですが, このような指定をするときは Maximaで利用されている変数も考慮しながら使用変数を指定する必要があることへの 注意も兼ねてそのまま残しました(言い訳?).
 「soukan_keisu」の値は,ちょっと紛らわしいかもしれません. 他の2桁の数と見比べると,3 と 2 の間には若干の空白があります. これは,\(\small 5/(3\cdot 2^{\frac32})\) を示しているので注意が必要です. 

 上図では,指定した出力結果から 相関係数を取り出して \(\small r\) に割り当てています. ただし,(%i61)のような指定をしてから「simple_linear_regression」を実行しても, 結果はデフォルトの形式である(%o43)の形で表示されます. 最初から自分の指定形式で表示させるには, 「simple_linear_regression」の実行と「inference_result」を指定して, 下図のように新たなコマンドを自分で作ることになると思われます.

 (%i63)では,上記の2つの処理を行うコマンドを「tankaiki」とし, block関数を利用して定義しています. 複数のコマンドを続けて実行したり,繰り返しや条件分岐を伴うような処理を 行わせる場合にはblock関数を利用します.stats.mac内の 「simple_linear_regression」でもblock関数が利用されています. (%i64)では,定義した「tankaiki」で「dat」を処理させています. 指定した書式で出力されているのが分かります. (%i63)では「float」をつけて小数で出力するようにしていますが, (%o29)や(%o43)を見ても分かるようにデフォルトの設定は回帰式以外は小数で 表示するようになっています.(%o61)が小数で表示されないのは, ちょっとおかしいです. 最新版では修正されており「float」をつける必要はありません. 小数第4位までになっているのは,(%i28)の指定がまだ生きているからです. (%i63)のコマンドは, 「inference_result」の指定に2行を追加しているだけです. 自分で書式を指定したいときは, 最初からこのような形で指定した方がよいでしょう.
 ただし,(%i63)の「inference_result」の第3引数で表示内容を指定しても, 私が使用している「Texmacs+Maxima」では第2引数の内容が全部出力されてしまいます. 最新版のwxMaximaで同じ処理をさせると,第3引数の指定通りに出力されました. バージョンが上がるごとに,このような細かい部分の修正が施されているようです.

▲戻る(トップメニュー)

重回帰分析

最小二乗法
 説明変数を増やして単回帰と同様の1次式での近似を考えるのが 重回帰(多重線形回帰)です.説明変数を \(\small x_1, x_2, \ldots, x_n\), 目的変数を \(\small y\) とするとき,目的変数を最も良く近似する1次式 \[\small y=b_0+b_1x_1+b_2x_2+\cdots+b_nx_n\] を最小二乗法により求めることになります. 簡単のため説明変数が2変数のときを考えて, 変数 \(\small y, x_1, x_2\) のとる具体的な値を, それぞれ \(\small y_i, x_{1i}, \ldots, x_{ni}\) とすると, 単回帰のときと同様に \[\small \sum_{i=1}^{n}\left\{y_i-(b_0+b_1x_{1i}+\cdots+b_nx_{ni})\right\}^2\] が最小になるような係数 \(\small b_0, b_1, \ldots, b_n\) を求めることに なります. 手法は単回帰のときと同様なのですが,変数が増えると式も複雑になります. 説明変数が2つの場合は,次のような計算を行うことになります. \[\small \sum_{i=1}^{n}\left\{y_i-(b_0+b_1x_{1i}+b_2x_{2i})\right\}^2\] を \(\small b_0, b_1, b_2\) の関数とみて \(\small F(b_0, b_1, b_2)\) とします. 最小になる箇所では極値を取ると考えると, \(\small b_0, b_1, b_2\) で偏微分すると 0 になるはずです. 実際に偏微分して \(\small -2\) で割ると次の式になります. なお,和の記号は,単に \(\small \Sigma\) で表しています. \[\begin{cases}\small \sum\left\{y_i-(b_0+b_1x_{1i}+b_2x_{2i})\right\}=0\\ \small \sum{x_{1i}}\left\{y_i-(b_0+b_1x_{1i}+b_2x_{2i})\right\}=0\\ \small \sum{x_{2i}}\left\{y_i-(b_0+b_1x_{1i}+b_2x_{2i})\right\}=0 \end{cases}\] これらを整理して,次の \(\small b_0, b_1, b_2\) に 関する連立方程式(1)が得られます. なお,左辺だけで改行しているのは, スマホでも適切に表示されるようにするためです. \[\left\{\begin{align*} &\small b_0n+b_1\sum{x_{1i}} + b_2\sum{x_{2i}}\\ &\small \hspace{38mm} =\sum{y_i}\\ &\small b_0\sum{x_{1i}}+b_1\sum{x_{1i}}^2 +b_2\sum{x_{1i}x_{2i}}\\ &\small \hspace{38mm} =\sum{x_{1i}}y_i\quad (1)\\ &\small b_0\sum{x_{2i}}+b_1\sum{x_{2i}x_{1i}} +b_2\sum{}x_{2i}^2\\ &\small \hspace{38mm} =\sum{x_{2i}}y_i\\ \end{align*}\right.\]  これにより,変量 \(\small x_1, x_2, y\) に関するデータから, 総和 \(\small \sum{x_{1i}}, \sum{x_{2i}}, \sum{y_i}\), 積和 \(\small \sum{x_{1i}x_{2i}}, \sum{x_{1i}y_i}, \sum{x_{2i}y_i}\) そして平方和 \(\sum{x_{1i}{}^2}, \sum{x_{2i}{}^2}\) を求めておけば, この連立方程式(1)を解くことで 係数 \(\small b_0, b_1, b_2\) の値を求めることができます. 第1式を \(\small n\) で割ると, \(\small \bar{y}=b_0+b_1\bar{x_1}+b_2\bar{x_2}\) であることが 分かります.
 ここで,一般論での議論も視野に入れて, この連立方程式(1)を行列を使って表してみます( 参照). \[\small X= \begin{pmatrix} 1&x_{11}&x_{21}\\ 1&x_{12}&x_{22}\\ \vdots&\vdots&\vdots\\ 1&x_{1n}&x_{2n}\end{pmatrix}\] \[\small \boldsymbol{y}=\begin{pmatrix}y_1\\y_2\\ \vdots\\ y_n\end{pmatrix},\quad \boldsymbol{b}=\begin{pmatrix}b_0\\b_1\\ b_2\end{pmatrix} \] とおくと,連立方程式(1)の右辺は \[\begin{align*} &\small {}^{t} X\boldsymbol{y}\\ &\small =\begin{pmatrix} 1&1&\cdots&1\\ x_{11}&x_{12}& \cdots&x_{1n}\\ x_{21}&x_{22}& \cdots&x_{2n} \end{pmatrix} \!\!\! \begin{pmatrix}y_1\\y_2\\ \vdots\\ y_n\end{pmatrix} \end{align*}\] と表すことができます.一方, (1)の左辺は行列とベクトルの積として \[\small \begin{pmatrix} n&\sum{x_{1i}}&\sum{x_{2i}}\\ \sum{x_{1i}}&\sum{x_{1i}^2}&\sum{x_{1i}x_{2i}}\\ \sum{x_{2i}}&\sum{x_{2i}x_{1i}}&\sum{x_{2i}^2} \end{pmatrix} \!\!\! \begin{pmatrix} b_0\\ b_1\\ b_2\end{pmatrix}\] と表すことができ,左側の行列は \(\small X\) の転置行列を利用すると \[\small \begin{pmatrix} 1&1&\cdots&1\\ x_{11}&x_{12}& \cdots&x_{1n}\\ x_{21}&x_{22}& \cdots&x_{2n} \end{pmatrix} \!\!\! \begin{pmatrix} 1&x_{11}&x_{21}\\ 1&x_{12}&x_{22}\\ \vdots&\vdots&\vdots\\ 1&x_{1n}&x_{2n} \end{pmatrix} \] と表すことができるので, 要するに(1)の左辺は次のように簡潔に表すことができます. \[\small {}^t X X \boldsymbol{b}\] したがって,連立方程式(1)は,行列を用いると 次のように表わされます. \[\small {}^t X X \boldsymbol{b}={}^{t} X\boldsymbol{y} \qquad (2)\] ここで,行列 \(\small X\) は説明変数のデータ行列に 1 だけからなる 列を追加した行列,列ベクトル \(\small \boldsymbol{y}\) は目的変数のデータ列,そして列ベクトル \(\small \boldsymbol{b}\) が求めるべき回帰係数の列です. 目的変数が増えると,\(\small X\) の列数と \(\small \boldsymbol{b}\) の 行数が増えるだけで, 式の形は(2)と全く同じです.
 式(2)において \(\small {}^t X X\) の逆行列が存在すれば, 回帰係数 \(\small \boldsymbol{b}\) は \[\small \boldsymbol{b}=({}^t X X)^{-1}\, {}^t X \boldsymbol{y} \quad (3)\] により求めることができます.したがって, 与えられたデータから \(\small X, \boldsymbol{y}\) を 作成して,(3)の右辺の計算をさせれば回帰係数が求められます.

▲戻る(トップメニュー)

基本統計量からの計算
 前述の式(3)を用いて,実際に計算してみましょう. 手計算で確かめやすいように,適当に作った次の例で考えます.

x1 5 4 3 1 2
x2 -1 5 -2 4 -3
y 12 -32 15 -21 36

 下図には掲載していませんが, 最初に「stats」パッケージを読み込んでいます. このパッケージを読み込むと,「descriptive」や「distrib」が 読み込まれていなくても,この2つのパッケージは自動で読み込まれます.

 (%i3)〜(%i5)では,サンプルデータをリストで入力しています. (%i6)では,作成したリストをまとめて行列にして, その転置行列を求めて 「linear_regression」に入力する行列「dat」を作っています. 実際の例の場合は, あらかじめcsvファイルを作成して, 「load」で読み込ませることになります.

 (%o6)のデータ行列「dat」をもとに, 以下ではstats.macの「linear_regression」の箇所に 書かれている計算過程に沿って計算します.
 (%i7)では入力されたデータ数を, (%i8)では「dat」の1行目の長さ (つまり,「dat」の列数)から 1 を引いて, 説明変数の個数を求めています.
 式(2)の行列 \(\small X\) を作るには, 入力されたデータ行列「dat」から最後尾の列である 目的変数のデータ列を除いて,さらに第1列に 1 だけからなる列を 追加する必要があります. 列を取り出すには「col」,追加するには「addcol」を利用します.

 (%i10)では,行列「dat」の p+1 列目を取り出して y に割り当てています. 「col」は列(column)のことです.p は(%i8)で求めた説明変数の個数です. (%i12)では,「dat」から p+1 列目を除いた行列を取り出しています. 「submatrix」の書式では, 「submatrix(i1,i2,M,j1,j2)」などとすると, 行列 M から i1行,i2行,そして j1列,j2列を除いた行列が返されます. i1, i2を省略すると列だけが除かれ,j1, j2を省略すると行だけが 除かれます.

 (%i13)では,1 だけからなる列を作っています. 最初に「makelist([1],i,n)」により, [1] だけからなるリストを作ります. 通常は「makelist(expr, i, 1, n, s)」などとすることで, 繰り返し変数「i」を 1 から n まで s 刻みに 変えながら「expr」からなるリストが作成されます. 初期値やステップ数を省略すると,それらは 1 として設定されます. 1 だけから成るのであれば,[1] ではなく単なる 1 でもよいように 思えるのですが,[1] としておかないと 「apply」コマンドの箇所でエラーが生じます. このコマンドは, 第2引数に対して第1引数で指定する処理を行わせるものです. (%i13)の書式により,第2引数で作成されたリストが, 1列だけからなる行列(matrix)として出力されます. (%i14)では,1 だけからなる列を行列「X」に追加しています.  この行列「X」に対して式(3)の右辺を計算すれば \(\small \boldsymbol{b}\) が求められます.

▲戻る(トップメニュー)

 ここでは,いきなり求めるのではなく, 入力したデータ行列「dat」から連立方程式(1)の係数となっている 各種の統計量を求めて, この連立方程式を地道に解いて求めてみましょう. それには,各データに対して 単純和,平方和,積和を求める必要があります.
 なお,以下の(%i27)〜(%i48)は,連立方程式(1)に対する理解や, Maximaのいろいろなコマンドの使用法の解説を兼ねた内容です. (%i14)の後に「stats.mac」で行われている計算は, (%i49)以降で解説します.

 (%i27)では,dat の平均に n を掛けることで \(\small x_1, x_2, y\) の 合計を求めています.(%i28)では,行列の2乗は各成分が2乗されることを 確認した上で,(%i29)では平均をもとに \(\small x_1, x_2, y\) の平方和を求めています. (%i30)(%i31)では,最初に定義したリストの積の平均を求めることで, \(\small x_1y, x_2y\) の積和を求めています. リストとリストの積では,対応する要素の積が計算されます.
 ここでは個別に計算しましたが,個々の列の平方和や積和の値は, 共分散行列を用いて求めることもできます. 「cov(dat)」により,共分散行列 \[\begin{pmatrix}\small \sigma_1^2& s_{12}& s_{1y}\\ s_{21} & \sigma_2^2 & s_{2y}\\ s_{y1} & s_{y2} & \sigma_y^2 \end{pmatrix}\] が得られます.たとえば,\(\small s_{2y}\) は, \[\begin{align*} \small s_{2y}&\small =\frac1{n}\sum_{i=1}^{n}(x_{2i}-\bar{x_2}) (y_i-\bar{y})\\ &\small =\frac1{n}\sum_{i=1}^{n}x_{2i}y_{i} -\bar{x_2}\bar{y} \end{align*}\] であるので,\(\small s_{2y}\) に2変量の平均の積 \(\small \bar{x_2}\bar{y}\) を加えて \(\small n\) を掛ければ 積和 \(\small \sum_{i=1}^{n}{x_{2i}y_i}\) が求められます.

 上図では,(%i34)により共分散行列が求めて「xycov」とします. (%i35)では,(2,3)成分に平均の積を加えています. (%i27)で「S」には総計を割り当てているので, 平均を求めるため \(\small n=5\) で割っています. (%o35)に 5 を掛けると,(%o31)と一致することが分かります.
 いずれにしろ,以上により連立方程式(1)で必要とされる係数が求まりました. 入力データ「dat」から作成される連立方程式は下記のようになります. \[\small\left\{\begin{align*} \small 5b_0+15b_1+3b_2&=10\\ \small 15b_0+55b_1+7b_2&=28\\ \small 3b_0+7b_1+55b_2&=-394 \end{align*}\right.\]  方程式を解くコマンド「solve」を利用して解くと, 下記のようになります.


▲戻る(トップメニュー)

行列を利用した計算
 今度は,逆行列を利用して求めてみましょう. この連立方程式を \(\small A\boldsymbol{b}=\boldsymbol{p}\) とすると, 解は \(\small \boldsymbol{b}=A^{-1}\boldsymbol{p}\) により求められます. Maximaは,連立方程式と変数のリストを指定すると, 係数行列や拡大係数行列を取り出すことができます.

 (%i39)では,拡大係数行列を取り出すコマンド「augcoefmatrix」を使用しています. ただし,定数項を移項して, 方程式を「\(\small \cdots=0\)」の形にしたときの係数が取り出されます. 「coefmatrix」を利用すると,係数行列だけが取り出されます. (%i40)では得られた拡大係数行列Abから「submatrix」を利用して4列目を除いた行列を, (%i41)では1〜3列目を除いたものを求めています. そして,(%i42)では「invert」により逆行列を求め, それを \(\small -\boldsymbol{p}\) に掛けることで 解を求めています. 当然ながら,(%o37)と同じ結果が得られます.

 最後に,行列 X を利用して,式(3) \[\small \boldsymbol{b}=({}^t X X)^{-1}\, {}^t X \boldsymbol{y}\] を直接計算して求めてみましょう. これは,実際に「stats.mac」の中で行われているもので, (%i14)からの続きの計算になります.

 (%i45)では行列 X の転置行列を XT とし,X との積 \(\small {}^t X X\) を求めています. 行列の積は \(\small *\) ではなく「.」を用います. (%i47)では \(\small {}^t X X\) の逆行列を求め, (%i48)で \(\small {}^t X\boldsymbol{y}\) との積を計算して \(\small \boldsymbol{b}\) を求めています. 行列とリストとの積「.」は, リストを行列に直さなくても1列だけからなる行列として計算します. (%o37)や(%o44)の結果と一致しています.
 ここでは,できるだけ個別に計算しましたが, stats.mac内では幾つかの計算がまとめて行われています. いずれにしろ,理論的な計算式に具体例をあてはめながら求めてみることも, 理解を深める上では必要と思われます. 「数式処理」を活用すると,数式を交えた文字計算も含めながら 確認していくことができます.

▲戻る(トップメニュー)

コマンドの表示結果
 多重線形回帰を行うコマンドは「linear_regression」です. 説明変数と目的変数の列データをまとめた行列を「dat」として 入力すると,分析得結果が表示されます. サンプルデータ「dat」を(%o6)として実行すると 下記が表示されます.

 この出力結果に関して,マニュアルでは下記のことが書かれています. 最新版のwxMaximaでは, 下記の項目のうち 1,4,5,6,7,8,9,11 の項目しか表示しません.
  1. b_estimation: 回帰係数推定
  2. b_covariances: 回帰係数推定の共分散行列
  3. b_conf_int:回帰係数の信頼区間
  4. b_statistics:係数テストの統計
  5. b_p_values:係数テストのp値
  6. b_distribution:係数テストの確率分布
  7. v_estimation:不偏分散推定量
  8. v_conf_int:分散信頼区間
  9. v_distribution:分散テストの確率分布
  10. residuals:残差
  11. adc:調整済み決定係数
  12. aic:赤池情報量基準
  13. bic:Bayes情報量基準
 とりあえず,回帰係数は「1. b_estimation」の箇所にあります. 当然ですが,これまでに求めたものと同じ値が小数で示されています. 実際の値との誤差 \(\small y_i-\hat{y}_i\) は「10. residuals」 の箇所です.調整済み決定係数は「11.adc」です. かなり近似の度合いが高いようです.
 それ以外の項目は, 求めた回帰係数や残差の推定に関する内容です. それをどのようにして計算しているのかが気になる方は, stats.macの後半にある「linear_regression」の箇所をみてください. ただし, その内容を理解するには確率分布や推定に関する知識が必要になります.
 wxMaximaでは,(%o52)のような多数の結果は表示されませんが, 同じ内容は保持しています. 下図は,wxMaximaの画面です.

 どのような項目が計算ずみであるかは,結果が表示された(%o9)に対して 「items_inference」を実行すると表示させることができます. (%o9)では画面の左半分しか表示されていませんが, マニュアルに書かれている全13項目の項目名が表示されます. 個々の項目の値は, 「take_inference」 を通して取り出すことができます. wxMaximaでの表示では取り出す項目名の書式を把握しにくいですが, マニュアルの該当箇所を参照しながら入力するとよいでしょう.

 この結果から予測値を計算するには次のようにします. 以下は,「TeXmacs+Maxima」の画面に戻ります.

 (%i54)では, 取り出した回帰係数を \(\small \boldsymbol{b}\) に割り当てて, (%i55)では \(\small X\boldsymbol{b}\) を計算しています. 「b」はリストですが,行列の積を計算するコマンド「.」は, \[\small X\boldsymbol{b} =\begin{pmatrix}1&5&-1\\ 1&4&5\\ 1&3&-2\\ 1&1&4\\ 1&2&-3\end{pmatrix} \begin{pmatrix}11.7\\ -1.717\\ -7.583\end{pmatrix}\] として計算します.  (%i56)は, (%i10)で定義済みの \(\small \boldsymbol{y}\) を表示させて, (%i58)で両者の差から残差を求めています. (%o52)と同じ結果であることが分かります.

▲戻る(トップメニュー)

いろいろな回帰

 変量間の関係を1次式で説明できないときは, 背景にある状況やデータ分布の様子などから予想される関係式を 自分で見いだす必要があります. そのようなとき,パッケージ「lsquares」に登録されている コマンド「lsquares_estimates」を利用すると, 自分で設定した式で回帰分析をすることができます.


lsquares_estimates
 最初に,このコマンドの利用例を見てみましょう. 具体例として,2010年式乗用車の 平均旅行速度(km/h)と燃料消費率(L/km)との関係をみてみます( 参照:表12).

 ここでは,入力番号をリセットしています. (%i1)で,最初に「lsquares」を読み込みます. 入力するデータ行列「dat」を(%i4)により作成します. そのグラフを(%i5)によりプロットすると,下図のようになります.

 グラフをみると両者の関係は2次関数のように見えるので, 回帰式を \(\small y=ax^2+bx+c\) に設定して 「lsquares_estimates」で分析してみましょう. このコマンドの書式は次のようになります.

lsquares_estimates(データ行列, 列の名前のリスト, 回帰式, パラメータのリスト)

 第1引数は分析するデータ列からなる行列です. 第2引数は,そのデータ列にリストで1列目から順に名前をつけます. 第3引数では,分析したい回帰式を第2引数の変数を使って記述します. したがって,必ずしも最後の列が目的変数である必要はありません. 第4引数は,第3引数で設定した式に使われているパラメータのリストです.
「lsquares_estimates」は,設定された回帰式に対して, 最小二乗法により誤差の平方和が最小となるようなパラメータの値を求めます. パラメータに関して得られる連立方程式から厳密解が求められるときは 厳密解を表示します. 厳密解を求めることができないときは, 近似解を求めるパッケージ「lbfgs」に必要な情報を引き渡して 近似計算を行います.回帰式は自分で自由に設定することができます. したがって,回帰係数を求めるだけであれば単回帰や重回帰も 「lsquares_estimates」で行うことができますが, 求めた回帰係数の信頼区間や決定係数のような内容には対応していません.

 (%i8)では,表示する小数の桁数を4桁に制限しています. (%i8)では,必要とされる4つの引数を指定しています. (%o9)は厳密解で,パラメータに関する連立方程式が解けたことを意味します. (%i10)では,得られた解を小数に直しています.

 連立方程式の解が複数個あるときは [[\(\small \cdots\)],[\(\small \cdots\)]] のような形で出力されます.(%i9)の解は1組だけであることから, 直前の式(%o9)にあるようにリストが重複して表示されます.そこで, (%i11)では内側のリストを取り出しています.  (%i12)では, 直前の結果(%)の個々の要素の右辺を取り出して \(\small a,b,c\) に値を割り当てています.しかし, 個別にいちいち「rhs」(等式の右辺:right-hand side)を指定するのは面倒です. そのようなときは,下図にあるような「map」関数を利用すると便利です.

 (%i16)では,(%o11)で出力されたリストの個々の要素の 右辺(rhs)を取り出すことを指定しています. (%o12)と同じ結果が得られます.

 (%i17)では,割り当てられた \(\small a, b, c\) を利用して 回帰式 \(\small f(x)\) を定義して,(%i18)でそのグラフを描画しています. わりと良く近似していることが分かります.

▲戻る(トップメニュー)

回帰式の残差
 「lsquares_estimates」は,回帰式を自由に設定して 与えられたデータを最も良く近似するパラメータを決定します. このコマンドの利用例として, マニュアルでは回帰式を \(\small (z+D)^2=Ax+By+C\) とする場合が紹介されています. しかし,単回帰や重回帰のような多岐にわたる情報は得られません. その後は,得られた係数の値をもとに自分で分析する必要があります. ただし,残差(誤差)を回帰式による左辺と右辺の差とすると, 残差と残差平方の平均を求めるコマンドがあります. 残差は「lsquares_residuals」,残差平方の平均は 「lsquares_residual_mse」により求めることができます. 前者は「residuals」,後者は「residual」なので気をつけてください.

 (%i19)は,求めた係数をもとに個々のデータの残差を求めるコマンドです. 第4引数は,(%o12)ではなく(%o11)のタイプのリストで引き渡します. 設定した回帰式の「(左辺)\(\small -\)(右辺)」の値がリストで表示されます. (%o9)[1]として引き渡すと結果は分数で表示されます. マニュアルでは,(%o9)を「reg」として割り当てて「first(reg)」と して取り出しています. 「first」は,リストの「最初の要素」を取り出すコマンドです. なお,(%o19)の右半分は省略しています.

 ここでは,実際に「(左辺)\(\small -\)(右辺)」を求めて確認しています. (%o19)と同じ値であることが分かります.

 (%i25)は,残差の平方の平均を求めるコマンドです. つまり,設定した式に対して \[\small \frac1{n}\sum_{i=1}^{n}\left\{(左辺)-(右辺)\right\}{}^2\] の値が返されます. (%i20)(%i21)で \(\small x, y\) にデータ列を割り当てたので, そのまま(%i25)を実行するとエラーになります.そこで, (%i24)では割り当てた変数の値を「kill」を利用して解除しています.
 以下では,(%o25)の値が(%o22)の二乗の平均(mean)と一致することを 確認しようとしています. 平均を求めるコマンド「mean」を利用するため,最初に パッケージ「descriptive」を読み込みます.

 (%i28)により(%o22)の二乗の平均を求めると, (%o25)と一致することが分かります. 「mean」は,行列が引数のときは各列の平均をリストで返します. (%i28)は引数が1列だけから成る行列なので, 結果は要素が1つだけのリストで返されます.

▲戻る(トップメニュー)

残差平方の平均
 「lsquares_estimates」は,設定された回帰式から「残差平方の平均」 が最小になるパラメータを決定しようとします. それを決定するための方程式は,単回帰や重回帰のときと同様に, 指定したパラメータで微分した式を 0 とする連立方程式です. コマンド「lsquares_mse」を利用すると, その連立方程式のもとになる「残差平方の平均」の式を出力させることができます.

 (%i30)で,いったん割り当てた \(\small a, b, c\) の値を解除します. (%i31)により,回帰方程式を \(\small y=ax^2+bx+c\) とするとき, 残差平方の平均の式が出力されます. 残差平方の平均は \[\small \frac1{n}\sum_{i=1}^{n}\left\{y_i-(ax_i^2+bx_i+c)\right\}^2\] であり,\(\small x_i\) は dat[i,1],\(\small y_i\) は dat[i,2] により 参照することができるので,その式は(%o31)のように出力されます. (%i30)で \(\small a,b,c\) への値の割り当てを削除しておかないと, (%o31)の式は割り当てた具体的な値を用いて表されます. (%o31)にデータ行列「dat」の値を代入するには, (%o31)に対して次のようにします.

 (%i32)のように,「''」と「nouns」をつけることで, (%o31)の式に「dat」の値を代入した式が出力されます. 画面の右半分は省略しています. なお,下図のように, 「'%o31」だけだと「%o31」が文字出力され, 「''%o31」とすると「%o31」と同じ式が出力されます. 「''%o31,nouns」とすることで「dat」の値が代入されます.


▲戻る(トップメニュー)

2次回帰の計算
 前項により2次回帰したときの残差平方の平均の式が出力できるので, それをパラメータで微分した式が 0 として得られる連立方程式を, 実際に解いて解が(%o9)となることを確認してみましょう. 微分のコマンド「diff」と, 方程式の解を求めるコマンド「solve」を利用します.

 出力された式から関数を定義するときは「define」を用います. (%i41)では,具体的なデータをあてはめた(%o32)の式を \(\small F(a,b,c)\) として定義して, それを \(\small a,b,c\) で微分した式を, それぞれ eqa, eqb, eqc に割り当てています. そして,(%i45)で,その式を \(\small a,b,c\) について解いています. 式が等式で与えられないとき,「solve」は 「\(\small \cdots=0\)」の式として解を求めます. (%o45)の解は(%o9)と一致します. なお,実際には, (%o45)の前にはもっと多数の 「rat: replaced \(\small \cdots\)」が出力されます. 「数式処理」を行う都合上, Maximaの内部では小数で与えられたデータを有理数に直して計算しています. どの小数をどのような有理数に直したかが表示されています.

 ここで,2次回帰 \(\small y=ax^2+bx+c\) のときの (平均ではなく)残差平方和を \(\small F(a,b,c)\) とすると,それは \[\small \sum_{i=1}^{n}\left\{y_i-(ax_i^2+bx_i+c)\right\}^2\] という式です.重回帰分析のときと同様にすると, \(\small a,b,c\) で微分して 0 であるとして得られる連立方程式は \[\left\{\!\! \begin{align*}\small \sum\left\{y_i-(ax_i^2+bx_i+c)\right\}x_i^2=0\\ \small \sum\left\{y_i-(ax_i^2+bx_i+c)\right\}x_i=0\\ \small \sum\left\{y_i-(ax_i^2+bx_i+c)\right\}\cdot 1=0 \end{align*}\right.\] であり,括弧をはずして整理すると次のようになります. \[\begin{align*}\small a\sum{x_i^4}+b\sum{x_i^3}+c\sum{x_i^2}&\small =\sum{x_i^2y_i}\\ \small a\sum{x_i^3}+b\sum{x_i^2}+c\sum{x_i}&\small =\sum{x_iy_i}\\ \small a\sum{x_i^2}+b\sum{x_i}+cn&\small =\sum{y_i} \end{align*}\]  今度は,この係数を求めることで解いてみましょう.

 使用するデータ「dat」には小数が含まれるので, (%i47)で小数の桁数をデフォルトの値に戻しておきます. (%i48)(%i49)では個々の列データを取り出して datx, daty とし, データ数を \(\small n\) としています. (%i51)(%i52)では datx, daty の総和を求めています. (%o52)は 0.491 のことです.

 (%i53)〜(%i57)では,べき乗和や積和の値を求めています.

 求めた値をもとに,個々の方程式を eq1, eq2, eq3 として (%i61)により解を求めると,(%o45)で求めたものと一致します. ここでは連立方程式を自分で入力しましたが, 入力間違いを防ぐには(%i50)〜(%i57)を適当な文字に割り当てて, その文字を利用してeq1,eq2,eq3を記述した方が良いでしょう.

 ここで,eq1, eq2, eq3 の式を,(%i42)〜(%i44)で定めた eqa, eqb, eqc と比較してみましょう. 係数の簡単な eq3 と eqc を比較すると, eq3 を 5 で割ったものが eqc になっています. この係数の違いは, eqa〜eqb は残差平方を n で割って微分したものなので, \(\small \frac1{10}\cdot 2=\frac15\) が付くことによります.

 いろいろな統計ソフトでは, コマンドにデータを入力すれば即座に結果が得られます. 理解を深めるためには,その結果を鵜呑みにする前に, 一度は理論式をもとに自分で求めてみることも必要と思います. その際, Maximaを「紙と鉛筆」代わりに用いて, 煩雑な基本統計量の計算や連立方程式の解法は 数式処理機能を利用すれば, あたかも自分で計算しているかのような感覚で結果を得ることができます.


▲戻る(トップメニュー)

指数回帰
 今度は,変量の間の関係を指数関数 \(\small y=ab^x\) で近似する場合を 考えてみましょう. これは,増加率や減少率が一定であるような変化です. 具体例としてタンチョウの個体数変化を取り上げます. 日本のタンチョウは渡り鳥ではなく, 北海道の一部の地域に生息する留鳥です. 一時は絶滅寸前にまで減少しましたが, 特別天然記念物に指定されて保護活動が本格化してから増加に転じています. 自然界の生物として,その個体数はかなり正確に把握されています. 下図では環境省のデータを利用します (参照).

 入力番号をリセットして,「descriptive」と「lsquares」の パッケージを改めて読み込んでおきます. (%i5)は昭和の年号,(%i6)は各年度に観測された個体数, そして,(%i7)で「lsquares_estimates」に入力するための行列データに しています.

 (%i8)によりグラフをプロットすると, グラフの形状は指数関数的であることが分かります. そこで,回帰式を \(\small y=ab^x\) に設定すると, 下図のような結果が得られます.

 (%i17)で(表示の見やすさのために)小数の桁数を4桁に制限し, (%i18)で回帰式を設定して係数を求めています. データ行列は整数だけなのに結果が小数で表示されたことは, 厳密解ではなく近似解が得られたことを意味します. なお,回帰式を \(\small y=ae^{bx}\) で設定するとエラーになりました.

 (%i19)はmap関数を利用しています. (%o18)の結果から「%[1]」により内部のリストを取り出し, 得られたリスト [\(\small a=0.1436, b=1.225]\) の右辺(rhs)の値を表示させています. (%i20)では,その値をもとに回帰式 \(\small f(x)\) を定義して, (%i22)でグラフを描画しています. かなり良く当てはまっています.

 タンチョウは,このまま指数関数的に増え続けていくわけではありません. この後は成長曲線(ロジスティック曲線)に当てはまるような変化をします. その変化は差分方程式で表されて教材化されており, 中学生や高校生にも教授可能です (参照).

 (%i23)では残差平方の平均の式を求めて, (%i24)では具体的なデータを当てはめています. \(\small a, b\) の値を求めるには, (%o24)(右側は省略)を \(\small a, b\) で微分すると 0 という連立方程式を 解くことになりますが,それは非線形連立方程式であり, 厳密解を求めることは困難です. そのようなとき,Maximaは自動的に近似解を求める手順に移ります. それを行うパッケージが「lbfgs」で, 必要な引数が自動的に引き渡されて計算が行われます.
 非線形連立方程式の近似解を求める方法に, ブロイデン・フレッチャー・ゴールドファーブ・シャノン法(BFGS法)という 方法があり,それを限られたメモリー(Limited memory)でも 実行できるように変更されたものが 「L-BFGS法」であるようです(参照1参照2).

▲戻る(トップメニュー)

対数変換
 いろいろな回帰分析を行う上で, 対数を取った値で分析する場合があります. 対数を取ることを「対数変換」といいます. なぜ対数を取るのか,試みにタンチョウの例で試してみましょう.

 (%oi25)では,kotaiの(自然)対数を取って小数に直しています. (%i26)では,グラフにプロットしています.前項のグラフと比較すると 直線的になっていることが分かります.
 対数を取ると,\(\small 10^n\) のような大きな値や \(\small 10^{-n}\) のような小さな値が取り扱いやすい値に変換されます. 対数関数の変化は極めてゆっくりです.たとえば, 常用対数 \(\small \log_{10}{x}\) の場合は, \(\small x\) の値が 1 から 100 まで変化しても 対数の値は 2 しか変化しません. 説明変数の変化が大きいときに対数を取ると, 変化の小さい変数に変換して考えることができるのです. タンチョウの個体数では, 33羽から172羽までの変化が, 自然対数を取ると 3.497 から 5.147 までの変化に変換されます.

 対数を取ると直線的な変化になることが多いので, 今度は単回帰を行ってみましょう. 下図では,(%i27)により 「simple_linear_regression」に引き渡す行列「dat2」を作成し, (%i28)では単回帰を行うコマンドが登録されているパッケージ「stats」 を読み込んでいます. そして,(%i29)により単回帰分析を行っています.

 Maximaのバージョン「5.38.1」ではすべての項目が表示されますが, 最新バージョン(5.42.2)が組み込まれているwxMaximaでは, このうちの一部の結果だけが表示されます. データ行列を引き渡すとき, (%i25)にあるような小数にしないとエラーになります. いずれ,回帰直線は \(\small 0.206x-2.046\) であることが分かります. ここに表示されている項目の詳細は 単回帰分析の項目を参照してください.

 (%i30)では直前の結果(%)からモデル式を取り出して, (%i31)で \(\small g(x)\) として定義しています. そして,(%i33)により直線を重ね描きしています

▲戻る(トップメニュー)


copyright