April 18, 2015

時系列解析入門 1章

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

March 28, 2015

python3による日付の扱い方メモ

twitterのAPIによるとタイムスタンプ Wed Dec 24 11:34:28 +0000 2014 python3で処理する from datetime import datetime, timezone, timedelta time = "Wed Dec 24 11:34:28 +0000 2014" d = datetime.strptime(time, '%a %b %d %H:%M:%S %z %Y') とすると、dは datetime.datetime(2014, 12, 24, 20, 34, 28, tzinfo=datetime.timezone(datetime.timedelta(0, 32400), 'JST')) となる。 unixtimeを得る d.timestamp() とすればよく、結果は 1419420868.0 となる。 日本時間にする JST = timezone(timedelta(hours=+9), 'JST') d = datetime.fromtimestamp(d.timestamp(), JST) とする。 文字列として出力 d.strftime("%Y/%m/%d %H:%M:%S") で出力は '2014/12/24 20:34:28' となる。 Read more

February 24, 2015

dplyrの暗黒面

> library(dplyr) > iris %>% + summarise_(heikin = "mean(Sepal.Length)") heikin 1 5.843333 と実行できるが、 > myfunc <- function(x){ + mean(x) + } > iris %>% + summarise_(heikin = "myfunc(Sepal.Length)") Error in summarise_impl(.data, dots) : could not find function "myfunc" は実行出来ない。 > iris %>% + summarise_(heikin = as.formula("~myfunc(Sepal.Length)")) heikin 1 5.843333 とすれば、実行できた。

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 31, 2014

2014年度のアクセス数

2014年がもうすぐ終わります。このブログを振り返るためにアクセスを載せます。 12月は記事の投稿が多くなったおかげで、多くの方にアクセスして頂きました。 流入元では、約40%がgoogleからの検索で、約25%がtwitterからでした。 一番人気の記事はhttp://blog.gepuro.net/archives/72で、検索エンジンからのアクセスが多いです。統計関連の記事でないのが寂しいですが、この調子で技術ネタを投稿していきたいです。また、最近増えている雑記カテゴリな記事もネタを思い付き次第書いてきます。 2015年もgepulogをよろしくお願いします。

December 23, 2014

新卒社会人のお金の計算をしてみる

ひとり暮らしの生活にかかるお金を計算してみた。 初任給 大卒の初任給は19万8000円!まさかこれより低賃金な奴はいないよなwwwwwwwと煽っている人がいますが、初任給の手取りを20万円として計算してみる。 食費 ひとり暮らしの生活で変動の自由度が最も高い食費を見る。 単身世帯では男性は女性より月1.2万円食費が多いによれば、単身男性の食費の割合は28.3%になっている。初任給20万円 * 0.283 = 56000円になる。飲み会や外食が多いのだろうか。食費の計算をしていた時は3万〜4万程度になることが多かった。社会人になるとランチが1000円を超える事も多いと思うので、朝100円 + ランチ1000円 + 夜400円=1500円/日 として、30日で4万5000円。土日はランチが安くなったり、週一で飲み会に参加すれば、5万6000円程度になるのだろうか。 住居費 その年収なら、どんな家に住める?(1)によれば住居費が収入の25%を超えると家計を圧迫するとの事。理想的な20%とすると家賃は20万円 * 0.2 = 4万円となる。都内ひとり暮らしで4万の家賃に住むならば、郊外に住む必要がある。感覚値であるが、都内で比較的に家賃が安い調布であっても5万円程度は必要になる。ギリギリである25%ならば家賃が5万円となる。実際にはもっと高い人が多いだろう。 初任給 - 食費 - 住居費 ここまでで残ったお金を計算すると、 20万円 - 5万6000円 - 5万 = 9万4000円になる。 生活費を残りのお金でやりくりする必要がある。 水道代や光熱費 6年間生活してる中での感覚値では、合計1万円程度になる。 参考に、家賃、食費、光熱費…気になる一人暮らしの出費平均を見ても、同じ程度である。 通信費 携帯料金については、 https://www.nttdocomo.co.jp/iphone/charge/index.htmlにある5GBプランならば、7200円になる。MVNOなどを駆使すれば、もう少し抑えられるだろう。 コンピュータのインターネット回線は、5000円ぐらいだろうか。合計すると1万2000円になる。 残りのお金を再計算する 9万4000円 - 1万円 - 1万2000円 = 7万2000円 奨学金の返済 苦学生にとっては忘れてはいけないお金である。 http://www.jasso.go.jp/henkan/henkanrei/daigaku/を見れば、国立で自宅外の学生の返還月額は14,400円となっている。大学院に進学すれば、もう少し増えるだろう。 年金など 国民年金保険料って、いくら?を見れば、1万5250円。 住民税の税率で確認して、10%。20万 * 0.1 / 12 = 約8000円です。住民税は一年目はかからないですが、一年で昇給すれば嬉しいですね。(汗 都道府県別 国民健康保険料ランキングを見ると東京の1人当たりの額は8万2936円になる。月額にすると、約7000円になる。 残ったお金は? 7万2000円 - 1万4000円 - 1万5000円 - 8000円 - 7000円 = 2万8000円 Read more

