gepuro.net
gepulog

データ分析エンジニアによる備忘録的ブログ

統計的学習の基礎4.4~を読みました

統計的学習の基礎読書会で発表してきました。以下がその時の資料です。

統計的学習の基礎 4.4~ from gepuro

確率変数の変換

何度か勉強しているのだけれど、毎回確認しているのでブログに書いておく。

数式転がし

確率変数 $X$ に対して、その関数の確率分布を求める。(例: $X^2$, $\log(X)$)

変換の関数を$Y=\phi(X)$とする。$\phi(x)$は単調増加。

$\phi(x)$は単調増加であるから、$x \leq X \leq x + \Delta x$と$Y \leq Y \leq y + \Delta y$は論理的に同値であるので、

\begin{align*} P(y \leq Y \leq y + \Delta y) = P(x \leq X \leq x + \Delta x) \end{align*} 確率変数$X$, $Y$の密度関数を$f(x)$, $g(y)$とおく。

ここで、$P(x \leq X \leq x + \Delta x) = f(x) \Delta x$であるから、

\begin{align*} g(y) = f(x) \Delta x / \Delta y \approx f(x) ( dx / dy ) \end{align*} となる。

関数$y=\phi(x)$をxについて解いた逆関数を$x=\psi(y)$とすると、

\begin{align*} g(y) = f(\psi(y)) * | d \psi(y) / d y | \end{align*}

として求められる。

確率変数$X \sim uni(0, 1)$に対して、 $y = x^2$(ただし$0 \leq x$)と変換する時の密度関数を求める。

$y = \phi(x) = x^2$の逆関数は、$x = \psi(y) = \sqrt{y}$となる。

よって、

\begin{align*} g(y) &= f(\psi(y)) * | d \psi(y) / d y | \\ & = 1 * d / dy \sqrt{y} \\ & = \frac{1}{2 \sqrt{y}} \end{align*}

となる。

参考

時系列解析の入門に良い資料

オハイオ州立大学のRobert de Jong先生が公開されている講義資料?がよさ気なので、紹介します。

誤植

部分的にしか読んでいないけど、誤植がありました。他にも誤植がありそう。

Lecture 7

式(3)

\begin{align*} \hat{\beta_n} = \left[\sum_{t=1}^n x_t x_t ' \right]^{-1} \left[ \sum_{t=1}^n x_t ' y_t \right] \end{align*}

だと思われます。(参考:自然科学の統計学)式(4)も転置記号が抜けてますね。

Proposition 1の1つ前の式

\begin{align*} Q^{-1} H_n^{-1}\left(\sum_{t=1}^n x_t u_t \right ) = Q^{-1} \begin{bmatrix} n^{-1/2} \sum_{t=1}^n u_t \\ n^{-3/2} \sum_{t=1}^n (t/n) u_t \end{bmatrix} \end{align*}

かな。そうじゃないと、辻褄が合わないかと。

時系列解析入門 1章

北川先生の時系列解析入門の読書記録です。

1.1 時系列データ

時間の経過とともに不規則に変動する現象の記録が時系列

1.2 時系列の分類

連続時間時系列と離散時間時系列

  • 連続時間時系列
    • 例:アナログレコーダ
  • 離散時間時系列
    • 例:一時間おきに計測された気圧

一変量時系列と多変量時系列

  • 一変量
    • 観測時点で一種類の情報
  • 多変量
    • 2つ以上の情報

定常時系列と非定常時系列

  • 時間的に変化しない一定の確率的モデルの実現値とみなせる
    • 定常時系列
  • 平均が時間とともに変動、平均のまわりの変動の仕方が時間的に変化
    • 非定常時系列

ガウス型時系列と非ガウス型時系列

  • 時系列の分布が正規分布に従う
    • ガウス型時系列
  • 時系列の分布が正規分布に従わない
  • 非ガウス型時系列

非ガウス型時系列のように上下非対称性が見られても、データを変換することによって、近似的にガウス時系列とみなせる場合がある

線形時系列と非線型時系列

  • 線形なモデルの出力として表現できるような時系列
    • 線形時系列
  • 非線形なモデルが必要な時系列
    • 非線型時系列

欠測値と異常値

  • 何らかの理由で観測値が記録されない
    • 欠測値
  • 現象の異常な振る舞い・観測装置の異常など
    • 異常値

時系列解析の目的

記述

図示、標本自己共分散関数、標本自己相関関数、ピリオドグラムなど

モデリング

  • 与えられた時系列の変動の仕方を表現する時系列モデルの構成
  • 時系列の確率的構造の解析

予測

