ブログ

統計検定準1級 問題解説 ~2017年6月実施 選択問題及び部分記述問題 問13~

過去問題

過去問題は統計検定公式が問題と解答例を公開しています。こちらを参照してください。

解答

[1]解答

$\boxed{ \ \mathsf{20}\ }$:①

結合されている終端ノードの距離は(ア),(イ)のいずれでも同じである.$A,B$の距離は$2$であり,表2から距離$2$なのは生徒$1$と生徒$2$の間だけなので,$A,B$が生徒$1$,$2$だとわかる.次に$D,E$の距離をみると約$1.5$であることが表2からわかる.生徒$3~5$の間で距離が$1.5$に近いのは生徒$4$と生徒$5$の$\sqrt{2}$であるからこれもわかる.最後に残ったのが$C$が生徒$3$である.

(ア)において,クラスター$\{A,B\}$と$C$がクラスターを形成しており,その距離は$2$より少し大きい.ここで,$A$と$C$の距離(つまり,生徒$1$と生徒$3$)及び,$B$と$C$(つまり生徒$2$と生徒$3$)の距離はそれぞれ,$\sqrt{17}$,$\sqrt{5}$となっている.$\sqrt{5} \approx 2.2$であるから(ア)が最近隣法であるとわかる.したがって,(イ)は最遠隣法となる.

[2]解答

(1) $\boxed{ \ \mathsf{21}\ }$:①

クラスター中心とそれに対応するクラスターは次のように変化する.

クラスター中心クラスター
初期点$(10,10)$,$(10,8)$$\{ 1 \}$,$\{2,3,4,5\}$
1回目の変化$(10,10)$,$(8,7.25)$$\{1,2\}$, $\{3,4,5\}$
2回目の変化$(10,9)$,$(7.3,7)$$\{1,2\}$, $\{3,4,5\}$

まず,初期点に対するクラスターを考える.下図から明らかに生徒$3$,$4$,$5$は生徒$1$よりも生徒$2$に近い.よって初期点に対して割当てられたクラスターは$\{1\}$,$\{2,3,4,5\}$である.

次に,このクラスターから新しいクラスター中心を求める(1回目の変化).クラスター$\{1\}$の中心は変化せず$(10,10)$である(赤い点).クラスター$\{2,3,4,5\}$に対応する中心は,

\[
\left( \dfrac{10+9+7+6}{4}, \dfrac{8+6+8+7}{4} \right)= (8,7.25)
\]

である(青い点).すると,生徒$1$,$2$が赤い点を中心とするクラスターに属し,生徒$3$,$4$,$5$が青い点を中心とするクラスターに属することがわかる.

したがって,出力されるクラスターは$\{1,2\}$,$\{3,4,5\}$である.

(2回目の変化)クラスター$\{1,2\}$,$\{3,4,5\}$に対応するクラスター中心はそれぞれ$(10,9)$,$(7.3,7)$である.これから出力されるクラスターは$\{1,2\}$,$\{3,4,5\}$であり変化しない.

以上から,出力されるクラスターは$\{1,2\}$,$\{3,4,5\}$である.

(2) $\boxed{ \ \mathsf{22}\ }$:④

下の図から初期点が$(a)$のときはクラスターは$\{1,2,3\}$,$\{4,5\}$となり,初期点が$(b)$のときは$\{1,4,5\}$,$\{2,3\}$となる.

Ch.12 「ベイズ推論」の章末問題の解答例 〜数理統計学(共立出版)〜

当記事は「数理統計学(共立出版)」の読解サポートにあたってChapter.$12$の「ベイズ推論」の章末問題の解答の作成を行いました。
基本的には書籍の購入者向けの解説なので、まだ入手されていない方は購入の上ご確認ください。また、解説はあくまでサイト運営者が独自に作成したものであり、書籍の公式ページではないことにご注意ください。

・解答まとめ
https://www.hello-statisticians.com/answer_textbook_math_stat#green

章末の演習問題について

問題12.1の解答例

$X \sim \mathrm{Bin}(n,p)$より、$p$の共役事前分布にベータ分布$\mathrm{Be}(\alpha,\beta)$を考えることができる。$[1]$、$[3]$もベータ分布の特殊な例と見ることができるので、実質的には$[2]$で$[1]$、$[3]$も内包することができる。

よって確率変数$Y$が$Y \sim \mathrm{Be}(\alpha,\beta)$であるときの$E[Y]$を先に導出し、その結果をそれぞれの設問で適用する。$E[Y]$は下記のように計算できる。
$$
\large
\begin{align}
E[Y] &= \int_{0}^{1} y \times \frac{1}{B(\alpha,\beta)} y^{\alpha-1} (1-y)^{\beta^1} dy \\
&= \frac{1}{B(\alpha,\beta)} \int_{0}^{1} y^{(\alpha+1)-1} (1-y)^{\beta^1} dy \\
&= \frac{B(\alpha+1,\beta)}{B(\alpha,\beta)} \\
&= \frac{\Gamma(\alpha+1)\cancel{\Gamma(\beta)}}{\Gamma(\alpha+\beta+1)} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\cancel{\Gamma(\beta)}} \\
&= \frac{\cancel{\Gamma(\alpha)}\alpha}{\cancel{\Gamma(\alpha+\beta)}(\alpha+\beta)} \frac{\cancel{\Gamma(\alpha+\beta)}}{\cancel{\Gamma(\alpha)}} \\
&= \frac{\alpha}{\alpha+\beta} \quad (1)
\end{align}
$$

$[1]$
ベータ分布$\mathrm{Be}(1,1)$は一様分布であるので、無情報事前分布と考えることができる。このとき$p$の事後分布を$P(p|x)$とおくと$P(p|x)$は下記のように考えることができる。
$$
\large
\begin{align}
P(p|x) & \propto P(x|p) P(p) \\
&= {}_{n} C_{x} p^{x}(1-p)^{n-x} \times \frac{1}{B(1,1)} p^{1-1}(1-p)^{1-1} \\
& \propto p^{(x+1)-1}(1-p)^{(n-x+1)-1}
\end{align}
$$

よって事後分布は$\mathrm{Be}(x+1,n-x+1)$であり、$p$のEAP推定量を$E[p|X]$とおくと$(1)$式より下記のように考えることができる。
$$
\large
\begin{align}
E[p|X] &= \frac{x+1}{(x+1)+(n-x+1)} \\
&= \frac{x+1}{n+2}
\end{align}
$$

$[2]$
$p$の事前分布にベータ分布$\mathrm{Be}(\alpha,\beta)$を考えるとき、事後分布$P(p|x)$は下記のように考えることができる。
$$
\large
\begin{align}
P(p|x) & \propto P(x|p) P(p) \\
&= {}_{n} C_{x} p^{x}(1-p)^{n-x} \times \frac{1}{B(\alpha,\beta)} p^{\alpha-1}(1-p)^{\beta-1} \\
& \propto p^{(x+\alpha)-1}(1-p)^{(n-x+\beta)-1}
\end{align}
$$

よって事後分布は$\mathrm{Be}(x+\alpha,n-x+\beta)$であり、$p$のEAP推定量を$E[p|X]$とおくと$(1)$式より下記のように考えることができる。
$$
\large
\begin{align}
E[p|X] &= \frac{x+\alpha}{(x+\alpha)+(n-x+\beta)} \\
&= \frac{x+\alpha}{n+\alpha+\beta}
\end{align}
$$

$[3]$
フィッシャー情報量$\displaystyle \sqrt{I(p)} = \sqrt{\frac{n}{p(1-p)}}$に比例するベータ分布$\displaystyle \mathrm{Be} \left( \frac{1}{2},\frac{1}{2} \right)$はジェフリーズの事前分布と考えられる。このとき$p$の事後分布を$P(p|x)$とおくと$P(p|x)$は下記のように考えることができる。
$$
\large
\begin{align}
P(p|x) & \propto P(x|p) P(p) \\
&= {}_{n} C_{x} p^{x}(1-p)^{n-x} \times \frac{1}{B \left( \frac{1}{2},\frac{1}{2} \right)} p^{\frac{1}{2}-1}(1-p)^{\frac{1}{2}-1} \\
& \propto p^{x+\frac{1}{2}-1}(1-p)^{n-x+\frac{1}{2}-1}
\end{align}
$$

