Type I/II/III error

まえがき

先日参加したニューヨークのカンファレンスのキーノートで、あるデータサイエンティストが以下のように話し始めた。

「みんなご存知のとおり、検定には3つのエラーがあります。Type I error、Type II error、Type III errorです。」

あれ、Type III errorって何だっけな？となったので調べてみた。

TLDR

Type I error

• 偽陽性
• 帰無仮説を誤って棄却してしまうこと

Type II error

• 偽陰性
• 帰無仮説を誤って採用してしまうこと

Type III error

• 帰無仮説を誤った理由によって正しく棄却すること

易しい解説

実験によってえられた結果が統計的に正しいかどうかを検証する作業を「検定」と言います。たとえば、「セミのオスとメスの平均寿命は異なる」という主張を証明するために実験をしたとします。

　このとき「セミのオスとメスの平均寿命は異なる」のように実験の結果として主張したい内容を「対立仮説」と呼びます。これに対し、主張したい内容と反対の仮説「セミのオスとメスの平均寿命に差はない」を「帰無仮説」と呼びます。

検定では「帰無仮説」が正しいとした場合に実験結果が観測される確率を計算し、その確率が十分に低い場合は「帰無仮説」が間違っており「対立仮説」が統計的に正しいだろうという結論付けをします。逆に実験結果が観測される確率が十分に低くない場合は、「帰無仮説」を棄却できず「対立仮説」が正しいとは言えないという結論付けをします。

セミの例で各種エラーを説明してみます。
本当はどうなのか分かりませんが、仮に「セミのオスとメスの平均寿命に差はない」というのが真実だったとします。

本当はオス・メスの差はないにも関わらず、たまたま実験結果に差が現れて帰無仮説が棄却されたとします。このような状況を「Type I error」と呼びます。

次に仮に「セミのオスの方がメスより平均寿命が長い」というのが真実だったとして話を進めます。

たまたま実験結果ではオス・メスの寿命に大きな差が見られず、本来棄却されるべき帰無仮説をたまたま棄却できなかったとします。この状況を「Type II error」と呼びます。

また、たまたま実験結果ではメスの方がオスよりも平均寿命が長く、帰無仮説を棄却したとします。帰無仮説を正しく棄却はできましたが、平均寿命の大小関係が真実と逆になっています。このような状況を「Type III error」と呼びます。

単位超球から一様サンプリング

　誰もが一度はモンテカルロ法を使って円周率を計算してみたことがあると思います。二次元の場合は、点が円の中に入ったり、入らなかったりするのですが、次元が大きくなると円の中に入る頻度が極端に少なくなります。
　
　じゃあ、もし高次元の単位超球内に一様分布する点をサンプリングしたくなったらどうすればよいでしょうか？

　以下自分がやってみたこと。

分散を計算するときの桁落ち対策

分散 = 二乗平均 - 期待値の二乗

という式を使って分散を計算すると、右辺の2つの項の値が非常に近い場合、桁落ちが発生する。

たとえば以下のようなデータを考える。

x₁ = 10¹⁰
x₂ = 10¹⁰ + 1

x₃ = 10¹⁰ + 2

x₄ = 10¹⁰ + 3

x₅ = 10¹⁰ + 4

この場合、二乗平均も期待値の二乗も10²⁰程度である。

両者の差は2だが、倍精度浮動小数点の有効桁数は15桁程度なので、この差の情報は失われてしまう。

分散は平行移動しても変わらないので、平均が0になるように中心化を行う。

x₁ = -2
x₂ = -1

x₃ = 0

x₄ = 1

x₅ = 2

すると、二乗平均は2、期待値の二乗は0となり、差の情報は失われない。

以下のプログラムを実行してみると、桁落ちの様子を確認できる。

#include <iostream>
#include <vector>
#include <iomanip>

using namespace std;

// データを中心化する
vector<double> centerize(const vector<double> &v) {
  int n = v.size();
  double avg = 0.0;
  for (auto &x : v)
    avg += x;
  avg /= n;
  vector<double> w(n);
  for (int i = 0; i < n; i++)
    w[i] = v[i] - avg;
  return w;
}

