2020-09-28

【統計学入門（東京大学出版会）】第10章練習問題解答

統計学

東京大学出版会から出版されている統計学入門（基礎統計学Ⅰ）について第10章の練習問題の解答を書いていく。

本章以外の解答
10.1
10.2
10.3
- i)
- ii)
10.4
10.5
10.6
10.7
- i)
- ii)
- iii)
- iv)
- v)
10.8
10.9
- i)
- ii)

本章以外の解答

本章以外の練習問題の解答は別の記事で公開している。
必要に応じて参照されたい。

10.1

誤差 $e$ の分布が $N(0, 0.1)$ であるから、その期待値は $E(e) = 0$ となる。したがって、測定値の標本分布 $\bar{X}$ の期待値は $E(\bar{X}) = 100 + E(e) = 100$ である。

一方で分散 $V(\bar{X})$ は、

$\begin{align} V(\bar{X}) &= \frac{V(e)}{n} \\ &= \frac{0.1}{10} \\ &= 0.01 \end{align}$

となる。

以上から、測定値の標本 $\bar{X}$ の分布は $N(100, 0.01)$ と求まる。

$N(100, 0.01)$ を $N(0, 1)$ に標準化する場合、

$\begin{align} z = \frac{\bar{X} - 100}{\sqrt{0.01}} \end{align}$

となるから、両側の確率を求めることに注意すると、

$\begin{align} P \left( |\bar{X} - 100|) \gt 0.3 \right) &= P \left( |z| \gt 3 \right) \\ &= 2 \cdot P \left( z \gt 3 \right) \end{align}$

を満たす確率を求めればよい。上記の確率を求めるPythonプログラムは次のとおりである。

from scipy.stats import norm

print(2 * norm.sf(x=3, loc=0, scale=1))

上記のプログラムを実行すると、次の結果が得られる。

0.0026997960632601866

10.2

測定回数を $n$ とすると10.1と同様、期待値は $E(\bar{X}) = 100$ 、分散 $V(\bar{X})$ は、

$\begin{align} V(\bar{X}) &= \frac{V(e)}{n} \\ &= \frac{0.1}{n} \end{align}$

であるから、標準化により、

$\begin{align} z = \frac{\bar{X} - 100}{\sqrt{\frac{0.1}{n}}} \end{align}$

を得る。

したがって、

$\begin{align} P \left( | \bar{X} - 100 | \lt 0.1 \right) &= P \left( \frac{\sqrt{0.1} z}{\sqrt{n}} \lt 0.1 \right) \\ &= 1 - P \left( |z| \ge \sqrt{\frac{n}{10}} \right) \\ &= 1 - 2 \cdot P \left( z \ge \sqrt{\frac{n}{10}} \right) \\ &\ge 0.9 \end{align}$

から、

$\begin{align} P \left( z \ge \sqrt{\frac{n}{10}} \right) &\lt 0.05 \end{align}$

を満たす $n$ を求めればよい。

なお、この式を満たす $z$ は、次のPythonプログラムによって求められる。

from scipy.stats import norm

print(norm.ppf(q=1-0.05, loc=0, scale=1))

結果は、

1.6448536269514722

である。したがって $n$ は、

$\begin{align} 1.6448536269514722 \lt \sqrt{\frac{n}{10}} \end{align}$

より $n \ge 27.06$ と求まるから、測定を28回以上繰り返す必要がある。

10.3

i)

正規母集団から標本を抽出するため、その標本平均の標本分布は $N( \mu, \sigma^{2} / n )$ に従う。したがってその標準化、

$\begin{align} Z = \frac{ \bar{X} - \mu }{ \sqrt{ \frac{\sigma^{2}}{n} } } \end{align}$

は標準正規分布 $N(0, 1)$ に従う。

したがって標準平均 $\bar{X}$ が3と6にある確率は、

$\begin{align} P \left( 3 \le \bar{X} \le 6 \right) &= P \left( \frac{3 - \mu }{ \sqrt{ \frac{\sigma^{2}}{n} } } \le \frac{\bar{X} - \mu }{ \sqrt{ \frac{\sigma^{2}}{n} } } \le \frac{6 - \mu }{ \sqrt{ \frac{\sigma^{2}}{n} } } \right) \\ &= P \left( -0.8165 \le z \le 1.633 \right) \\ & = P \left( 0 \le z \le 0.8165 \right) + P \left( 0 \le z \le 1.633 \right) \end{align}$

で求めることができる。これをPythonプログラムを使って求める。

from scipy.stats import norm

x1 = 0.5 - norm.sf(x=0.8165, loc=0, scale=1)
x2 = 0.5 - norm.sf(x=1.633, loc=0, scale=1)
print(x1 + x2)

上記のプログラムを実行すると、

0.7416583899007179

が得られる。

ii)

標本分散 $s^{2}$ は、 $\chi^{2}$ 分布と次のような関係がある。

$\begin{align} \chi^{2} = \frac{(n-1) s^{2}}{\sigma^{2}} \end{align}$

したがって $s^{2}$ が $a$ を超える確率が0.05であるという条件より、

$\begin{align} P \left( s^{2} \gt a \right) &= P \left( \frac{(n-1) s^{2}}{\sigma^{2}} \gt \frac{(n-1) a}{\sigma^{2}} \right) \\ &= P \left( \chi^{2} \gt \frac{(n-1) a}{\sigma^{2}} \right) \\ &= 0.05 \end{align}$

が得られるから、この式を満たす自由度 $10-1 = 9$ の $\chi^{2}$ の値を求める。

from scipy.stats import chi2

print(chi2.ppf(q=1-0.05, df=9))

上記のプログラムより、以下の結果が得られる。

16.918977604620448

したがって、

$\begin{align} \frac{(n-1) a}{\sigma^{2}} &= \frac{9a}{15} \\ &= 16.918977604620448 \end{align}$

を満たす $a$ は、28.198となる。

10.4

母集団の分散が未知であるから、スチューデントのt分布を利用する。

$\begin{align} t &= \frac{\bar{X} - \mu}{\sqrt{\frac{s^{2}}{n}}} \\ &= \frac{\bar{X} - 3}{\frac{s}{\sqrt{15}}} \end{align}$

これを用いて、

$\begin{align} P \left( \frac{\bar{X} - 3}{s} \gt a \right) &= P \left( \frac{1}{\sqrt{15}} t \gt a \right) \\ &= P\left( t \gt \sqrt{15} a \right) \\ &= 0.01 \end{align}$

この式を満たす自由度 $15 - 1 = 14$ の $t$ の値を求める。

from scipy.stats import t

print(t.ppf(q=1-0.01, df=14))

より、

2.624494067560231

と求まる。したがって、

$\begin{align} 2.624494067560231 = \sqrt{15} a \end{align}$

を満たす $a$ の値を求めればよく、その値は $a = 0.6776$ である。

10.5

2標本問題である。いずれの母集団の母平均と母分散が既知であるから、標本平均 $\bar{X}$ と $\bar{Y}$ の差は正規分布、

$\begin{align} N \left( \mu_{1} - \mu_{2}, \frac{\sigma_{1}^{2}}{m} + \frac{\sigma_{2}^{2}}{n} \right) \end{align}$

に従う。

各値を代入すると、

$\begin{align} N \left( \mu_{1} - \mu_{2}, \frac{\sigma_{1}^{2}}{m} + \frac{\sigma_{2}^{2}}{n} \right) &= N \left( -3, \frac{4}{5} \right) \end{align}$

が得られる。

10.6

母分散が既知のときに2つの標本分散の比を求めるときは、F分布

$\begin{align} F = \frac{\sigma_{2}^{2} s_{1}^{2}}{\sigma_{1}^{2} s_{2}^{2}} \end{align}$

を利用するのが便利である。

$\begin{align} P \left( \frac{s_{1}^{2}}{s_{2}^{2}} \gt c \right) &= P \left( F \gt \frac{\sigma_{2}^{2}}{\sigma_{1}^{2}} c \right) \\ &= P \left( F \gt \frac{4}{3} c \right) \\ &= 0.05 \end{align}$

この式を満たす自由度 $k_{1} = 10 - 1 = 9, k_{2} = 8 - 1 = 7$ の $F$ の値を求める。

from scipy.stats import f

print(f.ppf(q=1-0.05, dfn=9, dfd=7))

より、

3.6766746989395105

と求まる。

したがって、

$\begin{align} 3.6766746989395105 = \frac{4}{3} c \end{align}$

を満たす $c$ の値を求めればよく、その値は $c = 2.7575$ である。

10.7

i)

母平均 $\mu$ 、母分散 $\sigma^{2}$ とすると、標本平均の分布は $N(\mu, \sigma^{2} / 10)$ となる。したがって標準化、

$\begin{align} z = \frac{\bar{X} - \mu}{\sigma} \sqrt{10} \end{align}$

を考えると、

$\begin{align} P \left( | \bar{X} - \mu | \gt 0.8 \sigma \right) &= P \left( \frac{| \bar{X} - \mu |}{\sigma} \gt 0.8 \right) \\ &= P \left( |z| \gt 0.8 \right) \\ &= 2 \cdot P \left( z \gt 0.8 \right) \\ &= 0.4237 \end{align}$

なお、最後は次のPythonプログラムによって求めた。

from scipy.stats import norm

print(2 * norm.sf(x=0.8, loc=0, scale=1))

ii)

自由度が $10 - 1 = 9$ であることに注意してスチューデントのt分布

$\begin{align} t = \frac{\bar{X} - \mu}{\sqrt{\frac{s^{2}}{n}}} \end{align}$

を用いると、

$\begin{align} P \left( | \bar{X} - \mu | \gt 0.8 s \right) &= P \left( \frac{| \bar{X} - \mu |}{s} \gt 0.8 \right) \\ &= P \left( |t| \gt 0.8 \right) \\ &= 2 \cdot P \left( t \gt 0.8 \right) \\ &= 0.4443 \end{align}$

なお、最後は次のPythonプログラムによって求めた。

from scipy.stats import t

print(2 * t.sf(x=0.8, df=9))

iii)

自由度が $10 - 1 = 9$ であることに注意して $\chi^{2}$ 分布

$\begin{align} \chi^{2} = \frac{(n - 1) s^{2}}{\sigma^{2}} \end{align}$

を用いると、

$\begin{align} 1 - P \left( \frac{1}{2} \lt \frac{s^{2}}{\sigma^{2}} \lt 2 \right) &= 1 - P \left( \frac{9}{2} \lt 9 \cdot \frac{s^{2}}{\sigma^{2}} \lt 18 \right) \\ &= 1 - P \left( \frac{9}{2} \lt \chi^{2} \lt 18 \right) \\ &= 1 - \left( P \left( \chi^{2} \gt \frac{9}{2} \right) - P \left( \chi^{2} \gt 18 \right) \right) \\ &= 0.1596 \end{align}$

なお、最後は次のPythonプログラムによって求めた。

from scipy.stats import chi2

print(1- (chi2.sf(x=9/2, df=9) - chi2.sf(x=18, df=9)))

iv)

2標本問題である。

自由度が $10 - 1 = 9, 10 - 1 = 9$ のスチューデントのt分布

$\begin{align} t = \frac{ (\bar{X} - \bar{Y}) - (\mu_{1} - \mu_{2}) }{ s \sqrt{\frac{1}{m} + \frac{1}{n}} } \end{align}$

を用いると、母平均が等しい（すなわち、 $\mu_{1} = \mu_{2}$ ）ことに注意して、

$\begin{align} P \left( |\bar{X} - \bar{Y}| \ge 3s \right) &= P \left( \frac{|\bar{X} - \bar{Y} - (\mu_{1} - \mu_{2})|}{s \sqrt{\frac{1}{10} + \frac{1}{10}}} \ge \frac{3}{\sqrt{\frac{1}{10} + \frac{1}{10}}} \right) \\ &= P \left( |t| \ge 3 \sqrt{5} \right) \\ &= 0.0045718554422791826 \end{align}$

なお、最後は次のPythonプログラムによって求めた。

from scipy.stats import f
from math import sqrt

print(f.sf(x=3 * sqrt(5), dfn=9, dfd=9))

v)

2標本問題である。

自由度が $10 - 1 = 9, 10 - 1 = 9$ のF分布

$\begin{align} F &= \frac{\sigma_{2}^{2}}{\sigma_{1}^{2}} \cdot \frac{s_{1}^{2}}{s_{2}^{2}} \end{align}$

と母分散が等しい $\sigma_{1}^{2} = \sigma_{2}^{2}$ ことを用いて、

$\begin{align} P \left( \frac{s_{2}^{2}}{s_{1}^{2}} \ge 3 \right) + P \left( \frac{s_{1}^{2}}{s_{2}^{2}} \ge 3 \right) &= P \left( \frac{\sigma_{2}^{2}}{\sigma_{1}^{2}} \cdot \frac{s_{2}^{2}}{s_{1}^{2}} \ge 3 \right) + P \left( \frac{\sigma_{1}^{2}}{\sigma_{2}^{2}} \cdot \frac{s_{1}^{2}}{s_{2}^{2}} \ge 3 \right) \\ &= P \left( F \ge 3 \right) + P \left( F \le \frac{1}{3} \right) \\ &= P \left( F \ge 3 \right) + 1 - P \left( F \ge \frac{1}{3} \right) \\ &= 0.1173 \end{align}$

なお、最後は次のPythonプログラムによって求めた。

from scipy.stats import f

print(f.sf(x=3, dfn=9, dfd=9) + 1 - f.sf(x=1/3, dfn=9, dfd=9))

10.8

母集団分布が2次元の正規分布であるから、z変換を行った結果が標準正規分布になることを用いる。

$\begin{align} z = \frac{1}{2} \ln \left( \frac{1+r}{1-r} \right) \\ \eta = \frac{1}{2} \ln\left( \frac{1+\rho}{1-\rho} \right) \end{align}$

このとき、

$\begin{align} \sqrt{n-3} (z-\eta) = N(0, 1) \end{align}$

が確率0.95の範囲に収まるような $c_{1}, c_{2}$ を求めることになる。すなわち、

$\begin{align} -1.96 \le \sqrt{12} (z - \eta) \le 1.96 \end{align}$

を満たせばよい。なお標準正規分布で確率0.95（上側確率0.05）となる値は、次のPythonプログラムによって求めた。

from scipy.stats import norm

print(norm.ppf(q=(1-0.95)/2, loc=0, scale=1))

ここで、 $\rho = 15$ であるから、

$\begin{align} \eta &= \frac{1}{2} \ln \left( \frac{1+0.6}{1-0.6} \right) \\ &= 0.693147 \end{align}$

と求まる。したがって、

