読者です 読者をやめる 読者になる 読者になる

雑記: Python statsmodels をつかってみるだけ

以下の記事を全面的に参考にさせていただきました。
data.gunosy.io
上の記事では電力使用量のトレンドを抽出していますが、ここでは(電力使用量とかなり似通っていますが)R の組み込みデータである UKgas(1960年から1986年までの四半期ごとのイギリスのガス消費量)をつかいます。
Python から読み込む都合上、UKgas は以下のように csv に出力しておきます。

qtr, gas
1960-03, 160.1
1960-06, 129.7
1960-09, 84.8
1960-12, 120.1
1961-03, 160.1
(以下略)
UKgas をそのままプロットすると以下です。
  f:id:cookie-box:20170303235217p:plain:w600

トレンド成分、季節成分、残差に分解したのが以下です。
f:id:cookie-box:20170303235232p:plain:w560

季節成分は以下です。毎年1~3月がガスの消費量が大きい。

[ 175.13810096  -36.14122596 -168.96766827   29.97079327]

スクリプト

# -*- coding: utf-8 -*-
import numpy
import pandas
import matplotlib.pyplot as plt
import statsmodels.api as sm

if __name__ == "__main__":
  x = pandas.read_csv('UKgas.csv', skiprows=1, names=('qtr', 'gas'))
  x['qtr'] = pandas.to_datetime(x['qtr'])

  # plt.plot(x['qtr'], x['gas'])
  # plt.show()

  res = sm.tsa.seasonal_decompose(x['gas'].values, freq=4)
  res.plot()
  plt.show()

  print res.seasonal[0:4]

確率論セミナー(2/2, 2/23): 不参加メモ

Skype数学勉強会 確率論セミナー に参加できなかったメモ
読んでいる本(現在はサブテキスト): はじめての確率論 測度から確率へ : 佐藤 坦 : 本 : Amazon
前回: 確率論セミナー(51): 不参加メモ - クッキーの日記

参加できなかった回に進んだページ: 60~65ページ
以下雰囲気で書いたメモ

  • ルベーグ積分のつよいところ
    • いろんな関数が積分できる
    • lim と積分が多くの場合で交換する(リーマン積分はあまり交換しない)
  • じゃあどんな場合に lim と積分が交換するのか → 以下の3ステップで考える
    • まずは限定的な条件で lim と積分が交換することを示す―単調収束定理
      • 非負可測関数列がほとんどいたるところ単調増大していてかつ lim が有限でかつその積分値も有限なときは lim と積分が交換する
    • 一般の非負可測関数列で lim と積分を交換したときの不等式の関係をおさえる―ファトゥの補題
      •  E[\lim \inf _{n \to \infty} X_n ] \leq \lim \inf _{n \to \infty} E[ X_n ]\{X_n\} は非負可測関数列)
    • 特定の条件下で lim と積分が交換することを示す―ルベーグ収束定理(優収束定理)
      • 概収束する可測関数列 \{X_n\} の絶対値の sup がほとんどいたるところにおいて積分可能な非負可測関数でおさえられるとき lim と積分が交換する

確率論ではとにかく期待値をとりたい(積分をしたい)
どんな場合に期待値がとれるか知りたいので lim と積分がどんな場合に交換するかは重要

雑記: X-means法

X-means アルゴリズムを実装してみます。

参考文献

TLに以下の記事の引用ツイートが流れてきたのでやってみます。
スクリプトを全面的に参考にさせていただきました。
aaaazzzz036.hatenablog.com

使用データ

以前に以下の記事で、各都道府県の最低賃金合計特殊出生率、高齢化率を用意したので流用します。
入門 機械学習による異常検知―Rによる実践ガイド: ノート3 - クッキーの日記

やること