// 分散を計算する
double calc(const vector<double> &v) {
  int n = v.size();
  double avg_sq = 0.0;
  double avg = 0.0;
  for (auto &x : v) {
    avg += x;
    avg_sq += x * x;
  }
  avg_sq /= n;
  avg /= n;

  return avg_sq - avg * avg;
}

// データ生成
vector<double> gen() {
  vector<double> v;
  for (int i = 0; i < 5; i++)
    v.push_back(1e10 + i);
  return v;
}

int main(int argc, char *argv[])
{
  vector<double> v = gen();
  cout << fixed << setprecision(15);
  cout << calc(v) << endl;
  cout << calc(centerize(v)) << endl;
  return 0;
}

以下実行結果。

$ ./main       
0.000000000000000
2.000000000000000

PyStan入門（2）線形モデル

　PyStanで線形モデルのパラメータ推定をしてみた。
データyはNormal(w*x + b, σ²)から生成されるとしてモデルを組んで、パラメータw, b, σの分布を求めてみた。

分かったこと

• pystan.stan()メソッドの引数で以下を指定できる
• iter 各チェーンが行うサンプリング回数（warmup含む）
• warmup 初回の捨てるサンプリング数
• chain チェーンの数
• fit.extract()メソッドでサンプリング値を取れる
• permuted引数=Trueとすると、チェーンをマージしたサンプル列を取得できる
• permuted引数=Falseとすると、（iter, chain, parameters）の次元のnumpy.ndarrayを取得できる

コード

data {
    int<lower=0> N;
    vector[N] x;
    vector[N] y; 
} parameters {
    real w;
    real b;
    real<lower=0> sigma;
} model {
    y ~ normal(w * x + b, sigma);
}

import numpy as np
import pystan
import matplotlib.pyplot as plt
import pandas as pd


def gen():
  N = 30
  w = 10
  b = 3
  sigma = 0.5

  x = np.random.rand(N)
  err = sigma * np.random.randn(N)
  y = w * x + b + err

  return (x,y)


def train(x, y):
  data = {
    'N': len(x),
    'x': x,
    'y': y
  }
  fit = pystan.stan(file='linear_model.stan', data=data, iter=2000, warmup=1000, chains=4)
  return fit


def predict(param, x):
  return param['w'] * x + param['b']


if __name__ == '__main__':
  x, y = gen()
  fit = train(x, y)
  fit.plot()
  plt.show()

  xt = np.linspace(0, 1, 100)
  params = pd.DataFrame({
    'w': fit.extract(permuted=True).get('w'),
    'b': fit.extract(permuted=True).get('b')
  })

  median_params = params.median()
  yt = predict(median_params, xt)
  yt_lower = [np.percentile(predict(params, x), 2.5) for x in xt]
  yt_upper = [np.percentile(predict(params, x), 97.5) for x in xt]

  plt.plot(x, y, 'bo')
  plt.plot(xt, yt, 'k')
  plt.fill_between(xt, yt_lower, yt_upper, facecolor='lightgrey', edgecolor='none')
  plt.show()

結果
まずトレースプロットと、各パラメータの分布。それっぽい分布になっている。

次に学習結果を使って予測値のプロット。実線はパラメータの中央値を使ってyを予測した値。グレーの部分は得られたパラメータのサンプルを使って予測値を計算したときの2.5%パーセンタイルと97.5パーセンタイルの値以内の領域。

PyStan入門（1）とりあえず動かす

　緑本を読んで、MCMCで統計モデルのパラメータ推定をすることに興味がわいてきたので、PyStanをさわっていこうと思う。
　まずは、とりあえずインストールして動かしてみました。

PyStanのインストール
pipで入ります。

$ pip install pystan

モデルの記述
stanのモデルは.stanという拡張子のファイルに書くのが一般的なようです。他にもPythonのコードにベタ書きするという方法もあるようです。
ここではstanファイルにモデルを記述することにし、以下の内容を8schools.stanというファイルに記述します。