$\begin{align} 0.127344 \le z \le 1.25895 \end{align}$

一方で、

$\begin{align} r = \frac{e^{2z} - 1}{e^{2z} + 1} \end{align}$

であるから、

$\begin{align} 0.1267 \le r \le 0.8508 \end{align}$

が得られ、 $c_{1} = 0.1267, c_{2} = 0.8508$ と求まる。

10.9

i)

(a)

自由度nの $\chi_{(n)}^{2}$ 分布の定義は、 $z$ を標準正規分布に従う確率変数として、

$\begin{align} \chi_{(n)}^{2} = \sum_{k=1}^{n} z_{(k)}^{2} \end{align}$

である。したがって、自由度1の $\chi_{(1)}^{2}$ 分布は、

$\begin{align} \chi_{(1)}^{2} = z_{(1)}^{2} \end{align}$

となるが、 $\chi_{1}^{2}$ 分布は上側確率のみを与えるので、このための補正を行うと、

$\begin{align} \chi_{\alpha (1)}^{2} = z_{\frac{\alpha}{2} (1)}^{2} \end{align}$

が得られる。

(b)

本書より、自由度 $k$ のスチューデントのt分布の2乗 $t_{(k)}^{2}$ は、F分布 $F_{(1, k)}$ に従う。ただしt分布は両側確率を与えるものである一方、F分布は上側確率のみを与えるから、(a)の議論と同様に補正を行って、

$\begin{align} t_{\frac{\alpha}{2} (k)}^{2} = F_{\alpha (1, k)} \end{align}$

が得られる。

(c)

$k = 120$ であるから、 $k$ が十分に大きい値であると考えることができる。つまり、標本数を $\infty$ として考えてよく、このとき標本分散 $s^{2}$ は母分散 $\sigma^{2}$ に近づく。したがって、

$\begin{align} t_{\alpha (k)} &= \frac{z_{\alpha (k \ge 120)}}{\sqrt{ \frac{s^{2}}{\sigma^{2}} }} \\ &\rightarrow \frac{z_{\alpha (k \ge 120)}}{\sqrt{ \frac{\sigma^{2}}{\sigma^{2}} }} \\ &= z_{\alpha (k \ge 120)} \end{align}$

ii)

i)の各値を求めるPythonプログラムを次に示す。

from scipy.stats import norm, chi2, t, f

alpha = 0.05
k = 120

print("=== (a) ===")
print("Z2: {}".format(norm.ppf(q=1-alpha/2, loc=0, scale=1)**2))
print("chi2: {}".format(chi2.ppf(q=1-alpha, df=1)))

print("=== (b) ===")
print("t2: {}".format(t.ppf(q=1-alpha/2, df=k)**2))
print("F: {}".format(f.ppf(q=1-alpha, dfn=1, dfd=k)))

print("=== (c) ===")
print("t: {}".format(t.ppf(q=1-alpha, df=k)))
print("Z: {}".format(norm.ppf(q=1-alpha, loc=0, scale=1)))

上記のプログラムを実行すると、次のような結果が得られる。

=== (a) ===
Z2: 3.8414588206941254
chi2: 3.841458820694124
=== (b) ===
t2: 3.9201244088524523
F: 3.9201244089699054
=== (c) ===
t: 1.6576508993473795
Z: 1.6448536269514722

2020-09-24

【PyTorch】C++ APIを使ってResNet50を実装する

C++ Python 機械学習

はじめに

Experimentalな機能として提供されてきたPyTorchのC++ APIが、バージョン1.5よりStableになった¹。

C++ APIを利用することで、Pythonを利用できない環境や速度を求められる環境においてもPyTorchを使って深層学習のモデルを記述して学習できるようになる。 PyTorchのC++ APIは、既存のPython APIと可能な限り似た記述になるようにAPIが設計されている。このため、Python APIを使って深層学習のモデルをこれまで書いてきた人であれば、C++ APIを使いこなすのはそこまでハードルが高くないと思われる。ただしC++ APIは現在も開発中であり、Python APIに対応するAPIが存在しないことが多いため注意が必要である。例えば、Forward/BackwardのHook登録や一部のOptimizer は未実装となっている。

本記事では、ResNet50をC++ APIで実装するための解説を行う。いまさらResNet50かとツッコミどころはあるが、著名で実装しやすいこともあり、あえてこのモデルを選んだ。 Python APIとC++ APIを比較しやすくするために、Python APIとC++ APIの両方の実装を比較しながら説明していく。本記事を読むことにより、Python APIとC++ APIがよく似た仕様となっていることを理解できるだろう。

なお、PyTorchのC++ APIを使ってプログラムを書くための準備は、PyTorch公式のチュートリアルも参考になる。チュートリアルでは、入力層と出力層から構成される簡単なニューラルネットワークの作成から始まり、DCGANの実装までを丁寧に説明している。必要に応じてこちらも参照されたい。

本記事で紹介するソースコード一式は、GitHub にて公開している。

注意点

本記事で紹介するResNet50の実装は、以下の点で手を抜いた実装となっている。ソースコードを参考にする際は注意が必要である。

ImageNetの代わりにMNISTをデータセットとして使用
- 入力画像サイズ：[224, 224] → [28, 28]
- 入力画像Channel数：3(RGB) → 1(白黒)
- クラス数：1000 → 10
Learning Rate固定
- 本来であればepoch数にしたがってLearning Rateを減衰させる必要があるが、該当する機能がC++ APIとして提供されていないため実装をサボった

前提知識

本記事は、以下の読者を対象としている。

PyTorchのPython APIを使い、ニューラルネットワークを構築したことがある
C++を使ってプログラムを書いたことがある

前準備

PyTorchのC++ APIを利用するためには、libtorchと呼ばれるPyTorchのライブラリをダウンロードする必要がある。 CUDAなしの環境とCUDAありの環境とでライブラリが異なるため注意が必要である。

CUDAなし（CPU）

wget https://download.pytorch.org/libtorch/nightly/cpu/libtorch-shared-with-deps-latest.zip
unzip libtorch-shared-with-deps-latest.zip

CUDAあり（CPU+GPU）

wget https://download.pytorch.org/libtorch/cu102/libtorch-shared-with-deps-1.5.1.zip
unzip libtorch-shared-with-deps-1.5.1.zip

また、本プログラムでは前処理のためにOpenCVを利用しているため、こちらの記事を参考にしてインストールする。

プログラムのビルドにはCMakeを使うため、次のようにCMakeLists.txtを作成する。

対応するGitHub上のコード

cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
project(pytorch-cpp-example)

find_package(Torch REQUIRED)
find_package(OpenCV REQUIRED)

include_directories(
    ${PROJECT_SOURCE_DIR}
    ${OPENCV_INCLUDE_DIRS}
)
add_executable(train model.cpp train.cpp)
add_executable(predict model.cpp predict.cpp)
target_link_libraries(train ${TORCH_LIBRARIES})
target_link_libraries(predict ${TORCH_LIBRARIES} ${OpenCV_LIBRARIES})
set_property(TARGET train PROPERTY CXX_STANDARD 14)
set_property(TARGET predict PROPERTY CXX_STANDARD 14)

ResNet50のプログラムを書く

ResNet50の論文を参考に、プログラムを作成する。

1. Residual Blockの実装

ResNetは残差ブロック（Residual Block）を導入することにより、層の深いネットワークにおける勾配損失問題を解消したことで有名なニューラルネットワークモデルである。 Residual Blockの実装は、Python APIで記述すると次のようになる。

対応するGitHub上のコード

class ResidualBlock(torch.nn.Module):
    def __init__(self, in_channels, out_channels, stride=1):
        super(ResidualBlock, self).__init__()
        width = out_channels // 4

        # (1.a-1)
        self.conv1 = torch.nn.Conv2d(in_channels, width, kernel_size=(1, 1),
                                     stride=1, bias=False)
        self.bn1 = torch.nn.BatchNorm2d(width)
        self.relu1 = torch.nn.ReLU(inplace=True)

        self.conv2 = torch.nn.Conv2d(width, width, kernel_size=(3, 3),
                                     stride=stride, padding=1, groups=1,
                                     bias=False, dilation=1)
        self.bn2 = torch.nn.BatchNorm2d(width)
        self.relu2 = torch.nn.ReLU(inplace=True)

        self.conv3 = torch.nn.Conv2d(width, out_channels, kernel_size=(1, 1),
                                     stride=1, padding=0, bias=False)
        self.bn3 = torch.nn.BatchNorm2d(out_channels)

        # (1.a-2)
        def shortcut(in_, out):
            if in_ != out:
                return torch.nn.Sequential(
                    torch.nn.Conv2d(in_, out, kernel_size=(1, 1),
                                    stride=stride, padding=0, bias=False),
                    torch.nn.BatchNorm2d(out),
                )
            else:
                return lambda x: x
        self.shortcut = shortcut(in_channels, out_channels)

        self.relu3 = torch.nn.ReLU(inplace=True)

    def forward(self, x):
        identity = x

        out = self.conv1(x)
        out = self.bn1(out)
        out = self.relu1(out)

        out = self.conv2(out)
        out = self.bn2(out)
        out = self.relu2(out)

        out = self.conv3(out)
        out = self.bn3(out)

        shortcut = self.shortcut(x)

        out = self.relu3(out + shortcut)

        return out

torch.nn.Module を継承したResidualBlockクラスの中で、torch.nn モジュールで提供される各層に対応する演算APIを利用してResidual Blockを定義する (1.a-1)。入力channelと出力channelが異なる場合に、1x1のConvolutionによってUp-sampling (1.a-2) を行う shortcut 関数に注意しよう。

つづいて、C++ APIを用いてResidual Blockを実装した場合のソースコードを示す。

対応するGitHub上のコード

Python APIで記述した場合とC++ APIで記述した場合の違いは、おもに2つある。

1つ目は、Conv2dのStrideなどの設定について、PythonではAPIの引数に渡すことで実現しているのに対し、C++では Conv2dOptions 構造体などをAPIの引数に渡すことで実現している点である (1.b-1) 。

2つ目は、ResidualBlock 構造体のコンストラクタの最後の処理で register_module 関数を呼び出している点である (1.b-3) 。 register_module 関数は学習する層を登録するための関数で、この登録処理を忘れてしまうと学習時の誤差逆伝搬の処理が行えなくなってしまう。

ResidualBlockImpl::ResidualBlockImpl(int in_channels, int out_channels,
                                     int stride) {
    int width = out_channels / 4;

    // (1.b-1)
    conv1 = Conv2d(Conv2dOptions(in_channels, width, {1, 1})
                   .stride(1).bias(false));
    bn1 = BatchNorm2d(BatchNormOptions(width));
    relu1 = ReLU(ReLUOptions().inplace(true));

    conv2 = Conv2d(Conv2dOptions(width, width, {3, 3})
                   .stride(stride).padding(1).groups(1)
                   .bias(false).dilation(1));
    bn2 = BatchNorm2d(BatchNormOptions(width));
    relu2 = ReLU(ReLUOptions().inplace(true));

    conv3 = Conv2d(Conv2dOptions(width, out_channels, {1, 1})
                   .stride(1).padding(0).bias(false));
    bn3 = BatchNorm2d(BatchNormOptions(out_channels));

    // (1.b-2)
    Sequential shortcut(
        Conv2d(Conv2dOptions(in_channels, out_channels, {1, 1})
               .stride(stride).padding(0).bias(false)),
        BatchNorm2d(BatchNormOptions(out_channels))
    );
    this->shortcut = shortcut;
    relu3 = ReLU(ReLUOptions().inplace(true));

    this->in_channels = in_channels;
    this->out_channels = out_channels;

    // // (1.b-3)
    register_module("conv1", conv1);
    register_module("bn1", bn1);
    register_module("relu1", relu1);
    register_module("conv2", conv2);
    register_module("bn2", bn2);
    register_module("relu2", relu2);
    register_module("conv3", conv3);
    register_module("bn3", bn3);
    register_module("shortcut", shortcut);
    register_module("relu3", relu3);
}

torch::Tensor ResidualBlockImpl::forward(torch::Tensor input) {
    torch::Tensor out;
    torch::Tensor tmp;

    out = conv1->forward(input);
    out = bn1->forward(out);
    out = relu1->forward(out);

    out = conv2->forward(out);
    out = bn2->forward(out);
    out = relu2->forward(out);

    out = conv3->forward(out);
    out = bn3->forward(out);

    if (in_channels != out_channels) {
        tmp = shortcut->forward(input);
    } else {
        tmp = input;
    }
    out = relu3->forward(out + tmp);

    return out;
}

2. ResNet50の実装

Residual Blockを組み合わせて、ResNet50を実装する。最初にPython APIによるResNet50の実装を示す。

対応するGitHub上のコード

学習データとしてMNISTを利用するため、ResNet50モデルの入力データのサイズや出力クラス数に気をつけよう。 ResNet50モデルでは、最初の torch.nn.Conv2d (2.a-1) (2.b-1) と最後の torch.nn.Linear (2.a-2) (2.b-2) の引数に気をつける。

class ResNet50(torch.nn.Module):
    def __init__(self):
        super(ResNet50, self).__init__()
        # (2.a-1)
        self.conv1 = torch.nn.Conv2d(1, 64, kernel_size=(7, 7),
                                     stride=2, padding=3, bias=False)
        self.bn1 = torch.nn.BatchNorm2d(64)
        self.relu = torch.nn.ReLU(inplace=True)
        self.maxpool = torch.nn.MaxPool2d(kernel_size=(3, 3), stride=2,
                                          padding=1)

        self.layer1 = torch.nn.Sequential(
            ResidualBlock(64, 256),
            ResidualBlock(256, 256),
            ResidualBlock(256, 256),
        )

        self.layer2 = torch.nn.Sequential(
            ResidualBlock(256, 512, stride=2),
            ResidualBlock(512, 512),
            ResidualBlock(512, 512),
            ResidualBlock(512, 512),
        )

        self.layer3 = torch.nn.Sequential(
            ResidualBlock(512, 1024, stride=2),
            ResidualBlock(1024, 1024),
            ResidualBlock(1024, 1024),
            ResidualBlock(1024, 1024),
            ResidualBlock(1024, 1024),
            ResidualBlock(1024, 1024),
        )

        self.layer4 = torch.nn.Sequential(
            ResidualBlock(1024, 2048, stride=2),
            ResidualBlock(2048, 2048),
            ResidualBlock(2048, 2048),
        )

        self.avgpool = torch.nn.AdaptiveAvgPool2d((1, 1))
        self.flatten = torch.nn.Flatten(1)
        # (2.a-2)
        self.fc = torch.nn.Linear(2048, 10)

    def forward(self, x):
        x = self.conv1(x)
        x = self.bn1(x)
        x = self.relu(x)
        x = self.maxpool(x)

        x = self.layer1(x)
        x = self.layer2(x)
        x = self.layer3(x)
        x = self.layer4(x)

        x = self.avgpool(x)
        x = self.flatten(x)
        x = self.fc(x)

        return x