最低賃金合計特殊出生率、高齢化率の3変数に着目して都道府県をクラスタリングします。
k-means法はクラスタリングアルゴリズムとしてよく知られていますが、このアルゴリズムではクラスタ数を所与としなければなりません。でも、今回はクラスタ数をどうすべきかよくわかりません。
1つの考え方として、まず適当なクラスタ数でk-means法を一旦適用してみて、次はそのクラスタにさらにk-means法を適用して、もしさらにk-means法を適用した細かいクラスタリングの方がデータをよりよく説明するならそちらを採用した方がよいと考えられます。データが重心の周りに正規分布すると仮定し、ベイズ情報量基準が小さくなる限りk-means法を繰り返し適用するのがX-means法です。
Determining the number of clusters in a data set - Wikipedia

手順

以下のように最低賃金合計特殊出生率、高齢化率を data.txt にまとめてから読み込みます。

北海道,北海道,764,1.29,0.281
青森,東北,695,1.43,0.290
岩手,東北,695,1.50,0.296
宮城,東北,726,1.31,0.246
(以下略)

x <- read.table("data.txt", sep=",", header=FALSE)
colnames(x) <- c("pref", "region", "wage", "birth", "aging")

3次元プロットしてみると以下です。なんかよくわかりません。

library(scatterplot3d)
windowsFonts(MEI=windowsFont("Meiryo"))
par(family="MEI", las=2)
scatterplot3d(x[,3:5], pch=19)

f:id:cookie-box:20170301131830p:plain:w440
2次元に投影したプロットもしておきます。

par(mgp=c(2.7, 0.7, 0))
par(mar=c(4,4,1,1))
layout(matrix(c(1,0,2,3), 2, 2, byrow=TRUE))
plot(x$wage, x$birth, family="MEI", pch=19)
plot(x$wage, x$aging, family="MEI", pch=19)
plot(x$birth, x$aging, family="MEI", pch=19)

f:id:cookie-box:20170301164808p:plain:w690

それで、このデータをX-means法でクラスタリングすると以下です。ただし前提条件はコードを参照してください。

x[3:5] <- scale(x[3:5])   # 最低賃金、合計特殊出生率、高齢化率をスケーリング
source("xmeans.R")
xkm <- xmeans$new()
xkm$xmeans(x[3:5])
print(xkm$Cluster)
plot(x[3:5], col=ifelse(xkm$Cluster==1, "salmon", "royalblue"), pch=19)

f:id:cookie-box:20170302000449p:plain:w640

上の結果より、最低賃金合計特殊出生率、高齢化率をスケーリングしてX-means法でクラスタリングすると都道府県は2クラスタに分割されることがわかります。グラフからわかる各クラスタの特徴は以下です。
最低賃金合計特殊出生率高齢化率
クラスタ1高い低い低い
クラスタ2低い高い高い


追記
上の結果を日本地図にプロットします。以下の記事を全面的に参考にさせていただきました。
sudillap.hatenablog.com
地図データとレコードの順序を合わせるために、都道府県名をアルファベットにしておく必要があります。

x <- read.table("data.txt", sep=",", header=FALSE)
colnames(x) <- c("pref", "region", "wage", "birth", "aging")
x <- x[order(x$pref),] # 地図データとの整合性のため県名アルファベット順にソート
x[3:5] <- scale(x[3:5])

source("xmeans.R")
xkm <- xmeans$new()
xkm$xmeans(x[3:5])
print(xkm$Cluster)

##### 日本地図のプロット #####
library(maptools)
library(spsurvey)

jpn <- readShapePoly("JPN_adm_shp/JPN_adm1.shp")
jpn.dbf <- read.dbf("JPN_adm_shp/JPN_adm1.dbf")
xlim <- c(126, 146) # 経度
ylim <- c(26, 46)   # 緯度

par(mar=c(0,0,0,0))
library(RColorBrewer)
plot(jpn, xlim=xlim, ylim=ylim, col=ifelse(xkm$Cluster==1, "salmon", "royalblue"))

f:id:cookie-box:20170302124428p:plain:w500

上の散布図とクラスタが微妙に違います。元々k-means法が初期値依存なのと、レコード順序も変わったため。
上のスクリプトを何回も実行するとクラスタが1つになったり3つになったりもします。