時系列が互いに相関を持つことを利用して、現在までに得られた情報から今後の変動を予測すること

信号抽出

時系列から目的に応じて必要な信号や情報を取り出すこと

時系列データの前処理

非定常なデータを定常化する

変数変換

ロジット変換

$z_n = log \left( \frac{y_n}{1-y_n} \right)$ より一般的な方法として、Box-Cox変換がある

差分(階差)

$z_n = \Delta y_n = y_n - y_{n-1}$

$y_n = a + bn$であれば、 $z_n = \Delta y_n = b$となり直線の傾きが除去できる

y_n = a + bn + cn^2$ならば、

\begin{align*} $\Delta z_n & = z_n - z_{n-1} \\ & = \Delta y_n - \Delta y_{n-1} \\ & = (b + 2cn) - (b +2c(n-1)) \\ & = 2c \end{align*} となり2次成分及び1次成分を除去できる。

年周期であれば

\begin{align*} \Delta_p y_n = y_n = y_{n-p} \end{align*} として一周期前の値との差分を利用することがある。

前期比, 前年同期比

経済時系列などでは、減系列の前期比や前年同期比を分析することがある

$z_n = \frac{y_n}{y_n - 1}$

$x_n = \frac{y_n}{y_{n-p}}$

  • 時系列$y_n$が$y_n=T_n w_n$とトレンドとノイズの積
  • $T_n = (1+\alpha) T_{n-1}、$\alpha$は成長率

の場合、

\begin{align*} z_n & = \frac{y_n}{y_{n-1}} \\ & = \frac{T_n w_n}{T_{n-1} w_{n-1}} \\ (1 + \alpha) \frac{w_n}{w_{n-1}} \end{align*} で、成長率$\alpha$を検出出来る

  • 周期$p$, 周期関数$s_n$とノイズの積
  • $y_n = s_n \dot w_n$, $s_n = s_{n-p}$

の場合、

\begin{align*} x_n & = \frac{y_n}{y_{n-p}} \\ & = \frac{s_n w_n}{s_{n-p} w_{n-p}} \\ & = \frac{w_n}{w_{n-p}} \end{align*} となり周期分布を除去出来る

移動平均

$2K+1$項の移動平均は $ T_n = \frac{1}{2k+1} \sum_{j=-k}^k y_{n+k}$で定義

直線とノイズの和で表される場合

  • $y_n = t_n + w_n$
  • $t_n ~ a + bn$

$T_n$の平均は$t_n$と同じで分散は$w_n$の分散の$1/(2k+1)$となる。

  • 重み付き移動平均
    • $T_n = \sum_{j=-k}^{k} w_j y_{n-j}$

最強のハードディスクはどれだ?

