盆暗の学習記録

データサイエンスを中心として,日々学んだことの備忘録としていく予定です。初心者であり独学なので内容には誤りが含まれる可能性が大いにあります。

手元の本を電子書籍にした

就職を機に都心の狭くて安い物件に引っ越すので荷物を減らそうと努力しています。

私が持っている本棚の使い勝手がイマイチで,専有する場所の割には本があまり入らない(文庫本向けっぽくて専門書があまり入らない。未使用の空間が多く生まれる)ので,一つの策として本棚を捨てることを考えています。

そこで紙の本を「自炊」(自分で電子書籍化を行う)で電子書籍化しました。

そのあたりの意思決定に関して記しておきます。

なぜ自炊なのか。外注できないのか

自炊には時間的にも金銭的にもコストがかかります。

できれば電子書籍化代行サービス」のようなのを利用したいところですが,著作権法上そうしたサービスは違法なようです1

自炊して自分で読む分には適法のようなので,泣く泣く自炊用の設備を購入してきて自炊することにしました。

自炊のために揃えた設備

裁断機

買ったもの

個人用の裁断機のなかではかなりハイグレードと思われるこれにしました。

DURODEX 自炊裁断機 ブラック 200DX

DURODEX 自炊裁断機 ブラック 200DX

これに決めた理由:

  • 本のスライスと裁断が最も時間のかかる作業であるのでスペックが高いことは時短につながる
  • 一度に裁断できる枚数が多いほどスライス回数が少なく済んでよい
  • 小型の裁断機だと真っ直ぐ切れないものもあるらしいが,これは構造がしっかりしていて裁断精度が高い

例えば『マンキューミクロ経済学』は750ページ近くあるんですが,200DXなら一冊を2~3個のパートにスライスしてそれぞれ裁断できます。

もし,もっと安くてスペックの低いもの(例えば裁断枚数が50)なら8個のパートにスライスしてからそれぞれ裁断することになるでしょうから,それを考えるとだいぶ手間を省略できます。

大きくて置き場所に困りそうでちょっと怖いんですが,でも裁断の安定性は高いです。安全性もしっかりしています。

買わなかったもの

もう少し安くてコンパクトなものもあります。

買おうかちょっと迷ったのは次の2つです

コンパクトで,安くて,40枚を裁断できます。「うっかりしているとずれて裁断してしまう」といった裁断精度に関するレビューもあったので購入は見送りました。

プラス 裁断機 コンパクト 断裁機 PK-213 裁断幅 A4縦 26-366

プラス 裁断機 コンパクト 断裁機 PK-213 裁断幅 A4縦 26-366

コンパクトで,比較的まっすぐ切れそうですが,60枚しか切れない割に高い。

買えば良かったかも,と思ってるもの

今この記事を書いてて知ったのですが,裁断枚数400の

こういう裁断機もあるようです。非常にコスパが良さそうです。

これ買えばよかったかもしれません。

200DXよりも大きそうなので置き場所に困りそうですが…。

スキャナー

買ったもの

自炊の定番スキャナらしいです

富士通 ScanSnap iX500 (A4/両面)

富士通 ScanSnap iX500 (A4/両面)

これに決めた理由:

  • 評判がよかったから

他の選択肢を知らなかったので迷わずこれに決めました。

実際に使ってみても自炊用に作られていてクオリティが高いです。

例えば,裁断が甘くて糊がちょっと残ってて2枚の紙がくっついてる場合,スキャナ側が用紙を送ったときの超音波から異常を検知して中止・エラー通知してくれます。

また,スキャンした画像を検索可能なPDF(スキャンした画像の上にOCR検出したテキストデータが乗っていて文書内を検索できるPDF)にしてくれるソフトのライセンスなども付属していて,非常に役に立ちます。

買わなかったもの

もっとあたらしいバージョンが出ているんですが,さらに+1万円かかることになるので諦めました

今はその時とは価格が変わっているのでこちらもアリだとおもいます。

PFU ScanSnap iX1500 FI-IX1500

PFU ScanSnap iX1500 FI-IX1500

自炊に必要な費用

200DXとiX500で8万円近くかかりました。結構痛い出費です…。

裁断機を200DXにしなければもう少し安く抑えることができると思います。

成果

手元の本のうち「裁断してもいいや」とおもったものを9割方電子書籍化し終わった時点でこの程度の量でした

f:id:nigimitama:20190312140749j:plain

スキャナと裁断機の箱の体積分くらいにはなりましたかね

「結局,自炊用の設備が手元に残るからあんまり部屋の空きスペースが増えないんじゃないか」という不安もありますが,今後また

  • 電子書籍が出てない本
  • 1度目は紙で読みたい本
  • 新品の電子書籍だと高いけど紙の中古の本だと安い本

を買う時があれば少しずつ今回の投資を回収できるかなと思います。

ストレージ上に容量無限の本棚をもったようなものなので,置き場所の心配をすることなくガンガン本を買えます。

自炊用の設備にたっぷり投資した分,色んな本を沢山読んで回収していきたいと思います。

[R]R markdownでコードからmarkdownを書く

よくR markdownを使っているのですが,

  • 変数に格納した文字列を見出しに使う
  • forループを使って複数の節や小見出しを作る

といったコードからmarkdownを書くという操作がしたくなるときが時折あります。

ちょっと調べてみたのでメモ。

文字

rチャンクのオプションをresults='asis'にして,表示したいmarkdownテキストをcat()で表示すればいいみたいです。

見出し

例えば,results='asis'をオプションに指定したチャンクで次のように書けば見出しを書くことができます。

# ```{r, results='asis'}
library(tidyverse)

# 見出し ---------
midashi <- c("夏目漱石","草枕")

# h2見出し
paste0("## ", midashi[1],"\n") %>% cat()
# h3見出し
paste0("### ", midashi[2],"\n") %>% cat()
# ```

結果はこのようになります

f:id:nigimitama:20190309000443p:plain

文章

同様にテキストも("\n"で改行することを意識しつつ)cat()で表示できます

