TIM Labs

sato.yuichiroによるエントリー一覧

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第4章、GLMのモデル選択についてまとめます。

この章ではバイアスやAICについて説明がされています。 データを説明するモデルとして様々なモデルが考えられる時に、データへのあてはまりの良さだけを基準にモデルを選ぶことはできません。 それを理解するために、ポアソン回帰で推定されるパラメータがどれくらい真の値からずれるのかを見るためのコードを用意しました。 コードはRで書きました。

legends <- c("ideal", "regression")
ltys <- c("dashed", "solid")
for (i in 1:10) {
  y <- rpois(5, lambda = 4)
  data <- data.frame(y)

  fit <- glm(formula = y ~ 1, family = poisson, data = data)
  l <- exp(fit$coefficients[1])

  title <- paste("Poisson regression, estimated lambda = ", l)
  plot(1:10, dpois(1:10, lambda = 4), ylim = c(0.0, 0.25), type = "l", lty = "dashed", main = title, xlab = "", ylab = "")
  lines(1:10, dpois(1:10, lambda = l), type = "l")
  legend("topright", legend = legends, lty = ltys)
}

平均 4 のポアソン分布から5つのサンプルをサンプリングして、そのデータに対してポアソン回帰をしています。 リンク関数は定数です。 正確に推定されるなら、推定されるパラメータは 4 に近い値になるはずです。 実際に実行してみると、5個しかサンプルが無いわりには結構正確に推定されるのが分かります。 しかし、だいたい10回も推定すると、1回くらいは目に見えておかしな推定をするのが分かります。 例えば下図のような推定をします。

poisson_regresion_estimated4.png

真の値は 4 なのに、推定結果は 5.4 です。 このように、与えられたデータに対して最尤推定するしかない以上、推定は推定でしかなくどうしても真の値からずれてしまいます。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第3章、一般化線形モデル(GLM)についてのまとめの三回目です。

この章ではポアソン回帰について説明がされています。 ポアソン回帰によって得られるのは、説明変数と応答変数の分布の関係です。 変数が変化すると別の変数の分布が変化するという、結構複雑な関係が得られるわけです。 分布そのものにも興味はありますが、やはり一番興味があるのは、分布の中でどこが一番確率が高くなるかです。 そこで、本の中で使われているサンプルデータを使ってポアソン回帰をやってみて、得られた分布の中の一番確率が高い場所をプロットしてみました。 コードはRで書きました。

d <- read.csv("data3a.csv")

fit <- glm(y ~ x, data = d, family = poisson)

b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]

linkFunction <- function(x) {
  exp(b1 + b2 * x)
}

fitted <- function(x, y) {
  dpois(y, lambda = linkFunction(x))
}

cat("p(y | exp(a + bx)) where p is Poisson distribution\n")
cat("a =", b1, "b =", b2, "\n\n")

y <- 1:20
max <- numeric(length = 20)
for (x in 1:20) {
  at <- fitted(x, y)
  cat("when x = ", x, ", y = ", which.max(at), " is most probable, probability is ", max(at), "\n")
  max[x] <- which.max(at)
}
title <- paste0("p(y|exp(", b1, " + ", b2, "x))")
xl <- "x"
yl <- "the most probable y"
plot(1:20, max, type = "b", main = title, xlab = xl, ylab = yl)

実行すると以下のようなグラフが得られます。

most_probable.png

サンプルデータの説明変数 x と応答変数 y には正の相関があるので、 x の増加に従って y の分布は大きな値を取るようになります。 グラフを見ると一番高い確率で得られる y の値が x の増加に従って増えているのが分かります。 正の相関の一面が見られますね。

リンク関数は指数関数なので、もっときれいに指数関数的増加が見られないかとプロットする x の範囲を広げてみました。

most_probable2.png

細かい増加は分かりにくくなりましたが、全体的な指数関数的増加がよく分かるようになりました。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第3章、一般化線形モデル(GLM)についてのまとめの二回目です。

この章にはポアソン回帰の例が載っています。 本の中で使われているサンプルデータがここで公開されているので、解説されている通りにポアソン回帰をやってみました。 コードはRで書きました。

d <- read.csv("data3a.csv")
cat("summary of input data\n")
summary(d)

fit <- glm(y ~ x, data = d, family = poisson)

b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]

cat("\nlink function\n")
cat("exp(", b1, " + ", b2, "x)\n")

linkFunction <- function(x) {
  exp(b1 + b2 * x)
}

