2015年12月 – かものはしの分析ブログ

Rのleafletパッケージを使って蒙古タンメン中本をプロットしてみた

Rのleafletパッケージで蒙古タンメン中本の店舗をマッピングしました。
中本の全店舗地図は未だ観察できていないので、中本信者の方には嬉しい情報かもしれません。
元となるデータは中本のホームページから収集し、GoogleMapのAPIで住所から緯度・経度を求めました。

Rに読み込ませるデータはこんな感じになりました。

コードはこんな感じです。如何に簡単かがおわかりいただけると思います。

library(leaflet)
library(htmlwidgets)
setwd("任意")

nakamoto <- read.csv("nakamoto.csv",header=T)

m <- leaflet(nakamoto) %>%
  addTiles() %>%
  addMarkers(lng=~longitude,lat=~latitude) %>%
  addPopups(lng=~longitude,lat=~latitude,popup=~store_name)
m

saveWidget(m, file="m.html")

library(leaflet)

library(htmlwidgets)

setwd("任意")

nakamoto <- read.csv("nakamoto.csv",header=T)

m <- leaflet(nakamoto) %>%

addTiles() %>%

addMarkers(lng=~longitude,lat=~latitude) %>%

addPopups(lng=~longitude,lat=~latitude,popup=~store_name)

saveWidget(m, file="m.html")

参考情報

embed R leaflet map in wordpress

Japan.R 2015に参加して得た情報まとめ

Japan.Rで私にとって使えそうなアイデアや手法だけを厳選して載せたいと思います。（結局、長くなった気がしますが。。）
全てを見たい場合は、この方のブログがよいと思います。
Japan.R 2015に参加してきた

Spark Streaming

Spark Streamingストリーミングを使って、挙動に応じて自動的にメールやSlackに送るとかでサイトの監視に使えるそうな。

jupyter

RでいうところのRStudioのような分析しやすいPython用のGUI
　

rangerパッケージ

RFなどができる。C++で実装されているので速いらしい。
最近のRのランダムフォレストパッケージ -ranger/Rborist-

RcppRollパッケージ

移動平均とかの計算が楽になるらしい。
インターネットの予約のステップの伸び率をモデリングし、ある一定速度だとしたら表示を出し分けるなどの施策に活かしているそうな。

多腕バンディットアルゴリズム

サイト内のバナーなどの出し分けの最適化に使っているらしい。Rでは重くてできないそうな。

leafletパッケージ

これを使えば地理データの可視化が簡単にできる。
→デフォルトの地図がOpenStreetMapなので、サイトに載せる際も好都合。
→「ぐりぐり動かせる・レイヤー構造・多機能・Shinyでアプリケーションも作れる」
→Githubの方がバージョン上
→データテーブルで読み込んで緯度経度を指定する。
　サークルや線の太さも指定でき、ポップアップも立てれる。
　地点間でのラインも引ける。addPolylines
　プロットでエクスポートでWebページにすれば図が作れる。
leafletでMapbox（とかいろいろな地図）を使う

purrrパッケージ

→複数のモデルをリストにして、推定結果を一気に出すことができる。（都度、モデルを作る必要がない。）
{purrr} による非テーブルデータの処理
 purrr 0.1.0を使ってみる(1) map

ggfortifyパッケージ

purrrで実行した分析結果のプロットなどが簡単にできるらしい。
Concepts and Basics of ggfortify

Shiny

分析結果を可視化したりできるツールとしてはカジュアルにできていいかもしれない。
無料版では複数のアクセスに弱いので社内用のツールとしてはいい。

rApacheパッケージ

rとアパッチでwebアプリケーションを作れる。

deployR

APIを作るためのフレームワーク

Gradient Boosting Decision Tree

モデルの性能が良くなるように特徴量の選択ができる手法らしいです。これは進捗が良くなる、イケメンな手法ですね。
今度深掘ってみたいです。
Gradient Boosting Decision Treeでの特徴選択 in R

SparkRパッケージ

Rからスパークの処理を実行できる。
SparkRをEC2上で動かして分散処理してみる

Parquet（パーケ）

ビッグデータで一般的になってきているフォーマットのことらしい。
10GBのデータが300MBになった。→これはすごい！！
DMMゲームのログ解析~Parquetについて~