# ```{r, results='asis'}
# テキスト -------
# <br>タグで改行できる
text <- c("智に働けば角が立つ。情に棹させば流される。意地を通せば窮屈だ。<br>とかくに人の世は住みにくい。",
          "住みにくさが高じると、安い所へ引き越したくなる。どこへ越しても住みにくいと悟った時、詩が生れて、画が出来る。")

# 1つ目のパラグラフ
i = 1
paste0("\n", text[i], "\n") %>% cat()

# 2つ目のパラグラフ
i = i + 1
paste0("\n", text[i], "\n") %>% cat()
# ```

先程の見出しに続けて文章を出力するとこうなります。

f:id:nigimitama:20190309000752p:plain

グラフ

グラフでxに指定している変数名を節の見出しに使いたいときとかによくこの手法をつかいます

例えばirisデータセットのnumeric型の4変数をforループで順次参照して

  1. 参照している列名で見出しを表示
  2. 参照している列名でグラフを描画

という処理を繰り返す場合は,次のように書くことができます。

# ```{r, results='asis'}
# 変数名ベクトルを用意
VAR_NAMES = iris[,1:4] %>% colnames()

# forループしながらplot
for(VAR in VAR_NAMES){
  # h2見出し
  paste0("\n## ", VAR,"\n") %>% cat()
  
  # plot
  g <- ggplot(iris, aes_string(x = VAR, fill = "Species"))+
      geom_histogram(position = "identity", alpha = 0.5, bins = 20)
  print(g)
  
  # 改行
  cat("\n")
}
# ```

実行結果は次の「2.1 Sepal.Length」節以下の部分になります。

Sepal.LengthからPetal.Widthまで4つの変数について,見出しの表示とヒストグラムの作成を行うことができています。

f:id:nigimitama:20190309002116p:plain

おわりに

グラフを書く以外にも,例えば手元に複数の結果変数(被説明変数)があるときに結果変数ごとに節を分けて回帰分析を行ってみたりといったように,

markdownの節の作成」と「分析の処理」をforループでまとめて記述することができるのでとても便利です。おすすめです。

[R]ゼロから作る最尤法・ロジスティック回帰

「自分で数式をコードに落としていって動かす」という作業は非常に勉強になると思ったので,いろんなアルゴリズムをゼロから作っていくことで勉強してみたいと思います。

今回は最尤推定とロジスティック回帰についてです。最尤法のほうはおまけで,ロジスティック回帰がメインです。

恥ずかしながら,ロジスティック回帰モデルについては「パラメータ推定に最尤法を使うらしい」くらいしか理解しておらず,今回はじめて真面目に勉強しました…。

記事の前半で理論の要点を整理して,後半でコードを動かしていく構成になります。

最尤推定

まず,ロジスティック回帰でも使う推定方法について。

尤度

例えば,次のような問題について考えてみます。

あるコインを5回投げて2回が表になった。
このコインを投げたときに表が出る確率は?

この「表が出る確率」が推定したいパラメータになります。

表(head)が出ることを\(h=1\),裏が出ることを\(h = 0\)とし,\(h = 1\)となる確率を\(\theta\)\(h=0\)となる確率を\((1-\theta)\)とします。 \[ h = \begin{cases} 0 & \text{コインが裏のとき}\\ 1 & \text{コインが表のとき} \end{cases} \\ P(h = 1) = \theta, \hspace{2em} P(h=0) = (1-\theta) \] コインを5回投げて2回が表になる確率は,2項分布に従うと考えれば\({}_n\mathrm{C}_kp^k(1-p)^{n-k}\)なので

\[ P(h = {0, 0, 0, 1, 1}) = {}_{5}\mathrm{C}_{2} \theta^2 (1-\theta)^3 \]

となります。

この確率を尤度(ゆうど,likelihood)と呼び,多くの場合\(L(\theta)\)のように表記されます。 \[ L(\theta) = {}_{5}\mathrm{C}_{2} \theta^2 (1-\theta)^3 \] 尤度が最大になるパラメータ\(\theta\)が「最も尤もらしい(もっともらしい)」パラメータであるというわけです。

この尤度関数を描いてみると,0.4のあたりで最大になっていることが目視でもわかります。

library(tidyverse)

# 尤度関数
likelihood = function(θ, n = 5, k = 2){
  choose(n, k) * θ^k * (1 - θ)^(n-k) # nCk p^k (1-p)^{n-k}
}

# 描画用データ
θ = 1:100 * 0.01
L = likelihood(θ)
df = data_frame(θ, L)

# plot
ggplot(df, aes(x = θ, y = L))+
  geom_line(color = "dodgerblue")+
  labs(y = "L(θ)", title = expression(L(θ)==.[5]*C[2] * θ^2 * (1 - θ)^3))

f:id:nigimitama:20190217021011p:plain

対数尤度

尤度\(L(\theta)\)を最大にする\(\theta\)と対数尤度\(\log L(\theta)\)を最大にする\(\theta\)は同じなので,対数尤度を最大にする方法でも推定できます。

尤度は確率を沢山掛け合わせたものになって非常に小さな値を扱うことになりますが,対数尤度では対数にすることによって掛け算が足し算になります1。そのため,最尤推定では計算を簡単にするために対数尤度(log likelihood)が使われます。 \[ \begin{align} \log L(\theta) &= \log\{{}_{5}\mathrm{C}_{2} \theta^2 (1-\theta)^3\} \\ &= \log{}_{5}\mathrm{C}_{2} + 2\log \theta + 3\log (1-\theta)\\ \end{align} \]

df = data_frame(θ, logL=log(L))

# plot
ggplot(df, aes(x = θ, y = logL))+
  geom_line(color = "dodgerblue")+
  labs(y = "log L(θ)", title = expression(log(L(θ))==log(.[5]*C[2]) + 2*log(θ) + 3 * log(1 - θ)))

f:id:nigimitama:20190217021136p:plain