よって事後分布は$\displaystyle \mathrm{Be} \left( x+\frac{1}{2},n-x+\frac{1}{2} \right)$であり、$p$のEAP推定量を$E[p|X]$とおくと$(1)$式より下記のように考えることができる。
$$
\large
\begin{align}
E[p|X] &= \frac{x+1/2}{(x+1/2)+(n-x+1/2)} \\
&= \frac{x+1/2}{n+1} \\
&= \frac{2x+1}{2(n+1)}
\end{align}
$$

問題12.2の解答例

まとめ

ジェフリーズの事前分布(Jeffreys prior distribution)の概要と具体例

「尤度」と「事前確率」に関してベイズの定理を適用することでパラメータの事後分布が得られます。当記事では事前分布の選定にあたって用いられることのあるジェフリーズの事前分布(Jeffreys prior distribution)の概要と具体例に関して取り扱いました。
「数理統計学 統計的推論の基礎(共立出版)」の第$12$章「ベイズ推論」の内容などを参考に、作成を行いました。

・参考
共役事前分布(Conjugate Prior Distribution)まとめ

ジェフリーズの事前分布の概要

概要

ジェフリーズの事前分布(Jeffreys prior distribution)はフィッシャー情報量を$I(\theta)$とおく際に、$\sqrt{I(\theta)}$に比例するように定義される事前分布である。

フィッシャー情報量

パラメータ$\theta$の確率分布から標本$\mathbf{x} = x_1, \cdots , x_n$が得られたとき、尤度を$L(\theta|\mathbf{x})$とおく。このときフィッシャー情報量$I(\theta)$は下記のように定められる。
$$
\large
\begin{align}
I(\theta) = E \left[ \left( \frac{\partial \log{L(\theta|\mathbf{X})}}{\partial \theta} \right)^{2} \right] = – E \left[ \frac{\partial^2 \log{L(\theta|\mathbf{X})}}{\partial \theta^2} \right]
\end{align}
$$

ジェフリーズ事前分布の具体例

二項分布$\mathrm{Bin}(n,p)$

確率変数$X$に関して$X \sim \mathrm{Bin}(n,p)$が成立するとき、尤度を$L(p|X)$とおくと$L(p|X), \log{L(p|X)}$下記のように表せる。
$$
\large
\begin{align}
L(p|X) &= {}_{n} C_{X} p^{X} (1-p)^{n-X} \\
\log{L(p|X)} &= X \log{p} + (n-X) \log{(1-p)} + \mathrm{Const.}
\end{align}
$$

ここで$\log{L(p|X)}$を$p$で偏微分すると下記が得られる。
$$
\large
\begin{align}
\frac{\partial \log{L(p|X)}}{\partial p} &= \frac{X}{p} – \frac{n-X}{1-p} \\
&= \frac{X(1-p) – p(n-X)}{p(1-p)} = \frac{X – \cancel{Xp} – np + \cancel{Xp})}{p(1-p)} \\
&= \frac{X – np}{p(1-p)}
\end{align}
$$

このときフィッシャー情報量$I(p)$は下記のように得られる。
$$
\large
\begin{align}
I(p) &= E \left[ \left( \frac{\partial \log{L(p|\mathbf{x})}}{\partial p} \right)^{2} \right] \\
&= E \left[ \left( \frac{X – np}{p(1-p)} \right)^{2} \right] \\
&= \frac{1}{p^2(1-p)^2} E[X^2 – 2npX + n^2p^2] \\
&= \frac{1}{p^2(1-p)^2} (V[X] + E[X]^2 -2npE[X] + n^2p^2) \\
&= \frac{1}{p^2(1-p)^2} (np(1-p) + \cancel{n^2p^2} -2 \cancel{n^2p^2} + \cancel{n^2p^2}) \\
&= \frac{n}{p(1-p)}
\end{align}
$$

よって$\displaystyle \sqrt{I(p)} = \sqrt{\frac{n}{p(1-p)}}$に比例する事前分布を考えれば良いが、二項分布の共役事前分布であるベータ分布を元に考えると、$\displaystyle \mathrm{Be} \left( \frac{1}{2},\frac{1}{2} \right)$がジェフリーズ事前分布に該当する。

「数理統計学 統計的推論の基礎(共立出版)」の例$12.5$の内容を参考に上記は作成を行ったので合わせて参照すると良い。

ポアソン分布$\mathrm{Po}(\lambda)$

確率変数$X_i$に関して$X_1, X_2, \cdots , X_n \sim \mathrm{Po}(\lambda), \quad \mathrm{i.i.d}.,$が成立するとき、$\mathrm{X}=X_1, \cdots , X_n$に関して尤度を$L(\lambda|\mathrm{X})$とおくと$L(\lambda|\mathrm{X}), \log{L(\lambda|\mathrm{X})}$は下記のように表せる。
$$
\large
\begin{align}
L(\lambda|\mathrm{X}) &= \prod_{i=1}^{n} \frac{\lambda^{X_i e^{-\lambda}}}{X_i} \\
&= \prod_{i=1}^{n} \exp{(X_i \log{\lambda} – \lambda – \log{X_i})} \\
&= \exp{\left[ \sum_{i=1}^{n} \left( X_i \log{\lambda} – \lambda – \log{X_i} \right) \right]} \\
\log{L(\lambda|\mathrm{X})} &= \log{ \left( \exp{\left[ \sum_{i=1}^{n} \left( X_i \log{\lambda} – \lambda – \log{X_i} \right) \right]} \right) } \\
&= \sum_{i=1}^{n} \left( X_i \log{\lambda} – \lambda – \log{X_i} \right)
\end{align}
$$

ここで$\log{L(\lambda|\mathrm{X})}$を$\lambda$で偏微分すると下記が得られる。
$$
\large
\begin{align}
\frac{\partial \log{L(\lambda|\mathrm{X})}}{\partial \lambda} &= \sum_{i=1}^{n} \left( \frac{X_i}{\lambda} – 1 \right) \\
\frac{\partial^2 \log{L(\lambda|\mathrm{X})}}{\partial \lambda^2} &= \sum_{i=1}^{n} \left( -\frac{X_i}{\lambda^2} \right)
\end{align}
$$

このときフィッシャー情報量$I(\lambda)$は下記のように得られる。
$$
\large
\begin{align}
I(\lambda) &= – E \left[ \frac{\partial^2 \log{L(\lambda|\mathbf{x})}}{\partial \theta^2} \right] \\
&= – E \left[ \sum_{i=1}^{n} \left( -\frac{X_i}{\lambda^2} \right) \right] \\
&= \frac{n}{\lambda^2} E[X_i^2] \\
&= \frac{n}{\lambda^2} \times \lambda \\
&= \frac{n}{\lambda}
\end{align}
$$

よって$\displaystyle \sqrt{I(\lambda)} = \sqrt{\frac{n}{\lambda}}$に比例するジェフリーズ事前分布に$1/\sqrt{\lambda}$を考えることができる。

最急降下法(Gradient Descent)とステップ幅の直線探索

勾配ベクトルを用いた漸化式的に表される反復法を用いて最適解を計算する手法を最急降下法(Gradient Descent)といいます。当記事では最急降下法の数式と、ステップ幅の適応的計算にあたって用いられる直線探索について取りまとめを行いました。
「数理計画入門(朝倉書店)」の$4.4$節の「最急降下法」やWikipediaの内容を参考に作成を行いました。

・用語/公式解説
https://www.hello-statisticians.com/explain-terms

数式の確認

最急降下法

最急降下法の数式

最急降下法の漸化式は下記のように表される。
$$
\large
\begin{align}
\mathbf{x}_{k+1} = \mathbf{x}_{k} – \alpha_{k} \nabla f(\mathbf{x}_{k}) \quad (1)
\end{align}
$$

最急降下法の解釈