TinySegmenter

MeCabとかを入れてごにょごにょすることが許されていない方向け？
TinySegmenter:Javascriptだけで書かれたコンパクトな分かち書きソフトウェア

knitrパッケージ

Rの分析結果をhtmlでレポーティングができる。
R でレポートを作成するのに knitr が超便利

R Presentation

Studioでスライド作成できる。
Rとウェブの融合（５）ーslidify、R Presentationー

rvestパッケージ

RでWebスクレイピングができてしまう。（私はPythonでやっていますが。）
rvest でログインしてスクレイピング #rstatsj

pforeachパッケージ

デフォルトで並列計算するようにしたパッケージ
R で超簡単に並列処理を書けるパッケージ pforeach を作った

RstanでGoogleトレンドの推定

rstanを用いて、Googleトレンドデータの予測モデルを推定してみます。

ほとんど岩波データサイエンスのものですが、Googleトレンドのデータを月ごとの季節性を加味した状態空間モデルを用いて予測してみました。

今回の分析では、
・modelのstanコード(stan)
・Rでstanを動かすためのコード(R)
（・可視化のためのコード(R)）←必須ではない
を用意します。

データですが、GoogleTrendのサイトで任意のキーワードで検索して、
その時系列データをCSVでダウンロードすれば手に入ります。（ちょっと見つけにくい）

データの形式はシンプルで、
先頭にY
とおいて後はトレンドの値を行ごとに置いていけばいけます。

要はN行1列データをテキストファイルに保存すればOKです。（1行目はY）

まずstanのコードですが、岩波データサイエンスのサンプルコードの季節を4から12に変えています。（たったこれだけ）
Googleトレンドのデータは月単位でも結構値がふれることがあるので、月ごとに応じた潜在的な変数が必要だと思いました。

data {
   int<lower=1> T;
   int<lower=1> T_next;
   real Y[T];
}

parameters {
   real mu[T];
   real s[T];
   real<lower=0> s_mu;
   real<lower=0> s_s;
   real<lower=0> s_r;
}

model {
   for(t in 2:T)
      mu[t] ~ normal(mu[t-1], s_mu);
   for(t in 12:T)
      s[t] ~ normal(-s[t-1]-s[t-2]-s[t-3]-s[t-4]-s[t-5]-s[t-6]-s[t-7]-s[t-8]-s[t-9]-s[t-10]-s[t-11], s_s);
   for(t in 1:T)
      Y[t] ~ normal(mu[t]+s[t], s_r);
}

generated quantities {
   real mu_all[T+T_next];
   real s_all[T+T_next];
   real y_next[T_next];

   for (t in 1:T){
      mu_all[t] <- mu[t];
      s_all[t] <- s[t];
   }
   for (t in (T+1):(T+T_next)){
      mu_all[t] <- normal_rng(mu_all[t-1], s_mu);
      s_all[t] <- normal_rng(-s_all[t-1]-s_all[t-2]-s_all[t-3]-s_all[t-4]-s_all[t-5]-s_all[t-6]-s_all[t-7]-s_all[t-8]-s_all[t-9]-s_all[t-10]-s_all[t-11], s_s);
   }
   for (t in 1:T_next)
      y_next[t] <- normal_rng(mu_all[T+t]+s_all[T+t], s_r);
}

data {

int<lower=1> T;

int<lower=1> T_next;

real Y[T];

}

parameters {

real mu[T];

real s[T];

real<lower=0> s_mu;

real<lower=0> s_s;

real<lower=0> s_r;

}

model {

for(t in 2:T)

mu[t] ~ normal(mu[t-1], s_mu);

for(t in 12:T)

s[t] ~ normal(-s[t-1]-s[t-2]-s[t-3]-s[t-4]-s[t-5]-s[t-6]-s[t-7]-s[t-8]-s[t-9]-s[t-10]-s[t-11], s_s);

for(t in 1:T)

Y[t] ~ normal(mu[t]+s[t], s_r);

}