ps <- 0:20
xl <- "y"
yl <- "probability"
legends <- c("x = 0","x = 5","x = 10", "x = 15", "x = 20")
pchs <- c(0, 5, 10, 15, 20)
title <- paste("lambda = exp(", b1, "+", b2, "x)")
plot(0, 0, type = "n", xlim = c(0, 20), ylim = c(0.0, 0.25), main = title, xlab = xl, ylab = yl)
for (x in seq(0, 20, by = 5)) {
  d <- function(y) {
    dpois(y, lambda = linkFunction(x))
  }
  lines(ps, d(ps), type = "l")
  points(ps, d(ps), pch = x)
}
legend("topright", legend = legends, pch = pchs)

実行するとポアソン回帰を行い、得られる分布をプロットしてくれます。 実際にポアソン回帰をやっているコードは以下の部分です。

fit <- glm(y ~ x, data = d, family = poisson)

b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]

b1b2 が最尤推定されたパラメータです。 推定された結果は b1 = 1.291721b2 = 0.07566191 でした。 これらのパラメータが得られればリンク関数 λ = exp(b1 + b2 x) が得られるので、説明変数 x に対して y がどのような平均 λ でポアソン分布するかが分かります。 その分布は下図のようになります。

glm_result.png

グラフを見ると x が大きくなるにしたがって y の分布が大きな値をとるように移動しているのが分かります。 これは、推定されたパラメータのうち、 b2 が正だからです。 サンプルデータの xy には正の相関があるのですね。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第3章、一般化線形モデル(GLM)についてまとめます。

この章では主にポアソン回帰について解説されています。 ポアソン回帰というのは、ポアソン分布の平均がリンク関数という関数を通してデータと関連するモデルを使った回帰です。 モデルが複雑になったので、この回帰で得られる物が何なのかをグラフに描くコードを用意しました。 コードはRで書きました。

linkFunction <- function(b1, b2) {
  function(x) {
    exp(b1 + b2 * x)
  }
}

distribution <- function(x, linkFunction) {
  l <- match.fun(linkFunction)
  function(y) {
    dpois(y, lambda = l(x))
  }
}

ps <- 0:10
xl <- "y"
yl <- "probability"
legends <- c("x = 0","x = 5","x = 10", "x = 15", "x = 20")
pchs <- c(0, 5, 10, 15, 20)
for (b1 in seq(-0.1, 0.1, by  =  0.05)) {
  for (b2 in seq(-0.1, 0.1, by = 0.05)) {
    l <- linkFunction(b1, b2)
    title <- paste("lambda = exp(", b1, "+", b2, "x)")
    plot(0, 0, type = "n", xlim = c(0, 10), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl)
    for (x in seq(0, 20, by = 5)) {
      d <- distribution(x, l)
      lines(ps, d(ps), type = "l")
      points(ps, d(ps), pch = x)
    }
    legend("topright", legend = legends, pch = pchs)
  }
}

この例では、説明変数xに対してyが分布する場合について考えています。 リンク関数を作るコードは以下のようになります。 パラメータb1b2を与えると、そのパラメータを使ったリンク関数を返しています。

linkFunction <- function(b1, b2) {
  function(x) {
    exp(b1 + b2 * x)
  }
}

リンク関数は引数xを取って、パラメータと足したり掛けたりした後、累乗を計算する関数です。 ポアソン回帰ではリンク関数の戻り値を、ポアソン分布のパラメータとして使います。

リンク関数があれば、あるxに対するyの分布が得られます。 つまり、yを引数に取りyが起きる確率を返す関数が得られます。 以下の関数により得られる関数でyの起きる確率が計算されます。

distribution <- function(x, linkFunction) {
  l <- match.fun(linkFunction)
  function(y) {
    dpois(y, lambda = l(x))
  }
}

これがポアソン回帰のモデルです。 後はリンク関数で使われるパラメータを最尤推定で求めれば良いわけです。

今はとりあえずモデルの具体例が見たいだけなので、パラメータは自分で与えて、様々な値についてグラフをプロットしてみました。 コードを実行すると、いろいろな場合についてグラフがプロットされます。 具体例として、観ていて面白いのはb1 = -0.1b2 = 0.1の場合です。 以下に図を貼ります。

glm_b1_-0.1b2_0.1.png

リンク関数はexp(-0.1 + 0.1x)になります。 xがある数だった場合、yはこの値を平均に持つポアソン分布で分布するということです。 例えばx = 1の場合は1です。 x = 11の場合はe = 2.718...です。