data {
  int<lower=0> J; // number of schools 
  real y[J]; // estimated treatment effects
  real<lower=0> sigma[J]; // s.e. of effect estimates 
}
parameters {
  real mu; 
  real<lower=0> tau;
  real eta[J];
}
transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] <- mu + tau * eta[j];
}
model {
  eta ~ normal(0, 1);
  y ~ normal(theta, sigma);
}

上のモデルを以下のPythonコードから参照します。以下のコードをsample.pyという名前のファイルに記述します。

import pystan
import matplotlib.pyplot as plt

schools_dat = {
 'J': 8,
 'y': [28,  8, -3,  7, -1,  1, 18, 12],
 'sigma': [15, 10, 16, 11,  9, 11, 10, 18]
}

fit = pystan.stan(file='8schools.stan', data=schools_dat, iter=1000, chains=4)

print(fit)
fit.plot()
plt.show()

実行結果
実行してみます。

$ python sample.py
DIAGNOSTIC(S) FROM PARSER:
Warning (non-fatal): assignment operator <- deprecated in the Stan language; use = instead.

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_286b3180dfa752c4cfedaf0241add0e4 NOW.
Iteration:   1 / 1000 [  0%]  (Warmup) (Chain 1)
Iteration:   1 / 1000 [  0%]  (Warmup) (Chain 3)
Iteration:   1 / 1000 [  0%]  (Warmup) (Chain 2)
Iteration:   1 / 1000 [  0%]  (Warmup) (Chain 0)
Iteration: 100 / 1000 [ 10%]  (Warmup) (Chain 1)
Iteration: 100 / 1000 [ 10%]  (Warmup) (Chain 2)
Iteration: 100 / 1000 [ 10%]  (Warmup) (Chain 3)
Iteration: 100 / 1000 [ 10%]  (Warmup) (Chain 0)
Iteration: 200 / 1000 [ 20%]  (Warmup) (Chain 1)
Iteration: 200 / 1000 [ 20%]  (Warmup) (Chain 0)
Iteration: 200 / 1000 [ 20%]  (Warmup) (Chain 3)
Iteration: 300 / 1000 [ 30%]  (Warmup) (Chain 0)
Iteration: 300 / 1000 [ 30%]  (Warmup) (Chain 3)
Iteration: 300 / 1000 [ 30%]  (Warmup) (Chain 1)
Iteration: 200 / 1000 [ 20%]  (Warmup) (Chain 2)
Iteration: 400 / 1000 [ 40%]  (Warmup) (Chain 1)
Iteration: 400 / 1000 [ 40%]  (Warmup) (Chain 0)
Iteration: 400 / 1000 [ 40%]  (Warmup) (Chain 3)
Iteration: 500 / 1000 [ 50%]  (Warmup) (Chain 3)
Iteration: 501 / 1000 [ 50%]  (Sampling) (Chain 3)
Iteration: 300 / 1000 [ 30%]  (Warmup) (Chain 2)
Iteration: 500 / 1000 [ 50%]  (Warmup) (Chain 0)
Iteration: 500 / 1000 [ 50%]  (Warmup) (Chain 1)
Iteration: 501 / 1000 [ 50%]  (Sampling) (Chain 0)
Iteration: 501 / 1000 [ 50%]  (Sampling) (Chain 1)
Iteration: 600 / 1000 [ 60%]  (Sampling) (Chain 3)
Iteration: 600 / 1000 [ 60%]  (Sampling) (Chain 0)
Iteration: 400 / 1000 [ 40%]  (Warmup) (Chain 2)
Iteration: 600 / 1000 [ 60%]  (Sampling) (Chain 1)
Iteration: 700 / 1000 [ 70%]  (Sampling) (Chain 3)
Iteration: 700 / 1000 [ 70%]  (Sampling) (Chain 1)
Iteration: 500 / 1000 [ 50%]  (Warmup) (Chain 2)
Iteration: 501 / 1000 [ 50%]  (Sampling) (Chain 2)
Iteration: 700 / 1000 [ 70%]  (Sampling) (Chain 0)
Iteration: 800 / 1000 [ 80%]  (Sampling) (Chain 3)
Iteration: 800 / 1000 [ 80%]  (Sampling) (Chain 1)
Iteration: 900 / 1000 [ 90%]  (Sampling) (Chain 3)
Iteration: 600 / 1000 [ 60%]  (Sampling) (Chain 2)
Iteration: 800 / 1000 [ 80%]  (Sampling) (Chain 0)
Iteration: 1000 / 1000 [100%]  (Sampling) (Chain 3)
# 
#  Elapsed Time: 0.044725 seconds (Warm-up)
#                0.037462 seconds (Sampling)
#                0.082187 seconds (Total)
# 
Iteration: 900 / 1000 [ 90%]  (Sampling) (Chain 1)
Iteration: 900 / 1000 [ 90%]  (Sampling) (Chain 0)
Iteration: 700 / 1000 [ 70%]  (Sampling) (Chain 2)
Iteration: 1000 / 1000 [100%]  (Sampling) (Chain 1)
# 
#  Elapsed Time: 0.05144 seconds (Warm-up)
#                0.040057 seconds (Sampling)
#                0.091497 seconds (Total)
# 
Iteration: 1000 / 1000 [100%]  (Sampling) (Chain 0)
# 
#  Elapsed Time: 0.04517 seconds (Warm-up)
#                0.046124 seconds (Sampling)
#                0.091294 seconds (Total)
# 
Iteration: 800 / 1000 [ 80%]  (Sampling) (Chain 2)
Iteration: 900 / 1000 [ 90%]  (Sampling) (Chain 2)
Iteration: 1000 / 1000 [100%]  (Sampling) (Chain 2)
# 
#  Elapsed Time: 0.053818 seconds (Warm-up)
#                0.040396 seconds (Sampling)
#                0.094214 seconds (Total)
# 
Inference for Stan model: anon_model_286b3180dfa752c4cfedaf0241add0e4.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

           mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