December 22, 2014

回帰木+ヒストグラムの図を描く

回帰木の下にヒストグラムがあると格好良いと思ったので、メモしておく。 library("party") data(cars) model <- ctree(dist~., data=cars) t.style <- node_hist(model, ymax=0.06, xscale=c(0,150), col="red", fill=hsv(0.6, 0.5, 1)) plot(model, terminal_panel = t.style) 参考 [Rによるデータサイエンス13「樹木モデル」 ](http://www.slideshare.net/takemikami/r13-9821987)

December 22, 2014

集計時の順序を曜日で並べる

weekdays()を使って、Date型から曜日を取得できる。 > (today <- Sys.Date()) [1] "2014-12-22" > str(weekdays(today)) chr "月曜日" 戻り値が文字列になってるので、集計時に曜日順にならない。 > x <- c("月曜日","日曜日", "火曜日", "日曜日", "土曜日") > table(x) x 火曜日 月曜日 土曜日 日曜日 1 1 1 2 factor型にしてlevelsを指定することによって、順番を制御できる。 > yobi <- factor(x, levels=paste(c("日", "月", "火", "水", "木", "金", "土"), "曜日", sep="")) > table(yobi) yobi 日曜日 月曜日 火曜日 水曜日 木曜日 金曜日 土曜日 2 1 1 0 0 0 1 嬉しい! twitterで教えてもらった書き方は、曜日順での集計に関するメモ書き です。

December 19, 2014

Japan.R 2014を支える技術

R Advent Calendar 2014の18日目担当です。 12月6日にJapan.R 2014が開催されました。 総勢250名を超える方に参加して頂いて非常に嬉しいです。@0kayuさんを始め、会場提供およびスポンサーになって頂いたフリークアウトさん、 当日に手伝いをしてくれた方々、ありがとうございます! Japan.R 2014を開催するまでの話を書きます。 主催の一人に選ばれる ヤフー株式会社で開催されたJapan.R 2013での二次会後に、 主催の@holidayworkingさんに無茶振りされました。その時にいた@0kayuさんも巻き添え?になってしまいました。 2010年のJapan.Rに参加して以来、頻繁にTokyo.Rに参加していてRコミュニティに貢献したかったので、受けることにしました。 開催準備 6月頃に、12月6日にフリークアウトで開催することを決定 7月に参加して勉強会の懇親会で、パネルディスカッションのアイディアを頂く パネラー探しを始める。最終的に決まったのは9月下旬 30分枠の発表をお願い出来た 7月に オプトの@shsaixさん フリークアウトの@yanaokiさん 後に RCOの@TJO_datasciさん @teramonagiさんや@dichikaさんの紹介で@AriLamsteinさん 8月のTokyo.RでLTの募集を開始する(google docsのフォームを利用) 開催内容を紹介する。http://www.slideshare.net/yurieoka37/japanr 10月にフリークアウトさんがスポンサーになってくれることが決まった。ありがとうございます! 10月下旬にATNDをたてる Bar doradoraのお願いをする 11月上旬に@0kayuさんとピザを食いながら、パネルディスカッションの内容を打ち合わせ。ピザ生地が旨かった。 11月下旬に某社の会議室を貸して頂いて、パネラーの皆さんに集まってもらってパネルディスカッションの内容を打ち合わせ 12月始め レガラートで食べ物を注文 パネルディスカッション用のマイクをレンタル 二次会のエイトを予約 当日 11:15に当日スタッフの集合 会場の準備 Read more

December 15, 2014

オタクは未婚が多いのか?

2013年度のアイドル市場は前年度比約2割増/オタクの約7割は未婚者【オタク市場に関する調査】という見出しに衝撃を受けたので、オタク?の1人として、未婚率について少し調べてみた。統計的なツッコミを歓迎します。 オタクの年齢層 オタクの半数は関東在住、今年最も支持されたのは「東方Project」 SMS「おたくま」が発表の記事では、 という調査結果があり、21才〜25才が最も高い割合になっている。 年代別の未婚率 年齢階級別既婚・未婚の構成比に掲載されている国勢調査のまとめには、 と書かれており、20才〜24才では約90%が未婚で25才〜29才では約60%が未婚となっている。 まとめ オタクのボリュームゾーンの年代と、国勢調査の未婚率を並べて見てみると、オタク⇛未婚率が高いとは言えなさそうだ。 こんな事を気にしてるからモテないんだと言われると凹むので、言う時は優しめにお願いします。

December 14, 2014

Rstanでブログが検索結果に出ない話

で表示されないけど、 でなら表示される。 Rstan単独で検索してブログが表示される人もいるので、検索エンジンに嫌われてる可能性が高くて涙目です。 ちなみに、「stan ワイブル」で検索すれば、上位にヒットしました。

December 11, 2014

xargsまたはparalellを用いてRを並列処理させる

list.txtから入力を受け取り、hoge.Rで並列処理をさせる。 hoge.R argv <- commandArgs(TRUE) argv <- strsplit(argv, ",")[[1]] par1 <- as.numeric(argv[1]) par2 <- as.numeric(argv[2]) print(par1) print(par2) list.txt 1,1 1,2 2,1 2,2 xargsを用いる場合 cat list.txt | xargs -P8 -n1 Rscript hoge.R paralellを用いる場合 cat list.txt | parallel -j 8 Rscript hoge.R {}

December 10, 2014

Rstanで重回帰分析

ワイブルばかりだったので、メジャーな重回帰分析をしてみた。Rstanにハマりつつある。 linear_regression.stan data { int<lower=1> J; // number of data int<lower=1> K; // number of covariate real y[J]; matrix[J,K] x; } parameters { real<lower=0> sd0; vector[K] beta; real beta0; } model { for(j in 1:J){ increment_log_prob(normal_log(y[j], x[j] * beta + beta0, sd0)); } } linear_regression.R library(rstan) n <- 100 x1 <- rnorm(n, mean=1, sd=0.1) x2 <- rnorm(n, mean=1, sd=0.1) beta1 <- 1 beta2 <- 2 beta0 <- 3 y <- beta1 * x1 + beta2 * x2 + beta0 data. Read more

December 4, 2014

Rstanでワイブル回帰の推定(説明変数複数バージョン)

matrixやvectorの使い方を覚えた。 weibull_regression.stan data { int<lower=1> J; // number of data int<lower=1> K; // number of covariate int status[J]; real<lower=0> tt[J]; matrix[J,K] x; } parameters { real<lower=1, upper=3> m; vector[K] beta; real beta0; } transformed parameters { real<lower=0> eta[J]; for (j in 1:J) eta[j] <- exp(x[j]*beta + beta0); } model { for(j in 1:J){ if(status[j] == 1){ increment_log_prob(weibull_log(tt[j], m, eta[j])); } if(status[j] == 0){ increment_log_prob(weibull_ccdf_log(tt[j], m, eta[j])); } } } weibull_regression. Read more

December 4, 2014

Rstanでワイブル回帰の推定

weibull_regression.stan data { int<lower=1> J; // number of data int status[J]; real<lower=0> tt[J]; real x[J]; } parameters { real<lower=1, upper=3> m; real beta; real beta0; } transformed parameters { real<lower=0> eta[J]; for (j in 1:J) eta[j] <- exp(beta * x[j] + beta0); } model { for(j in 1:J){ if(status[j] == 1){ increment_log_prob(weibull_log(tt[j], m, eta[j])); } if(status[j] == 0){ increment_log_prob(weibull_ccdf_log(tt[j], m, eta[j])); } } } mを1から3ぐらいにしておかないと、推定値が遥か彼方に飛ばされることがある。ワイブル分布は厄介やね。 weibull_regression.R library(rstan) set.seed(0) n <- 1000 # 台数 m <- 2 # 形状パラメータ beta <- -1 # 係数 beta0 <- 10 # 切片 x <- rnorm(n) data. Read more

December 4, 2014

Rstanでワイブル分布のパラメータ推定(打切りのあるデータで)

weibull.stan data { int<lower=1> J; // number of data int status[J]; real<lower=0> tt[J]; } parameters { real<lower=1> m; real<lower=1> eta; } model { for(j in 1:J){ if(status[j] == 1){ increment_log_prob(weibull_log(tt[j], m, eta)); } if(status[j] == 0){ increment_log_prob(weibull_ccdf_log(tt[j], m, eta)); } } } weibull.R サンプルデータの生成 > set.seed(0) > data.df <- data.frame(tt = rweibull(1000, shape = 2, scale = 100), + status = rep(1,1000)) > data.df[data.df$tt > 30,]$status <- 0 > data. Read more

December 3, 2014

Rstanでワイブル分布のパラメータ推定

Rstanでワイブル分布のパラメータ推定をしてみる。打切りを考慮するにはどうしたら良いのだろうか? Stan Modeling Language User’s Guide and Reference Manual, v2.5.0の83ページ辺りを読めば出来そう。 weibull.stan data { int<lower=0> J; // number of data real<lower=0> tt[J]; // time } parameters { real<lower=0> m; real<lower=0> eta; } model { for(j in 1:J) tt[j] ~ weibull(m, eta); } weibull.R library(rstan) data.ls <- list(J = 1000, tt = rweibull(1000, shape = 2, scale = 100)) fit <- stan(file = 'weibull.stan', data = data.ls, iter = 1000, chains = 4) 結果 > fit Inference for Stan model: weibull. 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

© gepuro 2013