続いて、C++ APIによる実装を示す。

対応するGitHub上のコード

ResNet50Impl::ResNet50Impl() {
    // (2.b-1)
    conv1 = Conv2d(Conv2dOptions(1, 64, {7, 7})
                   .stride(2).padding(3).bias(false));
    bn1 = BatchNorm2d(BatchNormOptions(64));
    relu = ReLU(ReLUOptions().inplace(true));
    maxpool = MaxPool2d(MaxPoolOptions<2>({3, 3}).stride(2).padding(1));

    Sequential layer1(
        ResidualBlock(64, 256),
        ResidualBlock(256, 256),
        ResidualBlock(256, 256)
    );
    this->layer1 = layer1;

    Sequential layer2(
        ResidualBlock(256, 512, 2),
        ResidualBlock(512, 512),
        ResidualBlock(512, 512),
        ResidualBlock(512, 512)
    );
    this->layer2 = layer2;

    Sequential layer3(
        ResidualBlock(512, 1024, 2),
        ResidualBlock(1024, 1024),
        ResidualBlock(1024, 1024),
        ResidualBlock(1024, 1024),
        ResidualBlock(1024, 1024),
        ResidualBlock(1024, 1024)
    );
    this->layer3 = layer3;

    Sequential layer4(
        ResidualBlock(1024, 2048, 2),
        ResidualBlock(2048, 2048),
        ResidualBlock(2048, 2048)
    );
    this->layer4 = layer4;

    avgpool = AdaptiveAvgPool2d(AdaptiveAvgPool2dOptions({1, 1}));
    flatten = Flatten(FlattenOptions().start_dim(1));

    // (2.b-2)
    fc = Linear(2048, 10);

    register_module("conv1", conv1);
    register_module("bn1", bn1);
    register_module("relu", relu);
    register_module("maxpool", maxpool);
    register_module("layer1", this->layer1);
    register_module("layer2", this->layer2);
    register_module("layer3", this->layer3);
    register_module("layer4", this->layer4);
    register_module("avgpool", avgpool);
    register_module("flatten", flatten);
    register_module("fc", fc);
}

torch::Tensor ResNet50Impl::forward(torch::Tensor input) {
    torch::Tensor out;

    out = conv1->forward(input);
    out = bn1->forward(out);
    out = relu->forward(out);
    out = maxpool->forward(out);

    out = layer1->forward(out);
    out = layer2->forward(out);
    out = layer3->forward(out);
    out = layer4->forward(out);

    out = avgpool->forward(out);
    out = flatten->forward(out);
    out = fc->forward(out);

    return out;
}

ResidualBlock の実装と同じような要領で実装していけばよく、C++ APIでの実装に関して特筆すべきことはない。

ここで、実装が ResidualBlockImpl 構造体であるのに対して ResidualBlock としてアクセスできていることが気になるかもしれない。これはPyTorchをC++で書いた時の流儀で、実装を XXXXImpl という構造体名として TORCH_MODULE マクロでその構造体を登録することで、XXXX 構造体としてアクセスできるようになる。

3. 学習処理の実装

ResNet50を学習するための処理を実装する。実装は次の順番で行っていく。

i) コマンドライン引数の解析、乱数固定
ii) 学習を実行するデバイス（CPUかGPUか）の決定
iii) ResNet50モデルとOptimizerの定義
iv) MNISTデータセットの読み込み
v) 学習
vi) 評価
vii) 学習済みモデルの保存

i) コマンドライン引数の解析、乱数固定

最初にコマンドラインの解析と、再現性確保のための乱数固定を行う。乱数固定に関しては、こちらの記事を参照のこと。

対応するGitHub上のコード（Python API）

対応するGitHub上のコード（C++ API）

ii) 学習を実行するデバイス（CPUかGPUか）の決定

学習を実行するデバイスを決定する処理をPython APIで実装すると、次のようになる。

対応するGitHub上のコード

    # (3.a-1)
    # Parse arguments.
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "-m", dest="saved_model_path", type=str,
        help="Path to saved model", required=True)
    args = parser.parse_args()

    fix_randomness(1)

torch.cuda.is_available 関数を使い、CUDAが利用可能な場合はGPU、利用できない場合はCPUで学習を実行するようにデバイスを決定する (3.a-1)。

同様の処理をC++ APIで実装した場合は次のようになる。

対応するGitHub上のコード

    // (3.b-1)
    // Create device.
    torch::DeviceType device_type;
    if (torch::cuda::is_available()) {
        std::cout << "Train on GPU." << std::endl;
        device_type = torch::kCUDA;
    } else {
        std::cout << "Train on CPU." << std::endl;
        device_type = torch::kCPU;
    }
    torch::Device device(device_type);

対応するC++ APIに置き換えればよいだけなので、特別に注意する点はない。

iii) ResNet50モデルとOptimizerの定義

ResNet50モデルとOptimizerの定義は、Python APIでは次のようになる。

対応するGitHub上のコード

    # Build model.
    model = ResNet50()
    model.to(device)
    summary(model, (1, 28, 28))
    # (3.a-2)
    optimizer = torch.optim.Adam(model.parameters(), lr=0.01)

本記事の最初にも書いたが、OptimizerのLearning Rateはepochが進むにつれて減衰させるべきところを決め打ちで値を指定している (3.a-2) (3.b-2)。 Learning Rateを逐次変更したい場合は、Pythonでは torch.optim.lr_scheduler モジュールを利用できるが、C++ APIに同様のAPIが存在しないため自力で実装する必要がある。

C++ APIにおける実装を次に示す。

対応するGitHub上のコード

    // Build model.
    ResNet50 model;
    model->to(device);
    // (3.b-2)
    torch::optim::Adam optimizer(
        model->parameters(), torch::optim::AdamOptions(0.01));

ここに関しても特別に注意する点はないだろう。

iv) MNISTデータセットの読み込み

次に、MNISTのデータセットを読み込む。

対応するGitHub上のコード

    # Load dataset.
    train_loader = torch.utils.data.DataLoader(
        torchvision.datasets.MNIST(DATA_ROOT, train=True, download=True,
                                   transform=torchvision.transforms.Compose([
                                       torchvision.transforms.ToTensor(),
                                       torchvision.transforms.Normalize((0.1307,), (0.3081,)),
                                   ])
        ),
        batch_size=TRAIN_BATCH_SIZE,
    )
    test_loader = torch.utils.data.DataLoader(
        torchvision.datasets.MNIST(DATA_ROOT, train=False, download=True,
                                   transform=torchvision.transforms.Compose([
                                       torchvision.transforms.ToTensor(),
                                       torchvision.transforms.Normalize((0.1307,), (0.3081,)),
                                   ])
        ),
        batch_size=TEST_BATCH_SIZE,
    )

Python APIの場合は、torch.utils.data.DataLoader クラスを利用し、引数に torch.utils.data.Dataset から継承されたMNISTのデータセットを読み込むための便利クラス torchvision.datasets.MNIST のインスタンスを渡す。 torchvision.datasets.MNIST の引数 download に True を指定することで、DATA_ROOT にデータが存在しない場合は、MNISTのデータセットをダウンロードしてから読み込んでくれる。引数 shuffle が指定されていないことから、引数 shuffle はデフォルトで False となり、シーケンシャルにデータセットからサンプリングするようになる。

データ読み込み時は、transform 引数に入力データに対して適用する変形処理を渡す。ここでは、次のような変形処理を行った。

テンソルの各要素の値を[0, 255]から[0.0, 1.0]に正規化し、データレイアウトを(H, W, C)から(C, H, W)へ変換
各要素を平均0.1307、標準偏差0.3081に標準化

続いてC++ APIを用いた実装を示す。

対応するGitHub上のコード

    // Load dataset.
    auto train_dataset = torch::data::datasets::MNIST(kDataRoot)
        .map(torch::data::transforms::Normalize<>(0.1307, 0.3081))
        .map(torch::data::transforms::Stack<>());
    auto train_loader =
        torch::data::make_data_loader<torch::data::samplers::SequentialSampler>(
        std::move(train_dataset), kTrainBatchSize);

    auto test_dataset = torch::data::datasets::MNIST(
        kDataRoot, torch::data::datasets::MNIST::Mode::kTest)
        .map(torch::data::transforms::Normalize<>(0.1307, 0.3081))
        .map(torch::data::transforms::Stack<>());
    auto test_loader =
        torch::data::make_data_loader<torch::data::samplers::SequentialSampler>(
        std::move(test_dataset), kTestBatchSize);

Python APIの torch.utils.data.DataLoader に対応する DataLoader は、テンプレート関数 torch::data::make_data_loader を使って生成できる。テンプレート引数には、バッチを作るためのデータのサンプリング方法を指定可能である。今回はPython APIと同様にシーケンシャルにサンプリングするために、torch::data::samplers::SequentialSampler を利用した。

v) 学習

学習のループ処理を実装する。最初に、DataLoaderから学習用のデータを使ってResNet50のモデルを訓練するコードを示す。

対応するGitHub上のコード

    # Train loop.
    for epoch in range(NUMBER_OF_EPOCHS):
        print("Epoch {}:".format(epoch))

        # Train.
        print("Start train.")
        model.train()    # (3.a-3)
        for batch_idx, (data, target) in enumerate(train_loader):
            optimizer.zero_grad()    # (3.a-4)
            data = data.to(device)    # (3.a-5)
            target = target.to(device)

            output = model(data)    # (3.a-6)

            # (3.a-7)
            prob = F.log_softmax(output, dim=1)
            loss = F.nll_loss(prob, target)
            loss.backward()
            optimizer.step()

            if batch_idx % LOG_INTERVAL == 0:
                print("Batch: {}, Loss: {}".format(batch_idx, loss.item()))

ResNet50のモデル model をtraining modeに変更したあと (3.a-3) (3.b-3)、DataLoaderから学習用のデータを1バッチずつ取得する。 DataLoaderから読み込んだデータは、data に説明変数、target に目的変数がテンソル（torch.Tensor）として保存されている。 batch_idx は、一定間隔でloss値をログとして使用するときに利用する。

バッチを読み込んだあと、Backward Propagationの初期値をリセットするために optimizer.zero_grad を呼びだす (3.a-4) (3.b-4)。そして、dataとtargetのテンソルデータを学習を実行するデバイスのメモリに移動させ (3.a-5) (3.b-5)、ResNet50モデルのForward Propagationを実行する (3.a-6) (3.b-6)。最後に、Forward Propagationの結果を使ってロス値を計算（Log Softmax + NLL-Loss ＝ CrossEntropyLoss）したあと、Backward Propagationを行って各パラメータの勾配を求め、求めた勾配を使ってパラメータを更新する (3.a-7) (3.b-7)。

この処理をC++ APIで実装すると、次のようになる。

対応するGitHub上のコード

    // Train loop.
    for (size_t epoch = 0; epoch < kNumberOfEpochs; ++epoch) {
        std::cout << "Epoch " << epoch << ":" << std::endl;

        // Train.
        std::cout << "Start train." << std::endl;
        size_t batch_idx = 0;
        model->train();    // (3.b-3)
        for (auto& batch : *train_loader) {
            optimizer.zero_grad();    // (3.b-4)
            auto data = batch.data.to(device);    // (3.b-5)
            auto target = batch.target.to(device);

            auto output = model->forward(data);    // (3.b-6)

            // (3.b-7)
            auto prob = F::log_softmax(output, 1);
            auto loss = F::nll_loss(prob, target);
            AT_ASSERT(!std::isnan(loss.template item<float>()));
            loss.backward();
            optimizer.step();

            if ((batch_idx % kLogInterval) == 0) {
                std::cout << "Batch: " << batch_idx << ", Loss: "
                          << loss.template item<float>() << std::endl;
            }
            batch_idx++;
        }

vi) 評価

学習結果を評価する処理を実装する。

対応するGitHub上のコード

        # Evaluate.
        print("Start eval.")
        model.eval()    # (3.a-8)
        test_loss = 0
        correct = 0.0
        total = 0
        with torch.no_grad():    # (3.a-9)
            for data, target in test_loader:
                data = data.to(device)
                target = target.to(device)

                output = model(data)    # (3.a-10)

                prob = F.log_softmax(output, dim=1)
                test_loss += F.nll_loss(prob, target, reduction="sum").item()    # (3.a-11)
                pred = output.argmax(1, keepdim=True)    # (3.a-12)
                correct += pred.eq(target.view_as(pred)).sum().item()
                total += TEST_BATCH_SIZE

        print("Average loss: {}, Accuracy: {}"
              .format(test_loss / loss, correct / total))

ResNet50のモデル model をeval modeに変更したあと (3.a-8) (3.b.8)、DataLoaderから評価用のデータを1バッチずつ取得する。

評価時は勾配計算が行われないように、torch.no_grad コンテキスト内で評価のためのコードを実行させる (3.a-9) (3.b-9)。 torch.no_grad コンテキストにより、勾配のためにテンソルデータを保持する必要がなくなり、メモリの消費量を抑えることができる。

評価用の画像データに対してクラスを予測するため、ResNet50のForward Propagationを行う (3.a-10) (3.a-10)。その結果を利用し、loss値の計算 (3.a-11) (3.b-11) と評価用データのクラス予測 (3.a-12) (3.b-12) を行う。

上記の処理をC++ APIで実装すると、次のようになる。

対応するGitHub上のコード

        // Evaluate.
        std::cout << "Start eval." << std::endl;
        torch::NoGradGuard no_grad;    // (3.b-9)
        model->eval();    // (3.b-8)
        double test_loss = 0.0;
        size_t correct = 0;
        size_t total = 0;
        for (auto& batch : *test_loader) {
            auto data = batch.data.to(device);
            auto target = batch.target.to(device);
            auto output = model->forward(data);    // (3.b-10)

            auto prob = F::log_softmax(output, 1);
            test_loss += F::nll_loss(
                prob, target,
                F::NLLLossFuncOptions().reduction(torch::kSum)).template item<double>();    // (3.b-11)
            auto pred = output.argmax(1, true);    // (3.b-12)
            correct += pred.eq(target.view_as(pred)).sum().template item<int64_t>();
            total += kTestBatchSize;
        }

        std::cout << "Average loss: " << test_loss / total
                  << ", Accuracy: " << static_cast<double>(correct) / total
                  << std::endl;

torch.no_grad コンテキストに相当する torch::NoGradGuard の変数 no_grad が、変数の生存期間中有効になる点に注意が必要である。