スクリプト

ほぼ参考文献の記事( Rでk-means法とその拡張2 x-means編 - サボタージュ禁止のおさぼり日記 )のままです。
ただ、k-means法のクラスタ初期化はk-means++法で固定、細分したクラスタをマージするかどうかはするで固定、共分散行列の非対角成分(共分散)は無視で固定しています。元記事のスクリプトは引数で変更できるようになっています。

kmeans.R

library(dplyr)

kmeans <- setRefClass(
  Class="kmeans", 
  methods=list(
    initCluster_kmeanspp=function(dat, NUM_CLUSTER) { ##### クラスタの初期化 (k-means++法による) #####
      index        <- 1:nrow(dat) 
      center_index <- sample(index, 1) # データから無作為に1点を選ぶ (これが1クラスタ目の中心)
      center_      <- unlist(dat[center_index,])
      min_dist     <- rowSums((sweep(dat, 2, center_, "-"))**2) # 各データとcenter_index番目のデータとの2乗距離
      while (length(center_index) < NUM_CLUSTER) {
        # 目的のクラスタ数になるまでデータ中心を追加
        # 収集済みのどのデータ中心からも遠いデータほど (距離の2乗に比例した確率で) 選ばれやすいようにする
        prob         <- min_dist[-center_index] / sum(min_dist[-center_index])
        new_center   <- sample(index[-center_index], 1, prob=prob) # データから1点を選ぶ
        center_      <- unlist(dat[new_center,])
        new_dist     <- rowSums((sweep(dat, 2, center_, "-"))**2)
        min_dist     <- pmin(min_dist, new_dist) # 各データと収集済みデータ中心の中で1番近いものとの2乗距離
        center_index <- c(center_index, new_center)
      }
      cur_center <- list()
      for (i in seq_len(NUM_CLUSTER)) {
        cur_center[[i]] <- unlist(dat[center_index[i],])
      }
      dist <- sapply(seq_len(NUM_CLUSTER),
                     function(k) {
                       center_ <- cur_center[[k]]
                       rowSums((sweep(dat, 2, center_, "-"))**2)
                     })
      cur_clusters <- max.col(-dist)  
      return(list(cur_center=cur_center,
                  cur_clusters=cur_clusters))
    },
    updateCenter=function(dat, cur_clusters) { ##### クラスタ分けに従ってデータ中心を更新 #####
      cur_center <-
        cbind(dat, cc=cur_clusters) %>%
        group_by(cc) %>%
        do(rtn=colMeans (dplyr::select(., -cc))) %>%
        as.list() %>%
        "$"(rtn)
      return (cur_center)
    },
    kmeans=function(dat_, K_=2) { ##### メイン処理 #####
      dat          <- na.omit(dat_[, sapply(dat_, is.numeric)])
      NUM_CLUSTER  <- as.integer(K_)
      SIZE         <- as.integer(nrow(dat))
      tmp          <- .self$initCluster_kmeanspp(dat, NUM_CLUSTER) # クラスタ初期化
      cur_clusters <- tmp$cur_clusters
      cur_center   <- tmp$cur_center
      old_clusters <- cur_clusters
      convergence  <- FALSE
      iter         <- 1L
      while (!convergence) {
        iter <- iter + 1L
        if (iter > 100L) {
          cat("iter reached MAX_ITER\n")
          break
        }
        # 各データの所属クラスタを更新 (最短距離のデータ中心へ)
        distances <- sapply(seq_len(NUM_CLUSTER), 
                            function(k) {
                              center_ <- cur_center[[k]]
                              d_      <- sweep(dat, 2, center_, "-")
                              return(rowSums(d_*d_))
                            })
        cur_clusters  <- max.col(-distances)
        if (length(unique(cur_clusters)) != NUM_CLUSTER) { # クラスタが消滅してしまった場合やり直し
          cur_clusters <- .self$initCluster_kmeanspp(dat, NUM_CLUSTER)$cur_clusters
        }
        convergence  <- identical(old_clusters, cur_clusters) # 所属クラスタの更新がなければアルゴリズム終了
        old_clusters <- cur_clusters
        cur_center   <- .self$updateCenter(dat, cur_clusters)
      }
      return(as.integer(cur_clusters))
    }
  )
)