$(1)$式における$\nabla f(\mathbf{x}_{k})$は$\mathbf{x}=\mathbf{x}_{k}$における勾配ベクトルであり、最急降下法は勾配の計算結果に基づいてベクトル$\mathbf{x}$の位置を更新すると解釈できる。

ステップ幅の決定

前節$(1)$式の$\alpha_{k}$の決定にあたっては主に「①定数を用いる方法」と、「②適応的計算を用いる方法」がある。以下、「定数を用いる方法」と都度計算を行う「適応的計算」の$1$例である「直線探索」について取りまとめを行なった。

定数を用いる

$(1)$式の$\alpha_{k}$には定数$\alpha$を用いることが多い。定数$\alpha$の値が大きい場合は$\mathbf{x}$が発散し、小さい場合は収束までに時間がかかることは抑えておくと良い。

直線探索

直線探索は「勾配ベクトルの方向に沿って目的関数が最小となるようなステップ幅$\alpha_{k}$を用いる手法」である。

Pythonを用いた最急降下法の実行

問題設定

「数理計画入門(朝倉書店)」の$4.4$節の「最急降下法」と同様に下記の関数の最適化にあたって、最急降下法のプログラムの作成を行う。
$$
\large
\begin{align}
f(\mathbf{x}) = f(x_1,x_2) = (x_1-1)^2 + 10(x_1^2-x_2)^2
\end{align}
$$

勾配ベクトルの計算

勾配ベクトル$\nabla f(\mathbf{x})$は下記のように計算できる。
$$
\large
\begin{align}
\nabla f(\mathbf{x}) &= \left(\begin{array}{c} \displaystyle \frac{\partial f}{\partial x_1} \\ \displaystyle \frac{\partial f}{\partial x_2} \end{array} \right) \\
&= \left(\begin{array}{c} 2(x_1-1) + 20(x_1^2-x_2) \cdot 2x_1 \\ 20(x_1^2-x_2) \cdot (-1) \end{array} \right) \\
&= \left(\begin{array}{c} 40x_1^3-40x_1x_2+2x_1-2 \\ -20(x_1^2-x_2) \end{array} \right)
\end{align}
$$

最急降下法の実行

下記を実行することで最急降下法による最適化を行うことができる。

import numpy as np

def calc_func_value(x1,x2):
    return (x1-1)**2 + 10*(x1**2-x2)**2

def calc_gradient(x1,x2):
    return np.array([40*x1**3 - 40*x1*x2 + 2*x1 - 2, -20*(x1**2-x2)])

x = np.array([0., 1.])
gamma = np.arange(-0.5, 0.5, 0.0000001)

display_flag = [1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 100, 200, 300, 400, 500]

for k in range(500):
    grad = calc_gradient(x[0],x[1])
    x1 = x[0]+gamma*grad[0]
    x2 = x[1]+gamma*grad[1]
    f = calc_func_value(x1,x2)
    x[0] = x1[np.argmin(f)]
    x[1] = x2[np.argmin(f)]
    if k+1 in display_flag:
        print("Iter: {:.0f}, x: {}".format(k+1,x))

・実行結果

Iter: 1, x: [ 0.0998848  0.001152 ]
Iter: 2, x: [ 0.36070906  0.02723477]
Iter: 3, x: [ 0.35167097  0.11761506]
Iter: 4, x: [ 0.44424751  0.12687296]
Iter: 5, x: [ 0.43824528  0.18689394]
Iter: 10, x: [ 0.57217842  0.28642287]
Iter: 20, x: [ 0.6772309   0.43291189]
Iter: 30, x: [ 0.7394933   0.52795907]
Iter: 40, x: [ 0.78289396  0.59811139]
Iter: 50, x: [ 0.81554999  0.65307488]
Iter: 100, x: [ 0.90619007  0.81570209]
Iter: 200, x: [ 0.96896765  0.93720982]
Iter: 300, x: [ 0.98869242  0.97691039]
Iter: 400, x: [ 0.99575609  0.99130582]
Iter: 500, x: [ 0.9983905   0.99669874]

上記の結果は「数理計画入門(朝倉書店)」の表$4.1$の結果に基本的には一致する。「数理計画入門(朝倉書店)」に直線探索の設定がないので、厳密な再現ではないことに注意が必要である。

参考

matplotlibを利用したデータ可視化(data visualization)の基礎

データの可視化(data visualization)は、データを様々な角度から確認し、データ自体を理解する目的で行われることが多い印象です。探索的データ分析(EDA)と呼ばれることもあります。pythonを使ってデータを可視化する際に、手軽に静的なデータの可視化ができるライブラリとしてmatplotlibがよく使われていると思います。ここでは、matplotlibで手軽にデータを可視化するために覚えておくと良いと思われる基礎的な項目をまとめます。

matplotlib

matplotlibは、pythonにおけるデータ可視化のためのライブラリです。matplotlibを利用すると、簡潔に様々なグラフをプロットできます。主に静的なグラフのプロットに使われることが多い印象です。詳細なドキュメント以下の公式ドキュメントを参照してください。

https://matplotlib.org/

matplotlibの導入

インストール

インストールについては、公式ドキュメントで解説されています。特にこだわりがなければ、pipを使ってインストールできます。

> pip install matplotlib

matplotlibの始めの一歩

matplotlibを利用したシンプルな例を紹介します。

import matplotlib.pyplot as plt

x = [0, 1, 2, 3, 4, 5]
y = [0, 2, 1, 3, 2, 4]

plt.plot(x, y)
plt.show()

このコードを実行すると以下のグラフがプロットされます。

ここで重要なポイントは、以下で説明する2点です。

ライブラリのimport

import matplotlib.pyplot as plt

この文で、matplotlibのpyplotモジュールをpltという名前でインポートしています。pyplotモジュールというのは、matplotlibにおけるグラフ描画のためのインターフェースです。pyplotモジュールの配下に散布図やヒストグラムなどを描画するための関数が用意されています。

pyplotモジュールについては後述の「matplotlibでの描画方針」に追加で説明を記載します。

プロットの実行

plt.plot(x, y)

この文で折れ線グラフを描画しています。matplotlib.pyplotではこのように、描画したいグラフに対応する関数を呼び出すだけでグラフを描画することが出来ます。

matplotlibでの描画方針

matplotlibでグラフを描画するには以下の二つの方針があります。

  1. グラフ要素のオブジェクトを生成してグラフ描画していく方針(object-oriented interface)
  2. pyplotモジュールを利用して予め用意された関数を利用してグラフ描画をしていく方針(state-based interface)

1のobject-orientedな方針では、グラフ描画にあたって細かい指定ができるのですが、その分コードは複雑になります。2のstate-basedな方針では、matplotlibのpyplotモジュールで提供されている関数を呼び出すだけなので簡単に描画を実行できます。

本稿では、基本的には2の方針での解説を行います。ですが、複数のグラフを1枚にまとめる(subplot)などでオブジェクト思考的な利用をすることがあります。

matplotlibでの描画要素

matplotlibにおける描画要素(描画オブジェクト)は、pyplotモジュールを単純に利用するだけの場合はあまり意識する必要はありません。しかし、軸の調整や複数グラフをまとめる場合などで意識しておくとより理解が深まります。

このmatplotlibにおけるオブジェクトを理解する上で公式ドキュメントにある図(こちらを参照)がとてもわかりやすいです。ただ、基礎を理解するという視点では、少々複雑なので、必要な要素を抽出すると、以下の3点を覚えておくと良いです。

  • Figure: グラフをプロットする全体の領域。ここに下記のAxesオブジェクトを配置していくイメージ。
  • Axes: グラフ自体。Figureの中に複数のAxesを配置することが可能。
  • Axis: 軸。軸の大きさや表示範囲、スケールなどを定義する。AxesとAxis紛らわしいが別物なので注意。

基本的なグラフ描画

散布図

プロットの概要縦軸と横軸にそれぞれ対応する量を割り当て、データを点でプロットしたもの。
プロットから主に確認できるポイント縦軸と横軸に割り当てた量の相関関係。2次元のデータの分布。
向いているデータデータの順番に意味がないデータ(独立なデータ)

例として、y=2x+1に従ってランダムノイズが加わった線型モデルからランダムにサンプルしたデータを散布図で確認してみます。