vii) 学習済みモデルの保存

最後に学習済みのモデルを保存する。

対応するGitHub上のコード

    # Save trained model.
    os.makedirs(args.saved_model_path, exist_ok=True)
    model_path = "{}/{}".format(args.saved_model_path, SAVED_MODEL_NAME)
    torch.save(model.state_dict(), model_path)    # (3.a-13)
    print("Saved model to '{}'".format(model_path))

Python APIでは、torch.save の引数にResNet50モデルの state_dict を渡すことで、Optimizerの状態も含めてファイルに保存できる。

続いて、C++ APIを用いた実装を示す。

対応するGitHub上のコード

    // Save trained model.
    std::string model_path = args.saved_model_path + "/"
        + kSavedModelNamePrefix + "_model.pth";
    std::string optimizer_path = args.saved_model_path + "/"
        + kSavedModelNamePrefix + "_optimizer.pth";
    struct stat buf;
    if (stat(args.saved_model_path.c_str(), &buf)) {
        int rc;
        rc = mkdir(args.saved_model_path.c_str(), 0755);
        if (rc < 0) {
            std::cout << "Error: Failed to create diretory '"
                      << args.saved_model_path  << "' ("
                      << errno << ": " << strerror(errno) << ")" << std::endl;
            return 1;
        }
    }
    torch::save(model, model_path);
    torch::save(optimizer, optimizer_path);
    std::cout << "Saved model." << std::endl;
    std::cout << "  Model: " << model_path << std::endl;
    std::cout << "  Optimizer: " << model_path << std::endl;

ディレクトリ作成処理がC++ではやや煩雑になっているが、学習済みモデルの保存は torch::save を呼び出すだけで完了する。なお、C++ APIでは state_dict が提供されていないため、Optimizerの状態を別途 torch::save 関数に渡して保存しなければならないことに注意しよう。

以上で学習処理の実装が完了した。

ビルド

作成したソースコードをビルドする。

cd cpp
mkdir build
cd build
cmake -DCMAKE_PREFIX_PATH=/path/to/libtorch ..
cmake --build . --config Release

ビルドが完了したら、同ディレクトリに train と predict の実行プログラムが作られていることを確認しよう。

MNISTデータセットのダウンロード

MNISTデータセットをダウンロードし、ビルド時に作成した実行プログラムと同じディレクトリにMNISTデータセットを配置する。

mkdir mnist
wget http://yann.lecun.com/exdb/mnist/train-images-idx3-ubyte.gz
wget http://yann.lecun.com/exdb/mnist/train-labels-idx1-ubyte.gz
wget http://yann.lecun.com/exdb/mnist/t10k-images-idx3-ubyte.gz
wget http://yann.lecun.com/exdb/mnist/t10k-labels-idx1-ubyte.gz
gunzip *.gz

結果

ResNet50を学習するために次のコマンドを実行する。オプション -m は、学習済みモデルを保存するディレクトリを指定するものである。ディレクトリが存在しない場合は、自動的に指定されたディレクトリが作成される。

./train -m saved_model

プログラムの実行が完了すると、学習済みモデルがオプションで指定したディレクトリに保存される。

続いて、学習済みのモデルを使って、手書きの数字の画像が正しく推論できるか確認する。今回利用した画像は、新たに作成した手書きの数字の画像であり実行プログラム predict からの相対パス ../../data/digit.png に配置されているものとする。

次のコマンドを実行することで、saved_model に保存された学習済みモデルを利用して画像データ ../../data/digit.png について推論する。

./predict -i ../../data/digit.png -m saved_model

上記のコマンドを実行すると次のような出力が得られ、画像データが正しく推論できていることがわかる。

Predict: 7

なお、学習済みのモデルを使って推論するためのソースコードもGitHub上で公開している。 OpenCVを使って画像データの読み込み＆前処理を行っていることを除いて特に難しいところは無いため、具体的なソースコードの説明に関しては省略する。

Pythonによる実装

C++による実装

おわりに

CNNであるResNet50をPyTorchのC++ APIを使って実行し、MNISTデータセットを使って実際に学習させた。また学習したモデルを使用し、手書きの数字の画像が正しく推論できていることを確認した。

Python APIによる記述とC++ APIはどちらも似たようなAPI仕様となっているため、これまでPython APIを使ってPyTorchを使っていた人がC++ APIに移行するのはそれほど大変ではないだろう。一方で、C++ APIと同等の機能を持ったAPIが存在しないことも多く、C++ APIが安定版として提供され始めたとはいえ、まだまだ機能不足感は否めない。このため、ドキュメントやPyTorchのIssueを確認しながら、使いたいC++ APIがサポートされているか否かを事前に確認することが必要になるだろう。

C++ APIの開発は2019/9頃に GitHub にてContributorの募集があり、私もこの開発に参加してBatchNormなど12個のAPIを実装した。↩

2020-09-19

【統計学入門（東京大学出版会）】第9章練習問題解答

統計学

東京大学出版会から出版されている統計学入門（基礎統計学Ⅰ）について第9章の練習問題の解答を書いていく。

本章以外の解答
9.1
9.2
9.3
9.4
9.5
9.6
9.7
- i)
- ii)
9.8
- i)
- ii)
- iii)
9.9

本章以外の解答

本章以外の練習問題の解答は別の記事で公開している。
必要に応じて参照されたい。

9.1

$\bar{X}$ を標本平均とすると、その期待値 $E(\bar{X})$ は母集団の平均 $\mu$ になる。また、標本分散 $s^{2}$ の期待値 $E(s^{2})$ は母集団の分散 $\sigma^{2}$ になる。

この関係を用いることで、標本分布は母集団の分布を知るための手掛かりになると言える。

9.2

標本平均

$\begin{align} \bar{X} &= \frac{1}{n} \sum_{i=1}^{n} X_{i} \\ &= \frac{1.22+1.24+1.25+1.19+1.17+1.18}{6} \\ &\fallingdotseq 1.208333 \end{align}$

標本分散

$\begin{align} s^{2} &= \frac{1}{n-1} \sum_{i=1}^{n} (X_{i} - \bar{X})^{2} \\ &= \frac{1}{6-1} \left( (1.22 - 1.21)^{2} + (1.24 - 1.21)^{2} + (1.25 - 1.21)^{2} + (1.19 - 1.21)^{2} + (1.17 - 1.21)^{2} + (1.18 - 1.21)^{2} \right) \\ &\fallingdotseq 0.0011 \end{align}$

9.3

表に示されている系列に対して、 $s^{2}$ と $S^{2}$ をそれぞれ計算するPythonプログラムを次に示す。

import numpy as np

a = np.array([
    [0.3104913, 0.3304700, 0.0324358, 0.8283330, 0.1727581, 0.6306326, 0.7210595, 0.2451280, 0.7243750, 0.8197760],
    [0.2753351, 0.4359388, 0.7160295, 0.7775517, 0.3251019, 0.1736023, 0.0921532, 0.1318467, 0.0642188, 0.8002448],
    [0.3368585, 0.2513685, 0.2697405, 0.1164189, 0.3085003, 0.2234060, 0.9427391, 0.5800890, 0.7194922, 0.8344245],
    [0.4086511, 0.8016156, 0.3221239, 0.8498936, 0.4362011, 0.8559286, 0.9982964, 0.5540422, 0.3757575, 0.1312537],
    [0.4449823, 0.1457471, 0.9303545, 0.1033269, 0.4415264, 0.5430776, 0.8274743, 0.3946336, 0.8696082, 0.6028266],
])
ave = np.mean(a, axis=1, keepdims=True)
s2 = np.sum((a - ave)**2, axis=1) / 9
S2 = np.sum((a - ave)**2, axis=1) / 10

print("sigma^2: {}".format(1.0 / 12))
print(f"s^2: {s2}")
print(f"S^2: {S2}")

プログラムを実行すると、次の結果が得られる。

sigma^2: 0.08333333333333333
s^2: [0.08647192 0.08349042 0.08298067 0.08143628 0.08140686]
S^2: [0.07782473 0.07514138 0.0746826  0.07329265 0.07326618]

いずれの系列についても、 $s^{2}$ のほうが $\sigma^{2}$ に近い値であることがわかる。

9.4

$E(X_{i}^{2}) = V(X_{i}) + E(X_{i})^{2} = \sigma^{2} + \mu^{2} \\ \displaystyle E(\bar{X}^{2}) = V(\bar{X}) + E(\bar{X})^{2} = \frac{\sigma^{2}}{n} + \mu^{2} \\ \displaystyle \sum_{i=1}^{n} E(X_{i} \bar{X}) = E \left( \sum_{i=1}^{n} X_{i} \bar{X} \right) = E \left( n \bar{X} \sum_{i=1}^{n} \frac{X_{i}}{n} \right) = nE(\bar{X}^{2}) = n\mu^{2} + \sigma^{2}$

を利用すると、

$\begin{align} E(s^{2}) &= \frac{1}{n-1} \sum_{i=1}^{n} E(X_{i} - \bar{X})^{2} \\ &= \frac{1}{n-1} \sum_{i=1}^{n} \left( E(X_{i}^{2}) - 2E(X_{i}\bar{X}) + E(\bar{X}^{2}) \right) \\ &= \frac{1}{n-1} \left( \sum_{i=1}^{n} \left( \sigma^{2} + \mu^{2} + \frac{\sigma^{2}}{n} + \mu^{2} \right) - 2 (n\mu^{2} + \sigma^{2}) \right) \\ &= \sigma^{2} \end{align}$

となる。

9.5

$n=3$ のとき

A地点より信号を $n$ 回送信し、B地点で正しく受信する回数はそれぞれ次の確率で発生しうる。

正しく受信する回数	確率
0	$1 \cdot (0.1)^{3}$
1	${}_{3}C_{1} \cdot 0.9 \cdot (0.1)^{2}$
2	${}_{3}C_{2} \cdot (0.9)^{2} \cdot 0.1$
3	$1 \cdot (0.9)^{3}$

このうち、2回以上正しく受信できれば信号は正しく伝達されたと判断できるから、

$\begin{align} {}_{3}C_{1} \cdot (0.9)^{2} \cdot 0.1 + 1 \cdot (0.9)^{3} = 0.972 \end{align} となる。$

$n=5$ のとき

同様にして、

正しく受信する回数	確率
0	$1 \cdot (0.1)^{5}$
1	${}_{5}C_{1} \cdot 0.9 \cdot (0.1)^{4}$
2	${}_{5}C_{2} \cdot (0.9)^{2} \cdot (0.1)^{3}$
3	${}_{5}C_{3} \cdot (0.9)^{3} \cdot (0.1)^{2}$
4	${}_{5}C_{4} \cdot (0.9)^{4} \cdot 0.1$
5	$1 \cdot (0.9)^{5}$

このうち、3回以上正しく受信できれば信号は正しく伝達されたと判断できるから、

$\begin{align} {}_{5}C_{3} \cdot (0.9)^{3} \cdot (0.1)^{2} + {}_{5}C_{4} \cdot (0.9)^{4} \cdot 0.1 + 1 \cdot (0.9)^{5} = 0.99144 \end{align}$

となる。

9.6

1時間あたりの来客数を確率変数 $X_{i}$ とすると、午前3時間の来客数は $X = X_{1} + X_{2} + X_{3}$ で表せる。各確率変数 $X_{i}$ の確率分布が母集団 $P_{o} (\lambda)$ に従うとき、 $X = X_{1} + X_{2} + X_{3}$ の確率分布は $P_{o} (n \lambda) = P_{o} (3 \lambda)$ に従う。 $\lambda = 1.5$ であるから確率分布は、

$\begin{align} f_{4.5} (x) = e^{-4.5} \cdot \frac{(4.5)^{x}}{x!} \end{align}$

であるから、午前3時間の来客数が5人以上である確率は、

$\begin{align} 1 - \sum_{x=0}^{4} f_{4.5} (x) = 0.468 \end{align}$

と求まる。参考のため、この計算を行うPythonプログラムを次に示す。

from math import exp, pow, factorial

sum_ = 0.0
lambda_ = 4.5
for x in range(0, 5):
    sum_ += exp(-lambda_) * pow(lambda_, x) / factorial(x)

print(1 - sum_)

上記のプログラムを実行すると、次の結果が得られる。

0.4678964236252845

9.7

i)

交通事故のような稀に起こる事象に関しては、ポアソン分布が利用できる。表では1年あたりの交通事故死亡者数が記載されているから、 $\lambda$ を1年間あたりの交通事故死亡者数として確率分布が次のようになる。

$\begin{align} f_{\lambda} (x) = e^{-\lambda} \frac{\lambda^{x}}{x!} \end{align}$

したがって1年間あたりの交通事故死亡者数が10人未満である確率は、次の式から得られる。

$\begin{align} \sum_{x=0}^{9} f_{\lambda} (x) = \sum_{x=0}^{9} e^{-\lambda} \frac{\lambda^{x}}{x!} \end{align}$

この式を利用して各都道府県における1年間の交通事故死亡者数が10人未満である確率を求めるプログラムは、次のようになる。

from math import exp, pow, factorial

prefs = {
    "Hokkaido": 9.7,
    "Tokyo": 4.0,
    "Osaka": 5.7,
    "Fukuoka": 7.8,
}

for pref, lambda_ in prefs.items():
    sum_ = 0.0
    for x in range(0, 10):
        sum_ += exp(-lambda_) * pow(lambda_, x) / factorial(x)

    print(f"{pref}: {sum_}")

上記のプログラムを実行すると、次の結果が得られる。

Hokkaido: 0.49597884175020945
Tokyo: 0.991867757203066
Osaka: 0.9351825280291335
Fukuoka: 0.7411089165697634

ii)

1日の交通事故死傷者数は、1年の交通事故死傷者数を365で割ったものと考えてよい。したがって1日あたりの死者数が5人未満である確率は、1年の交通事故死傷者数を $\lambda$ として次の式から得られる。

$\begin{align} \sum_{x=0}^{5} f_{\lambda / 365} (x) = \sum_{x=0}^{5} \exp \left( \frac{-\lambda}{365} \right) \frac{(\lambda / 365)^{x}}{x!} \end{align}$

この式を利用して各都道府県における1年間の交通事故死傷者数が5人未満である確率を求めるプログラムは、次のようになる。

from math import exp, pow, factorial

prefs = {
    "Hokkaido": 526.6,
    "Tokyo": 508.7,
    "Osaka": 703.8,
    "Fukuoka": 867.2,
}

for pref, lambda_ in prefs.items():
    sum_ = 0.0
    for x in range(0, 5):
        sum_ += exp(-lambda_/365) * pow(lambda_/365, x) / factorial(x)

    print(f"{pref}: {sum_}")