第46回R勉強会@東京(#TokyoR)でLTをしてきました。 資料は、こちらです。

最強のハードディスクはどれだ? from Atsushi Hayakawa

また、https://github.com/gepuro/gweibullplotでワイブルプロットのパッケージを公開しました。本資料を作成する際に用いた関数とは別のものです。(研究室秘伝のソースでライセンスが謎なため)

ハードディスクの寿命分布を比較

Hard Drive Data Setsに4万台以上のハードディスクに関するデータが公開されている。これらの寿命分布を推定してみようと思う。

公開されているデータをSQLiteに展開し、前処理

create table total_days as
select serial_number, model, count(*) as total_days,
max(failure) as failure,
min(date(date)) as min_date, max(date(date)) as max_date
from drive_stats
group by serial_number
;

Rで読み込み、解析準備

library(RSQLite)
library(dplyr)

drv<-dbDriver("SQLite")
con<-dbConnect(drv,dbname="drive_stats.db")
totaldays <- dbGetQuery(con,"select * from total_days")
head(totaldays)

totaldays %>%
  group_by(model) %>%
  summarise(model_sum = n()) %>%
  arrange(model_sum)

hdd.df <- 
  totaldays %>%
  filter(model %in% c("ST4000DM000", "HGST HMS5C4040ALE640"))

weib.df <- data.frame(
  x = hdd.df$total_days,
  condition = ifelse(hdd.df$model == "ST4000DM000", 1, 2),
  status = hdd.df$failure
  )

ワイブルプロットでパラメータの推定

ワイブルプロット

  • 青色 : ST4000DM000, mhat = 0.7875, etahat = 3.550e+04
  • 赤色 : HGST HMS5C4040ALE640, mhat = 0.7628, etahat = 1.193e+05

です。

平均寿命で比較すると、 ST4000DM000は40685日で、HMS5C4040ALE640は140100日だった。

また、B1ライフ(1%壊れるまでの日数)は、 ST4000DM000は103.1日で、HMS5C4040ALE640は286.8日だった。

確率密度関数をプロット

pdf

こんな感じでHDDの寿命をモデル毎に比較出来る。(≧∇≦)/

今回対象としたモデルは、データ数が多い2つに注目したわけだが、こんなにも違いが出るなんて驚いた。

指数分布の確率プロット

2年ぐらい前にお勉強用に書いたコードを再利用するために発掘したので、共有しておきます。 "exp"の箇所を好きな分布に書きかえれば流用できるかと思います。

exponential.plot <- function(x){
  qq <- function(x){
    n <- length(x)
    qexp((1:n - 1/2)/n, rate=mean(x))[order(order(x))] # ミラーランク
  }

  probs<- c(0.001, 0.01, 0.1, 1, 5, 10, 20, 30, 50, 70, 80, 90, 95, 99, 99.9, 99.99, 99.999) / 100

  plot(c(min(x), max(x)),  qexp(c(probs[1], probs[length(probs)])), type="n", xaxt="n", yaxt="n", 
       ylab="Cumulative Percent", xlab="Observed Value")
  abline(h=qexp(probs), v=as.integer(seq(min(x), max(x))), col="grey")
  abline(h=qexp(0.5))
  axis(side=2, at=qexp(probs[c(4,8:17)]), las=1, labels=as.character(probs[c(4,8:17)]*100), cex.axis=0.8)
  axis(side=1, as.integer(seq(min(x), max(x))), cex.axis=0.8)
  points(x, qq(x))

  model <- lm(qq(x) ~ x - 1)
  abline(model)
  print(paste("r-squared",summary(model)$r.squared,sep=":"))
  # Maximum likelihood estimation
  print(paste("rate",mean(x),sep=":"))
}

rexp.data <- rexp(100,rate=1)
exponential.plot(rexp.data)

出力結果は、

exponential_plot

となります。

参考

条件付き推測木のサンプルデータ生成と推定

条件付き推測木(決定木)で生存時間解析

寿命分布が異なるように説明変数で分割する。

詳細はhttp://cran.r-project.org/web/packages/party/party.pdfで。

データ生成

set.seed(0)
n <- 1000
x <- rbinom(n, size=1, prob=1/3)
data.df <- data.frame(x = x,
                      tt = ifelse(x==0, rweibull(n, scale=200, shape=2), rweibull(n, scale=300, shape=3)),
                      status = rep(1,n)
                      )

推定

> library(party)
> (model <- ctree(Surv(tt, status)~x, data=data.df))

     Conditional inference tree with 2 terminal nodes

Response:  Surv(tt, status) 
Input:  x 
Number of observations:  1000 

1) x <= 0; criterion = 1, statistic = 150.161
  2)*  weights = 664 
1) x > 0
  3)*  weights = 336 
> data.df$node <- predict(model, type="node")

確認

> fitdistr(data.df[data.df$node==2,]$tt, "weibull")
      shape          scale    
    1.98999386   202.76429299 
 (  0.06034834) (  4.16490819)
Warning message:
In densfun(x, parm[1], parm[2], ...) : NaNs produced
> fitdistr(data.df[data.df$node==3,]$tt, "weibull")
     shape        scale   
    2.974397   306.660949 
 (  0.127542) (  5.919310)

ワイブル回帰と比例ハザードモデルのサンプルデータ生成と推定

ワイブル回帰

ワイブル分布を仮定して、形状パラメータ共通、尺度パラメータを説明変数で変わるモデル

\begin{align*} F(t|x) = 1 - \exp\left( - \left( \frac{t}{\eta(x)} \right)^m \right) \\ \eta(x) = \exp(\beta * x) \end{align*}

ワイブル回帰のサンプルデータ

set.seed(0)
n <- 1000 # 台数
m <- 2 # 形状パラメータ
beta <- -1 # 係数
beta0 <- 10 # 切片
x <- rnorm(n)
data.df <- data.frame(tt = rweibull(n, shape = m, scale=exp(beta * x + beta0)), x = x, status=rep(1,n))

weibregで推定

> library(eha)
> weibreg(Surv(tt, status)~x, data=data.df)
Call:
weibreg(formula = Surv(tt, status) ~ x, data = data.df)

Covariate           Mean       Coef Exp(Coef)  se(Coef)    Wald p
x                  -0.927     2.017     7.513     0.059     0.000 

log(scale)                   10.016 22386.800     0.016     0.000 
log(shape)                    0.710     2.034     0.025     0.000 

Events                    1000 
Total time at risk        32057517 
Max. log. likelihood      -10619 
LR test statistic         1495 
Degrees of freedom        1 
Overall p-value           0