N = 100
x = np.random.randn(N) * 1.0 + 3.0 # N(3, 1)からランダムに生成
y = 2. * x + 1 + np.random.randn(N)

plt.scatter(x, y)

上記コードを実行すると以下のような図が作成されます。

折れ線グラフ

プロットの概要散布図の一種で、データ間を線で結んだもの。
横軸に設定する量を基準にデータを並べることが多い。
プロットから主に確認できるポイント時間経過での変化など、データがどのように変わっていったかを確認できる。
向いているデータ時系列データに代表される、データ間の順番に意味のあるデータ

例として、横軸(時間)に対して指数関数に従って値が上昇していくようなデータをプロットして確認してみます。

t = np.arange(0, 10)
y = np.exp(0.1 * t)

plt.plot(t, y)

上記コードを実行すると以下のような図が作成されます。

折れ線グラフでは、データポイントを明示するように記号をつけることも多いです。下記のようにpyplt.plotのオプションを設定することで実現できます。

t = np.arange(0, 10)
y = np.exp(0.1 * t)

plt.plot(t, y, "-o")

コード例では、pyplot.plot の3つ目の引数として、フォーマットを指定しています。このフォーマットで線の色なども設定できます。詳しくはドキュメントを参照してください。

棒グラフ

プロットの概要データの大きさを棒(長方形)の高さで表現したグラフ。
横軸には特に既定はないが、カテゴリ(質的)データを使うことが多い印象(e.g. 営業員毎の売り上げなど)。
プロットから主に確認できるポイント横軸に設定する値ごとの量の比較
向いているデータ質的なデータとその量

例として、3つのカテゴリの量を比較する棒グラフを描いてみます。

x = np.array(["a", "b", "c"])
y = np.array([3.2, 1.1, 2.5])

plt.bar(x, y)

上記コードを実行すると以下のような図が作成されます。

ヒストグラム

プロットの概要度数分布表の度数を棒の高さで表現したグラフ。横軸には階級をとる。
プロットから主に確認できるポイントデータの分布状況を可視化する。
向いているデータ1次元の数量データ

例として、正規分布に従う乱数のヒストグラムを描いてみます。

x = np.random.randn(1000) * 1.0 + 0.0

plt.hist(x)

上記コードを実行すると以下のような図が作成されます。

ヒストグラムにとって、横軸(階級)は重要です。階級の幅を適切に設定しないとデータの解釈をミスリードする場合があります。例えば、上記のヒストグラムの階級幅を細かく設定してみます。matplotlibでは、pyploy.histbins に数値を設定するとデータの最大最小範囲を設定した値で均等に分割してくれます。

plt.hist(x, bins=100)

上記コードを実行すると以下のような図が作成されます。

単純な標準正規分布$\mathcal{N}\left(0, 1 \right)$から乱数を生成したのですが、上記の図ではいくつか山があるように見えないこともないですよね。適切な解像度でデータを確認しないとノイズに過度に適合してしまうデータ解釈をしてしまうことがあります。

なお、pyploy.histbinslistnumpy.arrayなどを設定すると、任意の階級を設定することができます。

グラフの修飾

色の変更

軸の変更

スケール、描画範囲、対数軸

日本語の利用

日本語などの文字を利用したい場合、フォントを指定する必要があります。matplotlibで日本語フォントを利用するには複数の方法があります(下記)。

  1. 設定ファイルを作成してフォントを指定する
  2. 日本語フォントに対応したFontPropertyを作成し、日本語を使う箇所で都度指定する
  3. japanize-matplotlibを利用する

1と2の方法は似ていますが、1の方法では予め設定ファイルを作成し、その中でデフォルトフォントを設定します。2の方法では、ソースコード単位でFont Propertyを作成します。3の方法ではmatplotlibを日本語化するライブラリが公開されており、そちらを利用する方法です。3の方法が最も手間が少なく簡単だと思います。

ここでは、2の方法について、フォントファイルの入手から手順を解説します。

まずフォントファイルのダウンロードと展開をします。フォントファイルはシステムにインストールされているフォントファイルを使っても良いですが、Linux環境などシステムフォントがインストールされていないケースもあると思います。

wget https://moji.or.jp/wp-content/ipafont/IPAexfont/ipaexg00401.zip
unzip ipaexg00401.zip

上記の例では、ターミナルでのコマンド実行を想定していますが、jupyterlabなどから実行する場合には、システムコマンドを利用するために各文の先頭に ! をつけます。

FontPropertyの作成と利用について、以下に例を示します。

import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
font_path = './ipaexg00401/ipaexg.ttf'
fp = fm.FontProperties(fname=font_path)

x = [0, 1, 2, 3, 4, 5]
y = [0, 2, 1, 3, 2, 4]
plt.plot(x, y, label="日本語描画の例")
plt.title("日本語タイトル", fontproperties=fp)
plt.legend(prop=fp)
plt.show()

上記のコードを実行すると、以下の図が出力されます。

上記のコードにおける重要なポイントを解説します。

import matplotlib.font_manager as fm
font_path = './ipaexg00401/ipaexg.ttf'
fp = fm.FontProperties(fname=font_path)

matplotlib.font_manager.FontPropertiesにフォントファイルを指定します。これで指定したフォントファイルを使うことが出来ます。なおこの例にあるように、フォントファイルはどこに展開したものでも構いません。フォントなどのインストールが自由に出来ない環境などに有効だと思います。

plt.plot(x, y, label="日本語描画の例")
plt.title("日本語タイトル", fontproperties=fp)
plt.legend(prop=fp)

plotのlabelとタイトルに日本語を利用しています。フォントを指定しない場合、日本語が表示できず四角で置き換わってしまうと思いますが、上図の通り、日本語表示できるようになっています。

pyplot.titleとpyplot.legendでフォントを指定するオプショナル引数が異なっているので注意してください。

[参考]

応用的な利用方法

3次元のデータの可視化

scatterplot, heatmap

複数のグラフをまとめる

1枚のFigureに複数のグラフをまとめたいケースはいくつか考えられると思います。このような場合、pyplotのsubplots関数を利用すると比較的容易に実現できます。複数のグラフをどのように配置するかという考え方については、上記「matplotlibでの描画要素」におけるFigure要素とAxes要素を意識すると良いです。

まず具体例を以下に示します。

x = np.random.randn(100) * 1.0 + 3.0 # N(3, 1)の分布から50個のデータを取得
y = 2. * x + 0.5
z = y + np.random.randn(100)


fig, axs = plt.subplots(1, 2, figsize=(10, 4))
axs[0].hist(x)
axs[1].scatter(x, z)
fig.show()

上記のコードを実行すると、以下の図が出力されます。なお乱数を利用しているため、実行毎に結果は変わります。

上記のコードにおける重要なポイントを解説します。

fig, axs = plt.subplots(1, 2, figsize=(10, 4))

まずここで、pyplot.subplotsを呼び出し、グラフを描画するFigure要素とAxes要素のオブジェクトを生成します。この時、subplotsの第1引数として行数(nrows)、第2引数に列数(ncols)を指定しているので、axsには二つのAxesオブジェクトのリストが返ってきます。

axs[0].hist(x)
axs[1].scatter(x, z)

二つのAxesオブジェクトそれぞれで、ヒストグラム(hist)と散布図(scatter)を描画しています。

[参考]

軸(Axis)を共有して複数のグラフをまとめる

一つのFigureに独立したグラフを並べるだけでは、複数のグラフをそれぞれ描画するのと変わりません。ここでは、複数のグラフの軸のを合わせて描画する方法について解説します。

まず、具体例を示します。

x1 = np.random.randn(500) * 1.0 + 0.0 # N(0, 1)の分布から500個のデータを取得
x2 = np.random.randn(500) * 2.0 + 5.0 # N(5, 2)の分布から500個のデータを取得

fig, axs = plt.subplots(2, 1, sharex=True)
axs[0].hist(x1)
axs[1].hist(x2)
fig.show()

上記のコードを実行すると、以下の図が出力されます。なお乱数を利用しているため、実行毎に結果は変わります。

上記のコードにおける重要なポイントを解説します。