k-means法を実行するときは例えば以下のようにする。

> source("kmeans.R")
> km <- kmeans$new()
> km$kmeans(x[3:5])

xmeans.R

library(matrixStats)
library(mvtnorm)
source("kmeans.R")

xmeans <- setRefClass(
  Class="xmeans",
  contains=c("kmeans"), 
  fields=c(
    "num_split_cluster" ="integer",
    "stack"             ="list",
    "Cluster"           ="integer"
    ),
  methods=list(
    initialize=function() {
      initFields(num_split_cluster=0L, stack=list())
    },
    xmeans=function(DATA_) { ##### メイン処理 #####
      data_ <- na.omit(DATA_[, sapply(DATA_, is.numeric)]) # 欠損レコード除去
      rownames(data_) <- 1:nrow(data_)
      .self$splitData(data_)   # まずデータを分割し切る
      .self$mergeSplitedData() # もしマージすべきだったらマージする
      .self$assignCluster()    # 各レコードのクラスタ番号を取得する
    },
    xBIC=function(dat_) { ##### ベイズ情報量基準 (分割前) #####
      dat <- dat_[, sapply (dat_, is.numeric)] 
      covar <- dat %>% 
               summarise_each(funs(var)) %>%
               diag() 
      df <- ncol(dat) * 2
      center <- colMeans(dat)
      bic <- -2 * sum(dmvnorm(dat, center, covar, log=TRUE)) + df * log(nrow(dat))
      return (bic)
    },
    xBICp=function(dat_, cluster_) { ##### ベイズ情報量基準 (分割後) #####
      dat <- dat_[, sapply(dat_, is.numeric)]
      cluster <- unlist(cluster_) 
      covar_eachCluster <- cbind(dat, cc=cluster) %>% 
                           group_by(cc) %>% 
                           summarise_each(funs(var)) %>%
                           dplyr::select(., -cc) %>%
                           rowwise() %>%
                           do(rtn=diag(.)) %>%
                           as.list() %>%
                           "$"(rtn)
      sumDeterminant <- sum(sapply(covar_eachCluster, function(mat) {prod(diag(mat))}))
      df <- 4 * ncol(dat)
      center_eachCluster <- cbind(dat, cc=cluster) %>%
                            group_by(cc) %>%
                            do(rtn=colMeans(dplyr::select(., -cc))) %>%
                            as.list() %>%
                            "$"(rtn)
      beta <- sqrt(c(crossprod(center_eachCluster[[1]] - center_eachCluster[[2]])) / sumDeterminant)
      logAlpha <- - log(2) - pnorm(beta, log.p=TRUE)
      logLikelihood_eachCluster <- sapply(seq_along(unique(cluster)), 
                                          function (k) {
                                            center_ <- center_eachCluster[[k]]
                                            covar_  <- covar_eachCluster[[k]] 
                                            dat_    <- dat[cluster==k, ]
                                            sum(dmvnorm(dat_, center_, covar_, log=TRUE))})
      bicp <- - 2 * (nrow(dat) * logAlpha + sum(logLikelihood_eachCluster)) + df * log(nrow(dat))
      return (bicp)
    },
    split2Cluster=function(dat_) { ##### データを2クラスタに分割する #####
      km       <- kmeans$new()
      cluster  <- km$kmeans(dat_, 2) # k-means で2クラスタにクラスタリング
      ccounts  <- table(cluster)     # 各クラスタに所属するデータ数を集計
      if (any(ccounts < 3)) {        # データ数3未満のクラスタができるようなら分割しない
        return (list(doSPLIT=FALSE))
      }
      bicp <- .self$xBICp(dat_, cluster=cluster)
      bic  <- .self$xBIC(dat_)
      if (bic > bicp) { # 分割後のベイズ情報量が小さくなるなら分割する
        return (list(doSPLIT=TRUE, d1=dat_[cluster==1,], d2=dat_[cluster==2,]))
      } else {
        return (list(doSPLIT=FALSE))
      }
    },
    merge2Cluster=function(dat_i, dat_j) { ##### 2クラスタをマージする #####
      par <- data.frame(len=c(nrow(dat_i), nrow(dat_j)), val=1:2)
      cluster <- par %>% 
                 rowwise() %>%
                 do(rval = rep(.$val, each=.$len)) %>%
                 "$"(rval) %>%
                 unlist()
      dat_ij  <- rbind(dat_i, dat_j)
      bicp    <- .self$xBICp(dat_ij, cluster)
      bic     <- .self$xBIC(dat_ij)
      if (bic < bicp) { # マージ後のベイズ情報量が小さくなるならマージする
        return(list(doMERGE=TRUE, D=dat_ij))
      } else {
        return(list(doMERGE=FALSE))
      }
    },
    splitData=function(data) { ##### データをクラスタに再帰的に分割する #####
      # simple x-means
      split_data <- split2Cluster(data)
      if (split_data$doSPLIT) {
        splitData(split_data$d1)
        splitData(split_data$d2)
      } else { # データがもうそれ以上分割されなかった場合は、そのデータを最終的なクラスタとしてメンバに書き込んで終わる
        num_split_cluster <<- num_split_cluster + 1L
        stack[[num_split_cluster]] <<- data
      }
    },
    mergeSplitedData=function() { ##### 分割し切ったクラスタにマージすべきものがあったらマージする #####
      ind_order           <- .self$getDataSize_order(.self$getDataSize_eachCluster())
      num_cluster         <- length(ind_order)
      isMerge_eachCluster <- numeric(num_cluster)
      for (i in seq_len (num_cluster - 1)) {
        i_ <- ind_order[sprintf ("%d", i)]
        if (isMerge_eachCluster[i_] == 1)
          next
        for (j in (i+1):num_cluster) {
          j_ <- ind_order[sprintf ("%d", j)]
           if (isMerge_eachCluster[j_] == 0) {
            merge_data <- merge2Cluster(stack[[i_]], stack[[j_]])
            if (merge_data$doMERGE) {
              stack[[i_]] <<- merge_data$D
              isMerge_eachCluster[j_] <- 1
              break
            }
          }
        }
      }
      stack <<- stack[isMerge_eachCluster==0]
      num_split_cluster <<- length(stack)
    },
    getDataSize_eachCluster=function() {
      out <- sapply(stack, nrow)
      return (out)
    },
    getDataSize_order=function(DataSize_eachCluster) {
      out <- seq_len(length(DataSize_eachCluster))
      out <- setNames(out, order(DataSize_eachCluster))
      return(out)
    },
    assignCluster=function() {
      out <- numeric(sum(.self$getDataSize_eachCluster()))
      ind <- lapply(stack, function(d){as.integer(rownames(d))})
      for (i in seq_along(ind)) {
        out[ind[[i]]] <- i
      }
      Cluster <<- as.integer(out)
    }
  )
)

