May 2, 2015

確率変数の変換

何度か勉強しているのだけれど、毎回確認しているのでブログに書いておく。 数式転がし 確率変数 $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} となる。 Read more

April 29, 2015

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

オハイオ州立大学のRobert de Jong先生が公開されている講義資料?がよさ気なので、紹介します。 Lecture 1: Stationary Time Series∗ Lecture 2: ARMA Models∗ Lecture 3: Spectral Analysis∗ Lecture 4: Asymptotic Distribution Theory∗ Lecture 5: Linear Regressions∗ Lecture 6: Vector Autoregression∗ Lecture 7: Processes with Deterministic Trends∗ Lecture 8: Univariate Processes with Unit Roots∗ [Lecture 9: Multivariate Unit Root Processes and Cointegration](http://www.econ.ohio-state.edu/dejong/note9.pdf) Lecture 10: Further Topics∗ 誤植 部分的にしか読んでいないけど、誤植がありました。他にも誤植がありそう。 Read more

April 18, 2015

時系列解析入門 1章

北川先生の時系列解析入門の読書記録です。 1.1 時系列データ 時間の経過とともに不規則に変動する現象の記録が時系列 1.2 時系列の分類 連続時間時系列と離散時間時系列 連続時間時系列 例:アナログレコーダ 離散時間時系列 例:一時間おきに計測された気圧 一変量時系列と多変量時系列 一変量 観測時点で一種類の情報 多変量 2つ以上の情報 定常時系列と非定常時系列 時間的に変化しない一定の確率的モデルの実現値とみなせる 定常時系列 平均が時間とともに変動、平均のまわりの変動の仕方が時間的に変化 非定常時系列 ガウス型時系列と非ガウス型時系列 時系列の分布が正規分布に従う ガウス型時系列 時系列の分布が正規分布に従わない 非ガウス型時系列 非ガウス型時系列のように上下非対称性が見られても、データを変換することによって、近似的にガウス時系列とみなせる場合がある 線形時系列と非線型時系列 線形なモデルの出力として表現できるような時系列 線形時系列 非線形なモデルが必要な時系列 非線型時系列 欠測値と異常値 何らかの理由で観測値が記録されない 欠測値 現象の異常な振る舞い・観測装置の異常など Read more

February 22, 2015

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

第46回R勉強会@東京(#TokyoR)でLTをしてきました。 資料は、こちらです。 最強のハードディスクはどれだ? from Atsushi Hayakawa また、https://github.com/gepuro/gweibullplotでワイブルプロットのパッケージを公開しました。本資料を作成する際に用いた関数とは別のものです。(研究室秘伝のソースでライセンスが謎なため)

February 7, 2015

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

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. Read more

January 26, 2015

指数分布の確率プロット

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. Read more

December 3, 2014

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

条件付き推測木(決定木)で生存時間解析 寿命分布が異なるように説明変数で分割する。 詳細は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. Read more

December 2, 2014

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

ワイブル回帰 ワイブル分布を仮定して、形状パラメータ共通、尺度パラメータを説明変数で変わるモデル \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. Read more

July 22, 2014

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

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 参考 http://cran.r-project.org/doc/contrib/Ricci-distributions-en.pdf http://lokad.jp/actuary/wp-content/uploads/2010/09/Gamma.pdf http://ja.wikipedia.org/wiki/%E3%82%AC%E3%83%B3%E3%83%9E%E5%88%86%E5%B8%83

August 8, 2013

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

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

August 7, 2013

チェビシェフの不等式で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 \\ Read more

© gepuro 2013