mu         7.81    0.17   5.09  -2.41   4.61   7.77  11.13  17.63    907    1.0
tau        6.63    0.24    5.6   0.29   2.45   5.33   9.17  21.19    531    1.0
eta[0]     0.42    0.02   0.94  -1.53  -0.17   0.44   1.05   2.17   1701    1.0
eta[1]     0.01    0.02   0.87  -1.67  -0.55 3.2e-3    0.6   1.74   1751    1.0
eta[2]    -0.23    0.02   0.89  -1.95  -0.81  -0.22   0.38   1.51   1866    1.0
eta[3]    -0.01    0.02    0.9  -1.82  -0.61-8.1e-3   0.58   1.79   1669    1.0
eta[4]    -0.33    0.02   0.87   -2.1  -0.87  -0.34   0.21   1.44   1603    1.0
eta[5]    -0.21    0.02    0.9   -1.9  -0.82  -0.23   0.39   1.55   1599    1.0
eta[6]     0.34    0.02    0.9  -1.42  -0.24   0.36   0.95   2.07   1537    1.0
eta[7]     0.06    0.02   0.94  -1.84  -0.56   0.07   0.66   1.87   1604    1.0
theta[0]  11.61    0.24   8.19  -1.43   6.15  10.55  15.65  31.15   1211    1.0
theta[1]   7.84    0.14   6.35  -5.56   4.05   7.76   11.7  20.83   1988    1.0
theta[2]   5.97     0.2    7.7  -11.7   1.85   6.38  10.91  20.08   1505    1.0
theta[3]   7.61    0.15   6.45  -5.51   3.46   7.67  11.86  20.35   1897    1.0
theta[4]   5.13    0.14   6.45  -8.65   0.99   5.58   9.56   17.0   2000    1.0
theta[5]   6.21    0.17   6.51  -7.34   2.29   6.46  10.55  18.69   1522    1.0
theta[6]  10.64    0.17   6.88  -1.89   6.09  10.06  14.68  26.06   1623    1.0
theta[7]   8.33     0.2   7.94  -7.88   3.77   8.07  12.54  25.84   1522    1.0
lp__      -4.82    0.11   2.66 -10.43  -6.54  -4.59  -2.84  -0.42    598    1.0

Samples were drawn using NUTS at Sat Apr  1 00:55:24 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

C++のコンパイラが走って、MCMCのサンプリングが走って、パラメータの分布が出たようです。続いて各パラメータの分布図とiterationごとの変化がプロットされます。