これからの強化学習: ノート4

読んでいる本(出典): これからの強化学習 | 牧野 貴樹, 澁谷 長史, 白川 真一 |本 | 通販 | Amazon

この本読みですが、スライドにまとめようとしたら2回で挫折したのでとりあえず感想ノートの方を続けることにします。
というかどこまで読んだかもよくわからなくなりましたが、とりあえず落ち着いて現在の進捗を整理すると以下です。

テキスト範囲 読書メモ スライド
1.1~1.3節(1~41ページ) ノート1 勉強会#1
1.4~1.5節(42~70ページ) まだ 勉強会#2
2.1節(71~111ページ) ノート3 まだ
2.2~2.4節(112~147ページ) この記事 まだ

以下、読書メモ。2章以降は1節1節が真面目に読むと重いのでふわっと雰囲気だけ。スライド作成時に補完したい。

2.2節
まず、強化学習にはさまざまな解法アルゴリズムがあるけど、アルゴリズムの性能を評価できないの?というお話。

  • 1つの考え方としては、最善の選択をできた場合の収益に届かなかった分(リグレット)でもって評価する。
    多腕バンディットタスクの場合、
    • UCB1は、リグレットを \varepsilon 減衰版 \varepsilon-greedy を長期的に適用した場合よりも抑えられる。このアルゴリズムはゲーム木の探索やその他の探索問題にも応用できる。
    • Thompson サンプリング(下図;各腕を引いたときにコインが出てくる確率 \mu の事後分布を更新していき、常にこの事後分布からのランダムサンプリング結果が最大の腕を選択する)でも「不確かなときは楽観的に」が実現でき、リグレットをUCB1と同等に抑えられる。
      f:id:cookie-box:20170219140944p:plain:w390
      ※ 確率分布の形は適当です。
  • 別の考え方としては、間違った行動(最善の選択をした場合より期待収益が劣る行動)を取ってしまった回数で評価する(サンプル複雑性)。