これをパラメータ\(\theta\)微分すると, \[ \begin{align} \frac{d \log L(\theta)}{d \theta} &= 2 \cdot \frac{1}{\theta} + 3 \cdot \frac{1}{1-\theta} \cdot (-1)\\ &= \frac{2}{\theta} - \frac{3}{1 - \theta} \end{align} \]

# 微分
logL = expression(log(c) + 2 * log(t) + 3 * log(1-t))
D(logL, "t")
> D(logL, "t")
2 * (1/t) - 3 * (1/(1 - t))

\(\log L(\theta)\)を最大にするパラメータ\(\theta^*\)は, \[ \frac{d \log L(\theta)}{d \theta} = \frac{2}{\theta} - \frac{3}{1 - \theta} = 0 \] となる\(\theta\)なので, \[ \theta^* = \frac{2}{5} = 0.4 \] となります。

このように尤度を最大にするパラメータを点推定する方法が最尤推定になります。

ロジスティック回帰

シグモイド関数

シグモイド関数(sigmoid function)は無限区間\((-\infty, \infty)\)の数を区間\((0,1)\)に変換する関数です2ロジスティック関数(logistic function)とも呼ばれます 3

ロジスティック回帰ではシグモイド関数を使うことで確率や二値変数を目的変数に使った回帰モデルを作ります。

# シグモイド関数
sigmoid = function(x) 1 / (1 + exp(-x))

# データ生成
x = -100:100 * 0.1
y = sigmoid(x)

# plot
ggplot(data_frame(x, y),
       aes(x = x, y = y))+
  geom_line(color = "dodgerblue")+
  labs(y = expression(σ(x)), x = expression(x))

f:id:nigimitama:20190217021952p:plain

ロジット関数

目的変数\(y\in \{0,1\}\)\(p\)次元説明変数\(\boldsymbol{x}=(x_1,...,x_p)^T\)で説明することを考えます。

\(\boldsymbol{x}\)の下で\(y=1\)になる条件付き確率と\(y=0\)になる条件付き確率を \[ P(y=1|\boldsymbol{x}) = \pi(\boldsymbol{x}), \hspace{2em} P(y=0|\boldsymbol{x}) = 1-\pi(\boldsymbol{x}) \]

とするとき,その確率の比をオッズ(odds)といいます。 \[ \text{odds}(\pi(\boldsymbol{x})) = \frac{\pi(\boldsymbol{x})}{1-\pi(\boldsymbol{x})} \] オッズの対数は対数オッズ(log odds)やロジット関数(logit function)と呼ばれます。 \[ \text{logit}(\pi(\boldsymbol{x})) = \log{\frac{\pi(\boldsymbol{x})}{1-\pi(\boldsymbol{x})}} \]

ロジスティック回帰モデル

一般化線形モデルの文脈では,リンク関数にロジット関数を用いて目的変数\(y\)と線形予測子\(\beta_0 + \beta_1 x_{1} + \cdots + \beta_p x_{p}=\boldsymbol{\beta}^{\top}\boldsymbol{x}\)を繋げた \[ \log{\frac{\pi(\boldsymbol{x})}{1-\pi(\boldsymbol{x})}} = \beta_0 + \beta_1 x_{1} + \cdots + \beta_p x_{p}=\boldsymbol{\beta}^{\top}\boldsymbol{x} \] をロジスティック回帰モデルと呼ぶようです(本橋 2015)。

これの両辺の指数をとると \[ \frac{\pi(\boldsymbol{x})}{1-\pi(\boldsymbol{x})} = \exp(\boldsymbol{\beta}^{\top}\boldsymbol{x}) \] であり,整理すると4 \[ \pi(\boldsymbol{x}) =\frac{\exp(\boldsymbol{\beta}^{\top}\boldsymbol{x})}{1 + \exp(\boldsymbol{\beta}^{\top}\boldsymbol{x})} \] と,シグモイド関数になります。

\(P(y=1|\boldsymbol{x}) = \pi(\boldsymbol{x})\)でしたので,\(\boldsymbol{x}\)の下で\(y=1\)になる確率\(P(y=1|\boldsymbol{x})\)を目的変数にして,ロジスティック関数を介して線形予測子\(\boldsymbol{\beta}^{\top}\boldsymbol{x}\)で説明している・・という感じになるようです。 \[ P(y=1|\boldsymbol{x}) = \pi(\boldsymbol{x}) =\frac{\exp(\boldsymbol{\beta}^{\top}\boldsymbol{x})}{1 + \exp(\boldsymbol{\beta}^{\top}\boldsymbol{x})} \]

パラメータの推定

ここまでで「確率や二値変数をどう\(\boldsymbol{\beta}^{\top}\boldsymbol{x}\)で説明するか」が定まったので,次はパラメータの推定について整理します。

誤差関数

最小二乗法では誤差関数に誤差ニ乗和(Sum of Squared Error:SSE)や平均二乗誤差(Mean Squared Error:MSE)を使いましたが,ロジスティック回帰では交差エントロピー誤差関数(cross-entropy error function)を使います。

モデルの目的変数\(y\in\{0,1\}\)が1をとる確率と0をとる確率をそれぞれ \[ P(y=1) = \pi, \hspace{1em} P(y=0) = 1- \pi \] と書くことにします。また確率変数\(y\)はパラメータ\(\pi\)のベルヌーイ分布に従う(\(y \sim Bern(\pi)\))とし, \[ f(y|\pi)=\pi^{y} (1 - \pi)^{1-y} \] とします。

\(N\)回の試行に基づく尤度関数 \[ L(\pi_1,...,\pi_N) = \prod_{i=1}^N f(y_i|\pi_i) = \prod_{i=1}^N \pi_i^{y_i} (1 - \pi_i)^{1-y_i} \] の対数をとって負の符号を付けたもの \[ \begin{align} \mathcal{L}(\pi_1,...,\pi_N) &= - \log L(\pi_1,...,\pi_N)\\ &= -\sum_{i=1}^N \{y_i \log \pi_i + (1-y_i) \log(1-\pi_i)\} \end{align} \]交差エントロピー誤差関数になります。