fig, axs = plt.subplots(2, 1, sharex=True)

pyplot.subplotsのオプショナル引数である sharex=True を指定しています。これで、二つのグラフにおけるx軸が共有されたことになります。同様に、y軸を共有するには sharey=True を指定します。

[参考]

まとめ

参考

数学検定2級 解説 〜過去問題② 解答例・解説 1次:計算技能検定・2次:数理技能検定〜

数学検定$2$級は数ⅡBまで相当の数学の基本トピックに関して取り扱った検定であり、統計学に必要な数学を身につける際の目安になります。当記事では「日本数学検定協会 監修」の「数学検定問題集 $2$級」の数学検定$2$級の内容に基づき、過去問題②の解答例と解説の作成を行いました。

・数学検定$2$級まとめ
https://www.hello-statisticians.com/math_certificate_2

$1$次:計算技能検定

問題$1 \,$ 展開の公式

下記のように変形できる。
$$
\large
\begin{align}
(x+2y)^{3} = x^3 + 6x^2 y + 12x y^2 + 8y^3
\end{align}
$$

問題$2 \,$ 因数分解

下記のように因数分解を行うことができる。
$$
\large
\begin{align}
a(b+c) – (b^2-c^2) &= a(b+c) – (b+c)(b-c) \\
&= (b+c)(a-b+c)
\end{align}
$$

問題$3 \,$ 分母の有理化

下記のように式変形できる。
$$
\large
\begin{align}
\frac{\sqrt{2}+3}{\sqrt{2}-1} – \frac{\sqrt{2}-3}{\sqrt{2}+1} &= \frac{((\sqrt{2}-1))(\sqrt{2}+3)}{(\sqrt{2}+1)(\sqrt{2}-1)} – \frac{(\sqrt{2}-3)((\sqrt{2}-1))}{(\sqrt{2}+1)((\sqrt{2}-1))} \\
&= \frac{2 + 4 \sqrt{2} + 3}{2-1} – \frac{2 – 4 \sqrt{2} + 3}{2-1} \\
&= 8 \sqrt{2}
\end{align}
$$

問題$4 \,$ 三角関数を用いた方程式

$\displaystyle \sin{\theta}-\cos{\theta} = \frac{1}{2}$より下記が得られる。
$$
\large
\begin{align}
(\sin{\theta}-\cos{\theta})^{2} &= \frac{1}{2} \\
\sin^{2}{\theta} + \cos^{2}{\theta} – 2 \sin{\theta} \cos{\theta} &= \frac{1}{4} \\
1 – 2 \sin{\theta} \cos{\theta} &= \frac{1}{4} \\
2 \sin{\theta} \cos{\theta} &= \frac{3}{4} \\
\sin{\theta} \cos{\theta} &= \frac{3}{8}
\end{align}
$$

問題$5 \,$ サイコロの目の確率

目の出方は$(1,3), (2,2), (3,1), (2,6), (3,5), (4,4), (5,3), (6,2), (6,6)$の$9$通りある。よって、確率は下記のように計算できる。
$$
\large
\begin{align}
\frac{9}{6^2} = \frac{1}{4}
\end{align}
$$

問題$6 \,$ 循環小数と数列の和

$x = 2.\dot{0}\dot{9}$とおくと、$100x = 209.\dot{0}\dot{9}$である。このとき$100x – x$は下記のように計算できる。
$$
\large
\begin{align}
100x – x &= 209.\dot{0}\dot{9} – 2.\dot{0}\dot{9} \\
99x &= 207 \\
x &= \frac{99}{207} = \frac{23}{11}
\end{align}
$$

問題$7 \,$ 二次方程式の判別式

二次方程式$2x^2 + 8x – 2k^2 – 9k + 13 = 0$が解を持つ$k$の範囲が得られれば良いので、判別式を$D$とおくと下記のような式が得られる。
$$
\large
\begin{align}
\frac{D}{4} = 4^2 – 2(-2k^2 – 9k + 13) & > 0 \\
2k^2 + 9k – 5 & > 0 \\
(2k – 1)(k + 5) & > 0
\end{align}
$$

上記より$\displaystyle k < -5, \frac{1}{2} < k$が得られる。

問題$8 \,$ 円の方程式

中心が$(2,0)$の円の方程式は下記のように表せる。
$$
\large
\begin{align}
(x-2)^2 + y^2 = r^2
\end{align}
$$

上記が原点を通るので下記が成立する。
$$
\large
\begin{align}
(0-2)^2 + 0^2 &= r^2 \\
r^2 &= 2^2
\end{align}
$$

よって、円の方程式は下記のように表せる。
$$
\large
\begin{align}
(x-2)^2 + y^2 = 2^2
\end{align}
$$

問題$9 \,$ 分母の複素数の実数化

下記の式変形が成立する。
$$
\large
\begin{align}
a + bi &= \frac{2-i}{1+i} \\
&= \frac{(2-i)(1-i)}{(1+i)(1-i)} \\
&= \frac{2-3i+i^2}{1-i^2} \\
&= \frac{2-1-3i}{1-(-1)} \\
&= \frac{1}{2} – \frac{3}{2}i
\end{align}
$$

上記より$\displaystyle a = \frac{1}{2}, b=-\frac{3}{2}$が成立する。

問題$10 \,$ 等比数列の一般項

一般項$a_n$は下記のように表せる。
$$
\large
\begin{align}
a_n = -64 \cdot \left( -\frac{1}{2} \right)^{n}
\end{align}
$$

よって$a_7$は下記のように計算できる。
$$
\large
\begin{align}
a_n &= -64 \cdot \left( -\frac{1}{2} \right)^{7} \\
&= 2^{6} \cdot \left( \frac{1}{2} \right)^{7} \\
&= \frac{1}{2}
\end{align}
$$

問題$11 \,$ 指数の計算

下記のように計算できる。
$$
\large
\begin{align}
3^{\frac{1}{6}} \div 3^{\frac{1}{2}} \times 3^{\frac{1}{3}} &= 3^{\frac{1}{6}-\frac{1}{2}+\frac{1}{3}} \\
&= 3^{0} \\
&= 1
\end{align}
$$

問題$12 \,$ 三角関数の二次関数

方程式は下記のように解ける。
$$
\large
\begin{align}
2 \sin^{2}{\theta} – 3 \sin{\theta} + 1 &= 0 \\
(2 \sin{\theta} – 1)(\sin{\theta} – 1) &= 0 \\
\sin{\theta} &= \frac{1}{2}, \, 1
\end{align}
$$

上記より$\displaystyle \theta = 30^{\circ}, 90^{\circ}, 150^{\circ}$である。

問題$13 \,$ 三次方程式の解

$$
\large
\begin{align}
x^3 + x^2 – 8x – 12 = 0
\end{align}
$$

上記に$x=3$を代入すると方程式が成立するので、左辺を$x-3$で割ることで下記のように変形できる。
$$
\large
\begin{align}
x^3 + x^2 – 8x – 12 &= 0 \\
(x-3)(x^2+4x+4) &= 0 \\
(x-3)(x+2)^2 &= 0
\end{align}
$$

上記より$x=-2, 3$である。

問題$14 \,$ 不定積分・定積分

$[1]$
$$
\large
\begin{align}
\int x(3x+2) dx &= \int (3x^2 + 2x) dx \\
&= x^3 + x^2 + C
\end{align}
$$

$[2]$
$$
\large
\begin{align}
\int_{-1}^{2} x(3x+2) dx &= \left[ x^3 + x^2 \right]_{-1}^{2} \\
&= (2^3+2^2) – ((-1)^3+(-1)^2) \\
&= 12
\end{align}
$$

問題$15 \,$ 空間ベクトルの成分表示

原点を$O$と定めると$\overrightarrow{OA}$と$\overrightarrow{OB}$はそれぞれ下記のように表せる。
$$
\large
\begin{align}
\overrightarrow{OA} = \left(\begin{array}{c} 2 \\ -2 \\ 0 \end{array} \right) , \, \overrightarrow{OB} = \left(\begin{array}{c} 0 \\ 2 \\ -5 \end{array} \right)
\end{align}
$$

