Weblog on mebius.tokaichiba.jp

国名しりとりの最長パターンを探した

子供に、国の名前でしりとりをしようと言われ、国名なんか少ないからすぐ行き詰まるだろうと思いながらやってみると、意外にも20ヶ国くらい繋がった。
それで、最長パターンはどれくらいなのだろう、と思って調べてみた。

■結果
最長しりとり (44ヶ国):
モルディブ → ブルンジ → ジンバブエ → エスワティニ → ニウエ → エジプト → トルコ → コモロ → ロシア → アルジェリア → アルバニア → アルメニア → アイルランド → ドイツ → ツバル → ルーマニア → アイスランド → ドミニカ → カーボヴェルデ → デンマーク → クウェート → トンガ → ガイアナ → ナミビア → アメリカ → カタール → ルクセンブルク → クック諸島 → ウクライナ → ナイジェリア → アラブ首長国連邦 → ウルグアイ → イスラエル → ルワンダ → 大韓民国 → クロアチア → アンドラ → ラオス → スイス → スリランカ → カンボジア → アンゴラ → ラトビア → アルゼンチン

ブログ記事を読む

中心極限定理の証明の理解でつまずいた

2年前に統計検定準1級の勉強をしていた時、どこかで特性関数を用いた中心極限定理の証明を読んで、なるほど、と思ってノートに書き残した。
先週、そのノートを読んだら、理解できない部分があって、結構な時間悩んでしまった。
同じように悩む人はあまり居ないと思うが、自分用にメモしておく。

中心極限定理とは、期待値がμ、分散がσ²である任意の独立同一分布(i.i.d)に従うn個の確率変数X_iの平均が、その分布がどういう形状であっても、nを大きくすると近似的に期待値μ、分散σ²/nの正規分布に従うというものである。

中心極限定理
$\lim_{n \to \infty} P\left(\frac{\overline{X}-\mu}{\sigma / \sqrt{n}} \le x \right) = \int_{-\infty}^{x} \frac{1}{\sqrt{2\pi}} e^{-\frac{t^2}{2}} dt$

証明方法としては、特性関数の連続定理により、ある確率変数X_nの特性関数φ_n(t)がn→∞とすると全てのtでXの特性関数φ(t)に一致するなら、X_nはXに分布収束するというのを使う。
つまり、X_iの平均を平均μ、分散σ²/nで標準化したZ_nの特性関数φ_Zn(t)が、標準正規分布の特性関数φ_Z(t)に一致することを示せば良い。

標準正規分布に従うZの特性関数は
$\begin{align} \phi_Z(t) & = E[e^{itZ}] = \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}} e^{-\frac{z^2}{2}} e^{itz} dz \\ &= \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{-\frac{1}{2}((z-it)^2 + t^2)} dz \\ &= \frac{1}{\sqrt{2\pi}} e^{-\frac{t^2}{2}} \int_{-\infty}^{\infty} e^{-\frac{(z-it)^2}{2}} dz \\ &= e^{-\frac{t^2}{2}} \end{align}$
である。最後は
$\int_{-\infty}^{\infty} e^{-a x^2} dx = \sqrt{\frac{\pi}{a}}$
であることを用いた。

本題の、Znの特性関数を計算する。
$\begin{align} \phi_{Z_n}(t) & = E[e^{itZ_n}] = E\left[\exp \left( it\frac{\sqrt{n}}{\sigma}\frac{1}{n}\sum_{j=1}^n (X_j-\mu) \right)\right] \\ &= E\left[\prod_{j=1}^n \exp\left(it\frac{(X_j-\mu)}{\sqrt{n} \sigma} \right)\right] \\ &= \left( E\left[ \exp\left(it\frac{(X_1-\mu)}{\sqrt{n} \sigma} \right)\right] \right)^n \\ \end{align}$

今回、この最後の変形が何故成り立つのかがわからず、調べまくってしまった。結局、九州大学の原先生の公開講座資料の第26頁の「最後のところでは，Xi の分布が同分布であることを用いました」を見て、X_iがi.i.dであることを思い出して解決した。
上の数式では敢えてX₁と書いた（原先生の資料でもそのようになっている）が、筆者のノートでは（おそらくその元資料でも）この部分がXと書かれていたのが、筆者にとって混乱の元だった。

計算を続ける。マクローリン展開により、
$e^x = 1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\cdots$
なので、2乗の項までを展開して、
$\begin{align} \phi_{Z_n}(t) & = \left( E\left[1 + \frac{it}{\sqrt{n}}\left(\frac{X_1-\mu}{ \sigma}\right) - \frac{t^2}{2n}\left(\frac{X_1-\mu}{ \sigma}\right)^2 - O(\frac{1}{n\sqrt{n}}) \right]\right)^n \\ &= \left(1 - \frac{t^2}{2n} - O(\frac{1}{n\sqrt{n}}) \right)^n \end{align}$
と変形できる。tの1乗の項はE[X_i-μ]=0、2乗の項はE[(X_i-μ)²]=σ²であることを用いた。O(n)は「ランダウの記号」というもので、nを大きくした時に高々f(n)の定数倍となることを表す。3乗の項以降はnを大きくした時に0に収束するので、このように省略している。