generated quantities {

real mu_all[T+T_next];

real s_all[T+T_next];

real y_next[T_next];

for (t in 1:T){

mu_all[t] <- mu[t];

s_all[t] <- s[t];

}

for (t in (T+1):(T+T_next)){

mu_all[t] <- normal_rng(mu_all[t-1], s_mu);

s_all[t] <- normal_rng(-s_all[t-1]-s_all[t-2]-s_all[t-3]-s_all[t-4]-s_all[t-5]-s_all[t-6]-s_all[t-7]-s_all[t-8]-s_all[t-9]-s_all[t-10]-s_all[t-11], s_s);

}

for (t in 1:T_next)

y_next[t] <- normal_rng(mu_all[T+t]+s_all[T+t], s_r);

}

Rでstanを動かすためのコードですが、ここはサンプルコードとほぼ一緒です。

library(rstan)
setwd("任意")

d <- read.csv('data-trend4.txt', header=TRUE)
T <- nrow(d)
T_next <- 8
data <- list(T=T, T_next=T_next, Y=d$Y)

stanmodel <- stan_model(file='model_custom.stan')
fit <- sampling(
   stanmodel, data=data, pars=c('mu_all','s_all','y_next','s_mu','s_s','s_r'),
   iter=10200, warmup=200, thin=10, chains=3,
   seed=123
)

library(rstan)

setwd("任意")

d <- read.csv('data-trend4.txt', header=TRUE)

T <- nrow(d)

T_next <- 8

data <- list(T=T, T_next=T_next, Y=d$Y)

stanmodel <- stan_model(file='model_custom.stan')

fit <- sampling(

stanmodel, data=data, pars=c('mu_all','s_all','y_next','s_mu','s_s','s_r'),

iter=10200, warmup=200, thin=10, chains=3,

seed=123

)

可視化のためのコードについてもサンプルコードとほぼ一緒です。

library(ggplot2)

# after estimation

d_obs <- data.frame(X=1:T, Y=d$Y)
p <- ggplot()
p <- p + theme_bw() + theme(text=element_text(size=18))
p <- p + geom_line(data=d_obs, aes(x=X, y=Y), color='black', alpha=0.8, size=2)
p <- p + labs(x='Time [month]', y='trend')
p <- p + coord_cartesian(xlim=c(0.9, 152.1))
ggsave(file='fig2-top-left.png', plot=p, dpi=300, width=6, height=4)


makeDataFrameQuantile <- function(x, y_smp){
  qua <- apply(y_smp, 2, quantile, prob=c(0.1, 0.25, 0.5, 0.75, 0.9))
  d_est <- data.frame(X=x, t(qua))
  colnames(d_est) <- c('X', 'p10', 'p25', 'p50', 'p75', 'p90')
  return(d_est)
}

plotTimecourse <- function(file, d_est, d_obs){
  p <- ggplot()
  p <- p + theme_bw() + theme(text=element_text(size=18))
  p <- p + geom_vline(xintercept=T, linetype='dashed')
  p <- p + geom_ribbon(data=d_est, aes(x=X, ymin=p10, ymax=p90), fill='black', alpha=0.25)
  p <- p + geom_ribbon(data=d_est, aes(x=X, ymin=p25, ymax=p75), fill='black', alpha=0.5)
  p <- p + geom_line(data=d_est, aes(x=X, y=p10), color='black', size=0.2)
  p <- p + geom_line(data=d_est, aes(x=X, y=p90), color='black', size=0.2)
  p <- p + geom_line(data=d_est, aes(x=X, y=p25), color='black', size=0.2)
  p <- p + geom_line(data=d_est, aes(x=X, y=p75), color='black', size=0.2)
  p <- p + geom_line(data=d_est, aes(x=X, y=p50), color='black', size=0.4)
  if (!is.null(d_obs)){
    p <- p + geom_line(data=d_obs, aes(x=X, y=Y), color='black', size=2, alpha=0.9)
  }
  p <- p + labs(x='month', y='trend')
  p <- p + coord_cartesian(xlim=c(0.9, 152.1))
  ggsave(file=file, plot=p, dpi=300, width=6, height=4)
}

la <- rstan::extract(fit)
d_est <- makeDataFrameQuantile(x=1:(T+T_next), y_smp=la$mu_all)
plotTimecourse(file='fig2-bottom-left.png', d_est=d_est, d_obs=d_obs)