$[1]$
$\overrightarrow{AB}$は下記のように計算できる。
$$
\large
\begin{align}
\overrightarrow{AB} &= \overrightarrow{OB} – \overrightarrow{OA} \\
&= \left(\begin{array}{c} 0 \\ 2 \\ -5 \end{array} \right) – \left(\begin{array}{c} 2 \\ -2 \\ 0 \end{array} \right) \\
&= \left(\begin{array}{c} -2 \\ 4 \\ -5 \end{array} \right)
\end{align}
$$

$[2]$
$|\overrightarrow{AB}|$は下記のように計算できる。
$$
\large
\begin{align}
|\overrightarrow{AB}| &= \sqrt{(-2)^2 + 4^2 + (-5)^2} \\
&= \sqrt{45} \\
&= 3 \sqrt{5}
\end{align}
$$

$2$次:数理技能検定

問題$1 \,$ 期待値の計算

コインの枚数の期待値を$E[X]$とおくと、$E[X]$は下記のように計算できる。
$$
\large
\begin{align}
E[X] &= 10 \cdot \frac{5}{100} + 4 \cdot \frac{25-5}{100} + 6 \cdot \frac{20-5}{100} \\
&= \frac{50+80+90}{100} \\
&= \frac{220}{100} = \frac{11}{5}
\end{align}
$$

問題$2 \,$

問題$3 \,$ ヘッセの標準形

直線$l$と平行なベクトルを$\vec{l}$とおくと、$\overrightarrow{OH}$と$\vec{l}$は下記のように表すことができる。
$$
\large
\begin{align}
\overrightarrow{OH} &= \left(\begin{array}{c} p \cos{\theta} \\ p \sin{\theta} \end{array} \right) \\
\vec{l} &= \left(\begin{array}{c} x-p \cos{\theta} \\ y-p \sin{\theta} \end{array} \right)
\end{align}
$$

このとき$\overrightarrow{OH}$と$\vec{l}$が垂直なので下記が成立する。
$$
\large
\begin{align}
\overrightarrow{OH} \cdot \vec{l} &= 0 \\
\left(\begin{array}{c} p \cos{\theta} \\ p \sin{\theta} \end{array} \right) \cdot \left(\begin{array}{c} x-p \cos{\theta} \\ y-p \sin{\theta} \end{array} \right) &= 0 \\
xp \cos{\theta} – p^{2} \cos^{2}{\theta} + yp \sin{\theta} – p^2 \sin^{2}{\theta} &= 0 \\
p(x \cos{\theta} + y \sin{\theta}) &= p^2(\cos^{2}{\theta}+\sin^{2}{\theta}) \\
x \cos{\theta} + y \sin{\theta} &= p \quad (1)
\end{align}
$$

$\overrightarrow{OH} \cdot \vec{l} = 0$はベクトルの大きさが$0$である場合も成立するので上記は必要条件であるが、$p>0$であれば$|\overrightarrow{OH}| \neq 0$、$|\vec{l}|=0$は直線上の点であることから十分条件であることも確認できる。よって$(1)$式が直線の方程式を表すことが確認できる。

問題$4 \,$

問題$5 \,$

問題$6 \,$ 二次方程式の解の公式の導出

下記のように導出を行うことができる。
$$
\large
\begin{align}
ax^2 + bx + c &= 0 \\
a \left( x^2 + \frac{b}{a}x \right) + c &= 0 \\
a \left( x + \frac{b}{2a} \right)^2 -\frac{ab^2}{(2a)^2} + c &= 0 \\
a \left( x + \frac{b}{2a} \right)^2 &= \frac{b^2}{4a} – c \\
a \left( x + \frac{b}{2a} \right)^2 &= \frac{b^2-4ac}{4a} \\
\left( x + \frac{b}{2a} \right)^2 &= \frac{b^2-4ac}{4a^2} \\
x + \frac{b}{2a} &= \frac{\pm \sqrt{b^2-4ac}}{2a} \\
x &= \frac{-b \pm \sqrt{b^2-4ac}}{2a}
\end{align}
$$

問題$7 \,$

【ゲーム×幾何分布】幾何分布の確率分布から考察するランダムエンカウント

幾何分布を理解するにあたって、初期のドラゴンクエストのようなRPGのエンカウントの調整の観点から確認することもできます。当記事では、初期のドラゴンクエストのエンカウントの設定とゲームの操作性について幾何分布の確率分布を考えることで考察を行います。

・ゲーム × 統計 まとめ
https://www.hello-statisticians.com/game_stat

前提の確認

エンカウントの仕組み

近年の$3$D化されたRPGではモンスターなどがフィールド上を動き回るシンボルエンカウントが増えましたが、従来の$2$Dのゲームではフィールドを動き回る際に予め設定された仕組みでエンカウントするようにプログラムされる場合が主流でした。

中でもファミコンの黎明期などは$1$歩ごとに一定の確率でエンカウント判定を行う「ランダムエンカウント」が多く、$0$歩エンカの映像なども探せば見つけることができます。

ファミコンのゲーム自体は旧世代のものではありますが、ゲームの作成側の視点でエンカウント調整を行うと考えるとなかなか難しく、「ランダムエンカウント」が黎明期に用いられたのは自然であるように考えられます。

一方で、$1$歩ごとに判定を行う「ランダムエンカウント」を用いる場合、次のエンカウントまでの歩数が「幾何分布」に従うことから、次のエンカウントまでの歩数のばらつきが大きいなどが問題になります。

以下では幾何分布の確率分布を元に「ランダムエンカウント」を用いた際の「次のエンカウントまでの歩数」がどのくらいばらつくかに関して計算や考察などを行います。

幾何分布

下記などに詳しくまとめました。

ドラゴンクエストにおけるランダムエンカウント

幾何分布の期待値

確率$p$でエンカウントする際に次のエンカウントまでの歩数を確率変数$X$で表すことを考えます。このとき、$X=x$歩でエンカウントする確率関数を$p(x)$とおくと、$p(x)$は下記のように表すことができます。
$$
\large
\begin{align}
p(x) = p(1-p)^{x-1}, \quad x=1,2,\cdots
\end{align}
$$

上記では簡易化にあたって、$0$歩目のエンカウントを考慮しないようにしました。この場合、「次のエンカウントまでの歩数」の期待値$E[X]$は下記のように計算することができます。
$$
\large
\begin{align}
E[X] &= \sum_{x=1}^{\infty} xP(X=x) \\
&= \sum_{x=1}^{\infty} xp(1-p)^{x-1} \\
&= p \sum_{x=1}^{\infty} x(1-p)^{x-1} \\
&= p \left( 1 + 2(1-p)^{2-1} + 3(1-p)^{3-1} \cdots \right) \\
&= p \frac{1}{(1-(1-p))^2} \\
&= p \frac{1}{p^2} \\
&= \frac{1}{p}
\end{align}
$$

よって、$p=0.1$であれば$E[X]=10$であり、エンカウント率を$10$%で設定することで、だいたい次のエンカウントまで$10$歩であろうと直感的には考えることができます。

ランダウの記号(Laudau symbol)・ランダウの漸近記法の定義と活用例の確認

テイラー展開による近似を行う際に「特定の次数までを取り扱い以降は誤差と見なす」ということはよく行われます。ここで誤差と見なす項を略記することで数式展開が行いやすくなりますが、この際に用いられるランダウの漸近記法について当記事では定義や活用例の確認を行います。
主に「大学教養 微分積分(数研出版)」の$3.4$節の「テイラーの定理」の内容を参考に作成を行いました。

・数学まとめ
https://www.hello-statisticians.com/math_basic

ランダウの漸近記法

ランダウの漸近記法の定義

ランダウの記号$\mathcal{o}$は下記のように定義される。
$$
\large
\begin{align}
\lim_{x \to a} \frac{f(x)}{g(x)} = 0 \iff f(x) = \mathcal{o}(g(x)), \, x \to a
\end{align}
$$

上記のような記法をランダウの漸近記法という。$\mathcal{o}(g(x)), \, x \to a$は「$x$が$a$に近づくとき、$f(x)$は$g(x)$に対しはるかに小さい」を表すと解釈すれば良い。