細かいことは置いといて、何やら楽しそうなことが出来るということは分かりました。
次回は自分で簡単な乱数生成モデルを作って、それをモデリングしてパラメータ推定ができるかどうかを試してみたいと思います。

参考
Getting started — PyStan 2.14.0.0 documentation

緑本 第11章 空間構造のある階層ベイズモデル

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。ようやく最終章までたどり着いた！第11章「空間構造のある階層ベイズモデル」を読んだので、理解度チェックリストを作っておく。

• 空間相関とは？
• どのようなときに考慮しないといけない？
• どのように表現できる？
• 条件つき自己回帰モデル（conditional auto regressive, CAR）とは？
• 確率場とは？
• マルコフ確率場とは？
• ガウス確率場とは？
• 観測データに欠測がある場合、空間相関を考慮する・しないで予測結果はどのように変わる？

緑本 第10章 階層ベイズモデル -GLMMのベイズモデル化-

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。今日は第10章「階層ベイズモデル -GLMMのベイズモデル化-」を読んだ。

　そもそもこの本を読み始めたのは、MCMCをどのようにベイズに適用するのかを理解するのが目的だったのだが、この章に来てようやくイメージができた気がする（9章からぼんやりわかり始めていたものが10章でしっかりと分かった）。この本を読む前までは、MCMCというのを最適化の手法のように勘違い（おそらくメトロポリス法が最適化手法の一つであるsimulated annealingに似ていたのでそう思ってしまったのかも）していた気がする。つまり最尤推定でやっているようにMCMCで対数尤度関数を最適化するパラメータを探すのだと勘違いしていて、それってどうやるのか全然わからない、っていうかわざわざそんなことする意味あるの？と悩んでいたのである。
　しかし、MCMCを「ベイズ」に適用するので、やろうとしていることはパラメータの点推定ではなく、パラメータの分布推定である。パラメータの分布が分かれば、特定パラメータの周辺分布や、応答変数の信頼区間や、ある説明変数が応答変数に影響するかどうかなどの統計検定ができる。こういうことをやるためにMCMCを使うというのが分かると、あぁなるほどな〜と使い方（ベイズのモデルをMCMCで計算するという意味）をイメージできるようになった。

　まえおきが長くなってしまったが、いつもの理解度をチェックするための確認リストを作っておく。

• 階層ベイズモデルとは？
• GLMMだと何が難しい？
• 階層事前分布とは？
• 超事前分布とは？
• ベイズ統計モデルで主観的な事前分布の使用を避けたい場合は、どのような分布を使えばよいか？
• 無情報事前分布と階層事前分布のどちらを使うべきかはどのように判断すればよいか？
• 大域的なパラメータと局所的なパラメータ

緑本 第9章 GLMのベイズモデル化と事後分布の推定

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。第9章「GLMのベイズモデル化と事後分布の推定」を読んだので理解度チェックリストを作っておく。

• 無情報事前分布とは？
• 一様分布以外にはどのようなものが使える？
• データの中央化とは？
• データの標準化とは？
• MCMCでサンプリングする際に設定すべき値は？
• どのようにチューニングするのがよいか？
• 周辺事後分布とは？
• ギブスサンプリングのアルゴリズムについて説明してください
• 全条件つき分布とは？
• メトロポリス法と比較したときの利点は？
• 共役事前分布とは？

緑本 第8章 マルコフ連鎖モンテカルロ（MCMC）法とベイズ統計モデル

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。8章「マルコフ連鎖モンテカルロ（MCMC）法とベイズ統計モデル」を読んだので、自分の理解度確認のためのチェックリストを作っておく。

• 統計モデルのパラメータqの最尤推定値を求めたいとします。qを解析的に求められないときにはどのような方法を使えばよいか？原始的な方法を説明してください。
• メトロポリス法のアルゴリズムを説明してください。
• MCMCとは？
• どの部分がマルコフ連鎖なのか？
• どの部分がモンテカルロ法なのか？
• MCMCの目的は？
• 何を出力する？
• どういうときに重宝する？
• MCMCの定常分布は何を表す？
• 頻度主義とベイズ統計学の違いについて説明せよ。
• モデリングにデータを当てはめた結果えられるものはそれぞれ何か？