b2が正なのでxが大きくなるほどリンク関数の指数の肩は大きな数になります。 つまり、大きなxではyは大きな平均を持つポアソン分布をします。 図を見ると、確かにxが大きくなるにしたがってyの分布は大きな値を取るようになっています。

このように、ポアソン回帰で得られるのはxの変化によってyの分布がどのように変化するかという傾向です。 ただ変化すると言うだけでなく、どれくらいの分布を伴って変化するかも分かります。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第2章、確率分布と統計モデルの最尤推定についてまとめます。

この章では主にポアソン分布と、ポアソン分布のパラメータの最尤推定について解説されています。 ポアソン分布とは、下図のような分布です。

pois5.png

ポアソン分布はパラメータを一つ持っていて、それが決まると分布が決まります。

最尤推定とは、あるデータが得られている時、そのデータが生成される確率が最大になるように、分布のパラメータを選ぶ推定方法です。 データは何でもいいですが、現実の世界で測定された数値です。 その数値をモデルを使って説明しようというのですから、実際に得られているデータに最もよく適合するモデルが欲しいと思うのは当然です。 最もよく適合するというのは、確率の言葉で言えば最大の確率を与えるということなので、このようにパラメータを推定するのが合理的ですね。

さて、問題が単純な場合には、解析的に確率を最大化するパラメータを求めることができるのですが、複雑な問題では普通は数値計算で求めます。 その様子を見やすくするデモを用意しました。 コードはRで書きました。

n <- runif(1, 1, 10)
data <- rpois(20, lambda=n[1])

cat("input data taken from Poisson distribution (lambda = ", n[1], ")\n\n")
data
summary(data)
cat("\n")

cat("likelihood of given data at\n")
res <- vector(length = 10)
for (i in 1:10) {
  prob <- function(x) {
    dpois(x, lambda=i)
  }
  res[i] <- sum(log(prob(data)))
  cat("lambda = ", i, " ", res[i], "\n")
}

cat("\nmaximum likelihood = ", max(res), "when lambda = ", which.max(res), "\n")

このコードが何をやっているか詳しく解説します。 まず最尤推定に使うデータを用意しています。 以下のコードで 1 から 10 までの乱数を一つ用意して、その乱数を平均に持つポアソン分布からサンプルを 20 用意しています。

n <- runif(1, 1, 10)
data <- rpois(20, lambda=n[1])

後は 1 から 10 までの整数 i について、このデータが平均 i を持つポアソン分布から得られる確率を計算しています。 確率を計算しているのは以下のコードです。

prob <- function(x) {
  dpois(x, lambda=i)
}
res[i] <- sum(log(prob(data)))

dpois(x, lambda=i)λ=i のポアソン分布から x が得られる確率です。 sum(log(prob(data)))data に含まれる x の全てについて、 dpois(x, lambda=i) を計算し対数を取って和を計算しています。確率の積をそのまま計算すると値が小さくなりすぎて計算できないので対数を取って和を計算しました。ちなみにこの確率には尤度という名前がついています。

後は結果をコンソールに出力しています。 実行結果は次のようになります。

input data taken from Poisson distribution (lambda =  2.923534 )

 [1] 5 3 4 6 5 4 0 2 6 3 6 0 3 4 7 2 1 1 2 2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     2.0     3.0     3.3     5.0     7.0 

likelihood of given data at
lambda =  1   -75.51993 
lambda =  2   -49.77221 
lambda =  3   -43.01152 
lambda =  4   -44.0245 
lambda =  5   -49.29702 
lambda =  6   -57.2638 
lambda =  7   -67.08986 
lambda =  8   -78.27679 
lambda =  9   -90.5031 
lambda =  10   -103.5493 

maximum likelihood =  -43.01152 when lambda =  3

データを生成するのに使われたパラメータは 2.923534 です。 なので、尤度はそれに近い 3 の時に最大になっています。 何度か実行してみると、データを生成するのに使ったポアソン分布の平均に近い時に尤度が一番大きくなっていることが分かると思います。 データの本当の分布は最初に引いた乱数を平均のポアソン分布なのだから、その乱数に近い推定を行った時に、最も当てはまりがよくなるのは順当な結果です。 このように最尤推定を行えば、統計処理したいデータに対し、順当な統計分布を見つけることができるのですね。

このアーカイブについて

このページには、sato.yuichiroが最近書いたブログ記事が含まれています。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。