ランダウの漸近記法の活用例

「大学教養 微分積分(数研出版)」の$3.4$節の「テイラーの定理」の例$3$の式は下記のように表せる。
$$
\large
\begin{align}
\lim_{x \to a} \frac{f(x)}{1} = \lim_{x \to a} f(x) = 0 & \iff f(x) = \mathcal{o}(1), \, x \to a \\
f(x) = 2x^2 + h(x), \, \lim_{x \to 0} \frac{h(x)}{g(x)} = 0 & \iff f(x) = 2x^2 + \mathcal{o}(x^2), \, x \to 0
\end{align}
$$

ランダウ記号の計算規則

$x \to 0$のとき下記が成立する。
$$
\large
\begin{align}
x^{m} \mathcal{o}(x^{n}) &= \mathcal{o}(x^{m+n}) \\
\mathcal{o}(x^{m}) \mathcal{o}(x^{n}) &= \mathcal{o}(x^{m+n}) \\
\mathcal{o}(x^{m}) + \mathcal{o}(x^{n}) &= \mathcal{o}(x^{l}), \, l = \min\{ m,n \}
\end{align}
$$

有限テイラー展開とランダウの漸近記法

式表記

$x \to a$における有限テイラー展開はランダウの記号を用いて下記のように表すことができる。
$$
\large
\begin{align}
f(x) = f(a) + f'(a)(x-a) + \cdots + \frac{f^{(n)}(a)}{n!}(x-a)^{n} + \mathcal{o}((x-a)^{n}), \, x \to a
\end{align}
$$

上記は「$(x-a)^{n+1}$以降の項を$\mathcal{o}((x-a)^{n})$でまとめた」と解釈すると良い。$x \to a$では$x-a$が$0$に近い値を取るので、次数が大きくなればなるほど$(x-a)^{n}$は$0$に近づくことに基づいて理解すると良い。下記を実行することで図を元に理解することができる。

import numpy as np
import matplotlib.pyplot as plt

n = np.arange(1, 11, 1)

a = 1.
x = np.array([0.9, 0.95, 0.99, 0.995, 0.999])

for i in range(x.shape[0]):
    delta = (x[i]-a)**n
    plt.plot(n, delta, label="x: {}".format(x[i]))

plt.legend(loc="lower right")
plt.show()

・実行結果

上図より、「次数が大きくなればなるほど$(x-a)^{n}$が$0$に近づくこと」が確認できる。

テイラー展開におけるランダウの漸近記法の活用例

$x=0$における$(1+x^2)e^{x}$の$5$次の漸近展開を以下行う。$e^{x}$のテイラー展開はランダウの記号を用いて下記のように表せる。
$$
\large
\begin{align}
e^{x} &= 1 + x + \frac{1}{2!} x^2 + \frac{1}{3!} x^3 + \frac{1}{4!} x^4 + \frac{1}{5!} x^5 + \mathcal{o}(x^5) \quad (1) \\
&= 1 + x + \frac{1}{2!} x^2 + \frac{1}{3!} x^3 + \mathcal{o}(x^3) \quad (2)
\end{align}
$$

$(2)$式より$x^2 e^{x}$は下記のように表せる。
$$
\large
\begin{align}
x^2 e^{x} &= x^2 (1 + x + \frac{1}{2!} x^2 + \frac{1}{3!} x^3 + \mathcal{o}(x^3)) \\
&= x^2 + x^3 + \frac{1}{2!} x^4 + \frac{1}{3!} x^5 + x^2 \mathcal{o}(x^3) \\
&= x^2 + x^3 + \frac{1}{2!} x^4 + \frac{1}{3!} x^5 + \mathcal{o}(x^5) \quad (3)
\end{align}
$$

$(1), (3)$式より$(1+x^2)e^{x}$の$5$次の漸近展開は下記のように表せる。
$$
\large
\begin{align}
(1+x^2) e^{x} &= 1 + x + \left( \frac{1}{2!} + 1 \right) x^2 + \left( \frac{1}{3!} + 1 \right) x^3 + \left( \frac{1}{4!} + \frac{1}{2!} \right) x^4 + \left( \frac{1}{5!} + \frac{1}{3!} \right) x^5 + \mathcal{o}(x^5) \\
&= 1 + x + \frac{1}{2!} x^2 + \frac{7}{3!} x^3 + \frac{13}{4!} x^4 + \frac{21}{5!} x^5 + \mathcal{o}(x^5)
\end{align}
$$

上記は「大学教養 微分積分(数研出版)」の$3.4$節の「テイラーの定理」の例題$2$を元に作成を行なった。

多変数の凸関数の制約なし最適化問題における勾配やヘッセ行列を用いた1次・2次の最適性条件

多変数の凸関数の制約なし最適化問題を取り扱うにあたっては多変数関数の勾配(gradient)やヘッセ行列(hessian matrix)を用いることで$1$次・$2$次の最適性条件を表すことができます。当記事ではそれぞれの式の表記や簡単な導出などの取りまとめを行いました。
「数理計画入門(朝倉書店)」の$4.3$節の「制約なし問題の最適性条件」の内容を参考に作成を行いました。

・用語/公式解説
https://www.hello-statisticians.com/explain-terms

前提の確認

問題設定

多変数関数$f(\mathbf{x})=f(x_1,\cdots,x_n)$について勾配$\nabla f(\mathbf{x})$とヘッセ行列$\nabla^2 f(\mathbf{x})$を下記のように定義する。
$$
\large
\begin{align}
\nabla f(\mathbf{x}) &= \left( \begin{array}{c} \displaystyle \frac{\partial f(\mathbf{x})}{\partial x_1} \\ \vdots \\ \displaystyle \frac{\partial f(\mathbf{x})}{\partial x_n} \end{array} \right) \\
\nabla^2 f(\mathbf{x}) &= \left( \begin{array}{c} \displaystyle \frac{\partial^2 f(\mathbf{x})}{\partial x_1^2} & \cdots & \displaystyle \frac{\partial^2 f(\mathbf{x})}{\partial x_1 \partial x_n} \\ \vdots & \ddots & \vdots \\ \displaystyle \frac{\partial^2 f(\mathbf{x})}{\partial x_n \partial x_1} & \cdots & \displaystyle \frac{\partial^2 f(\mathbf{x})}{\partial x_n^2} \end{array} \right)
\end{align}
$$

ここで上記の「多変数関数$f(\mathbf{x})=f(x_1,\cdots,x_n)$の最小化」を「制約なし最適化問題」と見なし、以下局所最適解$\mathbf{x}^{*}$に関して成立する最適性条件について確認を行う。

凸関数

関数$f$が下に凸の凸関数であるとき下記が成立する。
$$
\large
\begin{align}
\mathbf{x}, \mathbf{y} \in \mathbb{R}, \, 0 \leq \alpha \leq 1 \, \implies \, f(\alpha \mathbf{x} + (1-\alpha) \mathbf{y}) \leq \alpha f(\mathbf{x}) + (1-\alpha) f(\mathbf{y})
\end{align}
$$

$f(x) = x^2$のような二次関数のようなイメージで理解すると良く、たとえば$x=0, y=2, \alpha=0.5$のとき$f(\alpha x + (1-\alpha) y), \alpha f(x) + (1-\alpha) f(y)$はそれぞれ下記のように計算できる。
$$
\large
\begin{align}
f(\alpha x + (1-\alpha) y) &= f(0.5 \cdot 0 + (1 – 0.5) \cdot 2) \\
&= f(1) \\
&= 1^2 = 1 \\
\alpha f(x) + (1-\alpha) f(y) &= 0.5 f(0) + (1 – 0.5) f(2) \\
&= 0.5 \cdot 0^2 + 0.5 \cdot 2^2 \\
&= 2
\end{align}
$$

上記より$f(\alpha x + (1-\alpha) y) \leq \alpha f(x) + (1-\alpha) f(y)$が$x=0, y=2, \alpha=0.5$で成立することが確認できる。

上記の例は下記を実行することで図示できる。

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(0., 2.01, 0.01)
f_x = x**2
f_line = 2*x