緑本 第7章 一般化線形混合モデル -個体差のモデリング-

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。今日は、第7章「 一般化線形混合モデル -個体差のモデリング-」を読んだ。7章からは本の後半となり、より応用的な内容になる。
　以下は自分の理解度をチェックするために作った確認リスト。

• 過分散とは？
• 過分散が起こる原因は？
• 一般化線形混合モデル（GLMM）とは？
• どのようなときに使うモデル？
• 一般化線形混合モデルの最尤推定はどのように行う？
• 個体差の分布はどのように表すか？
• 混合分布とは？
• データの取り方において、反復、疑似反復の違いは？
• GLMMが必要なのはどっち？

緑本 第6章 GLMの応用範囲をひろげる -ロジスティック回帰など-

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。今日は第6章「GLMの応用範囲をひろげる -ロジスティック回帰など-」を読んだ。
　以下は自分の理解度をチェックするために作った確認リスト。

• 以下のようなデータをGLMで分析したい場合に、確率分布は何を用いるのがよさそうか？
• 上限のないカウントデータ
• 上限のあるカウントデータ
• 確率変数のとりうる範囲が0以上の連続値データ
• 負の値もとりうる場合の連続値データ
• 以下の文のX, Y, Zに当てはまる語句は何か？
• ロジスティック回帰は確率分布関数に(X)を、リンク関数に(Y)関数を用いたGLMである。
• (Y)関数は(Z)関数の逆関数である。
• オッズ比とは何か？
• ロジスティック回帰モデル（相互作用項無し）において、ある説明変数xが1増加した場合、オッズ比はどのように変化するか？
• 交互作用項とは何か？
• 交互作用項をいれる場合の注意点は？
• 観測データの割算値を使って回帰を行う場合の問題点は？
• オフセット項とは何か？そのメリットは？

緑本 第5章 GLMの尤度比検定と検定の非対称性

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。今日は第5章「GLMの尤度比検定と検定の非対称性」を読んだ。
　以下は自分の理解度をチェックするために作った確認リスト。

• 帰無仮説とは？
• 対立仮説とは？
• 有意水準とは？
• Neyman-Pearsonの検定の枠組みとは？
• 尤度比検定について説明してください。
• 第1種の過誤と第2種の過誤の違いは？
• Neyman-Pearsonの検定の枠組みではどちらに専念する？
• P値とは？
• 以下の尤度比検定を行うための手法について説明してください。
• パラメトリックブートストラップ法
• カイ二乗分布を使った近似
• 「帰無仮説が棄却できなかったので、帰無仮説は正しい」という主張を否定してください。
• Neyman-Pearsonの検定の枠組みの非対称性
• 検定力とは？
• 検定とモデル選択の目的の違いを説明してください。

緑本 第4章 GLMのモデル選択 -AICとモデルの予測の良さ-

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。今日は第4章「GLMのモデル選択 -AICとモデルの予測の良さ-」を読んだ。

　以下は自分の理解度をチェックするために作った確認リスト。

• モデル選択において、「線形モデルの場合は決定係数R²は高ければ高いほどよい。」という主張を否定してください。
• 逸脱度とは何でしょうか？
• 最大対数尤度を使って式で書いてみてください。
• フルモデル、ヌルモデルについて説明してください。
• それぞれのパラメータ数
• 最小逸脱度、最大逸脱度
• 残差逸脱度とは何でしょうか？
• AICとは何でしょうか？
• 逸脱度Dとパラメータ数kを使って式で書くとどうなる？
• 平均対数尤度とは何でしょうか？
• AICでモデル選択をしていい理由を以下の語句を使って説明してください。
• 最大対数尤度
• 平均対数尤度
• バイアス補正

緑本 第3章 一般化線形モデル（GLM）- ポアソン回帰 -

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。今日は第3章「一般化線形モデル（GLM）- ポアソン回帰 - 」を読んだ。

　以下は自分の理解度をチェックするために作った確認リスト。