d_est <- makeDataFrameQuantile(x=1:(T+T_next), y_smp=la$s_all)
plotTimecourse(file='fig2-bottom-right.png', d_est=d_est, d_obs=NULL)

d_est <- makeDataFrameQuantile(x=(T+1):(T+T_next), y_smp=la$y_next)
d_est <- rbind(data.frame(X=T, p10=d$Y[T], p25=d$Y[T], p50=d$Y[T], p75=d$Y[T], p90=d$Y[T]), d_est)
plotTimecourse(file='fig2-top-right.png', d_est=d_est, d_obs=d_obs)

library(ggplot2)

# after estimation

d_obs <- data.frame(X=1:T, Y=d$Y)

p <- ggplot()

p <- p + theme_bw() + theme(text=element_text(size=18))

p <- p + geom_line(data=d_obs, aes(x=X, y=Y), color='black', alpha=0.8, size=2)

p <- p + labs(x='Time [month]', y='trend')

p <- p + coord_cartesian(xlim=c(0.9, 152.1))

ggsave(file='fig2-top-left.png', plot=p, dpi=300, width=6, height=4)

makeDataFrameQuantile <- function(x, y_smp){

qua <- apply(y_smp, 2, quantile, prob=c(0.1, 0.25, 0.5, 0.75, 0.9))

d_est <- data.frame(X=x, t(qua))

colnames(d_est) <- c('X', 'p10', 'p25', 'p50', 'p75', 'p90')

return(d_est)

}

plotTimecourse <- function(file, d_est, d_obs){

p <- ggplot()

p <- p + theme_bw() + theme(text=element_text(size=18))

p <- p + geom_vline(xintercept=T, linetype='dashed')

p <- p + geom_ribbon(data=d_est, aes(x=X, ymin=p10, ymax=p90), fill='black', alpha=0.25)

p <- p + geom_ribbon(data=d_est, aes(x=X, ymin=p25, ymax=p75), fill='black', alpha=0.5)

p <- p + geom_line(data=d_est, aes(x=X, y=p10), color='black', size=0.2)

p <- p + geom_line(data=d_est, aes(x=X, y=p90), color='black', size=0.2)

p <- p + geom_line(data=d_est, aes(x=X, y=p25), color='black', size=0.2)

p <- p + geom_line(data=d_est, aes(x=X, y=p75), color='black', size=0.2)

p <- p + geom_line(data=d_est, aes(x=X, y=p50), color='black', size=0.4)

if (!is.null(d_obs)){

p <- p + geom_line(data=d_obs, aes(x=X, y=Y), color='black', size=2, alpha=0.9)

}

p <- p + labs(x='month', y='trend')

p <- p + coord_cartesian(xlim=c(0.9, 152.1))

ggsave(file=file, plot=p, dpi=300, width=6, height=4)

}

la <- rstan::extract(fit)

d_est <- makeDataFrameQuantile(x=1:(T+T_next), y_smp=la$mu_all)

plotTimecourse(file='fig2-bottom-left.png', d_est=d_est, d_obs=d_obs)

d_est <- makeDataFrameQuantile(x=1:(T+T_next), y_smp=la$s_all)

plotTimecourse(file='fig2-bottom-right.png', d_est=d_est, d_obs=NULL)

d_est <- makeDataFrameQuantile(x=(T+1):(T+T_next), y_smp=la$y_next)

d_est <- rbind(data.frame(X=T, p10=d$Y[T], p25=d$Y[T], p50=d$Y[T], p75=d$Y[T], p90=d$Y[T]), d_est)

plotTimecourse(file='fig2-top-right.png', d_est=d_est, d_obs=d_obs)

以上を実行した結果、以下のような図が出てきます。

こちらは実際の時系列データのプロットです。

8期先までの予測です。

8期先までの予測範囲です。信頼区間90%までの範囲となっています。

推定した潜在的な季節性のデータをプロットしています。

ついでに、4月までのデータを用いて、5~9月の予測を行い、その比較を行っています。

5月が大きく外れましたが、その後はある程度当てれているように見えます。
5月も当てれるようなモデルを作りたいものですね。

参考文献

<br />