上記のプログラムを実行すると、次の結果が得られる。

Hokkaido: 0.9839920186786263
Tokyo: 0.9859939779933685
Osaka: 0.9535909639906193
Fukuoka: 0.9071305858827781

9.8

i)

各データを $X_{i}$ とすると、母集団の平均は次の式で与えられる。

$\begin{align} X = \sum_{n=1}^{5} \frac{X_{n}}{5} \end{align}$

母集団の平均を求めるPythonプログラムを次に示す。

import numpy as np

a = np.array([171.0, 167.3, 170.6, 178.7, 162.3])

print(np.mean(a))

上記のプログラムを実行すると、次の結果が得られる。

169.97999999999996

ii)

3人の標本 $X_{1}, X_{2}, X_{3}$ を得たとして、標本平均 $\bar{X}$ と標本分散 $s^{2}$ は次の式で与えられる。

$\displaystyle \bar{X} = \sum_{n=1}^{3} \frac{X_{n}}{3} \\ \displaystyle s^{2} = \frac{1}{3-1} \sum_{n=1}^{3} (X_{n} - \bar{X})^{2}$

母集団から3人の標本を抽出する組み合わせと、その時の標本平均と標本分散を求めるPythonプログラムを次に示す。

import numpy as np
import itertools

a = [171.0, 167.3, 170.6, 178.7, 162.3]

for c in itertools.combinations(a, 3):
    ca = np.array(c)
    X = ca.mean()
    s2 = np.sum((ca - X)**2) / (ca.size - 1)
    print(f"{c}: X={X}, s^2={s2}")

上記のプログラムを実行すると、次の結果が得られる。

(171.0, 167.3, 170.6): X=169.63333333333333, s^2=4.123333333333301
(171.0, 167.3, 178.7): X=172.33333333333334, s^2=33.8233333333332
(171.0, 167.3, 162.3): X=166.86666666666667, s^2=19.06333333333329
(171.0, 170.6, 178.7): X=173.4333333333333, s^2=20.84333333333329
(171.0, 170.6, 162.3): X=167.96666666666667, s^2=24.123333333333253
(171.0, 178.7, 162.3): X=170.66666666666666, s^2=67.32333333333315
(167.3, 170.6, 178.7): X=172.19999999999996, s^2=34.40999999999988
(167.3, 170.6, 162.3): X=166.73333333333332, s^2=17.463333333333267
(167.3, 178.7, 162.3): X=169.43333333333334, s^2=70.65333333333312
(170.6, 178.7, 162.3): X=170.53333333333333, s^2=67.24333333333314

iii)

ii)の結果を用いる。ある組み合わせにおける標本平均を $\bar{X}_{i}$ とすると、標本平均の期待値 $E(\bar{X})$ と分散 $V(\bar{X})$ はそれぞれ、

$\begin{align} E(\bar{X}) = \sum_{i=1}^{n} \frac{\bar{X}_{i}}{n} \end{align}$

$\begin{align} V(\bar{X}) &= \sum_{i=1}^{n} \left( \bar{X}_{i} - E(\bar{X}) \right)^{2} \\ &= E(\bar{X}^{2}) - E(\bar{X})^{2} \\ &= \sum_{i=1}^{n} \frac{\bar{X}_{i}^{2}}{n} - \left( \sum_{i=1}^{n} \frac{\bar{X}_{i}}{n} \right)^{2} \end{align}$

となる。これらの式を用いて標本平均の期待値と分散を求めるPythonプログラムを次に示す。

import numpy as np
import itertools

a = [171.0, 167.3, 170.6, 178.7, 162.3]

Xi = []
for c in itertools.combinations(a, 3):
    ca = np.array(c)
    X = ca.mean()
    s2 = np.sum((ca - X)**2) / (ca.size - 1)
    
    Xi.append(X)

Xi_arr = np.array(Xi)
EX = np.sum(Xi_arr) / Xi_arr.size
VX = np.sum(np.mean(Xi_arr**2) - np.mean(EX)**2)
print(f"EX: {EX}, VX: {VX}")

上記のプログラムを実行した結果を次に示す。

EX: 169.98, VX: 4.7876000000032946

ここで母平均はi)より $169.97999999999996 \fallingdotseq 169.98$ となるから、標本平均の期待値は母平均に等しいことが示された。

一方、次の関係式、

$\begin{align} V(\bar{X}) = \frac{N-n}{N-1} \cdot \frac{\sigma^{2}}{n} \end{align}$

を求めるPythonプログラムを次に示す。

import numpy as np

a = np.array([171.0, 167.3, 170.6, 178.7, 162.3])
sigma2 = np.sum(np.mean(a**2) - np.mean(a)**2)

N = a.size
n = 3
VX = (N - n) * sigma2 / ((N - 1) * n)
print(VX)

上記のプログラムを実行すると、次の結果が得られる。

4.787600000002082

これは標本平均の分散と等しい。

9.9

Pythonの random.random を利用して[0, 1]の正規分布を作って、これに2500をかけることで乱数表の場所を決めてそこから4桁の数値を取得する。重複が生じた場合は再度同様の処理を行う。 4桁の数値が1000を超える場合（★）についても再度同様の処理を行う。

5000人の場合も同様の手続きで行えるが、上記の★に相当する可能性が少なくなるのでより容易に乱数を取得できる。

2020-09-15

【統計学入門（東京大学出版会）】第8章練習問題解答

統計学

東京大学出版会から出版されている統計学入門（基礎統計学Ⅰ）について第8章の練習問題の解答を書いていく。

本章以外の解答
8.1
8.2
- i)
- ii)
8.3

本章以外の解答

本章以外の練習問題の解答は別の記事で公開している。
必要に応じて参照されたい。

8.1

母集団分布が何であれ、確率変数の和 $\bar{X} = X_{1} + X_{2} + ... + X_{n}$ は中心極限定理により正規分布 $N(n \mu, n \sigma^{2})$ に従う。このため、

$\begin{align} z = \frac{\bar{X} - n \mu}{n \sigma^{2}} \end{align}$

の標準化を行えば、

$\begin{align} P(L \le X_{1} + X_{2} + ... + X_{n} \le U) &= P \left( L - n \mu \le \bar{X} - n \mu \le U - n \mu \right) \\ &= P \left( \frac{L - n \mu}{\sqrt{n} \sigma} \le z \le \frac{U - n \mu}{\sqrt{n} \sigma} \right) \end{align}$

が得られる。 $z$ は標準正規分布に従うから、数値表より上側確率 $(1 - 0.95) / 2 = 0.025$ となる値を探すと、 $1.96$ であることがわかる。標準正規分布は偶関数であるから、

$\displaystyle \frac{L - n \mu}{\sqrt{n} \sigma} = -1.96 \\ \displaystyle \frac{U - n \mu}{\sqrt{n} \sigma} = 1.96$

を満たす $L, U$ を求めればよい。ベルヌーイ分布 $Bi(1, p)$ の期待値 $E(X) = p$ と分散 $V(X) = p(1 - p)$ を利用して式変形すると、

$L = n \mu -1.95 \sqrt{n} \sigma = np - 1.96 \sqrt{np(1-p)} \\ U = n \mu + 1.95 \sqrt{n} \sigma = np + 1.96 \sqrt{np(1-p)}$

が得られる。 $n = 700, p = 0.4$ を代入すると、 $L = 254.596, U= 305.405$ が得られる。

8.2

i)

確率分布 $f(x)$ が、

$\begin{equation} f(x) = \begin{cases} p & (X_{i} = 1) \\ q & (X_{i} = -1) \\ 0 & (X_{i} \neq 1, -1) \end{cases} \end{equation}$

であるから $p + q = 1$ を利用して、期待値 $E(X)$ は、

$\begin{align} E(X) &= \sum_{x} x f(x) \\ &= p - q = 2p - 1 \end{align}$

分散 $V(X)$ は、

$\begin{align} E(X^{2}) &= \sum_{x} x^{2} f(x) \\ &= p + q = 1 \end{align}$

より、

$\begin{align} V(X) &= E(X^{2}) - E(X)^{2} \\ &= 1 - (2p - 1)^{2} \\ &= - 4p^{2} + 4p \\ &= 4p (1 - p) =4pq \end{align}$

となる。

確率変数 $S_{n}$ は、 $n$ が大きいとき中心極限定理により正規分布 $N(n \mu, n \sigma^{2})$ に従うから、近似的確率分布は $N( n (2p - 1), 4npq)$ に従うことがわかる。

ii)

近似的確率分布 $S_{10}, S_{20}$ を描画するPythonプログラムを次に示す。

from math import sqrt, pi
import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(6, 12))
fig.patch.set_alpha(1)

def plot(idx, n):
    p = 0.4
    q = 1 - p

    mu = n * (2 * p - 1)
    sigma = sqrt(4 * n * p * q)

    x = np.arange(-20, 20, 0.1)
    y = np.exp(- (x - mu)**2 / (2 * sigma**2)) / (sqrt(2 * pi) * sigma)

    plt.subplot(2, 1, idx)
    plt.plot(x, y)
    plt.xlabel("x")
    plt.ylabel("f(x)")
    plt.title(f"3.1 ii) n={n}")
    
plot(1, 10)
plot(2, 20)

上記のプログラムを実行すると、次のグラフが描画される。

f:id:nuka137:20200915081618p:plain

8.3

450打数のときに3割のバッターになれる確率

打者がヒットを打つ確率変数を $X_{i}$ としその和を $S_{n} = X_{1} + X_{2} + ... + X_{n}$ とすると、今回の問題では $n$ が十分大きいことから中心極限定理を適用でき、 $S_{n}$ は正規分布 $N( n \mu, n \sigma^{2})$ に従う。

また今回の問題では、確率変数 $X_{i}$ はベルヌーイ分布に従うから、今シーズンにおける打者のヒット数の期待値 $E(X)$ と分散 $V(X)$ は、

$\begin{align} \mu = E(X) &= np \\ &= 450 \cdot 0.28 = 126 \end{align}$

$\begin{align} \sigma^{2} = V(X) &= np(1-p) \\ &= 450 \cdot 0.28 \cdot (1 - 0.28) = 90.72 \end{align}$

となる。

今シーズンで打率が3割以上となる確率は、 $z = \frac{S_{n} - \mu}{\sigma}$ による標準化により、

$\begin{align} P(S_{n} \ge 0.3 n) &= P(z \ge \frac{0.3n - \mu}{\sigma}) \\ &= P(z \ge 0.944911) \\ &= 0.17361 \end{align}$

と求まる。なお、最後は付表の値を参照した。

3割のバッターになれる確率が0.2以上となる必要打数

打数を $m$ とすると、

$\begin{align} \mu = E(X) &= mp \end{align}$

$\begin{align} \sigma^{2} = V(X) &= mp(1-p) \end{align}$

であるから、 $z = \frac{S_{n} - \mu}{\sigma}$ による標準化により、

$\begin{align} P(S_{n} \ge 0.3 m) &= P(z \ge \frac{0.3m - mp}{\sqrt{mp(1-p)}}) \\ &\ge 0.2 \end{align}$

を満たす $m$ を求めればよいことになる。付表を参照すると、

$\begin{align} \frac{0.3m - mp}{\sqrt{mp(1-p)}} \le 0.84 \end{align}$

より $m \le 355.62$ 、すなわち355打数以下でなければならない。

2020-09-14

【Python】cuMLの性能を評価してみた

機械学習 Python

KaggleでRAPIDSを使う方法が公開されていたので、RAPIDSに含まれている機械学習ライブラリcuMLについて、簡単な性能評価を行ってみた。
cuMLを利用することで、GPU上でランダムフォレストなどの機械学習のアルゴリズムを利用することができ、学習の高速化が期待できる。

準備

今回はKaggle上でcuMLの性能評価を行うため、Kaggleで新しいNotebookを作成した上で次の手順に従い、RAPIDSを利用できる環境を整える必要がある。

RAPIDSのDataset をNotebookに追加
NotebookのAcceleratorをGPUに変更
Notebookの最初に次のコマンド一式をコピペ

import sys
!cp ../input/rapids/rapids.0.14.0 /opt/conda/envs/rapids.tar.gz
!cd /opt/conda/envs/ && tar -xzvf rapids.tar.gz > /dev/null
sys.path = ["/opt/conda/envs/rapids/lib/python3.7/site-packages"] + sys.path
sys.path = ["/opt/conda/envs/rapids/lib/python3.7"] + sys.path
sys.path = ["/opt/conda/envs/rapids/lib"] + sys.path 
!cp /opt/conda/envs/rapids/lib/libxgboost.so /opt/conda/lib/

※ Kaggleのコンテナに搭載されているPythonのバージョン等が異なる場合は上記のコマンドがエラーとなるため、適宜修正が必要

今回はお試し利用であるため、TitanicのDatasetを利用する。ベンチマークに使用したソースコード一式は、Kaggle上で公開している。

cuMLのベンチマークプログラム

今回は、RAPIDSに含まれているDataFrameライブラリ cuDF と機械学習ライブラリ cuML を利用し、scikit-learn と性能比較した。

cuDFとcuMLを利用するためには、以下のパッケージをインポートする必要がある。

import cudf
import cuml

pandasのDataFrameからcuDFのDataFrameへは、cudf.from_pandas で変換できる。

train_df = pd.read_csv("../input/titanic/train.csv")
train_df_gpu = cudf.from_pandas(train_df)

cuMLはscikit-learnと似たようなインターフェースを提供しているため、過去にscikit-learnを使ったことがある人なら比較的簡単にcuMLを使うことができる。しかし一部のAPIの引数が異なっていることと、入力データとして受け付けるデータ型に制約があることに注意しなければならない。特にデータ型については、scikit-learnに入力したデータ型と同じデータ型でcuMLのAPIを利用しようとして何度もエラーで悩まされたため、結局全てのデータ型をfloat32に型変換した。

今回は、以下の3つのアルゴリズムを性能評価対象とした。

k近傍法
サポートベクタ―マシン
ランダムフォレスト

scikit-learnとcuMLを用いて性能評価で使用したベンチマークプログラムを次に示す。

import time

import pandas as pd

import sklearn.neighbors
import sklearn.svm
import sklearn.ensemble
from sklearn.model_selection import KFold

import cudf
import cuml

import matplotlib.pyplot as plt
import numpy as np

NFOLDS = 5
ITERATION = 300

# Load data.
train_orig_df = pd.read_csv("../input/titanic/train.csv")
test_orig_df = pd.read_csv("../input/titanic/test.csv")

