ねこすたっと

ねこの気持ちと統計について悩む筆者の備忘録的ページ。

高校生のためのデータ分析入門 (23):離散型アウトカムの回帰モデル(前編)

数学が苦手なうちのJKに、将来必要となるかもしれないデータ分析への抵抗感をなくしてもらう目的で記事を書くことにしました。

前回:高校生のためのデータ分析入門 (22):稀なイベントの回数ならポアソン分布! - ねこすたっと

おさらい:線形回帰モデルでは上手くいかないケース

以前、線形回帰モデルを説明しました。

necostat.hatenablog.jp

線形回帰モデルは、

  • 観測値 = 期待値+誤差
  • 期待値は説明変数の関数
  • 誤差は正規分布に従う

というモデルでした(線形性や誤差に関する他の仮定など、他にも大事な前提条件があるので、忘れていたら前の記事を振り返ってください)。

誤差が、正規分布という連続型確率分布に従うので、応答変数も連続変数でなければなりませんでしたね。次の例に挙げるような、誤差の分布が正規分布とみなせない離散型の応答変数の場合は、線形回帰モデルはうまく当てはまりません。

  • 試験に合格するかどうか(2値変数*1
  • 学校の売店でたまにしか売れない商品の月販売個数(平均値が大きくない個数・回数)

2値変数をモデル化する

例えば、模擬試験の結果をもとに、3か月後の国家試験に合格するかどうかを予測したいとしましょう。この国家試験は得点開示がなく、合格・不合格しか分かりません。これまでに模擬試験を受けて、本番の結果を教えてくれた1000人のデータを使います。

  • 説明変数X:模擬試験の点数
  • 応答変数Y:1=合格, 0=不合格

期待値を説明変数の関数で表す

線形回帰モデルでは、期待値を説明変数の関数で表しました。今回も、

 
\begin{aligned}
E[Y_i|X_i] &= \beta_0 + \beta_1 X_i
\end{aligned}

としたいところですが、ここで一つ問題が...。

応答変数が2値変数の場合、期待値はY=1のカテゴリーを取る割合です。 つまり、左辺の E[Y|X] は「模試でX点だった人が合格する割合」なので、0から1の範囲しか取りません。 これに対して、右辺  \beta_0 + \beta_1 X_i は、βやXの値によって-∞から+∞まであらゆる値を取る可能性があります。 0から1の範囲から外れた値は、合格割合として解釈できません(下図)。

もちろん、応答変数が連続変数であっても、データによっては期待値がYの取りうる範囲外になることはあります。 下の図では、XもYも全て0以上の値ですが、回帰直線の一部は(Xの取りうる範囲内にも関わらず)0未満になっています。

ただ、このような状況は2値変数のときより遥かに少ないです(2値変数ではほぼ毎回起こる!)。

ロジット変換でつなぐ!

そこで、期待値と  \beta_0 + \beta_1 X を結びつける間に、一工夫はさみます。 見やすさのため、E[Y|X]=pと書き直しました。

 
\begin{aligned}
\log \left( \frac{p_i}{1-p_i} \right) &= \beta_0 + \beta_1 X_i
\end{aligned}

左辺の変換、つまり確率を「あることが起こる確率を起こらない確率で割ったものの自然対数*2」にする変換ロジット変換(logit transformation)と呼びます。下は、 f(p) = \log \left( \frac{p_i}{1-p_i} \right)(ロジット関数) のグラフです。確率pが0に近づくと、ロジットlog(p/1-p)はp=0の線に近付きながら-∞に向かい、pが1に近づくと、ロジットはp=1に漸近して+∞に向かいます。単調増加関数なので、ロジットが1つに決まれば、対応する確率も1つです。

ロジット変換を挟めば、  \beta_0 + \beta_1 X がどんな値になっても、確率として適切な範囲にある値が1つに決まりますね(すばらしい!)

さっきのデータに当てはめると、下のグラフのようになります(確率の軸が縦方向なので、90°向きが変わっていることに注意)。

上の図のS字カーブは、ロジット関数の逆関数(pについて解いたもの)で、ロジスティック関数(logistic function)と呼びます。期待値の式を書き換えると、

 
\begin{aligned}
E[Y_i|X_i]  &= \frac{e^{\beta_0 + \beta_1 X_i}}{e^{\beta_0 + \beta_1 X_i}+1} \\
&= \frac{1}{1+e^{-(\beta_0 + \beta_1 X_i)}} \\
\end{aligned}

となります。

この例のように、「2値変数の期待値にロジスティック関数を当てはめるモデル」をロジスティック回帰モデル(logistic regression model)と呼びます。

誤差は?

「期待値」に関しては説明しましたが、「誤差」についてはどうでしょうか。 ここで、二項分布の特殊形、ベルヌーイ分布の登場です。

例えば、模試の点数が60点だった人の期待値(=合格確率)が90%と推定されたとしましよまう。模試の点数が60点だった人のなかには、無事合格する人もいれば、残念な結果になる人もいます。模試の点数が同じなのに、結果が違うのは全くの運です。 あたかも神様が、当りが90%のクジを引いて、受験者の合否を決めているみたいですね。

線形回帰モデルでは、期待値と実際に観察された値の違いは正規分布で説明されました。 ロジスティック回帰モデルでは、その違いを二項分布(ベルヌーイ分布)で説明しています。

ロジスティック回帰モデルのまとめ

まず、ロジスティック回帰モデルを式で表してみます。

 
\begin{aligned}
\log \left( \frac{p_i}{1-p_i} \right) &= \beta_0 + \beta_1 X_i \\
Y_i &\sim Binom(1, p_i)
\end{aligned}

これを日本語で噛み砕くと、

  • 説明変数  X_i の線形関数として確率  p_i が決まる。つまり、それぞれの人は  X_i に応じて  p_i が与えられる。
  • 説明変数と確率を直接はつなげないので、ロジット変換を介してつなぐ。
  • 応答変数  Y_i は、確率  p_i をパラメータとするベルヌーイ分布に従う。

図で表すとこんな感じでしょうか。

係数βはどう解釈したらいいの?

線形回帰モデルでは、ある説明変数の係数βは「その説明変数が1単位変化したときの期待値の変化量」として解釈できました。

necostat.hatenablog.jp

ロジスティック回帰モデルではどうでしょうか。線形回帰モデルのときと同じように、E[Y|X=x] と E[Y|X=x+1] の差を計算してみましょう。読みやすくなるように、E[Y|X=x+1] = P1, E[Y|X=x] = P0と書くことにします。

 
\begin{aligned}
P_1  &= \frac{1}{1+e^{-(\beta_0 + \beta_1 (x+1))}} \\
P_0 &= \frac{1}{1+e^{-(\beta_0 + \beta_1 x)}} \\
\end{aligned}

数式は整理できましたか?全然キレイな形にできませんね...。

ロジスティック回帰モデルでは、係数βと期待値の変化量を結びつけることが簡単ではありません。代わりに、下のように期待値のロジットを計算してみましょう。

 
\begin{aligned}
\log \left( \frac{P_1}{1-P_1} \right) &= \beta_0 + \beta_1 (x+1) \\
\log \left( \frac{P_0}{1-P_0} \right) &= \beta_0 + \beta_1 x \\
\end{aligned}

第1式から第2式を引くと、下のようになります。

 
\begin{aligned}
\log \left( \frac{P_1}{1-P_1} \right) - \log \left( \frac{P_0}{1-P_0} \right)&= \beta_1 \\
\end{aligned}

右辺はすっきりしましたね。左辺は、log A-log B= log(A/B) を使って変形します。

 
\begin{aligned}
\log \left( \frac{P_1/(1-P_1)}{P_0/(1-P_0)} \right) &= \beta_1 \\
\frac{P_1/(1-P_1)}{P_0/(1-P_0)} &= e^{\beta_1}
\end{aligned}

さて、ここで左辺の分母と分子に登場した P/(1-P)は「あることが起こる確率を、起こらない確率で割ったもの」で、オッズ(odds)と呼ばれます。オッズ起こる確率Pが小さいとき(目安:<10%)は、起こる確率Pとオッズは近い値になります*3

左辺は「X=x+1のときにY=1が起こる確率」と「X=xのときにY=1が起こる確率」の比で、オッズ比(odds ratio)と呼ばれます。 オッズ比が1であれば、2つの発生確率は同じであり、オッズ比が1から離れるほど発生確率に違いがある、ということになります。

係数βを指数変換したもの  e^{\beta} を使えば、「説明変数が1単位変化したときに、オッズが何倍になるか」、すごく平たく言えば「何倍起こりやすくなるか*4」を推定できるわけです。

おわりに

*1:取る値が2種類だけのカテゴリー変数

*2:底eは省略しました。

*3:Pが小さければ、1-P≒1なので。

*4:オッズは発生確率ではないので、この表現は正確ではありませんが、このように言い換えずに「オッズとは何か」を考えるのは難しいです。

*5:読み:もっともらしい

高校生のためのデータ分析入門 (22):稀なイベントの回数ならポアソン分布!

数学が苦手なうちのJKに、将来必要となるかもしれないデータ分析への抵抗感をなくしてもらう目的で記事を書くことにしました。

前回:高校生のためのデータ分析入門 (21):成功?失敗?それなら二項分布! - ねこすたっと

売店でたまにしか売れない商品、いくつ仕入れたらいい?

学校の売店はすごく繁盛してるんですが、月に2個くらいしか売れないキーホルダーが置いてあります。買いたい人が多くて品切れになったら困りますが、そんなに売れないものを置いておくのも無駄なので、毎月4個仕入れているそうです。

月に売れる個数を変数Xとして、確率分布を考えてみましょう。

  • n:店に来た人
  • p:店に来た人がキーホルダーを買う確率

として、前回説明した二項分布を使おうかと思いましたが、「店に来た人」が正確に何人か分かりません。 防犯カメラがないので数えることができないんですが、そもそも次のような人たちをnに含むべきかどうか判断できません。

  • 店の前で立ち止まった人
  • 店の方を覗き込んだ人

このケースだとnが分からないので、二項分布を当てはめるのは難しそうです。

稀なイベントの回数にはポアソン分布!

冒頭の例のように、

  • 対象となるイベントはとても稀なこと(pが非常に小さい)
  • 対象となる集団は大きい(nが非常に大きい)

ようなイベントの回数を考えるときは、ポアソン分布(Poisson distribution)を使って二項分布を近似することができます。

ポアソン分布の確率質量関数は、次の式で表されます。

 
\begin{aligned}
P(X=x) = \frac{\lambda^x e^{-\lambda}}{x!} \ \ \ (\lambda>0, x = 0, 1, 2, ...)
\end{aligned}

ここで、eはネイピア数(自然対数の底)で約2.718です。 二項分布の確率関数と比べると、だいぶ難しいですね...。どうしてこのような式の関数が登場したのかは置いといて、どんな性質の関数なのかを確認していきましょう。

ポアソン分布と二項分布の共通点は、両方とも離散変数の確率分布であることですね。 一方、ポアソン分布と二項分布の違いは、

  • 分布の形状を決めるパラメータ*1 \lambda(ラムダ)1つだけ(二項分布はnとpの2つ)
  • Xが取りうる値の上限がない(二項分布は最大でもn)

ということです。

変数Xがポアソン分布に従うことを、

 
\begin{aligned}
X \sim Poisson(\lambda)
\end{aligned}

とか

 
\begin{aligned}
X \sim Pois(\lambda)
\end{aligned}

とか、単に、

X \sim Po(\lambda)

と書くことがあります(Pois(λ)もPo(λ)もメジャーな書き方ではないかもしれません)。

ポアソン分布の平均と分散

Xがパラメータλのポアソン分布に従うとき、Xの期待値を計算してみます。

 
\begin{aligned}
E[X] &= \sum_{k=0}^{\infty} k \frac{\lambda^k e^{-\lambda}}{k!} \\
&= \sum_{k=1}^{\infty} k \frac{\lambda^k e^{-\lambda}}{k!} \\
&= \sum_{k=1}^{\infty} \lambda \frac{\lambda^{k-1} e^{-\lambda}}{(k-1)!} \\
&= \lambda  \sum_{k=1}^{\infty} \frac{\lambda^{k-1} e^{-\lambda}}{(k-1)!} \\
&= \lambda  \sum_{k'=0}^{\infty} \frac{\lambda^{k'} e^{-\lambda}}{k'!} \\
&= \lambda
\end{aligned}

式変形のヒント:

  • 右辺第1式 → 第2式は、k=0のときは足すものも0なので、k=1から足しても同じだからです。
  • 第5式の  \sum_{k'=0}^{\infty} \frac{\lambda^{k'} e^{-\lambda}}{k'!} は、別のポアソン分布の確率の総和なので1になります*2

次に分散を計算します。まずは E[X2] を求めます。

 
\begin{aligned}
E[X^2] &= \sum_{k=0}^{\infty} k^2 \frac{\lambda^k e^{-\lambda}}{k!} \\
&= \sum_{k=1}^{\infty} \{k(k-1) + k\} \frac{\lambda^k e^{-\lambda}}{k!} \\
&= \sum_{k=1}^{\infty} k(k-1) \frac{\lambda^k e^{-\lambda}}{k!}  + \sum_{k=1}^{\infty} k \frac{\lambda^k e^{-\lambda}}{k!}\\
&= \sum_{k=2}^{\infty} k(k-1) \frac{\lambda^k e^{-\lambda}}{k!}  + \lambda\\
&= \lambda^2 \sum_{k=2}^{\infty} \frac{\lambda^{k-2} e^{-\lambda}}{(k-2)!}  + \lambda\\
&= \lambda^2 \sum_{k'=0}^{\infty} \frac{\lambda^{k'} e^{-\lambda}}{k'!}  + \lambda\\
&= \lambda^2  + \lambda\\
\end{aligned}

式変形のヒント:

  • 第3式で2つの項に分かれます。後ろ側は期待値を求めたときと同じなのでλになります。
  • 前側も同じように \sum_{k'=0}^{\infty} \frac{\lambda^{k'} e^{-\lambda}}{k'!} = 1 を使っています。

これを使うと、分散は

 
\begin{aligned}
Var(X) &= E[X^2] - (E[X])^2 \\
&= (\lambda^2  + \lambda) - \lambda^2 \\
&= \lambda
\end{aligned}

となります。

まとめると、ポアソン分布では平均も分散もλになります。

売れないキーホルダーが在庫切れになる確率

1か月間で売れるキーホルダーの個数をXとします。「月平均で2個売れる」とのことなので、Xがλ=2のポアソン分布*3に従うとしましょう。

 
\begin{aligned}
X \sim Poisson(2)
\end{aligned}

前述の確率質量関数の式にx=0,1,2,...を代入して、売上個数の確率分布を計算してみます。

この表からは、売り上げが5個を超える確率が約5.3%と計算できる*4ので、仕入れ個数4個だとだいたい1年半〜2年に1回くらいは在庫切れになる月がありそうです*5(それくらいなら許容派範囲ですね、多分)。

ちなみに、表をグラフにするとこんな感じです。

おわりに

*1:「母数」とも言う

*2:本当はこれも示さないといけないんですが、高校で習う数学の範囲は超えてしまいます。気になる人はこちら → ポアソン分布の定義と例と性質まとめ | 数学の景色

*3:λ=2のポアソン分布の平均は2だから。

*4:1-Pr(X≦4)より。

*5:1/24=0.042, 1/18=0.056だから。

高校生のためのデータ分析入門 (21):成功?失敗?それなら二項分布!

数学が苦手なうちのJKに、将来必要となるかもしれないデータ分析への抵抗感をなくしてもらう目的で記事を書くことにしました。

前回:

高校生のためのデータ分析入門 (20):期待値と分散・共分散の計算シャワー - ねこすたっと

残り試合でヒットは何本出そう?

確率変数が取りうる値と、その確率を示した関数を確率分布といい、正規分布は連続型確率分布の代表格である、と以前説明しました。

necostat.hatenablog.jp

今回は、離散型確率変数の代表である二項分布について説明します。

例えば、真の打率が3割ちょうどの選手がいるとします。残り5試合(15打数)で3本ヒットを打つことができると、チームの歴代最多安打記録を更新することができます。 「15打数×3割 = 4.5本だから、3本でOKなら確実じゃん!」と思ったあなた、本当にそうでしょうか?

「順位や試合展開によっては敬遠されて…」とか考え出すとキリがないので、問題をシンプルにするために、打数は15で固定しましょう。 ヒットの本数を変数Xとすると、Xの取りうる値xは0〜15ですね*1。ヒットの本数がxである確率  Pr(X=x) は、以下のようになります。

 
\begin{aligned}
P(X=x) = \ _{15} C_x (0.3)^x (0.7)^{15-x}
\end{aligned}

xを0から15まで変えて、  Pr(X=x) を計算してみます。

この表から、3本以上ヒットを打てる確率は約87%と分かりますね。意外に確実でもなさそうです。

ちなみに、さっきの表をグラフにするとこんな感じ。

二項分布とは?

二項分布で扱うのは、

  • コインを投げて表裏のどちらが出るか
  • 引き分けのない試合で勝つか負けるか
  • 1年間で入院することがあるかどうか
  • 事業に成功するか失敗するか

のような「2種類の事象のどちらかを取る試行(ベルヌーイ試行, Bernoulli trial)」です。2種類の事象は必ずどちらかが起こるので、一方の事象(例:成功)が起こる確率をpとすると、他方(例:失敗)が起こる確率は1-pです。「確率pがどちらの事象に注目しているか」を必ず明確にしなくてはいけません。

互いに独立な*2ベルヌーイ試行をn回繰り返したときに観察される成功の数Xが従う確率分布が、二項分布(binomial distribution)です。 変数Xが二項分布を従うことを、

 
\begin{aligned}
X \sim Binom(n, p)
\end{aligned}

あるいは単に、

 
\begin{aligned}
X \sim B(n, p)
\end{aligned}

と書くことがあります。

ベルヌーイ分布:二項分布の特殊形

ベルヌーイ試行を1回だけ行ったときの確率分布ベルヌーイ分布(Bernoulli distribution)と言います。 二項分布でn=1とした場合( B(1,p) )です。

ベルヌーイ分布に従う変数Xは、X=1になる確率がp、X=0になる確率が1-pなので、 期待値は、

 
\begin{aligned}
E[X] &= p \times 1 + (1-p) \times 0\\
&= p 
\end{aligned}

となります。これは直感でも分かりますね。

分散は、

 
\begin{aligned}
Var(X) &= E[X^2] - (E[X])^2 \\
&= (p \times 1^2 + (1-p) \times 0^2) - p^2\\
&= p - p^2 \\
&= p(1-p)
\end{aligned}

となります。これはpの2次関数で、p=0.5のとき最大で、p=0または1のときに最小となりますね。 「分散が大きい」ということは、「Xの値がバラついていて、X=1になるかX=0になるか見当もつかない」ということなので、どっちに転ぶかわからないとき、つまりp=0.5のときに分散が最大になるのも直感的に納得できます。 逆にpが0や1に近いときは、「Xは0あるいは1ばっかりになって、結果はそんなにバラつかない」ということです。

二項分布の平均と分散

前回、期待値と分散の計算を色々確認しました。

necostat.hatenablog.jp

次のように、二項分布はベルヌーイ分布の集合とみなすことができます。 例えば、15打数でのヒット本数Xは、

  •  X_1:1打数目でヒットかどうか
  •  X_2:2打数目でヒットかどうか
    ...
  •  X_{15}:15打数目でヒットかどうか

の合計ですよね。 X_i はそれぞれ平均p, 分散p(1-p)のベルヌーイ分布に従うので、 「和の期待値」を使えば、

 
\begin{aligned}
E[X] &= E[X_1 + X_2 + X_3 + ... + X_{15} ] \\
&= E[X_1] + E[X_2] + E[X_3] + ... E[X_{15}] \\
&= p + p + p + ... + p\\
&= np
\end{aligned}

となります。「15打数×3割 = 4.5本」と見積もったので、これも直感的にわかりますね。

分散もベルヌーイ分布から計算します。 足し合わせる変数が互いに独立ならば、「和の分散は分散の和」なので、

 
\begin{aligned}
Var(X) &= Var(X_1 + X_2 + X_3 + ... + X_{15} ) \\
&= Var(X_1) + Var(X_2) + Var(X_3) + ... Var(X_{15}) \\
&= p(1-p) + p(1-p) + p(1-p) + ... + p(1-p)\\
&= np(1-p)
\end{aligned}

となります。

おわりに

*1:これまでは適当に書いてきましたが、変数の名前は大文字、実際に観察される値は小文字で書いて区別するのが一般的です。

*2:「試行同士が影響しあわず、試行の結果が他の試行の結果と関連しない」ということです。「独立かどうか」という話と「確率pが変わらないかどうか」という話がごちゃ混ぜになりそうなので、続きはまた今度。

高校生のためのデータ分析入門 (20):期待値と分散・共分散の計算シャワー

数学が苦手なうちのJKに、将来必要となるかもしれないデータ分析への抵抗感をなくしてもらう目的で記事を書くことにしました。

前回:高校生のためのデータ分析入門 (19):カテゴリー変数を説明変数に使う - ねこすたっと

はじめに

期待値や分散についての計算が色々出てきたので、一度まとめておきます。特に必要がなければ、結果だけ確認して読み飛ばしてもOKだけど、新しく導入される教科の「情報」では出てくることがあるかもしれないみたいなので、知っておいて損はないかも。

期待値の計算

期待値は、「変数の実現値を、その値を取る確率で重みづけて平均を取ったもの」で、 次のように計算されるんでしたね。

 
\begin{aligned}
E[X] &= \sum_i^n x_i Pr(X=x_i)
\end{aligned}

離散型変数を想定して書きましたが、連続型変数のときは

  •  \sum_i^n → 積分  \int_{-\infty}^{\infty}
  • 確率質量関数  Pr(X=x_i) → 確率密度関数  f(x)

に変換してください。

定数の期待値

値が固定されているので、期待値もその値です。

 
\begin{aligned}
E[c] &= c
\end{aligned}

定数倍の期待値

Xをc倍したものの期待値は、E[X]をc倍したものになります。

 
\begin{aligned}
E[cX] &= \sum_i^n c x_i Pr(X=x_i)\\
&= c \sum_i^n x_i Pr(X=x_i)\\
&= cE[X]
\end{aligned}

和の期待値

和の期待値」は「期待値の和」です。これはXとYが独立でなくても成り立ちます。

 
\begin{aligned}
E[X + Y ] &= \sum_i^n \sum_j^m (x_i+y_j) Pr(X=x_i, Y=y_j)\\
&= \sum_i^n \sum_j^m x_i Pr(X=x_i, Y=y_j) + \sum_i^n \sum_j^m y_i Pr(X=x_i, Y=y_j) \\
&= \sum_i^n x_i \sum_j^m  Pr(X=x_i, Y=y_j) +  \sum_j^m y_i \sum_i^n Pr(X=x_i, Y=y_j) \\
&= \sum_i^n x_i  Pr(X=x_i) +  \sum_j^m y_i Pr(Y=y_j) \\
&= E[X] + E[Y]
\end{aligned}

第3式から第4式の変形で、以下のような計算が出てきます。

 
\begin{aligned}
\sum_j^m  Pr(X=x_i, Y=y_j) &= Pr(X=x_i) \\
\sum_i^n Pr(X=x_i, Y=y_j) &= Pr(Y=y_j)\\
\end{aligned}

これは、下のようなXとYの組み合わせとそれに対応している確率が書かれている表で、  X = x_i (赤)あるいは  Y = y_j (青)に該当する全ての確率を足し合わせることに相当します。

この表で、それぞれのセルには  X = x_i かつ  Y = y_j が同時に起こる確率(= 同時確率, joint probability)が収められています。 例えば、セルをY方向に足し合わせると、Xがある値のときの全てのパターンの確率がもとまります。 これは、表の周辺部分(= 一番右あるいは下)に書いてある、確率の小計を求めることになるので、周辺化(marginalization)と言います。

積の期待値

変数が互いに独立であるときに限り、「積の期待値」は「期待値の積」です。

 
\begin{aligned}
E[XY ] &= \sum_i^n \sum_j^m x_i y_j Pr(X=x_i, Y=y_j)\\
&=  \sum_i^n \sum_j^m x_i y_j Pr(X=x_i) Pr(Y=y_j)\\
&=  \sum_i^n x_i Pr(X=x_i) \sum_j^m  y_j  Pr(Y=y_j)\\
&= E[X] E[Y]
\end{aligned}

「互いに独立である」という条件がないと、 Pr(X=x_i, Y=y_j) = Pr(X=x_i) Pr(Y=y_j)と変形できません。

分散・共分散の計算

分散の計算

下のように、分散は「2乗の期待値 - 期待値の2乗」で計算できます。

 
\begin{aligned}
V(X) &= E[(X - E[X])^2] \\
&= E[X^2] - (E[X])^2
\end{aligned}

前にも取り上げたことがあるので、中間の式変形が気になる人は下の記事を参考にしてください。

necostat.hatenablog.jp

共分散の計算

共分散は「積の期待値 - 期待値の積」で計算できます。

 
\begin{aligned}
Cov(X, Y) &= E[(X - E[X])(Y-E[Y])] \\
&= E[XY - Y E[X] - X E[Y] + E[X]E[Y] ] \\
&= E[XY] - E[Y] E[X] - E[X]E[Y] + E[X]E[Y] \\
&= E[XY] - E[X]E[Y]
\end{aligned}

YをXに置き換えると、「分散 = 2乗の期待値 - 期待値の2乗」になります。 分散は共分散の特殊形とみなせるので、共分散の計算を覚えておけば導出できます。

定数の分散・共分散

定数は揺らぐことがないので、分散は0です。

 
\begin{aligned}
V(c) &= E[c^2] - (E[c])^2 \\
&= c^2 - c^2 \\
&= 0
\end{aligned}

変数と定数との共分散は0です。もちろん定数同士の共分散も0です。

 
\begin{aligned}
Cov(aX, b) &= E[abX] - E[aX]E[b] \\
&= ab E[X] - a E[X]b \\
&= 0
\end{aligned}

定数倍の分散・共分散

係数を2乗して前に出します

 
\begin{aligned}
V(cX) &= E[(cX)^2] - (E[cX])^2 \\
&= E[c^2 X^2] - (c E[X])^2 \\
&= c^2 E[ X^2] - c^2 (E[X])^2 \\
&= c^2 \left( E[X^2] - (E[X])^2 \right) \\
&= c^2 V(X)
\end{aligned}

定数倍した変数同士の共分散は、変数の係数の積を前に出します。

 
\begin{aligned}
Cov(aX, bY) &= E[aXbY] - E[aX]E[bY] \\
&= ab E[XY] - a E[X] b E[Y] \\
&= ab (E[XY] - E[X]E[Y]) \\
&= ab Cov(X,Y)
\end{aligned}

和の分散

和の分散は「分散の和 + 2×共分散」です。何となく(a+b)2の展開の式みたいですね。

 
\begin{aligned}
V(X+Y) &= E[(X+Y)^2] - (E[X+Y])^2\\
&= E[ X^2 + 2XY + Y^2] - (E[X]^2 + 2E[X]E[Y] + E[Y]^2) \\
&= (E[X^2] - E[X]^2) + (E[Y^2] - E[Y]^2) + 2(E[XY] - E[X]E[Y]) \\
&= V(X) + V(Y) +2Cov(X, Y)
\end{aligned}

XとYが独立のときはCov(X,Y)=0なので、シンプルに「和の分散」は「分散の和」になります。

 
\begin{aligned}
V(X+Y) &= V(X) + V(Y)
\end{aligned}

和の共分散も同様に式変形すれば求めることができますが、公式として覚えておくほど登場機会がなさそうなので割愛します。

積の分散

変数が互いに独立であるときに限り、以下のように変形できます。これも同じ要領で展開して変形すると証明できます。

 
\begin{aligned}
V(XY) &= V(X)V(Y) + E[Y]^2 V(X) + E[X]^2 V(Y) \\
\end{aligned}

あんまり使う機会がなさそうなので、詳細は書きませんでした。証明が気になる人は下のリンクを参照してください。

mathwords.net

互いに独立であることが前提なので、共分散についての計算はありません。

おわりに

高校生のためのデータ分析入門 (19):カテゴリー変数を説明変数に使う

数学が苦手なうちのJKに、将来必要となるかもしれないデータ分析への抵抗感をなくしてもらう目的で記事を書くことにしました。

前回:高校生のためのデータ分析入門 (18):曲線を当てはめる - ねこすたっと

カテゴリー変数も数字で扱う

以前、変数には量的変数と質的変数があると説明しました。

necostat.hatenablog.jp

質的変数は対象が属するグループを示す変数で、カテゴリー変数(categorical variable)とも言うんでしたね。 例えば、性別というカテゴリー変数は、

  • 男, 男, 女, 男, 女, ...
  • M, M, F, M, F, ...

などのように、変数の取る値が「グループ名を表す文字・記号」を使って表されます。 しかし、回帰モデルで扱うためには数式に組み込めないといけませんので、カテゴリー変数も数値として扱う必要があります

そこで、男=0, 女=1として、

  • 0, 0, 1, 0, 1, …

のようにします。もちろん、男女は逆でも構いません。 どの数字がどのグループを意味しているのか明確にしていれば、0,1の代わりに、1,2を使っても構いません。 しかし、一般的には0,1を使います。グループを0,1で定義しておけば、変数の平均値が「1が表すカテゴリーの割合」になって便利だからです。

係数はどう解釈したらいいの?

例えば、性別と英語の点数のデータに回帰モデルを当てはめるとします。

  • 説明変数(X):性別(0=男, 1=女)
  • 応答変数(Y):英語の点数

データの分布を下に示します。Xは0か1しか取らないので、普通にプロットすると点が重なってよく分からなくなってしまいます。そこで、ランダムに左右にずらして描き込みました*1。ちなみに、緑の四角はそれぞれのグループの平均値です。

下の式のような線形回帰モデルを当てはめたとしましょう。

 
\begin{aligned}
E[Y|X] &= \beta_0 + \beta_1 X
\end{aligned}

回帰直線を重ねると、下のようになります。

傾き \beta_1 は、Xが1単位増えたときのYの期待値の変化量でしたね。 ここではXは0か1しか取らないので、「X=0のときに比べて、X=1のときには期待値がどれくらい違うのか」を表している、と解釈します。 カテゴリー変数はXが取りうる値が飛び飛びなので、「回帰直線の傾き」というよりも、「期待値の増分(下図)」という方が適切かもしれませんね(見方の違いだけです)。

係数の意味するものを整理すると、

  •  \beta_0:男子(X=0)の平均値
  •  \beta_0 + \beta_1:女子(X=1)の平均値
  •  \beta_1:男子と女子の平均値の差(女子 - 男子)

カテゴリーが3つ以上あるとき

例えば、英語の先生が授業についてアンケートを取ったとしましょう。 アンケートでは、授業のわかりやすさを「普通にわかりやすい・とてもわかりやすい・めちゃめちゃわかりやすい」の中から1つ選んで答えてもらいます。

  • 説明変数(X):アンケートの回答
  • 応答変数(Y):英語の点数

アンケートの回答はカテゴリー変数ですが、次のように数値を割り振ってみましょうか。

  • 0 = 普通に〜
  • 1 = とても〜
  • 2 = めちゃめちゃ〜

データの分布は下のようになっています。

下の式で表される、さっきと同じモデルを当てはめてみましょう。

 
\begin{aligned}
E[Y|X] &= \beta_0 + \beta_1 X
\end{aligned}

各グループの平均値は、

  • 「普通に(X=0)」:  \beta_0
  • 「とても(X=1)」: \beta_0 +  \beta_1
  • 「めっちゃ(X=2)」: \beta_0 +  2 \beta_1

となります。グラフにすると下のとおり。

計算された傾き  \beta_1 は、Xが1単位増えたときのYの期待値の変化量ですが、これはXが0→1のときも、1→2のときも同じになっています。

でも実際のデータは、一定の傾きではないようです。0→1よりも1→2の方が、英語の点数はグッと高くなっていますね。 先程のモデル式では、どう頑張っても0→1と1→2で増分が同じになってしまい、増分の違いは表現できません。 さらに、回帰直線はそれぞれのグループの平均値を通っていません。

「2次関数を当てはめる」というのも1つの作戦かもしれませんが、0→1と1→2で別々の増分を設定したモデルを使う作戦が一般的です。

ダミー変数を使う

そこで次のように、「それぞれのカテゴリーに当てはまるかどうか」を示した2値変数を導入します。

  •  X_1:「とても〜」なら1、それ以外は0を取る
  •  X_2:「めちゃめちゃ〜」なら1、それ以外は0を取る

 X_1 X_2 も0であれば、「普通に〜」を表すことになるので、 X_0 を導入する必要はありません。

このように、k個のカテゴリーを持つ変数は (k-1)個の2値変数の組み合わせで表すことができます。 このときに用いる2値変数をダミー変数(dummy variable)と呼びます。

ダミー変数に変換したデータに対して、下のようなモデルを当てはめます。元のカテゴリー変数Xではなく、2つのダミー変数を説明変数にしたモデルです。

 
\begin{aligned}
E[Y|X] &= \beta_0 + \beta_1 X_1 + \beta_2 X_2
\end{aligned}

このモデルのもとで、それぞれのカテゴリーの平均点は、次のようになります。

  • 「普通に」:  \beta_0
  • 「とても」: \beta_0 +  \beta_1
  • 「めっちゃ」: \beta_0 +  \beta_2

グラフにすると下のとおり。

こうすれば、X=0→1の増分は \beta_1、X=0→2の増分は \beta_2 になるので、0→1と0→2で別々の増分を設定することができました。

例の「普通に〜」のように、全てのダミー変数=0で表されるカテゴリーは、他のカテゴリーの比較対照になります。このように比較の基準となるカテゴリーを参照カテゴリー(reference category)と呼びます。

おわりに

*1:jitteringと言います

高校生のためのデータ分析入門 (18):曲線を当てはめる

数学が苦手なうちのJKに、将来必要となるかもしれないデータ分析への抵抗感をなくしてもらう目的で記事を書くことにしました。

前回:高校生のためのデータ分析入門 (17):複数の説明変数を使った回帰モデル - ねこすたっと

曲線的な関係を当てはめる

これまで、回帰モデルでは直線や平面など、曲がっていない関係性を当てはめてきました。 今回は、「曲線的な関係性」を当てはめてみましょう。 次のような、新興感染症の患者数のデータを例に考えていきます。

  • 説明変数(X):初めて確認されてから経過した日数
  • 応答変数(Y):累積感染者数*1

直線を当てはめてみると、下のようになりました。

「単調増加している」という大まかな傾向はつかめていますが、ズレているところも目立ちますね。20日目は実際よりも回帰直線が大きく上を行っていますし、30日目は逆に回帰直線がだいぶ下に来ています。 このデータには、直線よりも曲線を当てはめる方がズレが少なくて済みそうです。

曲線を使いたいときは、2次関数が一番お手軽です。 2次関数の一般式は y = ax2 + bx + c ですから、期待値のモデル式は次のようになります。

 
\begin{aligned}
E[Y|X] &= \beta_0 + \beta_1 X + \beta_2 X^2
\end{aligned}

直線を当てはめたときに比べると、ズレの具合はだいぶ小さくなりましたね!

線形回帰モデルとは?

これまで紹介してきたモデルは、全て「線形回帰モデル(linear regression model)」と呼ばれるものです。どういうモデルのことなのか、まとめてみます*2

変数の種類

まず、線形回帰モデルは、応答変数が連続変数のときに用いることができるモデルです。 なぜなら、誤差が正規分布に従うことを想定しているからです。正規分布は連続変数の分布の代表格でしたね。

一方、説明変数は連続変数でなくても構いません。カテゴリー変数を説明変数としたモデルについては、次回説明します。

変数の個数

応答変数は1つでなければなりません。応答変数が2つあったら、そのYを決める「期待値+誤差」の式がもう1つ必要になってしまいます。

これに対し、説明変数は(理論的には)いくつあっても構いません(複数の説明変数を含むモデルは前回取り上げましたね)。ただし、サンプル数に比して説明変数の個数が多すぎたり、類似した情報を持った説明変数(例:身長と座高)が含まれていると、推定が上手くいかないことがあります

期待値が「線形式」で表される

「線形」は英語にするとlinear(リニア、直線)なので、線形回帰モデルは「期待値が直線で表されるモデル」と考えたくなりますが、そういうわけではありません。実際、さっき曲線を当てはめた例を出しました。

線形回帰モデルでは、期待値を決める式に係数βがいくつか含まれていますが、それぞれの係数βについては1次式になっています(つまり、係数同士の積はない)。1次式は直線的な関係を表す式ですから、期待値が係数βの1次式で表されているモデルを「線形回帰モデル」と呼ぶのです。

例えば、次の式で表されるモデルは線形回帰モデルではありません(非線形モデル(non-linear model)といいます)。

  • 例1: E[Y|X] = \beta_0 e^{\beta_1 X}
  • 例2: E[Y|X] = \frac{\beta_0 X}{\beta_1 + X}

おまけ:指数関数を当てはめる

累積感染者数のモデル化では、2次関数を使うことで残差は小さくなりましたが、逆に問題が出てきてしまいました(わかりますか?)。

累積感染者数というのは、新しく発生した感染者の数を足していきますが、治った人や亡くなった人の数を引くことはありません。つまり、単調増加する数値なのです。 応答変数の期待値は単調増加するはずなのに、2次関数は単調増加ではありません。実際、さっき当てはめた曲線も前半のところで減少している部分が出てしまっています。

では、単調増加する曲線を当てはめればいい、ということになります。単調増加する曲線の代表は指数関数です。また、指数関数 y = ax は常に正なので、累積感染者数の期待値として負の値を返してしまうことがない、というのもイイところです(1次式、2次式では、期待値がマイナスになってしまっている部分がありました)。

そこで、期待値のモデル式を次のようにしました。

 
\begin{aligned}
E[Y|X] &= c a^X \ \ \ (a>0, c>0)
\end{aligned}

ここで、aはXが増えたときに、Yがどれくらい急峻に増加するかを決めます。 cは、X=0のときのYの期待値、つまりy軸の切片を決めます。

これは先程説明した非線形モデルの例1と同じような式ですね。このままだと線形回帰モデルとして考えることができません。 そこで、両辺の対数を取ってみます。

 
\begin{aligned}
E[\log(Y)|X] &= \log (c a^X) \\
&= \log c + X \log a
\end{aligned}

ここで、 Y^* = \log(Y) で定義される変数Y*を考えます。Yがデータとして収集されていれば、 Y^* もデータとして作ることができますよね。 Y^* のことを「Yを対数変換した変数」などと言ったりします。 見やすいように係数も  \log (c) = \beta_0,  \log (a) = \beta_1 として書き直すと、

 
\begin{aligned}
E[Y^*|X] &= \beta_0 + \beta_1 X
\end{aligned}

となります。こうすれば、見慣れた線形回帰のモデル式にすることができました。 下のグラフに示すように、ピッタリ当てはまっています(データを指数関数をもとにして作成したので当たり前ですが)。

非線形モデルの例1は、対数変換をすることで線形回帰モデルにすることができましたが、例2は線形回帰モデルにできません。

おわりに

*1:「累積」なので、治った人や亡くなった人も数に含まれています。

*2:誤差についての仮定は前回説明したので割愛

高校生のためのデータ分析入門 (17):複数の説明変数を使った回帰モデル

数学が苦手なうちのJKに、将来必要となるかもしれないデータ分析への抵抗感をなくしてもらう目的で記事を書くことにしました。

前回:高校生のためのデータ分析入門 (16):回帰モデルを当てはめてみよう - ねこすたっと

モデルに複数の説明変数を含める

前回は気温とおでんの売り上げの関係を例にして、回帰モデルを当てはめました。このときは、1つの説明変数(=気温)から1つの応答変数(=おでん売り上げ)の期待値を推定しましたが、現実ではもっと色々な要因が影響しているはずです。

例えば下のように、気温のほかに来場者数のデータもあったとしましょう。同じ気温でも、来場者が多い方がおでんは売れそうですよね。

  • 説明変数1( X_1):最高気温(℃)
  • 説明変数2( X_2):来場者数(人)
  • 応答変数( Y):おでんの売り上げ(個)

 Y の期待値が X_1 X_2 によって決まるようにしたいんですが、どうしたらいいでしょうか。 1番簡単な方法は、下のように X_2 の項を足してやる方法です。

 
\begin{aligned}
E[Y|X_1, X_2] &= \beta_0 + \beta_1 X_1 + \beta_2 X_2
\end{aligned}

説明変数が1つのときは直線を表していましたが、説明変数が2つになると平面を表します(下図)。3つ以上のときは図にできませんが、超平面(hyperplane)と呼びます。

係数はどう解釈したらいい?

説明変数が増えても、「期待値がどのように決まるか」を表していることには変わりません。 例えば、気温10℃で来場者が100人のときの売り上げの期待値は、

 
\begin{aligned}
E[Y|X_1=10, X_2=100] &= \beta_0 + 10 \beta_1 + 100 \beta_2
\end{aligned}

であることが推定されます。

傾きは他の変数が固定されている状態で解釈する

では、傾き  \beta_1, \beta_2 そのものはどのように解釈したらいいのでしょうか?

「傾き」は、その説明変数が1単位増えたときに、応答変数がどれくらい変化すると期待されるか、を表しています。説明変数が1つのときは、回帰直線からイメージしやすいですね。

説明変数が2つになると、軸がもう1つ追加されて、直線ではなく平面になります。でも平面の傾きと言われてもピンと来ませんよね。まずは図ではなく数式で考えてみましょう。

例えば、気温が10℃のときに来場者数が100人、あるいは200人だった場合、売り上げの期待値はそれぞれ、

 
\begin{aligned}
E[Y|X_1=10, X_2=100]  &= \beta_0 + 10 \beta_1 + 100 \beta_2  \\
E[Y|X_1=10, X_2=200] &= \beta_0 + 10 \beta_1 + 200 \beta_2
\end{aligned}

となります。2番目の式から1番目の式を引くと、来場者が100人増えると期待値は  100 \beta_2 だけ増えることが分かります。つまり、 \beta_2来場者が1人増加するのに対して、どれくらいおでんの売り上げが増えると期待されるか、を表しています。 気温の話は無視して、来場者の変化の影響だけを考えられるので、とてもシンプルで分かりやすいですね。

ただし、このように解釈できるのは、気温が同じときだけです。気温が違うと、 \beta_1 の項が打ち消し合わずに残ってしまいますもんね。 つまり、傾き  \beta_2 を「売り上げに対する来場者の変化の影響度」として解釈することができるのは、「他の説明変数(ここでは気温)は同じ値に固定されているとしたら」という条件付きなのです。

傾きは一定

さっきは、「気温が10℃のとき、来場者が100人増えると売り上げは  100 \beta_2 増える」と言いました。では、気温が20℃のときはどうでしょうか。

 
\begin{aligned}
E[Y|X_1=20, X_2=100]  &= \beta_0 + 20 \beta_1 + 100 \beta_2  \\
E[Y|X_1=20, X_2=200] &= \beta_0 + 20 \beta_1 + 200 \beta_2
\end{aligned}

書くまでもなかったかもしれませんが、第2式から第1式を引くと、来場者が100人増えると期待値は  100 \beta_2 増える(つまり、 \beta_2 / 人)ので、気温10℃のときと同じになります。 このモデルでは、気温が何℃であれ、来場者1人あたりの売り上げ増加分の期待値は同じになります。

回帰分析の結果は空想上の世界の話

ところで、「他の説明変数が一定だとしたら」という条件は、現実世界で成立するんでしょうか?例えば、気温が低すぎたら来場者は少なくなるかもしれないし、そしたらおでんの売り上げが逆に減るかもしれません。

あるいは、気温が10℃のときも20℃のときも、来場者が100人増えたら同じくらい売り上げが伸びるんでしょうか?20℃のときは来場者がいくら増えても、おでんを買う人は少数派かもしれません。

回帰分析の結果は、

  • 気温が同じのまま、来場者が増えたとすると...
  • 気温が変わっても、来場者増加の影響が一定だとすると...

のような「仮想の世界」の中でデータを眺めた結果に過ぎません。 そもそも、上記のような仮定が妥当かどうかは、どうやって判断したらいいのでしょうか。

1つ目は、経験や知識にもとづいて、分析者が主観的に判断する方法です。「気温が上がっても、来場者増加の効果は一定だろう」と考えれば、先ほどのモデルを使えばいいですが、結果を聞いて解釈する人がその仮定を受け入れてくれなければ、結論に納得してくれないでしょう。

2つ目は、データがモデルの仮定をどれくらい満たしているかを確認する方法です。

データからモデルの妥当性を確認する方法

誤差は直接観察できません*1が、モデルから推測される値(期待値・予測値)と実際に観測された値の差(=残差)を使えば、誤差の様子を推定できます。 そこで、残差が誤差に求められる仮定を満たしているかどうかを確認します。

誤差に求められる仮定について忘れた人は、下の記事を復習してください。

necostat.hatenablog.jp

残差が正規分布しているか

1つ目は、「誤差は正規分布に従う」という仮定です(正規性の仮定)。これは残差のヒストグラムを描くことで確認します。

例えば、左下のようになれば正規分布といって良さそうですが、右下のような場合は仮定が満たされていないと考えます。

残差のバラツキ具合が一定か

2つ目は、「誤差の分散は一定である」という仮定です(等分散性の仮定)。 これは、説明変数Xに対して残差をプロットした散布図で確認します。

例えば、左下はXが変化しても残差のバラツキ具合は一定に見えますが、右下ではXが大きくなるにつれて残差のバラツキも大きくなっているのがわかります。

残差に一定の傾向がないか

3つ目は、「誤差は互いに無関係である」という仮定です(独立性の仮定)。 例えば、観察された順に残差を並べてみて一定の傾向があるか、などで確認します。

例えば、左下は残差の大きさや向きに一定の傾向がなさそうに見えますが、右下では周期的な規則性があるように見えます。

2つ目に説明した「説明変数に対するプロット」でも確認できます。 例えば下の図では、説明変数が小さいうちは残差が負、大きくなると正になっていて、一定の傾向があるように見えます*2

おわりに

  • 回帰分析とは、モデルを使ってデータを解釈することです。常に「モデルが正しいとすれば」という前提条件が付くことを忘れてはいけません。
  • どのようなモデルが正しいかは、神様しか分かりません。さらに言えば、「全てのモデルは間違っている」のです。
  • なぜなら、真実は複雑すぎるので、人間が理解するには簡素化(が必要だからです。データにモデルを当てはめるということは、データを見て「要はこういうことじゃない?」と考えることなのです。
  • モデルの正しさとは別の話ですが、来場者数で売り上げを予測する場合、来場者数はいつ、どうやって知るのでしょう?予想最高気温と違って、前日に知ることは難しそうですね。使うシチュエーションも考慮して回帰モデルを考える必要があります。
  • 次回:高校生のためのデータ分析入門 (18):曲線を当てはめる - ねこすたっと

*1:真の値が分からないので。

*2:正確には、これは残差同士の独立性ではなく、残差と説明変数の独立性の話ですが。