確率論セミナー(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変数あたりのマハラノビス距離をこの閾値と比較する。

数理論理学: ノート4

読んでいる本(出典): 数理論理学 | 戸次 大介 |本 | 通販 | Amazon

前回: ノート3 / 次回: まだ
目次: 数理論理学

以下、第5章後半の自分の理解。
最後の節(標準形)がまだ読めていないので、後でこの記事に加筆するか、新しい記事にまとめる。

前回までのあらすじ(一階述語論理の統語論

  • 命題論理では扱えない「すべての x について~」「ある x が存在して~」のような知見を表現すべく、論理式を名前と述語に分けます。
  •  \forall \exists といった量化子を導入します。
  • 論理式がすごい増えてしまいそうですが、自然数で番号付けできます(ゲーデル数)。
    • 但し、論理式のパーツになる記号の集合  \mathcal{A} は高々可算集合とする。

一階述語論理の意味論

  • 命題論理のときと同様、「 I: 論理式の集合  \rightarrow 真偽の領域 D_t \equiv \{1, 0\}」であるような解釈  I を考える。
    命題論理のときは、 I は「原子命題の真偽の組合せの全パターン」ということで  2原子命題の数 個あった。
    一方、一階述語論理では原子命題がさらに項と述語に分かれるので、「各項が指し示しているもの」「各述語について、項が指し示しているものに対応する、論理式の真偽」として考えるのがよい(構造  M)。
    あと変項もあるので、変項に何を割り当てたかも論理式の真偽にかかわる(割り当て  g)。
    つまり、ある1つの解釈  I = \langle M, g \rangle を与えるには、以下を完全に決めればよい。
    • そもそもどんなものたちについて命題を述べるのか(存在物)の集合:  D_M
    • 名前記号が指す存在物は何か:  F_M : \{ a_1, a_2, \cdots \} \to D_M
    • n個の存在物に n演算子を適用したときに指す存在物は何か:  F_M : \{ o_{n,1}, o_{n,2}, \cdots \} \to {D_M}^{{D_M}^n}
    • n個の存在物に n項述語を適用したときの真偽はどうか:  F_M : \{ \theta_{n,1}, \theta_{n,2}, \cdots \} \to {D_t}^{{D_M}^n}
    • 変項に何を割り当てたか:  g : \{\xi_1 , \xi_2, \cdots \} \to D_M

Ex.あるゲームのキャラクターたちについて命題を述べるとする(以下、ネタバレな上にしようもない例)。
    M_1ダンガンロンパのキャラクターの構造とすると、 D_{M_1}=\{苗木誠, 舞園さやか, 桑田怜恩, \cdots \}
   名前記号から存在物への対応付けは例えば、 F_{M_1}(a_1)= 苗木誠 ,  F_{M_1}(a_2)= 舞園さやか , \cdots
   例えば  o_{1,1}(\tau_1) が「 \tau_1 を殺した犯人」のような1項演算子なら、 F_{M_1}(o_{1,1})( 舞園さやか ) = 桑田怜恩 。
   例えば  \theta_{2,1}(\tau_1, \tau_2) が「 \tau_1 \tau_2 を殺した」のような2項述語なら、 F_{M_1}(\theta_{2,1})( 苗木誠, 舞園さやか ) = 0
   このように、各名前記号が指す存在物、各n演算子が指す存在物(項が指す存在物による)、各n項述語の
   真偽(項が指す存在物による)を完全に決めれば、この論理体系のすべての基本述語の真偽は完全に決まる。
   ※ ここで、「桑田怜恩が舞園さやかを殺したから  F_{M_1}(o_{1,1})( 舞園さやか ) = 桑田怜恩」ということではなく、
     あくまで  o_{1,1}(\tau_1) という1項演算子をこの写像で定義したということに注意。
   ※ ただ、  o_{1,1}(\tau_1) が本当に「 \tau_1 を殺した犯人」になるように定義されていて、  \theta_{2,1}(\tau_1, \tau_2) が本当に「 \tau_1
      \tau_2 を殺した」かどうかを表すように定義されているなら、 \theta_{2,1}(o_{1,1}(\tau_1), \tau_1) はどの解釈でも真になりそう。
     殺されていないキャラクターを代入したときの取り扱いは適切にする必要があるけど。
   他方、別の解釈の下では、論理体系にはスーパーダンガンロンパ2のキャラクターの構造  M_2 が当てはめられて
   いるかもしれない。 D_{M_2}=\{日向創, 七海千秋, \cdots \}
   この場合の名前記号の対応付けは例えば、 F_{M_2}(a_1)= 日向創 ,  F_{M_2}(a_2)= 七海千秋 , \cdots かもしれない。

  • 量化論理式 \forall \xi \varphi, \; \exists \xi \varphi の真偽を解釈するには、 \xi への割り当てをすべての存在物に変異させる(というより、これが \forall \xi \varphi, \; \exists \xi \varphi という論理式の定義そのものである)。
    •  \langle M, g \rangle (\forall \xi \varphi)=1 \; \Longleftrightarrow \; すべての  a \in D_M について  \langle M, g \rangle [ \xi \mapsto a] (\varphi) = 1 である.
    •  \langle M, g \rangle (\exists \xi \varphi)=1 \; \Longleftrightarrow \; \langle M, g \rangle [ \xi \mapsto a] (\varphi) = 1 となる  a \in D_M が存在する.
  • 解釈について、以下が成り立つ。
    • その変項が項/論理式の自由変項でないなら、その変項を変異させても項が指す存在物/論理式の真偽は変わらない。
    • 項/論理式の変項にある項を代入した上で解釈しても、解釈の割り当ての変異によって変項をある項に変えても、どちらの場合も項が指す存在物/論理式の真偽は変わらない。
    •  \tau_1 \tau_2 が指す存在物が同じなら、項/論理式の変項を  \tau_1 に変異させても  \tau_2 に変異させても項が指す存在物/論理式の真偽は変わらない(外延性)。

一階述語論理における妥当な推論

  • 意味論的含意の概念は一階命題論理のときと同じ。
  • 量化論理式が推論の前提 or 帰結になる場合に色々な定理が成り立つ。



練習問題5.48
n 項述語に渡せる項のセットのパターン数は、 |D_M|^n パターンある。
それぞれのパターンが1か0を取りうるので n 項述語は  2^{|D_M|^n} 種類ありうる。
例えば、 D_M =\{A_1, A_2\} で、 n=2 だったら、2項述語  \theta_{2, i} は16種類ある。
\tau_1 \tau_2 \theta_{2,1} \theta_{2,2} \theta_{2,3} \theta_{2,4} \theta_{2,5} \theta_{2,6} \theta_{2,7} \theta_{2,8} \theta_{2,9} \theta_{2,10} \theta_{2,11} \theta_{2,12} \theta_{2,13} \theta_{2,14} \theta_{2,15} \theta_{2,16}
 A_1  A_1 11111111 00000000
 A_1  A_2 11110000 11110000
 A_2  A_1 11001100 11001100
 A_2  A_2 10101010 10101010
練習問題5.55

  • 「奇数の二乗は奇数である」は真だが、 \forall x ( 奇数 (x) \; \land 奇数  (x \times x) \, ) は、 x に例えば 2 を割り当てたとき偽。
  • 「二乗が奇数であるような奇数が存在する」は真で、 \exists x ( 奇数 (x) \to 奇数  (x \times x) \, ) も真だが、前者は 2 を割り当てると偽で、後者は 2 を割り当てても真(というか後者はどんな数を割り当てても真)。

練習問題5.58 面倒なので略。
練習問題5.60 面倒なので略。
練習問題5.62 面倒なので略。
練習問題5.65 面倒なので略。
練習問題5.82 面倒なので略。
練習問題5.84 面倒なので略。
練習問題5.86 面倒なので略。
練習問題5.88 面倒なので略。
練習問題5.90 面倒なので略。
練習問題5.91
 \forall x (F(x)) \land G(a) = \! \! \! | | \! \! \! = \forall x (\lnot \lnot F(x)) \land G(a) を証明するには、二重否定律より  F(x) = \! \! \! | | \! \! \! = \lnot \lnot F(x) なので、これに量化論理式の置き換えを用いて  \forall x (F(x)) = \! \! \! | | \! \! \! = \forall x (\lnot \lnot F(x)) 。これと  G(a) = \! \! \! | | \! \! \! = G(a) に対し複合論理式の置き換えを適用。


所感

  • x は哺乳類である」「x は卵生である」の変項 x に「友達の鈴木君」を割り当てる(98ページ): 確かに友達の鈴木君は哺乳類だし、卵生ではないけど、鈴木君的には友達からの扱いが「哺乳類の一個体」でいいのだろうか…。
  • 集合における外延性は、「集合はそれが含む要素によって一意に定まる」。外延性の公理 - Wikipedia
    一階述語論理の解釈の外延性(定理5.64)は、「論理式の真偽は項が指す存在物によって一意に定まる」。当たり前に感じられすぎてよくわからない。「定義されたこと以外は知らないふりをする《知らないふりゲーム》」(数学ガール ゲーデル不完全性定理の31ページ)は、こと「ある論理式が真か偽か」については経験がありすぎるので難しい。
  • ラムダ計算(114ページ)って何。