# Pre-process
train_df = train_orig_df.copy()
train_df.drop(["Cabin", "Ticket", "Name"], axis=1, inplace=True)
train_df = pd.get_dummies(train_df.iloc[:, 1:], columns=["Pclass", "Sex", "Embarked"])
train_df.dropna(inplace=True)

X_all = train_df.drop(["Survived"], axis=1).astype("float32")
y_all = train_df["Survived"].astype("int32")

X_all_gpu = cudf.from_pandas(X_all)
y_all_gpu = cudf.from_pandas(y_all)

# Benchmark code.
def bench(X, y, classifiers, params):
    elapsed = {}
    for name, clf_class in classifiers.items():
        elapsed_list = []

        for _ in range(ITERATION):
            kf = KFold(n_splits=NFOLDS)
            clf = clf_class()
            clf.set_params(**params[name])

            elapsed_sum = 0
            for i, (train_idx, val_idx) in enumerate(kf.split(X, y)):
                X_train = X_all.iloc[train_idx]
                y_train = y_all.iloc[train_idx]
                X_val = X_all.iloc[val_idx]
                y_val = y_all.iloc[val_idx]

                start = time.time()
                clf.fit(X_train, y_train)
                elapsed_sum += time.time() - start

            elapsed_list.append(elapsed_sum)

        elapsed[name] = pd.Series(elapsed_list).mean()
    return elapsed

# scikit-learn
classifiers = {
    "KNN": sklearn.neighbors.KNeighborsClassifier,
    "SVM": sklearn.svm.SVC,
    "RandomForest": sklearn.ensemble.RandomForestClassifier
}

params = {
    "KNN": {},
    "SVM": {
        "random_state": 47
    },
    "RandomForest": {
        "n_estimators": 100,
        "random_state": 47
    }
}

elapsed_sklearn = bench(X_all, y_all, classifiers, params)

# cuML
classifiers = {
    "KNN": cuml.neighbors.KNeighborsClassifier,
    "SVM": cuml.svm.SVC,
    "RandomForest": cuml.ensemble.RandomForestClassifier
}

params = {
    "KNN": {},
    "SVM": {},
    "RandomForest": {
        "n_estimators": 100
    }
}

elapsed_cuml = bench(X_all_gpu, y_all_gpu, classifiers, params)

# Results.
left = np.arange(len(elapsed_sklearn.keys()))
width = 0.3

fig = plt.figure(figsize=(6, 6))
fig.patch.set_alpha(1)

plt.subplot(1, 1, 1)

plt.bar(left, elapsed_sklearn.values(), color='b', width=width, label="Scikit Learn", align="center")
plt.bar(left + width, elapsed_cuml.values(), color="g", width=width, label="cuML", align="center")

plt.xticks(left + width / 2, elapsed_sklearn.keys())
plt.legend(loc=2)
plt.ylabel("sec / iter")
plt.title("fit() performance")
plt.show()

性能評価

ベンチマークプログラムを実行すると、次のような結果が得られた。

f:id:nuka137:20200914132442p:plain

いずれのアルゴリズムについてもcuMLの方が優れているが、今回性能評価に利用したTitanicのデータ量が小さいこともあり、ランダムフォレスト除くとその性能差はわずかしかない。より大きなデータを用いた学習など、計算量が多いものではcuMLのほうが優勢になると考えられる。

データ型の制約によりcuMLの使い勝手はscikit-learnと比較して劣るものの、計算量の多いアルゴリズムや非常に大きなデータを扱うときにはcuMLによる高速化が活きてくると考える。 cuMLなどを開発しているRAPIDSには KaggleのGrandmasterの人も在席しているため、使い勝手の面でも性能面でも改善が期待できるcuMLの発展が楽しみである。

【2020.12.11 追記】

Home Credit Default Risk のデータセットを利用し、KNNとランダムフォレストについて性能評価した。

f:id:nuka137:20201211082109p:plain

Titanicと比較してデータサイズが大きいこともあり、scikit-learnと比較してよりcuMLのほうが性能が高い結果が得られた。このことから、データサイズが大きい場合はcuMLを利用したときの効果が非常に高くなると言える。一方で、データサイズが小さいとGPUの処理量に対してCPUの処理量の割合が多くなるため、cuMLがあまり活きてこないと考えられる。

なお今回評価に使用したソースコードは、KaggleのNotebook として公開している。

2020-09-13

【統計学入門（東京大学出版会）】第7章練習問題解答

統計学

東京大学出版会から出版されている統計学入門（基礎統計学Ⅰ）について第7章の練習問題の解答を書いていく。

本章以外の解答
7.1
- i)
- ii)
7.2
- i)
- ii)
- iii)
7.3
- 7.4
7.5
7.6
- i)
- ii)
- iii)
7.7
- i)
- ii)
7.8
- i)
- ii)
- iii)
7.9
- i)
- ii)
- iii)

本章以外の解答

本章以外の練習問題の解答は別の記事で公開している。
必要に応じて参照されたい。

7.1

i)

$E(X) = \mu_{X}$ かつ $E(Y) = \mu_{Y}$ とすると、共分散 ${\rm Cov} (X, Y)$ の定義より、

$\begin{align} {\rm Cov} (X, Y) &= E( (X - \mu_{X})(Y - \mu_{Y}) ) \\ &= E(XY) - \mu_{X} E(Y) - \mu_{Y} E(X) + \mu_{X} \mu_{Y} \\ &= E(XY) - \mu_{X}\mu_{Y} \end{align}$

が得られ、

$\begin{align} V(X, Y) &= E( ( (X+Y) - (\mu_{X} + \mu_{Y}) )^{2} ) \\ &= E( (X + Y)^{2} ) - 2 ( \mu_{X} + \mu_{Y} ) E( X + Y ) + (\mu_{X} + \mu_{Y})^{2} ) \\ &= E( X^{2} + 2XY + Y^{2} ) - (\mu_{X} + \mu_{Y})^{2} \\ &= E( X^{2} ) - \mu_{X}^{2} + E( Y^{2} ) - \mu_{Y}^{2} + 2 ( E( XY ) - \mu_{X} \mu_{Y} ) \\ &= V(X) + V(Y) + 2 {\rm Cov} (X, Y) \end{align}$

が得られる。

ii)

i)と同様にして、

$\begin{align} V(X, Y) &= E( ( (aX+bY) - (a \mu_{X} + b \mu_{Y}) )^{2} ) \\ &= E( (aX + bY)^{2} ) - 2 (a \mu_{X} + b \mu_{Y} ) E( aX + bY ) + (a \mu_{X} + b \mu_{Y})^{2} ) \\ &= E( a^{2}X^{2} + 2abXY + b^{2}Y^{2} ) - (a^{2} \mu_{X}^{2} + 2ab \mu_{X} \mu_{Y} + b^{2} \mu_{Y}^{2}) \\ &= a^{2} (E(X^{2}) - \mu_{X}^{2}) + b^{2} ( E(Y^{2}) - \mu_{Y}^{2} ) + 2ab ( E(XY) - \mu_{X} \mu_{Y} ) \\ &= a^{2} V(X) + b^{2} V(Y) + 2ab {\rm Cov} (X, Y) \end{align}$

7.2

i)

期待値

$\begin{align} E(R_{p}) &= E( xR_{1} + (1 - x)R_{2} ) \\ &= xE(R_{1}) + (1 - x)E(R_{2}) \\ &= x e_{1} + (1 - x) e_{2} \end{align}$

分散

$\begin{align} V(R_{p}) &= V( xR_{1} + (1 - x)R_{2} ) \\ &= x^{2}V(R_{1}) + (1 - x)^{2}V(R_{2}) + 2 {\rm Cov} (xR_{1}, (1 - x)R_{2}) \\ &= x^{2} \sigma_{1}^{2} + (1 - x)^{2} \sigma_{2}^{2} + 2 \rho \sqrt{ x^{2} \sigma_{1}^{2}} \sqrt{(1 - x)^{2} \sigma_{2}^{2}} \\ &= (\sigma_{1}^{2} - 2 \rho \sigma_{1} \sigma_{2} + \sigma_{2}^{2} ) x^{2} + 2 \sigma_{2} ( \rho \sigma_{1} - \sigma_{2} ) x + \sigma_{2}^{2} \end{align}$

ii)

i)で求めた分散 $V(R_{p})$ を式変形すると、

$\begin{align} V(R_{p}) &= (\sigma_{1}^{2} - 2 \rho \sigma_{1} \sigma_{2} + \sigma_{2}^{2} ) x^{2} + 2 \sigma_{2} ( \rho \sigma_{1} - \sigma_{2} ) x + \sigma_{2}^{2} \\ &= (\sigma_{1}^{2} - 2 \rho \sigma_{1} \sigma_{2} + \sigma_{2}^{2} ) \left( x + \frac{\sigma_{2} ( \rho \sigma_{1} - \sigma_{2} )}{(\sigma_{1}^{2} - 2 \rho \sigma_{1} \sigma_{2} + \sigma_{2}^{2} )} \right)^{2} - \left( \frac{\sigma_{2} ( \rho \sigma_{1} - \sigma_{2} )}{\sigma_{1}^{2} - 2 \rho \sigma_{1} \sigma_{2} + \sigma_{2}^{2}} \right)^{2} + \sigma_{2}^{2} \\ &\equiv \alpha ( x - \gamma )^{2} + \beta \end{align}$

が得られる。

なお、 $\alpha$ は $-1 \le \rho \le 1$ より、

$\begin{align} \alpha &\equiv \sigma_{1}^{2} - 2 \rho \sigma_{1} \sigma_{2} + \sigma_{2}^{2} \\ &= (\sigma_{1} - \sigma_{2})^{2} + 2\sigma_{1}\sigma_{2} (1 - \rho) \\ & \ge 0 \end{align}$

となり、 $V(R_{p})$ は下に凸の関数であることがわかる。

ここで $0 \le x \le 1$ を考慮すると、最小値をとる $x$ の値は $\gamma$ の値によって変わることに注意する。

	$\gamma$ の値	$V(R_{p})$ が最小となる $x$
①	$\gamma \le 0$	$x = 0$
②	$0 \lt \gamma \lt 1$	$x = \gamma$
③	$\gamma \ge 1$	$x = 1$

① $\gamma \le 0$

$\gamma \le 0$ 、すなわち $\rho \ge \frac{\sigma_{2}}{\sigma_{1}}$ のとき、 $x = 0$ で $V(R_{p})$ が最小値をとる。

したがって $x = 0$ を $V(R_{p})$ に代入すると、 $min(V(R_{p})) = \sigma_{2}^{2}$ が得られる。

② $0 \lt \gamma \lt 1$

$0 \lt \gamma \lt 1$ のときは、 $x = \gamma$ で $V(R_{p})$ が最小値をとる。

なお $\gamma$ は、

$\displaystyle \gamma \equiv - \frac{\sigma_{2} ( \rho \sigma_{1} - \sigma_{2} )}{\sigma_{1}^{2} - 2 \rho \sigma_{1} \sigma_{2} + \sigma_{2}^{2} }$

である。

したがって $x = \gamma$ を $V(R_{p})$ に代入すると、

$\begin{align} min( V(R_{p}) ) &= (\sigma_{1}^{2} - 2 \rho \sigma_{1} \sigma_{2} + \sigma_{2}^{2} )\gamma^{2} + 2 \sigma_{2} ( \rho \sigma_{1} - \sigma_{2} ) \gamma + \sigma_{2}^{2} \\ &= \frac{\sigma_{1}^{2} \sigma_{2}^{2} ( 1 - \rho^{2})}{\sigma_{1}^{2} - 2 \rho \sigma_{1} \sigma_{2} + \sigma_{2}^{2} } \end{align}$

が得られる。

③ $\gamma \ge 1$

$\gamma \ge 1$ 、すなわち $\rho \ge \frac{\sigma_{1}}{\sigma_{2}}$ のとき、 $x = 1$ で $V(R_{p})$ が最小値をとる。

したがって $x = 1$ を $V(R_{p})$ に代入すると、 $min(V(R_{p})) = \sigma_{1}^{2}$ が得られる。

iii)

$E(R_{p})$ と $V(R_{p})$ を描画するPythonプログラムを次に示す。

import numpy as np
import matplotlib.pyplot as plt


e1 = 0.198
e2 = 0.055
sigma1 = 0.357
sigma2 = 0.203
rho = 0.18

x = np.arange(0, 5, 0.01)
E = x * e1 + (1 - x) * e2
V = (sigma1**2 - 2 * rho * sigma1 * sigma2 + sigma2**2) * x**2 + 2 * sigma2 * (rho * sigma1 - sigma2) * x + sigma2**2

fig = plt.figure(figsize=(6, 12))
fig.patch.set_alpha(1)

plt.subplot(2, 1, 1)
plt.plot(x, E)
plt.xlabel("x")
plt.ylabel("E")
plt.title("7.2 iii) E")

plt.subplot(2, 1, 2)
plt.plot(x, V)
plt.xlabel("x")
plt.ylabel("V")
plt.title("7.2 iii) V")

プログラムを実行すると次のようなグラフが描画される。

f:id:nuka137:20200908085637p:plain

7.3

A, Bそれぞれのつぼにボールが入る確率を1/2とする。

$X$ の確率分布 $g(x)$ を示す。

$X$	確率
0	$\displaystyle {}_{3}C_{0} \left( \frac{1}{2} \right)^{3} = \frac{1}{8}$
1	$\displaystyle {}_{3}C_{1} \left( \frac{1}{2} \right)^{2} \cdot \left( \frac{1}{2} \right) = \frac{3}{8}$
2	$\displaystyle {}_{3}C_{2} \left( \frac{1}{2} \right) \cdot \left( \frac{1}{2} \right)^{2} = \frac{3}{8}$
3	$\displaystyle {}_{3}C_{3} \left( \frac{1}{2} \right)^{3} = \frac{1}{8}$

3つのボールが1つのつぼに入るのは $X = 0, X = 3$ のとき、2つのつぼに入るのは $X = 1, X = 2$ のときである。このことから、 $Y$ の確率分布 $h(y)$ は次のように得られる。

$Y$	確率
1	$\displaystyle \frac{1}{8} + \frac{1}{8} = \frac{1}{4}$
2	$\displaystyle \frac{3}{8} + \frac{3}{8} = \frac{3}{4}$

これを同時確率分布表にすると、次のようになる。


		Y
X		1	2
	0	1/8	0
	1	0	3/8
	2	0	3/8
	3	1/8	0

独立の条件は、 $f(x, y) = g(x) \cdot h(x)$ が成り立つことが条件であるが、

$\begin{align} f(x = 0, y = 1) &= \frac{1}{8} \cdot \frac{1}{4} = \frac{1}{32} \\ &\neq \frac{1}{8} \end{align}$