• ポアソン回帰について説明せよ。
• どのようなデータに対して適用するのが自然か？
• モデルを式で書くとどうなる？
• 一般化線形モデルを用いる際に指定する3つの項目を答えよ。
• GLMにおいて因子型の説明変数はどのように扱えばよいか説明せよ。
• 対数リンク関数と恒等リンク関数の違いについて説明せよ。
• 各説明変数が応答変数に与える影響はどうなるか？
• 直線回帰と一般化線形モデルの関係について説明せよ。
• 直線回帰による限界を具体例を用いて説明せよ。
• 直線回帰が前提としている2つの仮定は何か？

緑本 第2章 確率分布と統計モデルの最尤推定

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。今日は第2章「確率分布と統計モデルの最尤推定」を読んだ。

　以下は自分の理解度をチェックするために作った確認リスト。

• 確率分布とは何か？
• 以下の確率分布はどのような場面で使うのが適切か？
• ポワソン分布
• 二項分布
• 正規分布
• ガンマ分布
• 尤度とは何か？
• 以下のデータはポワソン分布にしたがってランダムに抽出されたものとする。
• 4, 3, 4, 2, 2, 2, 1, 3, 6, 3
• 対数尤度を最大化するパラメータと標本平均が一致することを確認せよ。（ポワソン分布の式をそらで書けないようだとダメ）
• 標準誤差とは何か？
• 推定と予測の違いについて説明せよ。

緑本 第1章 データを理解するために統計モデルを作る

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読み出した。自分自身の理解度チェックのために各章ごとに確認リスト的なものを作ることにした。確認リストが質問形式なのは何となく詰め寄られてる気がして緊張感が出るから。

　今日は第1章「データを理解するために統計モデルを作る」を読んだ。

• 統計モデルとは何でしょうか？
• 応答変数と説明変数について説明してください。
• リンク関数について説明してください。
• 統計モデルのどのあたりが「統計」なモデルなのか説明してください。
• 統計学の誤用あるあるネタを挙げてください。

EMアルゴリズムの分かりやすい資料

　EMアルゴリズムの分かりやすい資料を見つけた。半年くらい前に混合ガウス分布を例にした資料をもとに勉強してだいたい理解したつもりでいたが、いまいちしっくりきていない部分があったのでもう一度別の資料で勉強しなおした。

What is the expectation maximization algorithm?

これを読むと以下のことが理解できる。EMアルゴリズムの復習をしたいときは読み直すとよい。

• やりたいことのイメージ
• コイントスを用いたシンプルな例
• Eステップ、Mステップの名前の由来
• Eステップでは、潜在変数の確率分布を求める。実際には分布そのものではなく、十分な情報をもつ統計量の期待値を求める。
• Mステップでは、対数尤度を最大化するようにモデルを再構築する。
• EMアルゴリズムの応用例
• バイオインフォマティクス領域でよく使われるらしい
• DNAシーケンスのクラスタリング（多次元正規分布）
• よくあるEMアルゴリズム更新のイメージ図の理解（*注）
• Eステップで出した凹関数が求めたい真の関数の下限になっているのは何故か？
• Eステップで出した凹関数と真の関数が接するのは何故か？
• EMアルゴリズムでなぜ局所解に到達するのか？

注) 数学的に厳密な説明は載ってないが、イメージ的に分かる。

関数の括弧の中のセミコロン

　関数の定義で f(x; θ) のような表記があった。

xとθの関数なので f(x, θ) でいいのではと思ったが、セミコロンで区切られている。
セミコロンで区切ると、セミコロンの前は変数、セミコロンの後はパラメータという意味らしい。

f(x; θ)はパラメータθによって特徴付けられる、xを変数として受け取る関数fみたいなイメージ。
θは変数ではなく、fを定義するときに決定する固定値というふうに考える。

信頼区間付きの線形回帰予測

　pythonで線形回帰の信頼区間付きpredictionをしたかったけど、sklearnには該当する機能が無かった。Gradient Boosting Regressionだと出来るのだけど[2]、単純なLeast Squareでやりたかったのでどうやるのか調べていたら、stasmodelsを使えばいい[1]らしい。

サンプルソース

import numpy as np
import statsmodels.api as sm
from statsmodels.sandbox.regression.predstd import wls_prediction_std
import matplotlib.pyplot as plt