plt.plot(x, f_x, color="limegreen", label="f(x)=x^2")
plt.plot(x, f_line, color="orange", label="chord")

plt.scatter(0.5*0+0.5*2., (0.5*0+0.5*2.)**2, color="green")
plt.scatter(0.5*0+0.5*2., 0.5*0**2+0.5*2.**2, color="red")

plt.legend()
plt.show()

・実行結果

上図のように下に凸の凸関数の数式は「関数の上に弦が描画できる」と解釈することもできる。

最適性条件

最適性の1次の必要条件

$$
\large
\begin{align}
\nabla f(\mathbf{x}^{*}) &= \mathbf{0} \\
\left( \begin{array}{c} \displaystyle \frac{\partial f(\mathbf{x})}{\partial x_1} \\ \vdots \\ \displaystyle \frac{\partial f(\mathbf{x})}{\partial x_n} \end{array} \middle) \right|_{x_1=x_1^{*}, \cdots , x_n=x_n^{*}} &= \left( \begin{array}{c} 0 \\ \vdots \\ 0 \end{array} \right)
\end{align}
$$

最適性の$1$必要条件は上記のように表される。上記の式は「関数の最大値・最小値では各方向の勾配が$0$」であることに対応し、直感的には「山の頂上は平らである」と解釈すれば良い。

最適性の2次の必要条件

「ヘッセ行列$\nabla^2 f(\mathbf{x}^{*})$が半正定値行列」が$\mathbf{x}^{*}$が局所最適解である必要条件である。

最適性の2次の十分条件

$\mathbf{x}=\mathbf{x}^{*}$で「①勾配$\nabla f(\mathbf{x}^{*})$が$0$ベクトル」、「②ヘッセ行列$\nabla^2 f(\mathbf{x}^{*})$が正定値行列」であることが$\mathbf{x}^{*}$が局所最適解である十分条件である。このことを以下に示す。

目的関数$f$は最適解$\mathbf{x}^{*}$の周辺で下記のように表せる。
$$
\large
\begin{align}
f(\mathbf{x}^{*} + \alpha \mathbf{d}) = f(\mathbf{x}^{*}) + \alpha \nabla f(\mathbf{x}^{*})^{\mathrm{T}} \mathbf{d} + \frac{\alpha^2}{2} \mathbf{d}^{\mathrm{T}} \nabla^2 f(\mathbf{x}^{*}) \mathbf{d} + \mathcal{o}(||\mathbf{d}||^2)
\end{align}
$$

ここで$\nabla^2 f(\mathbf{x}^{*})$が正定値行列であれば、$\alpha$が十分小さいとき下記が成立する。
$$
\large
\begin{align}
\frac{\alpha^2}{2} \mathbf{d}^{\mathrm{T}} \nabla^2 f(\mathbf{x}^{}) \mathbf{d} + \mathcal{o}(||\mathbf{d}||^2) > 0
\end{align}
$$

よって、「①勾配$\nabla f(\mathbf{x}^{*})$が$0$ベクトル」、「②ヘッセ行列$\nabla^2 f(\mathbf{x}^{*})$が正定値行列」のとき下記が成立する。
$$
\large
\begin{align}
f(\mathbf{x}^{*} + \alpha \mathbf{d}) – f(\mathbf{x}^{*}) &= \alpha \nabla f(\mathbf{x}^{*})^{\mathrm{T}} \mathbf{d} + \frac{\alpha^2}{2} \mathbf{d}^{\mathrm{T}} \nabla^2 f(\mathbf{x}^{*}) \mathbf{d} + \mathcal{o}(||\mathbf{d}||^2) > 0 \\
f(\mathbf{x}^{*} + \alpha \mathbf{d}) – f(\mathbf{x}^{*}) & > 0 \\
f(\mathbf{x}^{*} + \alpha \mathbf{d}) & > f(\mathbf{x}^{*})
\end{align}
$$

上記より、「①勾配$\nabla f(\mathbf{x}^{*})$が$0$ベクトル」、「②ヘッセ行列$\nabla^2 f(\mathbf{x}^{*})$が正定値行列」が$\mathbf{x}^{*}$が最適解であることの十分条件であることが確認できる。

参考

・数理計画法 講義資料

Ch.17 「行列の基本変形」の演習問題の解答例 〜統計学のための数学入門30講(朝倉書店)〜

当記事は「統計学のための数学入門$30$講(朝倉書店)」の読解サポートにあたってChapter.$17$の「行列の基本変形」の章末問題の解答の作成を行いました。
基本的には書籍の購入者向けの解説なので、まだ入手されていない方は購入の上ご確認ください。また、解説はあくまでサイト運営者が独自に作成したものであり、書籍の公式ページではないことにご注意ください。

・書籍解答まとめ
https://www.hello-statisticians.com/answer_textbook_math#math_stat

本章のまとめ

行列の基本変形の概要

基本行列

行列の基本変形と逆行列

演習問題解答

問題$17.1$

$$
\large
\begin{align}
A = \left(\begin{array}{cccc} 1 & 1 & 0 & 1 \\ 0 & 1 & -1 & 0 \\ 0 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 \end{array} \right)
\end{align}
$$

上記は下記のように基本変形を行うことができる。
$$
\large
\begin{align}
A &= \left(\begin{array}{cccc} 1 & 1 & 0 & 1 \\ 0 & 1 & -1 & 0 \\ 0 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 \end{array} \right) \\
& \implies \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & -1 & 0 \\ 0 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 \end{array} \right) \\
& \implies \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 1 & 2 & 1 \\ 0 & 0 & 1 & 1 \end{array} \right) \implies \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 1/2 & 1 & 1/2 \\ 0 & 0 & 1 & 1 \end{array} \right) \\
& \implies \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & -1/2 & 1 & 1/2 \end{array} \right) \implies \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & -1 & 2 & 1 \end{array} \right) \\
& \implies \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array} \right)
\end{align}
$$

よって逆行列は存在することが確認できる。以下、基本変形を行うことによって逆行列の導出を行う。
$$
\large
\begin{align}
(A|I_4) &= \left(\begin{array}{cccc} 1 & 1 & 0 & 1 \\ 0 & 1 & -1 & 0 \\ 0 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 \end{array} \,\, \middle\vert \,\, \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array} \right) \\
& \implies \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & -1 & 0 \\ 0 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 \end{array} \,\, \middle\vert \,\, \begin{array}{cccc} 1 & -1 & 0 & -1 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array} \right) \\
& \implies \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 1 & 2 & 1 \\ 0 & 0 & 1 & 1 \end{array} \,\, \middle\vert \,\, \begin{array}{cccc} 1 & -1 & -1 & -1 \\ 0 & 1 & 1 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array} \right) \\
& \implies \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 1 & 1 & 1 \\ 0 & 0 & 1/2 & 1 \end{array} \,\, \middle\vert \,\, \begin{array}{cccc} 1 & -1 & -1/2 & -1 \\ 0 & 1 & 1/2 & 0 \\ 0 & 0 & 1/2 & 0 \\ 0 & 0 & 0 & 1 \end{array} \right) \\
& \implies \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & -1/2 & 1/2 & 1/2 \end{array} \,\, \middle\vert \,\, \begin{array}{cccc} 1 & -1/2 & -1/2 & -1/2 \\ 0 & 1/2 & 1/2 & -1/2 \\ 0 & -1/2 & 1/2 & -1/2 \\ 0 & 0 & 0 & 1 \end{array} \right) \\
& \implies \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1/2 \end{array} \,\, \middle\vert \,\, \begin{array}{cccc} 1 & -1 & 0 & -1/2 \\ 0 & 0 & 1 & -1/2 \\ 0 & -1 & 1 & -1/2 \\ 0 & 1 & -1 & 1 \end{array} \right) \\
& \implies \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array} \,\, \middle\vert \,\, \begin{array}{cccc} 1 & -1 & 0 & -1 \\ 0 & 0 & 1 & -1 \\ 0 & -1 & 1 & -1 \\ 0 & 1 & -1 & 2 \end{array} \right) = (I_4|A^{-1})
\end{align}
$$

問題$17.2$