であるため、独立ではないと言える。

一方で無相関の条件は ${\rm Cov} (X, Y) = 0$ です。

$\begin{align} E(X) &= 0 \cdot \left( \frac{1}{8} \right) + 1 \cdot \left( \frac{3}{8} \right) + 2 \cdot \left( \frac{3}{8} \right) + 3 \cdot \left( \frac{1}{8} \right) \\ &= \frac{3}{2} \end{align}$

$\begin{align} E(Y) &= 1 \cdot \left( \frac{1}{4} \right) + 2 \cdot \left( \frac{3}{4} \right) \\ &= \frac{7}{4} \end{align}$

$\begin{align} E(XY) &= 0 \cdot 1 \cdot \left( \frac{1}{8} \right) + 1 \cdot 2 \cdot \left( \frac{3}{8} \right) + 2 \cdot 2 \cdot \left( \frac{3}{8} \right) + 3 \cdot 1 \cdot \left( \frac{1}{8} \right) \\ &= \frac{21}{8} \end{align}$

より、

$\begin{align} {\rm Cov} (X, Y) &= E(XY) - E(X) \cdot E(Y) \\ &= \frac{21}{8} - \frac{3}{2} \cdot \frac{7}{4} \\ &= 0 \end{align}$

となるため、無相関であると言える。

7.4

(Ⅰ)の方法で物体AとBを測定したときの分散は、それぞれ $V(m_{A}) = \sigma^{2}$ , $V(m_{B}) = \sigma^{2}$ となる。

一方で(Ⅱ)の方法で測定したときは、それぞれ $V(m_{X}) = \sigma^{2}$ , $V(m_{Y}) = \sigma^{2}$ となる。

$\displaystyle m_{X} = m_{A} + m_{B} \\ \displaystyle m_{Y} = m_{A} - m_{B}$

であるため、物体AとBがそれぞれ相関がないこと7.1 ii)の結果を利用して、

$\begin{align} V(m_{A}) &= V \left( \frac{m_{X} + m_{Y}}{2} \right) \\ &= \frac{1}{4} V(m_{X}) + \frac{1}{4} V(m_{Y}) \\ & = \frac{\sigma^{2}}{2} \end{align}$

同様にして、

$\begin{align} V(m_{B}) &= V \left( \frac{m_{X} - m_{Y}}{2} \right) \\ &= \frac{1}{4} V(m_{X}) + \frac{1}{4} V(m_{Y}) \\ & = \frac{\sigma^{2}}{2} \end{align}$

が得られる。

したがって、(Ⅰ)の方法で測定するよりも(Ⅱ)の方法で測定した方が分散が少なくなる（測定によるばらつきが小さい）ため、より優れた方法であると言える。

7.5

$\begin{align} \rho_{UV} &= \frac{{\rm Cov} (U, V)}{\sqrt{V(U)} \sqrt{V(V)}} \\ &= \frac{{\rm Cov} (aX+b, cY+d)}{\sqrt{V(aX+b)} \sqrt{V(cY+d)}} \end{align}$

ここで、

$\begin{align} {\rm Cov} (aX+b, cY+d) &= E( (aX + b) (cY + d) ) - E( aX + b ) E( cY + d ) \\ &= E( (aX + b) (cY + d) ) - (a E(X) + b) (c E(Y) + d) \\ &= E( acXY + adX + bcY + bd ) - (a E(X) + b) (c E(Y) + d) \\ &= ac( E(XY) - E(X) E(Y) ) \end{align}$

$V(aX+b) = a^{2} V(X) \\ V(cY+d) = c^{2} V(Y)$

より、

$\begin{align} \rho_{UV} &= \frac{{\rm Cov} (aX+b, cY+d)}{\sqrt{V(aX+b)} \sqrt{V(cY+d)}} \\ &= \frac{ac( E(XY) - E(X) E(Y) )}{a \sqrt{V(X)} \cdot c \sqrt{V(Y)}} \\ &= \frac{{\rm Cov} (X, Y)}{\sqrt{V(X)} \sqrt{V(Y)}} \\ &= \rho_{XY} \end{align}$

7.6

i)

$X, Y$ はともに標準正規分布に従い、かつ互いに独立であるから、

$E(X) = E(Y) = 0 \\ V(X) = V(Y) = 1 \\ E(X^{2}) = E(Y^{2}) = 1 \\ E(XY) = E(X)E(Y) = 0$

となる。

このため、

$\begin{equation} E(cX + Y) = cE(X) + E(Y) = 0 \\ V(cX + Y) = c^{2}V(X) + V(Y) = c^{2} + 1 \\ E(X(cX+Y)) = cE(X^{2}) + E(XY) = c \\ {\rm Cov} (X, cX+Y) = E(X(cX+Y)) - E(X)E(Y) = c \end{equation}$

が得られる。

これらを用いると、相関係数 $\rho$ は、

$\begin{align} \rho &= \frac{{\rm Cov} (X, cX+Y)}{\sqrt{V(X)} \sqrt{V(cX+Y)}} \\ &= \frac{c}{\sqrt{c^{2} + 1}} \\ &= 0.5 \end{align}$

の関係式が得られるため、 $c$ について解くと、

$\displaystyle c = \frac{1}{\sqrt{3}}$

が得られる。

ii)

i)より、

$\displaystyle \rho = \frac{c}{ \sqrt{ c^{2} + 1 } }$

であるから、 $c$ について解くと、

$\displaystyle c = \frac{\rho}{ \sqrt{ 1 - \rho^{2} } }$

が得られる。

iii)

$U = aX + bY \\ V = cX + dY$

とすると、本書の (7.30) より、

$\sigma_{1}^{2} = a^{2} + b^{2} \\ \sigma_{2}^{2} = c^{2} + d^{2} \\ \sigma_{12} = ac + bd \\ \displaystyle \rho = \frac{\sigma_{12}}{\sigma_{1} \sigma_{2}}$

となる。

$U, V$ は2次元正規分布 ${\rm N} ( (0, 0), (\sigma_{1}^{2}, \sigma_{2}^{2}, \rho ) )$ に従うから、

$\sigma_{12} = ac + bd = \rho \\ \sigma_{1} \sigma_{2} = 1$

が得られる。

ここで $b = 0$ とおくと、

$\sigma_{1}^{2} = a^{2} + b^{2}$

より、

$a = \sigma_{1}$

が得られる。

また、

$\sigma_{12} = ac + bd = \rho$

より、

$\displaystyle c = \frac{\rho}{a} = \frac{\rho}{\sigma_{1}} = \sigma_{2} \rho$

が得られる。

最後に、

$\sigma_{2}^{2} = c^{2} + d^{2}$

から、

$d = \sigma_{2} \sqrt{1 - \rho^{2}}$

が得られる。

したがって、

$U = \sigma_{1} X \\ V = \sigma_{2} \rho X + \sigma_{2} \sqrt{1 - \rho^{2}} Y$

と求まる。

7.7

i)

互いに独立なシステム $S_{1}, S_{2}$ が並列に結合されているとき、全体の寿命 $Y$ は2つのシステムが寿命を迎えたタイミングであるから、 $max( X_{1}, X_{2} )$ が全体の寿命となる。

システム $S_{1}, S_{2}$ の寿命は指数分布、

$\begin{equation} Ex(\lambda) = \begin{cases} \displaystyle \lambda e^{- \lambda x} & (x \ge 0) \\ 0 & (x \lt 0) \end{cases} \end{equation}$

で与えられ、累積分布関数 $F(x)$ は、

$\begin{align} F(x) &= \int_{-\infty}^{x} \\ &= \int_{0}^{x} \lambda e^{- \lambda x} dx \\ &= 1 - e^{- \lambda x} \end{align}$

である。

$Y$ の寿命を $y$ とすると、その累積分布関数は $S_{1}, S_{2}$ のどちらか一方が $y$ まで寿命を迎えていなければよいため、

$\begin{align} P(Y \le y) &= P(X_{1} \le y, X_{2} \le y) \\ &= P(X_{1} \le y) P(X_{2} \le y) \\ &= (1 - e^{- \lambda y} )^{2} \end{align}$

と求まる。

したがってその確率密度関数は、

$\begin{align} \frac{d}{dy} P(Y \le y) &= \frac{d}{dy} (1 - e^{- \lambda y} )^{2} \\ &= 2 \lambda e^{- \lambda y} (1 - e^{- \lambda y} ) \end{align}$

ii)

互いに独立なシステム $S_{1}, S_{2}$ が直列に結合されているとき、全体の寿命 $Y$ はどちらか一方のシステムが寿命を迎えたタイミングであるから、 $min( X_{1}, X_{2} )$ が全体の寿命となる。

したがって、全確率から $S_{1}, S_{2}$ 両方のシステムが寿命を迎えていない確率を引けばよいため、

$\begin{align} P(Y \le y) &= 1 - P(X_{1} \gt y, X_{2} \gt y) \\ &= 1 - P(X_{1} \gt y) P(X_{2} \gt y) \\ &= 1 - (1 - (1 - e^{- \lambda y} ) )^{2} \\ &= 1 - e^{- 2 \lambda y} \end{align}$

と求まる。

確率密度関数は、

$\begin{align} \frac{d}{dy} P(Y \le y) &= \frac{d}{dy} (1 - e^{- 2 \lambda y}) \\ &= 2 \lambda e^{-2 \lambda y} \end{align}$

7.8

i)

確率密度関数が、

$\begin{equation} f(x) = \begin{cases} \displaystyle 1 & ( 0 \le x \le 1) \\ 0 & (x \lt 0, 1 \lt x) \end{cases} \end{equation}$

で与えられるから、累積分布関数 $F(x)$ は、

$\begin{equation} F(x) = \begin{cases} \displaystyle x & ( 0 \le x \le 1 ) \\ 0 & (x \lt 0, 1 \lt x) \end{cases} \end{equation}$

となる。

したがって、UとVの累積分布関数 $U(u), V(v)$ はそれぞれ、

$\begin{align} U(u) &= P(max (X_{1}, X_{2}, ..., X_{n}) \le u) \\ &= P(X_{1} \le u) P(X_{2} \le u) ... P(X_{n} \le u) \\ &= (F(u))^{n} \\ &= u^{n} \quad ( 0 \le u \le 1 ) \end{align}$

$\begin{align} V(v) &= P(min (X_{1}, X_{2}, ..., X_{n}) \le v) \\ &= 1 - P(X_{1} \gt v) P(X_{2} \gt v) ... P(X_{n} \gt v) \\ &= 1 - ( 1 - v )^{n} \quad ( 0 \le v \le 1 ) \end{align}$

となる。このため、UとVの確率密度関数 $g(u), h(v)$ はそれぞれ、

$\begin{align} g(u) &= \frac{dU}{du} \\ &= nu^{n-1} \quad ( 0 \le u \le 1 ) \end{align}$

$\begin{align} h(v) &= \frac{dV}{dv} \\ &= n(1-v)^{n-1} \quad ( 0 \le v \le 1 ) \end{align}$

となる。

ii)

確率密度関数が、

$\begin{equation} f(x) = \begin{cases} \displaystyle \lambda e^{- \lambda x} & ( x \ge 0) \\ 0 & (x \lt 0) \end{cases} \end{equation}$

で与えられるから、累積分布関数 $F(x)$ は、

$\begin{equation} F(x) = \begin{cases} \displaystyle 1 - e^{- \lambda x} & ( x \ge 0 ) \\ 0 & (x \lt 0) \end{cases} \end{equation}$

となる。

したがって、UとVの累積分布関数 $U(u), V(v)$ はそれぞれ、

$\begin{align} U(u) &= P(max (X_{1}, X_{2}, ..., X_{n}) \le u) \\ &= P(X_{1} \le u) P(X_{2} \le u) ... P(X_{n} \le u) \\ &= (1 - e^{- \lambda u})^{n} \quad ( u \ge 0 ) \end{align}$

$\begin{align} V(v) &= P(min (X_{1}, X_{2}, ..., X_{n}) \le v) \\ &= 1 - P(X_{1} \gt v) P(X_{2} \gt v) ... P(X_{n} \gt v) \\ &= 1 - e^{- \lambda nv} \quad ( v \ge 0 ) \end{align}$

となる。このため、UとVの確率密度関数 $g(u), h(v)$ はそれぞれ、

$\begin{align} g(u) &= \frac{dU}{du} \\ &= n \lambda e^{- \lambda u} ( 1 - e^{- \lambda u } )^{n-1} \quad ( u \ge 0 ) \end{align}$

$\begin{align} h(v) &= \frac{dV}{dv} \\ &= \lambda n e^{- \lambda n v} \quad ( v \ge 0 ) \end{align}$

となる。

iii)

各 $X_{i}$ の確率密度関数が $f(x)$ 、累積分布関数が $F(x)$ で与えられるから、UとVの累積分布関数 $U(u), V(v)$ はそれぞれ、

$\begin{align} U(u) &= P(max (X_{1}, X_{2}, ..., X_{n}) \le u) \\ &= P(X_{1} \le u) P(X_{2} \le u) ... P(X_{n} \le u) \\ &= F(u)^{n} \end{align}$

$\begin{align} V(v) &= P(min (X_{1}, X_{2}, ..., X_{n}) \le v) \\ &= 1 - P(X_{1} \gt v) P(X_{2} \gt v) ... P(X_{n} \gt v) \\ &= 1 - (1 - F(v))^{n} \end{align}$

で与えられる。

このため、UとVの確率密度関数 $g(u), h(v)$ はそれぞれ、

$\begin{align} g(u) &= \frac{dU}{du} \\ &= n \cdot \frac{dF}{du} \cdot F(u)^{n-1} \\ &= n \cdot f(u) \cdot F(u)^{n-1} \end{align}$

$\begin{align} h(v) &= \frac{dV}{dv} \\ &= n \cdot \frac{dF}{dv} \cdot (1 - F(v))^{n-1} \\ &= n \cdot f(v) \cdot (1 - F(v))^{n-1} \end{align}$

となる。

7.9

たたみこみの結果元と同じ確率分布となる場合、その確率分布は再生性も持つと言える。

i)

二項分布の確率分布は、

$f(x) = Bi(n, p) = {}_{n} C_{x} p^{x} (1 - p)^{n-x}$

であるから、再生性が持つことを証明するためには、

$Bi(n, p) * Bi(m, p) = Bi(n + m, p)$

が成り立つことを確認すればよい。

$\sum_{x=0}^{z} {}_{n} C_{x} \cdot {}_{m} C_{z-x} = {}_{n+m} C_{z}$

に注意すると、