np.random.seed(1234)
nsample = 30
x1 = np.linspace(0, 20, nsample)
x = np.column_stack((x1, x1**0.5))
x = sm.add_constant(x)
beta = [1, 0.5, 0.5]
y = np.dot(x, beta) + np.random.normal(size=nsample)
res = sm.OLS(y, x).fit()
prstd, iv_l, iv_u = wls_prediction_std(res)

plt.plot(x1, y, 'o', label="Data")
plt.plot(x1, res.fittedvalues, 'r--')
plt.plot(x1, iv_u, 'r--', label="Ordinary Least Squares")
plt.plot(x1, iv_l, 'r--')
plt.legend(loc="best");
plt.show()

wls_prediction_stdは、学習データに対する予測値の有意水準5%信頼区間を返してくれる。有意水準は引数alphaで自由に設定できる。また、学習データ以外のデータに対する予測を行いたい場合は、引数exogにデータを渡せばよい[3]。

結果

参考
[1] https://www.reddit.com/r/MachineLearning/comments/3raivl/code_to_calculate_confidence_interval_for_linear/
[2] http://scikit-learn.org/stable/auto_examples/ensemble/plot_gradient_boosting_quantile.html
[3] https://github.com/statsmodels/statsmodels/blob/master/statsmodels/sandbox/regression/predstd.py

ウィルコクソンの順位和検定

　2群のデータの差を検定するノンパラメトリック手法。以下のような場面で使う。

• データが順位データの場合（5:とても満足、4:満足、3:普通、2:不満、1:とても不満のようなリッカート尺度など）
• データが量的データだが、正規性を仮定できない場合

例題1

東京と大阪の住人を無作為に8人ずつ選んで、内閣支持に関するアンケートを行った。

5: 支持する、4:やや支持する、3:どちらとも言えない、2:やや支持しない、1:支持しない

以下の結果から、東京と大阪で内閣に対する意見は異なっているか検定せよ。

東京	大阪
4	3
3	3
3	1
2	5
1	2
3	5
4	1

データがリッカート尺度で与えられるため、ウィルコクソンの順位和検定を用いる。

> wilcox.test(c(4,3,3,2,1,3,4), c(3,3,1,5,2,5,1))

 Wilcoxon rank sum test with continuity correction

data:  c(4, 3, 3, 2, 1, 3, 4) and c(3, 3, 1, 5, 2, 5, 1)
W = 25.5, p-value = 0.9475
alternative hypothesis: true location shift is not equal to 0

p値 = 0.9475なので、 帰無仮説は棄却できない。
よって東京・大阪で有意な差があるとは言えない。

例題2
ある新薬の効果を試すために被験者8名を集めた。
まず8名のγ-GTP値を計測し、その後2ヶ月新薬を服用してもらい、2ヶ月後γ-GTP値を再計測した。
結果は以下のとおりである。
新薬に効果があるかどうかを検定せよ。

被験者番号	投薬前	投薬後
1	100	95
2	110	100
3	120	105
4	130	110
5	135	120
6	140	145
7	220	170
8	250	195

まず検定対象の分布に正規性があるかどうか確認する。

> shapiro.test(c(100, 110, 120, 130, 135, 140, 220, 250))

 Shapiro-Wilk normality test

data:  c(100, 110, 120, 130, 135, 140, 220, 250)
W = 0.81118, p-value = 0.03766

すると有意水準5%で投薬前のデータは正規性を持たないことが分かる。

正規性が仮定できないため、パラメトリック手法であるt検定は使うことができない。
よって、ノンパラメトリック手法であるウィルコクソンの順位和検定で検定する。
今回は同じ被験者の投薬前と投薬後の比較をするため、pairedオプションをTrue（群間の対応あり）にする。

> x_before <- c=""> x_after <- c=""> wilcox.test(x_before, x_after, paired=T)

 Wilcoxon signed rank test with continuity correction

data:  x_before and x_after
V = 34.5, p-value = 0.02471
alternative hypothesis: true location shift is not equal to 0

上記の結果から、有意水準5%で新薬に効果ありと言える。