survregで推定

> library(survival)
> survreg(Surv(tt, status) ~ x, data=data.df, dist="weibull")
Call:
survreg(formula = Surv(tt, status) ~ x, data = data.df, dist = "weibull")

Coefficients:
(Intercept)           x 
 10.0162268  -0.9916226 

Scale= 0.491728 

Loglik(model)= -10618.8   Loglik(intercept only)= -11366.3
    Chisq= 1494.89 on 1 degrees of freedom, p= 0 
n= 1000 

比例ハザードモデル

ベースラインのハザード(微小時間における故障率)は共通で、ハザードが説明変数によって変動するモデル \begin{align*} \lambda(t|x) = \lambda_0(t) exp(\beta*x) \end{align*}

比例ハザードモデルのサンプルデータ

set.seed(0)
n <- 1000 # 台数
x <- rbinom(n, size=1, prob = 1/3)
tt <- runif(n , min=1, max=1000)
beta <- -1
# とりあえずベースラインはワイブル分布にしておいた
# ベースラインはノンパラなので、適当に
baseline <- pweibull(1:500, scale=100, shape=2)
hazard0 <- baseline * exp(beta * 0)
hazard1 <- baseline * exp(beta * 1)
data.df <- data.frame(tt = tt, status=ifelse(runif(tt, min=0, max=1) <= ifelse(x==0,hazard0[tt], hazard1[tt]), 1, 0), x = x)

coxphで推定

> library(survival)
> coxph(Surv(tt, status) ~ x, data=data.df)
Call:
coxph(formula = Surv(tt, status) ~ x, data = data.df)


    coef exp(coef) se(coef)     z       p
x -0.979     0.376     0.15 -6.51 7.7e-11

Likelihood ratio test=51.9  on 1 df, p=5.78e-13  n= 520, number of events= 333 
   (480 observations deleted due to missingness)

第六回「データ解析のための統計モデリング入門」で発表してきました

以下が発表資料です。

第六回「データ解析のための統計モデリング入門」前半 from Atsushi Hayakawa

熟読した方は気づくと思いますが、ヘッセ行列を求める微分をした元の関数の一部を省略して定義しています。微分したら消えるから良いよね?

ガンマ分布のパラメータ推定をする

R言語でガンマ分布のパラメータを推定する方法を2通り示す。

データの生成

> x <- rgamma(1000, shape=2, rate=0.5)

最尤法による推定

> fitdistr(x, "gamma")
     shape         rate   
  1.81070179   0.45850801 
 (0.07471446) (0.02177244)

モーメント法による推定

> med.gam <- mean(x)
> var.gam <- var(x)
> (shape <- med.gam/var.gam)
[1] 0.45973
> (rate <- ((med.gam)^2)/var.gam)
[1] 1.815528

参考

第14回 「はじめてのパターン認識」 読書会で発表しました

その時の資料です。

はじパタ11章 後半 from Atsushi Hayakawa

データサイエンティスト養成読本を書きました

データサイエンティスト養成読本 [ビッグデータ時代のビジネスを支えるデータ分析力が身につく! ] という本を書かせて頂きました。 共に書いた人は、普段からお世話になっている人ばかりで、安心感がありました。修士1年という時期に、素晴らしい経験が出来たと思います。今後も精進を続けて行きますので、よろしくお願いします。

チェビシェフの不等式でMathJaxの練習

MathJax使いやすいですね。

\begin{align*} \sigma^2 & = \int_{\infty}^{\infty} (x-\mu)^2 f(x) dx \\ & = \int_{\infty}^{\mu - k\sigma} (x-\mu)^2 f(x) dx + \int_{\mu-k\sigma}^{\mu+k\sigma} (x-\mu) ^2 f(x) dx + \int_{\mu+k\sigma}^{\infty} (x-\mu)^2 f(x) dx \\ & \geq \int_{-infty}^{\mu-k\sigma} (x-\mu)^2 f(x) dx + \int_{\mu+k\sigma}^{\infty} (x-\mu)^2 f(x) dx \\ & = \int_{|x-\mu| \geq k \sigma} (x-\mu)^2 f(x) dx \\ & \geq \int_{|x-\mu| \geq k\sigma} (x-\mu)^2 f(x) dx \\ & \geq k^2 \sigma^2 \int_{|x-\mu| \geq k \sigma} f(x) dx \\ & = k^2 \sigma^2 P(|X-\mu| \geq k \sigma) \end{align*} よって、 \begin{align*} P(|X - \mu| \geq k \sigma) \leq \frac{1}{k^2} \end{align*} となる。

参考: チェビシェフの定理