$\begin{align} Bi(n, p) * Bi(m, p) &= \sum_{x=0}^{z} {}_{n} C_{x} p^{x} (1 - p)^{n -x} \dot {}_{m} C_{z-x} p^{z - x} (1 - p)^{m - (z - x)} \\ &= p^{z} (1-p)^{n+m-z} \cdot \sum_{x=0}^{z} {}_{n} C_{x} \cdot {}_{m} C_{z-x} \\ &= {}_{n+m} C_{z} \cdot p^{z} (1-p)^{n+m-z} \\ &= Bi(n+m, p) \end{align}$

となり、再生性を持つことが証明された。

ii)

二項分布の確率分布は、

$\begin{align} f(x) = P_{o}(\lambda) = \sum_{x} e^{-\lambda} \frac{\lambda^{x}}{x!} \end{align}$

であるから、再生性が持つことを証明するためには、

$\begin{align} P_{o}(\lambda) * P_{o}(\mu) = P_{o}(\lambda + \mu) \end{align}$

が成り立つことを確認すればよい。

二項定理

$\begin{align} (a + b)^{n} &= \sum_{x=0}^{n} {}_{n} C_{x} a^{x} b^{n-x} \\ &= \sum_{x=0}^{n} \frac{n!}{x! (n-x)!} a^{x} b^{n-x} \end{align}$

を利用すると、

$\begin{align} P_{o}(\lambda) * P_{o}(\mu) &= \sum_{x} e^{-\lambda} \frac{\lambda^{x}}{x!} \cdot e^{-\mu} \frac{\mu^{z-x}}{(z-x)!} \\ &= \frac{e^{-(\lambda + \mu)}}{z!} \sum_{x=0}^{z} \frac{z!}{x! (z-x)!} \lambda^{x} \mu^{z-x} \\ &= \frac{e^{-(\lambda + \mu)}}{z!} (\lambda + \mu)^{z} \\ &= P_{o}(\lambda + \mu) \end{align}$

となり、再生性を持つことが証明された。

iii)

二項分布の確率分布は、

$\begin{align} f(x) = \frac{1}{\sqrt{2 \pi} \sigma} \exp \left( - \frac{(x-\mu)^{2}}{2 \sigma^{2}} \right) \end{align}$

であるから、再生性が持つことを証明するためには、

$\begin{align} f_{\mu_{1}, \sigma_{1}^{2}}(x) * f_{\mu_{2}, \sigma_{2}^{2}}(x) = f_{\mu_{1} + \mu_{2}, \sigma_{1}^{2} + \sigma_{2}^{2}} (x) \end{align}$

が成り立つことを確認すればよい。

ガウス積分の公式を利用すると、

$\begin{align} f_{\mu_{1}, \sigma_{1}^{2}}(x) * f_{\mu_{2}, \sigma_{2}^{2}}(x) &= \int_{-\infty}^{\infty} \frac{1}{\sqrt{2 \pi} \sigma_{1}} \exp \left( - \frac{(x - \mu_{1})^{2}}{2 \sigma_{1}^{2}} \right) \cdot \frac{1}{\sqrt{2 \pi} \sigma_{2}} \exp \left( - \frac{( (z - x) - \mu_{2})^{2}}{2 \sigma_{2}^{2}} \right) dx \\ &= \frac{1}{2 \pi \sigma_{1} \sigma_{2}} \int_{-\infty}^{\infty} \exp \left( - \frac{(x - \mu_{1})^{2}}{2 \sigma_{1}^{2}} - \frac{( (z - x) - \mu_{2})^{2}}{2 \sigma_{2}^{2}} \right) dx \\ &= \frac{1}{2 \pi \sigma_{1} \sigma_{2}} \int_{-\infty}^{\infty} \exp \left( - \frac{1}{2} \cdot \frac{ \sigma_{1}^{2} + \sigma_{2}^{2} }{ \sigma_{1}^{2} \sigma_{2}^{2} } \left( x - \frac{1}{\sigma_{1}^{2} + \sigma_{2}^{2}} (\mu_{1} \sigma_{2}^{2} + (z - \mu_{2}) \sigma_{1}^{2} ) \right)^{2} \right) dx - \frac{1}{2 (\sigma_{1}^{2} + \sigma_{2}^{2})} (z - \mu_{1} - \mu_{2})^{2} \\ &\equiv \frac{1}{2 \pi \sigma_{1} \sigma_{2}} \int_{-\infty}^{\infty} \exp \left( - \frac{1}{2} B ( x - A )^{2} + C \right) dx \\ &= \frac{e^{C}}{2 \pi \sigma_{1} \sigma_{2}} \int_{-\infty}^{\infty} \exp \left( - \frac{1}{2} B X^{2} \right) dX \\ &= \frac{e^{C}}{2 \pi \sigma_{1} \sigma_{2}} \sqrt{\frac{2 \pi}{B}} \\ &= \frac{1}{\sqrt{2 \pi}} \frac{1}{ \sqrt{ \sigma_{1}^{2} + \sigma_{2}^{2} } } \exp \left( - \frac{(z - \mu_{1} - \mu_{2})^{2}}{2 (\sigma_{1}^{2} + \sigma_{2}^{2})} \right) \\ &= f_{\mu_{1} + \mu_{2}, \sigma_{1}^{2} + \sigma_{2}^{2}} (x) \end{align}$

となり、再生性を持つことが証明された。

なお途中の式変形で、

$\begin{align} A = \frac{-\sigma_{1}^{2} z - \sigma_{2}^{2} \mu_{1} + \sigma_{1}^{2} \mu_{2} }{ \sigma_{1}^{2} + \sigma_{2}^{2} } \end{align}$

$\begin{align} B = \frac{\sigma_{1}^{2} + \sigma_{2}^{2} }{ \sigma_{1}^{2} \sigma_{2}^{2}} \end{align}$

$\begin{align} C = - \frac{1}{2(\sigma_{1}^{2} + \sigma_{2}^{2})} (z - \mu_{1} - \mu_{2})^{2} \end{align}$

$\begin{align} X = x - A \end{align}$

とおいた。

2020-09-09

【論文読み】EfficientNet

機械学習画像処理論文

arxiv.org

Kaggleの画像コンペで多用されているEfficientNetについて、中身を知らずに使うのもどうかと思って元の論文を読んでみたので、理解したことをまとめてみる。

EfficientNetとは？

従来のCNNにおいては、次の3点の中から1点を選んでモデルをスケールアップさせることが、精度を向上させるための主な施策であった。

ネットワークの層
ネットワークの幅（Channel数）
入力画像の解像度

もちろん過去にも、上記の中から2点以上を選んでモデルをスケールアップさせる試みも行われていなかったわけではない。しかし2点以上を選んでモデルをスケールアップしようとすると、選択自由度が増える分探索空間が広くなるため、チューニングに手間がかかる問題があった。

論文ではこの問題に対して、CNNにおいてモデルをスケールアップする方式 Compound Model Scaling を提案している。 Compound Model Scalingでは、先の3点についてGrid Searchによりスケーリング係数を決め、モデルをスケールアップするための関係式を提供する。モデルをスケールアウトする際にこの関係式を利用することで、チューニングの手間を省きつつ限られた計算資源の中で精度のよいモデルを学習できる。もちろん環境ごとにチューニングしてスケーリング係数を決めればよりよい精度を持つモデルを作れるが、本論文はチューニングの手間を削減しつつ精度のよいモデルを作ることが目的であることに注意したい。

f:id:nuka137:20200904095539p:plain

Compound Model Scalingは、探索したスケーリング係数をもとにBaseline Modelをスケールアップすることであるため、Baseline Modelの選択も重要である。従来のCNN（ResNetやMobileNet）をBaseline Modelとし、Compound Model Scalingを適用しても精度の改善はみられた。しかし本論文では、新たに設計した EfficientNet をBaseline Modelとして採用している。

EfficientNetをBaseline Modelに採用し、Compound Model Scalingを適用したモデルEfficientNet-B7では、ImageNetにおいてTop1で84.4%、Top5で97.1%の精度となりSOTAを達成した。また、最高性能を誇っていたCNN（GPipe）と比較してパラメータ数が8.4倍程度少なく、推論では6.1倍高速であった。

EfficientNetの実装

TensorFlow実装（論文に記載の実装）
- https://github.com/tensorflow/tpu/tree/master/models/official/efficientnet
PyTorch実装
- https://github.com/lukemelas/EfficientNet-PyTorch

Compound Model Scaling

限られた計算資源（メモリや演算量）の中で精度を最大にするCNNを探索するためには、次の最適化問題を解くことと同じと言える。

f:id:nuka137:20200908105522p:plain:w400

$N (d, w, r)$ : ネットワーク
$\hat{F_{i}}$ : 演算
$\hat{L_{i}}$ : Baseline Modelで定義したネットワークの層数
$X$ : Baseline Modelに入力されるShapeが $( \hat{H_{i}}, \hat{W_{i}}, \hat{C_{i}}$ である入力テンソル
$d$ : ネットワークの層数
$w$ : ネットワークの幅（Channel数）
$r$ : 入力画像の解像度

各スケーリングパラメータの考察

f:id:nuka137:20200904134009p:plain

ネットワークの層数
- CNNにおいてネットワークの層を増やすことは、より複雑な特徴量を捉えることができて汎化性能の向上につながる一方で、勾配消失問題が発生しやすくなる
- 勾配消失問題を解消するため、ResNetのSkip ConnectionやBatch Normlizationなどの手法が提案されている
ネットワークの幅
- ネットワークの幅を広げることでより細かい特徴を捉えられ、より学習がうまくいくようになる
- ネットワークの幅が広くてネットワークの層数が少ない場合は、複雑な特徴を捉えきれないという問題がある
入力画像の解像度
- 入力画像の解像度を高めることで、より細かい特徴を捉えることができる

いずれのスケーリングパラメータについてもパラメータを増加させることで精度は向上するが、パラメータ増加に伴って精度の向上速度が鈍化する傾向が見られた。

f:id:nuka137:20200904135046p:plain:w400

また、これらのスケーリングパラメータは相互に関連しているものである。例えば、大きな画像を扱う場合はネットワークの層と幅のスケーリング量も大きくする必要がある。実際にスケーリングパラメータ1つだけを固定してモデルをスケールアップさせるよりも、バランスよく各スケーリングパラメータを調整してモデルをスケールアップしたほうが精度が良くなる傾向が見られた。

Compound Model Scalingの導入

ここまでの議論で、各スケーリングパラメータをバランスよくスケーリングすることが重要だとわかった。そこで本論文では、各スケーリングパラメータを程よい具合にスケールアップさせるための式を提案した。これを、Compound Model Scalingと呼んでいる。

f:id:nuka137:20200908105552p:plain:w350

$\phi$ はユーザが指定するスケーリング係数、 $\alpha , \beta , \gamma$ はそれぞれBaseline Modelにおける「ネットワークの層数」「ネットワークの幅」「入力画像の解像度」である。制約条件において $\beta , \gamma$ が2乗されているのは、畳み込み演算が主となるCNNにおいて「ネットワークの幅」と「入力画像の解像度」が演算量に2乗のオーダーで効いてくるためである。

EfficientNet

Compound Model Scalingの効果を正しく評価するためには、Baseline Modelの選択が重要である。本論文ではEfficientNetと呼ばれる新たなモデルを設計し、Baseline Modelとして採用した。

EfficientNetは、NASを利用して精度と演算量を目的関数に設定し、最適化して設計した。最適化の目的関数は、モデル $m$ 、目的演算量 $T$ 、精度と演算量を調整するハイパーパラメータ $w$ として、次のように設定した。

f:id:nuka137:20200908105821p:plain:w200

NASによって最適化されたネットワークをEfficientNet-B0と名付けた。 EfficientNet-B0は、MobileNetV2 で導入されたMobile Inverted Bottleneckを複数段重ねた構成となっている。 Mobile Inverted Bottleneckでは、SENet によるSE Blockも導入している。

f:id:nuka137:20200908132617p:plain:w400

なお、Grid Searchによってスケーリング係数 $\alpha = 1.2, \beta = 1.1, \gamma = 1.15$ が得られた。 EfficientNet-B0のスケーリング係数を固定した状態で $\phi$ の値を変更してBaseline Modelをスケールアップさせ、EfficientNet-B1からEfficientNet-B7までを得た。

実験

Baseline ModelにMobileNetやResNetを利用しCompound Model Scalingを適用すると、従来のように1点（例えばネットワークの深さのみ）を選んでモデルをスケールアップしたときと比較して高い精度が得られることが確認できた。これは、Compound Model Scalingが一般的なCNNに対しても有効であることを裏付けている。

EfficientNetの学習は次の条件のもとで行い、EfficientNet-B7ではTop1で84.4%、Top5で97.1%の精度を達成した。

Optimizer: RMSProp (Decay 0.9, Momentum 0.9)
Batch Normalization: Momentum 0.99
Dropout: $\phi$ に比例して増加（EfficientNet-B0は0.2、EfficientNet-B7は0.5）
Weight Decay: 1e-5
Learning Rate: 初期 0.256, 2.4 epoch毎に0.97減衰
Swish Activation
Fixed AutoAugument Policy
Stochastic Depth: Survival Probability 0.8

f:id:nuka137:20200908141057p:plain:w400 f:id:nuka137:20200908141014p:plain:w400

EfficientNet-B7のパラメータ数は66M、計算量は37BFLOPSである。過去に最高精度を達成したGPipeよりも8.4倍パラメータが少なく、かつ精度が優れている。これらのことから、EfficientNetはパラメータ数と計算量が他のCNNよりも小さいにもかかわらず、精度がよいモデルであると言える。転移学習においてもEfficientNetは、他のCNNを超える精度を達成している。

考察

f:id:nuka137:20200908143507p:plain:w400

EfficientNetにCompound Model Scalingを適用してモデルをスケールアップした場合、従来のように1点（例えばネットワークの深さのみ）を選んでモデルをスケールアップしたときと比較すると、ImageNetのTop-1について2.5%の精度向上が見られた。このことから、精度向上にCompound Model Scalingが効果的であることがわかる。

f:id:nuka137:20200908142419p:plain

Class Activation Map を見てみると、Compound Model Scalingを適用してスケールアップしたモデルは、きちんと関連する領域に焦点を当ててオブジェクトの詳細な特徴まで捉えられていることがわかる。一方で従来のスケールアップ方式では、オブジェクトの詳細まで捉えられていない、もしくはオブジェクトの全体を捉えられていないことがわかる。

【2020.12.11 追記】

MNISTのデータセットを使って、PyTorchにてPre-trainモデルからEfficientNetを学習するソースコードを KaggleのNotebook として公開した。 Notebookの最後では学習済みのEfficientNetのモデルを使用し、手書き数字の画像が正しく分類できていることも示している。