あと、ちょうどUCB1や Thompson サンプリングで評価値の確率分布を更新していくように、探索と利用のトレードオフベイズ的に取り扱えるよね、というお話。

2.3節
報酬ってどう定義すればいいんだろう?というお話。

  • 例えばボードゲームの学習で勝ったときにのみ +1 という報酬の定義はあり得るが、ゲームの性質が以下のような場合、このように終端状態のみで報酬を定義するのは得策ではない。
    • 初期方策=ランダムな手で勝てる見込みが薄い(そもそも学習が進まない)。
    • 手数が多く、終端状態が初手の行動価値を修正する幅が小さい(学習が遅い)。
    • 自分が勝ったとき、自分の行動と相手の行動のどちらが寄与したのか明らかにはわからない。
  • 学習の結果を最善手に導いてくれるような報酬を定義したいんだけど、それなら最善手がわかっているときに、最善手が学習されるように(?)各状態に報酬を割り当てればいい(逆強化学習)。状態集合が有限集合なら、報酬関数の推定は線形計画問題になる。

2.4節
学習を速くしたいというお話。

  • 経験強化型学習では、報酬を得たときに、前回報酬を得て以降に現れた状態行動対の評価値を更新する。このとき、エピソード内の各ステップの状態行動対に報酬からどれだけ過去かに応じてポイントを割り当て、そのエピソードにおけるある状態行動対の評価は、その状態行動対に割り当てられたポイントの和とする(その状態行動対を4回通れば4回分のポイントの和を取る:138ページの図2.4.1)。過去にさかのぼるほど等比減少関数でポイントを割り当てれば、迂回路となる状態行動対を優先してしまうことがない(ということか?)。
    • Sutton本での単純なモンテカルロ法との違いは、単純なモンテカルロ法においては、ある状態行動対の評価はエピソード内の初回発生後の収益となるので、迂回路も平等に評価されてしまう。だからモンテカルロ法だと学習が遅い。ということだと思う。


  • 「リグレットは『探索コスト』を表す指標に見えるが(113ページ)」: 現実に「最善の選択をできた場合の収益」が出せないのは、現実には探索にコストをかける必要があるので、そのコストの分が差し引かれているように見える、という見方と考えられる。ただ、本当は、この箇所の直後に書かれているように、「探索と利用」から「探索」だけ切り離したようなものではない。
  • 「探索と利用を同時に実現する(113ページ)」: UCB1は、どれだけ探索できているか(=どれだけの精度で評価できているか)を加味して行動を選択するアルゴリズムだった(「これからの強化学習」勉強会#1 - クッキーの日記 7ページ)。その行動が未探索なほど(=評価の精度が粗いほど)より高く評価するようにしておくことで、利用のために取った行動で探索(精度向上)も実現している。
  • Thompson サンプリング(114ページ): フリーの原論文見つからず。以下のPDFが参考になる?
    https://bandits.wikischolars.columbia.edu/file/view/Lecture+4.pdf
  • 「(略)強化学習が優位な点としては、(略)たとえば、動的計画法のように環境のダイナミクスを表す状態遷移行列は不要であるし(127ページ)」動的計画法が強化学習の範疇ではないような書き方なんだけど、動的計画法は強化学習の一解法だと思っていたんだけど、モンテカルロ法やTD学習による解き方のみを強化学習とみなすような線引きもあるということ?

入門 機械学習による異常検知―Rによる実践ガイド: ノート3

読んでいる本(出典): 入門 機械学習による異常検知―Rによる実践ガイド | 井手 剛 |本 | 通販 | Amazon

前回: ノート2 / 次回: まだ
目次: 入門 機械学習による異常検知―Rによる実践ガイド

読んだページ: 44~52ページ
以下、メモと雑談。

前回までのあらすじ

  • データが多次元正規分布にしたがう場合、 a(x') = (x' - \hat{\mu}) ^{\rm T} \hat{\Sigma} ^{-1} (x' - \hat{\mu}) で異常度を定義できる。
    観測データが i.i.d. に多次元正規分布にしたがうなら、この異常度は F 分布にしたがう。
    • 例.(ノート2と同じ例)2017年1月1日~18日の東京の最低/最高気温は、1.63±1.99℃/10.47±2.77℃だった。
      気象庁|過去の気象データ検索
      このうち、2017年1月8日の東京の最低気温は1.6℃、最高気温は6.0℃だった。2017年1月8日は異常な日だったのだろうか。但し、ここで恣意的に、「14%も起こらないこと」を異常の閾値にする。
      → A. 最低気温と最高気温を1変数ずつみると異常ではないが、2変数の組合せは異常となる(以下)。
      • 2017年1月8日の最低気温だけみると1.6℃で、これはこの18日間の平均値1.63℃にかなり近いので異常値ではないだろう(証明略)。
      • 2017年1月8日の最高気温だけみると6.0℃とこの18日間で2番目に低いが、最高気温が独立同一正規分布にしたがうなら、「14%も起こらないこと」ではない(異常度 2.331211 < F分布の上側14%点 2.396771)。
        x1 <- c(2.0, 3.8, 3.5, 3.6, 3.7, 1.5, 0.1, 1.6, 3.8, 3.5, 3.9, 0.7, 1.4, -1.3, -2.3, -2.0, 0.7, 1.1)
        x2 <- c(13.8, 13.3, 13.7, 14.0, 10.4, 8.8, 8.7, 6.0, 11.1, 12.7, 11.0, 12.1, 12.7, 6.3, 4.7, 8.0, 10.9, 10.3)
        > a <- (18-1)/(18+1) * (x2[8]-mean(x2))^2 / (var(x2)*(18-1)/18)
        > a
        [1] 2.331211
        > F_0.14 <- qf(df1=1, df2=18-1, 0.86)
        > F_0.14
        [1] 2.396771
      • 2017年1月8日の最低/最高気温の組合せだと、最低/最高気温が独立同一2次元正規分布にしたがうなら、「14%も起こらないこと」といえる(F分布の上側14%点 2.228783 < 異常度 2.248451)。
        > mu <- colMeans(cbind(x1, x2))
        > xc <- cbind(x1, x2) - matrix(1, 18, 1) %*% mu
        > sigma <- t(xc) %*% xc / 18
        > a <- rowSums( (xc %*% solve(sigma) ) * xc)
        > a <- a * (18 - 2) / ((18 + 1) * 2)
        > a[8]
        [1] 2.248451
        > F2_0.14 <- qf(df1=2, df2=18-2, 0.86)
        > F2_0.14
        [1] 2.228783

  • なぜ異常度が F 分布にしたがうかという以前に、i.i.d. に M 次元正規分布にしたがう観測データからの、平均ベクトルと共分散行列の最尤推定値の導出の復習…。
    • 対数尤度は  \displaystyle L(\mu, \Sigma | x_{1:N}) = - \frac{MN}{2} \ln (2\pi) -\frac{N}{2} \ln |\Sigma| - \frac{1}{2} \sum_{n=1}^{N}(x_n - \mu)^{\rm T} \Sigma^{-1}(x_n - \mu) 。これを  \mu \Sigma について偏微分しないといけない。
      •  \mu について偏微分するには、行列の積のトレースの微分公式  \displaystyle \frac{\partial}{\partial A} {\rm Tr} (ABA^{\rm T}) = (B + B^{\rm T})A をつかうと  \displaystyle \frac{\partial L }{\partial \mu} = - \frac{1}{2} \sum_{n=1}^{N} \frac{\partial}{\partial \mu} {\rm Tr} \Bigl( (x_n - \mu)^{\rm T} \Sigma^{-1}(x_n - \mu) \Bigr) = - \Sigma^{-1} \sum_{n=1}^{N} (x_n - \mu)
        (∵ 対称行列  \Sigma逆行列もまた対称なので、 \Sigma^{-1} + (\Sigma^{-1})^{\rm T} = 2\Sigma^{-1}
      •  \Sigma について偏微分するには、正則な正方行列に関する微分公式  \displaystyle \frac{\partial \ln |A|}{\partial A} = (A^{-1})^{\rm T} と、
        ベクトル a と対称行列 B について成り立つ a^{\rm T} B a = \sum_i \sum_j b_{i,j} a_i a_j = {\rm Tr}(B a a^{\rm T}) と、
        行列の積のトレースの微分公式  \displaystyle \frac{\partial}{\partial B} {\rm Tr} (BA) = A^{\rm T} をつかうと
         \displaystyle \frac{\partial L}{\partial (\Sigma^{-1})} = \frac{N}{2} \frac{\partial}{\partial (\Sigma^{-1})} \ln |\Sigma ^{-1}| - \frac{1}{2} \sum_{n=1}^{N} \frac{\partial}{\partial (\Sigma^{-1})} (x_n - \mu)^{\rm T} \Sigma^{-1}(x_n - \mu)
              \displaystyle \, = \frac{N}{2} \Sigma - \frac{1}{2} \sum_{n=1}^{N} (x_n - \mu) (x_n - \mu)^{\rm T}
  • それで、なぜ異常度  a(x') = (x' - \hat{\mu}) ^{\rm T} \hat{\Sigma} ^{-1} (x' - \hat{\mu}) が F 分布にしたがうのか(44~48ページ)。
    • 証明の仕方は技巧的(47ページ)で、 \displaystyle a(x') = \frac{(x' - \hat{\mu}) ^{\rm T} \Sigma ^{-1} (x' - \hat{\mu})}{(x' - \hat{\mu}) ^{\rm T} \Sigma ^{-1} (x' - \hat{\mu}) / (x' - \hat{\mu}) ^{\rm T} \hat{\Sigma} ^{-1} (x' - \hat{\mu})} とすると分子も分母もカイ2乗分布にしたがうので、全体はF分布にしたがう、という仕組み。
    • 分母がカイ2乗分布にしたがう説明: 「分母は。もともと自由度 N-1 のウィシャート分布にしたがう  M \times M 行列を、ブロック分割行列の逆行列の公式(付録の定理 A.4 参照)を使って  1 \times 1 行列に縮約したものに対応している(48ページ)」とあるけど、定理 A.4 を参照してもここの意味がよくわからない…。なぜ逆行列の公式で縮約ができるのか…。線型代数は大学で習ったけどもう覚えていなかった…。
    • ウィキペディアには証明はなかった…原論文へのリンクは有。

今回のお話

  • 多次元のホテリング理論で「最高/最低気温の組合せが異常」というのは検知できるが、もっといえば、最高気温と最低気温のどちらが/あるいは両方が異常だったのかも知りたい。→ 正常であることが期待されるデータ群  \mathcal{D} から、全データの1変数あたりのマハラノビス距離を計算し、ここまでが正常という閾値を適当に決め、異常なデータの1変数あたりのマハラノビス距離をこの閾値と比較する。