$\lim_{n\to\infty} \left(1 + \frac{x}{n} \right)^n = e^x$
なので、
$\begin{align} \lim_{n\to\infty} \phi_{Z_n}(t) & = \lim_{n\to\infty} \left(1 - \frac{t^2}{2n} - O(\frac{1}{n\sqrt{n}}) \right)^n \\ &= \lim_{n\to\infty} \left(1 + \frac{- \frac{t^2}{2} - O(\frac{1}{\sqrt{n}})}{n} \right)^n \\ &= e^{- \frac{t^2}{2}} \end{align}$
となり、標準正規分布と特性関数が一致するので、Znは標準正規分布に収束する。

ブログ記事を読む

「統計学が最強の学問である［数学編］」のロジスティック回帰をやってみた

何年も前に買って読んでいなかった「統計学が最強の学問である」シリーズを、今更ながら読んでいる。その［数学編］はほとんど知っている内容だったので、半分くらいは流し読みしたのだが、所々、特に微分・積分以降はわかっているようでわかっていなかったことを思い知らされて勉強になったし、大変面白かった。

さて、p.490のロジスティック回帰の最急降下法の説明の所に、「気になる方はエクセルや何らかのプログラミング言語を使って挑戦していただければと思いますが、...」と書かれている。気にならない訳がないので、やってみた。

問題は、次のようなデータについて、トラブルの有無がどの程度リピート率に影響するかをロジスティック回帰分析で調べるというものである。

利用時間	トラブル	リピート	人数
ランチ	なし	なし	207
ランチ	なし	あり	23
ランチ	あり	なし	18
ランチ	あり	あり	2
ディナー	なし	なし	435
ディナー	なし	あり	290
ディナー	あり	なし	15
ディナー	あり	あり	10

●コード

import numpy as np

# 元データ(p.483)
data = np.array([
    # ディナーダミー、トラブルダミー、リピートダミー、人数
    [0, 0, 0, 207],
    [0, 0, 1, 23],
    [0, 1, 0, 18],
    [0, 1, 1, 2],
    [1, 0, 0, 435],
    [1, 0, 1, 290],
    [1, 1, 0, 15],
    [1, 1, 1, 10]])

# 各行を人数分に展開したもの
data_flat = np.repeat(data[:, :3], data[:, 3], axis=0)

X = data_flat[:, :2] # ディナーダミー、トラブルダミー
X = np.hstack((np.ones_like(X[:,:1]), X)) # 1列目に全て1の列を追加
y = data_flat[:, 2]  # リピートダミー

# コスト関数 C(β) = -log L(β) = - \sum { (y - 1) x^T β - log(1 + exp(-x^T β)) }
# L(β) = \Pi p^y (1-p)^(1-y), p = 1 / (1 + exp(-x^T β))
def cost_function(X, y, beta):
    return -np.sum((y - 1) * (X @ beta) - np.log(1 + np.exp(-X @ beta)))

# 勾配 ∂C/∂β のテスト(p.490)
beta = np.array([1.0, 1.0, 1.0])
delta = 0.0001
gradient = [
    (cost_function(X, y, beta + [delta, 0, 0]) - cost_function(X, y, beta)) / delta,
    (cost_function(X, y, beta + [0, delta, 0]) - cost_function(X, y, beta)) / delta,
    (cost_function(X, y, beta + [0, 0, delta]) - cost_function(X, y, beta)) / delta,
]

print("first gradient =", gradient)

# 勾配が最小になるβを求める
# βから勾配に学習率を掛けたものを引くのを繰り返す
learning_rate = 0.005

for i in range(120):
    gradient = np.array([
        (cost_function(X, y, beta + [delta, 0, 0]) - cost_function(X, y, beta)) / delta,
        (cost_function(X, y, beta + [0, delta, 0]) - cost_function(X, y, beta)) / delta,
        (cost_function(X, y, beta + [0, 0, delta]) - cost_function(X, y, beta)) / delta,
    ])
    beta -= learning_rate * gradient

print("beta =", beta)

# オッズ比(p.492)
print("odds of beta1, beta2 =", np.exp(beta[1:]))

# scikit-learnのLogisticRegressionで検算する
from sklearn.linear_model import LogisticRegression

X = data_flat[:, :2] # ディナーダミー、トラブルダミー
y = data_flat[:, 2]  # リピートダミー

lr = LogisticRegression().fit(X, y)
print("sklearn beta0 =", lr.intercept_, ", (beta1, beta2) =", lr.coef_[0])
print("sklearn odds of beta1, beta2 =", np.exp(lr.coef_[0]))

●実行結果

first gradient = [523.157877723861, 362.3960970844564, 29.430456188492826] beta = [-2.19686992 1.79119201 0.00347654] odds of beta1, beta2 = [5.99659623 1.00348259] sklearn beta0 = [-2.12298445] , (beta1, beta2) = [ 1.70835681 -0.00918478] sklearn odds of beta1, beta2 = [5.51988382 0.99085727]

betaの1つ目が切片、2つ目と3つ目が回帰係数で、それぞれ小数点第2位を四捨五入すると[-2.2, 1.8, 0.0]、回帰係数をオッズ比にすると[6.0, 1.0]である。p.493の結果と一致した。トラブル有無のオッズ比が1.0なので、リピート率への影響は無しという結果である。
検算の為にscikit-learnのLogisticRegressionでもやってみたら、切片が-2.1、回帰係数の1つ目が1.7と少し違う結果になり、LogisticRegressionのハイパーバラメーターを色々変えてみても本の記載と同じ結果にはできなかった。そういうものなのだろうか。コスト関数の勾配は本の結果の方が小さかったので、本の結果の方が良いと思われる。

ブログ記事を読む