先程の \[ \pi_i = \frac{\exp(\boldsymbol{\beta}^{\top}\boldsymbol{x}_i)}{1 + \exp(\boldsymbol{\beta}^{\top}\boldsymbol{x}_i)} \] を代入して整理すれば, \[ \mathcal{L}(\pi_1,...,\pi_N) = \mathcal{L}(\boldsymbol{\beta}) = -\sum_{i=1}^N \{y_i \boldsymbol{\beta}^{\top}\boldsymbol{x}_i - \log (1 + \exp(\boldsymbol{\beta}^{\top}\boldsymbol{x}_i))\} \] になります(ただし,\(\boldsymbol{x}_i\)は入力ベクトル…レコード\(i\)の説明変数達です)

また,これの勾配(微分したもの)は \[ \frac{\partial \mathcal{L}(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = -\sum_{i=1}^N \left\{y_i \boldsymbol{x}_i - \frac{\boldsymbol{x}_i\exp(\boldsymbol{\beta}^{\top}\boldsymbol{x}_i)}{1 + \exp(\boldsymbol{\beta}^{\top}\boldsymbol{x}_i)}\right\} = \sum_{i=1}^N \boldsymbol{x}_i (\pi_i - y_i) \] になります。

尤度関数・交差エントロピー誤差は最小二乗法のMSEなどとは異なり解析的に解が出せないので,勾配降下法ニュートン・ラフソン法を使って,コンピュータでゴリゴリと数値計算をして解を出します。

交差エントロピー誤差を\(N\)で割った平均交差エントロピー誤差関数 \[ -\frac{1}{N}\sum_{i=1}^N \{y_i \boldsymbol{\beta}^{\top}\boldsymbol{x}_i - \log (1 + \exp(\boldsymbol{\beta}^{\top}\boldsymbol{x}_i))\} \] というのもあるようです。こちらのほうが誤差がデータ数の影響を受けにくいのだとか(伊藤 2018)。

勾配降下法

勾配降下法(gradient descent method,あるいは最急降下法steepest descent method)5

  1. 係数(重み)の初期値を決める
  2. 誤差を評価する
  3. 係数を誤差関数の勾配に基づいて更新する
  4. 2~3を繰り返す

という感じの方法です。

係数の更新には次の式が使われます。 ここで学習率(learning rate)と呼ばれるもので,係数を更新するときに値を大きく変えすぎないように0.01とか小さな値を設定します。

Rで実装

irisデータセットのうちSpecies列が2値のデータを使います。

手元の本に具体的なアルゴリズムの例がないので一部想像で補完してコードを書きます(係数の初期値とか)。

# データ
df = iris[51:150,] # 51 ~ 150行
X = cbind(Intercept = 1, df[,-5]) %>% as.matrix() # 切片項用の列を追加
y = ifelse(df[,5] == "versicolor", 0, 1) # versicolor = 0, virginica = 1

# シグモイド関数
sigmoid = function(x) exp(x) / (1 + exp(x))

# 勾配降下法の関数を定義
GradientDescent = function(X, y, K = 10, eta = 0.001) { # K:重み更新の反復回数, eta: 学習率
  n = dim(X)[1] # サンプルサイズ
  p = dim(X)[2] # 説明変数の数
  
  beta = runif(p)
  
  k = 1  # カウントをリセット
  while (k != K) {
    # 線形予測子Xβをシグモイド関数に通す
    pi = sigmoid(X %*% beta)
    
    # 交差エントロピー誤差関数の勾配
    gradient = 0 # 初期化
    for(i in 1:n){ # 総和
      gradient = gradient + X[i,] * (pi[i] - y[i])
    }
    
    # 係数の更新
    ## β := β - eta * 誤差関数の勾配
    beta = beta - eta * gradient
    
    # (係数の更新の様子を記録)
    ## 誤差
    MSE = mean((pi - y)^2)
    ## データフレームにまとめる
    df_k = data.frame(k, t(beta), MSE)
    ## 結合
    if(k == 1) df = df_k
    if(k != 1) df = rbind(df, df_k)
    
    # カウントを更新
    k = 1 + k
    }
  
  # 結果
  result = list(coefficients = beta, # データにフィットさせた係数
                MSE = MSE, # 最終的なMSE
                log = df)  # 係数の更新の様子
  return(result)
}

推定したパラメータは次のようになりました。

# パラメータの推定
result = GradientDescent(X, y, K = 20000)

# 推定した係数
result$coefficients
##    Intercept Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
##    -7.597939    -4.262963    -4.924182     6.707644     9.197969

20000回の重みの更新とMSEの推移が次のグラフです。

20000回やってもまだ推定値が収束せず,反復回数を増やすたびにじわじわ動いていて,MSEもじわじわ動いているので,反復回数が足りないようです。

(MSEを予測精度の評価に使うのは雑な気がしますが,まぁ参考程度ということで…。)

# plot
temp = result$log %>% gather(variable, value, 2:6)
ggplot() + 
  geom_line(aes(x = temp$k, y = temp$value, color = temp$variable))+
  geom_line(aes(x = temp$k, y = temp$MSE * 150), color = "red", linetype = 2) +
  scale_y_continuous(sec.axis = sec_axis(~ . /150, name = "MSE")) +
  annotate("text", x = temp$k[1], y = temp$MSE[1] * 150 + 2, label = "MSE", color = "red")+
  labs(y = "Estimate", x = "iteration", color = "variable")

f:id:nigimitama:20190217121712p:plain

最終的なMSEはこうなりました。

str_c("MSE: ", result$MSE) # MSE
## [1] "MSE: 0.0234984687192942"

ちなみに,ニュートン法を使用しているglm()とはだいぶ結果が違います…。

# glmと比較
res_glm = glm(Species ~ . , data = iris[51:150,], family = binomial(link = "logit"))
res_glm
## 
## Call:  glm(formula = Species ~ ., family = binomial(link = "logit"), 
##     data = iris[51:150, ])
## 
## Coefficients:
##  (Intercept)  Sepal.Length   Sepal.Width  Petal.Length   Petal.Width  
##      -42.638        -2.465        -6.681         9.429        18.286  
## 
## Degrees of Freedom: 99 Total (i.e. Null);  95 Residual
## Null Deviance:       138.6 
## Residual Deviance: 11.9  AIC: 21.9
y_hat = res_glm$fitted.values
str_c("glmのMSE: ",(y - y_hat)^2 %>% mean())
## [1] "glmのMSE: 0.0188203802832239"

glmのほうがMSEを小さくできていますね

ニュートン法

ニュートン法(Newton's method)あるいはニュートン・ラフソン法(Newton-Raphson method)は,誤差関数をテイラー展開で近似したものを使って解く方法です。

次のように係数(重み)\(\boldsymbol{\beta}\)を更新していきます。 \[ \boldsymbol{\beta}^{(new)} = \boldsymbol{\beta}^{(old)} - \boldsymbol{H}^{-1} \frac{\partial \mathcal{L}(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} \] ここで\(\boldsymbol{H}\)ヘッセ行列(Hessian matrix)という,\(\boldsymbol{\beta}\)の2次の導関数を要素に持つ行列で,その\((i,j)\)要素は \[ (\boldsymbol{H})_{ij} = \frac{\partial^2 \mathcal{L}(\boldsymbol{\beta})}{\partial \beta_i \partial \beta_j} \] つまり \[ \boldsymbol{H}= \left( \begin{array}{cccc} \frac{\partial^2 \mathcal{L}(\boldsymbol{\beta})}{\partial \beta_0^2} & \frac{\partial^2 \mathcal{L}(\boldsymbol{\beta})}{\partial \beta_0 \partial \beta_1} & \cdots & \frac{\partial^2 \mathcal{L}(\boldsymbol{\beta})}{\partial \beta_0 \partial \beta_p} \\ \frac{\partial^2 \mathcal{L}(\boldsymbol{\beta})}{\partial \beta_1 \partial \beta_0} & \frac{\partial^2 \mathcal{L}(\boldsymbol{\beta})}{\partial \beta_1^2} & \cdots & \frac{\partial^2 \mathcal{L}(\boldsymbol{\beta})}{\partial \beta_1 \partial \beta_p} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2 \mathcal{L}(\boldsymbol{\beta})}{\partial \beta_p \partial \beta_0} & \frac{\partial^2 \mathcal{L}(\boldsymbol{\beta})}{\partial \beta_p \partial \beta_1} & \cdots & \frac{\partial^2 \mathcal{L}(\boldsymbol{\beta})}{\partial \beta_p^2} \end{array} \right) \] のような行列になります。

反復再重み付け最小二乗法(IRLS)

ニュートン・ラフソン法を交差エントロピー誤差関数に適用したものは反復再重み付け最小二乗法(iteratively reweighted least squares:IRLS)と呼ばれるようで,Rのglm.fit()でもそう呼ばれています。

The default method "glm.fit" uses iteratively reweighted least squares (IWLS)
glm function | R Documentationより)

ニュートン・ラフソン法を交差エントロピー誤差関数に適用する場合,ニュートン・ラフソン法のヘッセ行列は \[ \boldsymbol{H} = \frac{\partial^2 \mathcal{L}(\boldsymbol{\beta})}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^{\top}} = \sum_{i=1}^N \pi_i (1 - \pi_i) \boldsymbol{x}_i \boldsymbol{x}_i^{\top} = \boldsymbol{X}^{\top} \boldsymbol{W}\boldsymbol{X} \] になります。ここで\(\boldsymbol{X}\)はデータの行列で,\(\boldsymbol{W}\)\(N\times N\)対角行列で,その要素は \[ W_{ii} = \pi_i(1 - \pi_i) \] になります。

そして,重みの更新は \[ \boldsymbol{\beta}^{(new)} = \boldsymbol{\beta}^{(old)} - (\boldsymbol{X}^{\top} \boldsymbol{W}\boldsymbol{X})^{-1} \frac{\partial \mathcal{L}(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} \] となり, \[ \frac{\partial \mathcal{L}(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = \sum_{i=1}^N \boldsymbol{x}_i (\pi_i - y_i) = \boldsymbol{X}^{\top}(\boldsymbol{\pi} - \boldsymbol{y}) \] なので \[ \boldsymbol{\beta}^{(new)} = \boldsymbol{\beta}^{(old)} - (\boldsymbol{X}^{\top} \boldsymbol{W}\boldsymbol{X})^{-1} \boldsymbol{X}^{\top}(\boldsymbol{\pi} - \boldsymbol{y}) \] となります。

また, \[ \boldsymbol{z} = \boldsymbol{X} \boldsymbol{\beta} - \boldsymbol{W}^{-1}(\boldsymbol{\pi} - \boldsymbol{y}) \] とおけば \[ \begin{align} \boldsymbol{\beta}^{(new)} &= \boldsymbol{\beta}^{(old)} - (\boldsymbol{X}^{\top} \boldsymbol{W}\boldsymbol{X})^{-1} \boldsymbol{X}^{\top}(\boldsymbol{\pi} - \boldsymbol{y}) \\ &=(\boldsymbol{X}^{\top} \boldsymbol{W}\boldsymbol{X})^{-1} \{\boldsymbol{\beta} - \boldsymbol{X}^{\top}(\boldsymbol{\pi} - \boldsymbol{y})\}\\ &=(\boldsymbol{X}^{\top} \boldsymbol{W}\boldsymbol{X})^{-1} \boldsymbol{X}^{\top} \boldsymbol{W} \boldsymbol{z} \end{align} \] と書くこともできます(最終的に加重最小二乗法の正規方程式みたいな形になります)。

Rで実装

手元の本に具体的なアルゴリズムの例がないので一部想像で補完してコードを書きます。

試してみたところ,係数の初期値は乱数でも固定値でもよさそうですが,0.1とか微小な値にするのが良さそうです。0.7とかだと収束しない場合が…

# IRLS関数を定義
IRLS = function(X, y, K = 20) { # K:重みを更新する反復回数
  n = dim(X)[1] # サンプルサイズ
  p = dim(X)[2] # 説明変数の数
  
  # 係数βの初期値をランダムに決める
  beta = runif(n = p, max = 0.1)
  
  k = 1  # カウントをリセット
  while (k != K) {
    # 線形予測子Xβをシグモイド関数に通す
    pi = sigmoid(X %*% beta)
    
    # 重み行列 W の作成
    W = matrix(0, ncol = n, nrow = n) # 値を入れるための空行列
    for(i in 1:n) {
      W[i,i] = pi[i] * (1 - pi[i])
    }
 
    # 係数の更新
    ## β := β - (X'WX)^{-1}X'(pi-y)
    beta = beta - solve(t(X) %*% W %*% X) %*% t(X) %*% (pi - y)
    
    # (係数の更新の様子を記録)
    ## 誤差
    MSE = mean((pi - y)^2)
    ## データフレームにまとめる
    df_k = data.frame(k, t(beta), MSE)
    ## 結合
    if(k == 1) df = df_k
    if(k != 1) df = rbind(df, df_k)
    
    # カウントを更新
    k = 1 + k
    }
  
  # 結果
  result = list(coefficients = beta, # データにフィットさせた係数
                MSE = MSE, # 最終的なMSE
                log = df) # 係数の更新の様子
  return(result)
}

20回の重みの更新を行うようにして推定したところ,推定したパラメータは次のようになりました。

# パラメータの推定
result = IRLS(X, y, K = 20)

# 推定した係数
result$coefficients
##                    [,1]
## Intercept    -42.637804
## Sepal.Length  -2.465220
## Sepal.Width   -6.680887
## Petal.Length   9.429385
## Petal.Width   18.286137
# glmと比較
glm(Species ~ . , data = iris[51:150,], family = binomial(link = "logit"))
## 
## Call:  glm(formula = Species ~ ., family = binomial(link = "logit"), 
##     data = iris[51:150, ])
## 
## Coefficients:
##  (Intercept)  Sepal.Length   Sepal.Width  Petal.Length   Petal.Width  
##      -42.638        -2.465        -6.681         9.429        18.286  
## 
## Degrees of Freedom: 99 Total (i.e. Null);  95 Residual
## Null Deviance:       138.6 
## Residual Deviance: 11.9  AIC: 21.9

IRLSで解いているglm()の係数と一致しました。コードに誤りはなさそうです。

20回の重みの更新の推移が次のグラフです。10回以前にすでに推定の収束が済んでいますね。なかなか優秀です。

# plot
temp = result$log %>% gather(variable, value, 2:6)
ggplot() + 
  geom_line(aes(x = temp$k, y = temp$value, color = temp$variable))+
  geom_line(aes(x = temp$k, y = temp$MSE * 150), color = "red", linetype = 2) +
  scale_y_continuous(sec.axis = sec_axis(~ . /150, name = "MSE")) +
  annotate("text", x = temp$k[1], y = temp$MSE[1] * 150 + 2, label = "MSE", color = "red")+
  labs(y = "Estimate", x = "iteration", color = "variable")

f:id:nigimitama:20190217122337p:plain

最終的なMSEは

str_c("MSE: ", result$MSE) # MSE
## [1] "MSE: 0.0188203802832454"

になりました。勾配降下法よりもいい感じです。

まとめ

  • ロジスティック回帰は,対数尤度を負にした交差エントロピーを誤差関数に使う
  • 交差エントロピーは(最小二乗法の平均二乗誤差とは違い)解析的に解を得られないので,勾配降下法ニュートン・ラフソン法でゴリゴリ計算していって数値的に解を出す
  • Rのglm()関数で採用されているのはニュートン法のほう。少ない反復回数でうまく収束してくれる。

参考


  1. 対数法則:(1) \(\log (xy) = \log x + \log y\), (2) \(\log (x^p) = p \log x\)

  2. \(\exp(-x)\)\(e^{-x}\)と表記されることもあります。

  3. この関数はロジスティック関数と呼ばれていたり,シグモイド関数と呼ばれていたり,ロジスティック・シグモイド関数(logistic sigmoid function)と呼ばれたりします。

  4. \(\begin{align}\pi(\boldsymbol{x}) &= \exp(\boldsymbol{\beta}^{\top}\boldsymbol{x})(1-\pi(\boldsymbol{x})) \\&= \exp(\boldsymbol{\beta}^{\top}\boldsymbol{x})- \exp(\boldsymbol{\beta}^{\top}\boldsymbol{x}) \pi(\boldsymbol{x})\\ \pi(\boldsymbol{x}) + \exp(\boldsymbol{\beta}^{\top}\boldsymbol{x}) \pi(\boldsymbol{x})&= \exp(\boldsymbol{\beta}^{\top}\boldsymbol{x}) \\ \pi(\boldsymbol{x}) (1 + \exp(\boldsymbol{\beta}^{\top}\boldsymbol{x})) &=\exp(\boldsymbol{\beta}^{\top}\boldsymbol{x}) \\ \pi(\boldsymbol{x}) &=\frac{\exp(\boldsymbol{\beta}^{\top}\boldsymbol{x})}{1 + \exp(\boldsymbol{\beta}^{\top}\boldsymbol{x})} \end{align}\)

  5. バッチ学習と呼ばれることもあります。ただし,名前が少し似ている確率的勾配降下法(stochastic gradient descent method: SGDオンライン学習とも)はこれとは少し異なるアルゴリズムです。

[読書メモ]『ビジネスフレームワーク図鑑』

ビジネスフレームワーク図鑑 すぐ使える問題解決・アイデア発想ツール70

ビジネスフレームワーク図鑑 すぐ使える問題解決・アイデア発想ツール70

半額セールを利用して買ってみました。

フレームワーク”と呼ばれるものは就職活動を始めた頃から多く目にするようになり,結構いろんなものがあるので,「世の中にはいったいいくつのフレームワークがあるんだろう…」と思っていて,図鑑的なものが欲しいとちょうど思っていました。

この本は私が求めていた要素が十分に揃っていて,買ってよかったなと思うものでした。

具体的には

  1. 網羅されていること
    • 有名なものも,細かいものも載っています
  2. 整理されていること
    • 「問題・課題を発見したいときはこれらのフレームワーク」など目的のジャンル別にまとまっています
  3. 読みやすいこと
    • 図解と程よい分量の解説文
  4. 使い方がわかること
    • 「使い方」という節がすべてのフレームワークについていて,そこを読めばだいたいわかりそうです

です。

「広く浅くざっと知る」のを目的にする分には十分なもので,入門には最適そうです。

[読書メモ]話し方・会話に関する本

話し方に関する本3冊に目を通したのでざっくりメモ。

『「話し方」の授業』

ビジネスでの話し方や雑談での話し方に関する本です。

具体的な問題に対する細かなアドバイスがまとまっています。

例えば「スピーチ中に頭が真っ白になった時」→「①自分がどこまで話したか聴衆に訊いてしまう,②『ここまでで何か質問はありますか』と一旦区切る,etc.」など。

解決策が複数提示してある場合が多いのも情報が多くていいなと思います。

細かいことまで書いてある分,文字数はやや多めなので,全体を通しで一気に読む本というよりは自分が困ってるシチュエーションの部分だけピックアップして読む方法が合いそうな感じの本です。

『コミュ障で損しない方法38』

コミュ障で損しない方法38

コミュ障で損しない方法38

雑談での話し方の本。

「会話は『気まずさ』を駆逐する(協力型の)ゲームである」という考え方をベースに書かれています。

コミュ障向けにかかれているためか会話に対する姿勢・考え方についての記述が多めな印象ですが,「質問の仕方」など悩みへの解決策の話も書かれています。

『マンガでわかる!会話がとぎれない!話し方』

雑談での話し方の本。

Kindle Unlimitedだったので読んでみたところ,かなり良かったです。

「会話は気持ちのキャッチボール」という考え方をベースにして書かれています。

話し方,聞き方,尋ね方,答え方など,基礎的な部分・高頻度で直面する問題に対するアドバイスが主な内容なのが私には合ってました。

マンガなので読むハードルが低く,また絵があると会話のシチュエーションがイメージしやすいのもあり,非常に読みやすかったです。

[R]ゼロから作る最小二乗法2:重回帰

「自分で数式をコードに落としていって動かす」という作業は非常に勉強になると思ったので,いろんなアルゴリズムをゼロから作っていきたいと思います。

理論の要点を整理してから実装する構成で述べていきます。

モデル

線形回帰(linear regression)は,予測の目的変数\(Y_i\)と説明変数\(X_{i1}, X_{i2}, ..., X_{ip}\)の間に次のような線形関係を仮定したモデルを置いて予測する手法です。 \[ Y_i = \beta_0 + \beta_1X_{i1} + \beta_2 X_{i2} + \cdots + \beta_p X_{ip} + u_i \]

このように説明変数が複数ある線形回帰は重回帰モデル(multiple regression model)とも呼ばれます。

行列を使って表記すると,次のように表現できます

\[ \boldsymbol{y}=\boldsymbol{X} \boldsymbol{\beta}+\boldsymbol{u} \]

ただし,

\[ \boldsymbol{y} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} , \ \boldsymbol{X} = \begin{bmatrix} 1 & X_{11} & \cdots & X_{1p} \\ 1 & X_{21} & \cdots & X_{2p} \\ \vdots &\vdots & &\vdots &\\ 1 & X_{n1} & \cdots & X_{np} \\ \end{bmatrix} ,\ \boldsymbol{\beta} = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_p \end{bmatrix} , \ \boldsymbol{u} = \begin{bmatrix} u_1 \\ u_2 \\ \vdots \\ u_n \end{bmatrix} \]

誤差関数

\(\boldsymbol{\beta}\)の最小2乗推定量(ordinary least square’s estimator: OLSE)を\(\hat{\boldsymbol{\beta}}\)とすると,モデルからの\(\boldsymbol{y}\)の予測値は\(\hat{\boldsymbol{y}} = \boldsymbol{X}\hat{\boldsymbol{\beta}}\),実測値と予測値の誤差(残差)は\(\boldsymbol{e}=\boldsymbol{y}-\hat{\boldsymbol{y}}\)になります。

パラメータ\(\boldsymbol{\beta}\)を推定するために多く使われる方法は,実測値\(\boldsymbol{y}\)と予測値\(\hat{\boldsymbol{y}}\)誤差2乗和(Sum of Squared Error)1

\[ SSE = \sum_{i=1}^n (y_i - \hat{y}_i)^2 = \sum_{i=1}^n e_i^2=\boldsymbol{e}'\boldsymbol{e} \]

を最小にするパラメータを採用するという最小二乗法(least squares method)です。

パラメータの推定

最小二乗推定量は「誤差2乗和を微分してゼロになる点」という次の条件

\[ \frac{\partial \boldsymbol{e}'\boldsymbol{e}}{\partial \hat{\boldsymbol{\beta}}}=\boldsymbol{0} \]

の解となるものです。

誤差2乗和は

\[ \begin{align} \boldsymbol{e}'\boldsymbol{e} &= (\boldsymbol{y} - \hat{\boldsymbol{y}})'(\boldsymbol{y} -\hat{\boldsymbol{y}})\\ &= (\boldsymbol{y} - \boldsymbol{X}\hat{\boldsymbol{\beta}})'(\boldsymbol{y} - \boldsymbol{X}\hat{\boldsymbol{\beta}})\\ &= \boldsymbol{y}'\boldsymbol{y} - \boldsymbol{y}'\boldsymbol{X}\hat{\boldsymbol{\beta}} - \hat{\boldsymbol{\beta}}'\boldsymbol{X}'\boldsymbol{y} + \hat{\boldsymbol{\beta}}'\boldsymbol{X}'\boldsymbol{X}\hat{\boldsymbol{\beta}}\\ &= \boldsymbol{y}'\boldsymbol{y} - 2\hat{\boldsymbol{\beta}}'\boldsymbol{X}'\boldsymbol{y} + \hat{\boldsymbol{\beta}}'\boldsymbol{X}'\boldsymbol{X}\hat{\boldsymbol{\beta}} \end{align} \]

であるから23

\[ \frac{\partial \boldsymbol{e}'\boldsymbol{e}}{\partial \hat{\boldsymbol{\beta}}} =-2\boldsymbol{X}'\boldsymbol{y}+2(\boldsymbol{X}'\boldsymbol{X})\hat{\boldsymbol{\beta}} =\boldsymbol{0} \] となり,これを書き直すと,最小二乗法の正規方程式と言われる次の式を得ます \[ (\boldsymbol{X}'\boldsymbol{X})\hat{\boldsymbol{\beta}}=\boldsymbol{X}'\boldsymbol{y} \] この正規方程式を\(\hat{\boldsymbol{\beta}}\)について解けば \[ \hat{\boldsymbol{\beta}}=(\boldsymbol{X}'\boldsymbol{X})^{-1}\boldsymbol{X}'\boldsymbol{y} \] が得られます。

Rで実装

Rでの行列の掛け算は%*%という演算子で行います。

また,転置はt()逆行列solve()になります。

# 関数を定義
OLS <- function(X, y) {
  # 切片項の追加
  X = data.frame(Intercept = rep(1, nrow(X)), X)
  
  # 入力されたデータフレームをmatrixに変える
  if(is.data.frame(X)) X = as.matrix(X)
  if(is.data.frame(y)) y = as.matrix(y)
  
  # パラメータの推定: β = (X'X)^{-1} X'y
  beta = solve(t(X) %*% X) %*% t(X) %*% y
  return(beta)
}

# treesデータセットで試す
X = trees[c("Girth","Height")]
y = trees["Volume"]
OLS(X, y)
##                Volume
## Intercept -57.9876589
## Girth       4.7081605
## Height      0.3392512
# lm()と比較
lm(Volume ~ Girth + Height, trees)
## 
## Call:
## lm(formula = Volume ~ Girth + Height, data = trees)
## 
## Coefficients:
## (Intercept)        Girth       Height  
##    -57.9877       4.7082       0.3393

ちゃんと同じ値が出ました。

参考

経済分析のための統計的方法

経済分析のための統計的方法

計量経済学大全

計量経済学大全


  1. 残差平方和(Sum of Squared Residuals: SSR,あるいはResidual Sum of Squares: RSS)とも呼ばれます(統計学の分野だと残差,機械学習の分野だと誤差と呼ぶことが多いように思います)

  2. ※転置の基本公式から,\((\boldsymbol{X}\hat{\boldsymbol{\beta}})'=\hat{\boldsymbol{\beta}}'\boldsymbol{X}'\)

  3. \(\boldsymbol{y}'\boldsymbol{X}\hat{\boldsymbol{\beta}}=(\boldsymbol{X}\hat{\boldsymbol{\beta}})'\boldsymbol{y}=\hat{\boldsymbol{\beta}}'\boldsymbol{X}'\boldsymbol{y}\)

本の安売りセール情報を知りたい

最近SNS経由で

翔泳社の50%OFFセール

「まんがで読破」シリーズの99%OFFセール

を偶然知り,たいへん得しました。

今回はたまたまTwitter等で見かけて知ることができたものの,こういうお得な情報を見逃すことは減らしたいものです。

どうしたらお得な情報が得られるのでしょうか?

軽く調べたところ,本(特にKindle本)の安売り情報を知りたい場合は,以下のような方法がありそうでした。

(1) 自分から調べに行く

自分からWebサイトを訪れることを厭わないのであれば,次のような方法があるかと思います

Kindleのトップページで知る

Kindleのトップページには出版社側のセール情報の告知バナーが出ているのでここで確認できそうです

f:id:nigimitama:20190211133727p:plain

Kindle日替わり・月替りセールのページで知る

出版社側のセールではなくAmazon側のセールということだと思うのですが,日替わり・月替りのセールがあるようです。

Amazon.co.jp: Kindle日替わりセール: Kindleストア

Amazon.co.jp: Kindle月替わりセール: Kindleストア

③セール情報をまとめているサイトを使う

電子書籍の更地というブログを見つけました。

主に出版社側のセール情報を掲載されているようです(Amazonの月替りセール情報も載っているが,日替わりセールについては載って無いかも)

また,Kindle以外の電子書籍のセール情報も掲載されているようです。

(2) セール情報を通知してもらう

自分から調べに行くのはめんどうなので,できれば通知してもらいたいものです。

Amazon.co.jpデリバーズ(メール通知)

f:id:nigimitama:20190211153714p:plain

Kindle日替わりセール月替わりセールのページに

Amazon.co.jpデリバーズ Kindle本の登録はこちら

という文言があります。

これを登録することでEメールでセール情報を通知してくれるようです。

②セール情報をまとめているサイトを使う

先程紹介した電子書籍の更地さんは,はてなブログなので,はてなの会員の方は「読者になる」ボタンを押せばメール等で通知がくるようになりますし,Twitterもされているようなのでフォローすればそちらからでも情報を受けることができそうです。

この方法が一番ラクかもしれません。

まとめ

Kindle本のセール情報を知る方法について,「Amazon側のセール / 出版社側のセール」の軸,「自分から調べに行く / 通知してもらう」の軸との2軸で整理すると

自分で調べに行く 通知してもらう
Amazon側のセール Kindle日替わりセール月替わりセールのページで知る Amazonデリバーズ(購読Eメール)を使う
出版社側のセール Kindleのトップページで知る
電子書籍の更地で知る
電子書籍の更地さんを使う

という感じですかね。

私はとりあえずAmazonデリバーズの登録と,「電子書籍の更地」さんの購読をしてしばらく使い心地を確認してみようと思います