Project Euler Problem 1:10

久しぶりにProject Eulerをやってみようと思ったが、せっかくなので進捗をリセットして最初からにした。
とりあえず10問やってみたが大分脳がヤバくなってる感あったのでちゃんと継続してボケ防止に努めたい。
解答コードはgithubで管理することにして、ここには要点をメモっていく。

Problem 1: Multiples of 3 and 5

簡易版FizzBuzz。添字操作と条件式の組み合わせの基本確認。

Problem 2: Even Fibonacci numbers

RでProject Eulerをやる場合、gmpパッケージを縛るか縛らないかで難易度に雲泥の差が生ずるが、基本的には使っていく方針で進めていきたい。便利そうなパッケージがあれば積極的に利用していく。
フィボナッチ数はgmp::fibnum()で求められる。

Problem 3: Largest prime factor

素因数分解はgmp::factorize()で実行できる。

Problem 4: Largest palindrome product

文字列の反転はstringi::stri_reverse()で実行できる。
組み合わせ毎に調べたほうが効率が良いが、探索範囲が狭いので999*999からスタートして下っていく方針でもすぐに終わる。
クリアすると読めるPDFファイルで勉強するための問題だろう。

Problem 5: Smallest multiple

大きな階乗はgmp::factorialZ()を使うと良い。
Big Integerに演算を施す場合は専用の関数が必要になる場合がしばしばある。剰余を求めるにはgmp::mod.bigz()を使う。

Problem 6: Sum square difference

ベクトル演算が出来るRでは基本操作レベル。

Problem 7: 10001st prime

gmpにはgmp::nextprime()という便利な関数がある。

Problem 8: Largest product in a series

文字列の操作方法が問題となる。

  • 置換: sub(最初の1つだけ)、gsub(全て)
  • 切り出し: substring
  • 分割: strsplit(結果はlistなのでunlistする)

Problem 9: Special Pythagorean triplet

多重ループを抜けるにはlabeledLoopパッケージによるラベル付きループが便利だ。
探索範囲が狭いのでbrute forceでもすぐに終わってしまうのだが、これもProblem 4と同様、クリア後にPDFファイルで学習する用に簡単になっているのだろう。

Problem 10: Summation of primes

エラトステネスの篩は覚えるまで何度か自分で実装してみたほうが良いが、numbersというパッケージにエラトステネスの篩が入っている。numbers::Primes()は2つの引数を与えるとその数以下の、2つの値を与えると2数間の素数を列挙する。2つの値を与える場合に2数が20万くらい以上離れているとエラーが出るようだが、一つだけなら2^53-1まで指定できるようだ(あまり大きいと終わらなくなるが)

「データ解析のための統計モデリング入門」第11章メモ

p.254の図11.7をstanを使って再現したかった。

stanコードはこんな感じ。

data {
   int<lower=1> N;     // 欠測値を含むデータ数
   int<lower=1> N_obs; // 欠測値以外のデータの数
   int<lower=1, upper=N> Idx_obs[N_obs]; // 欠測値以外のデータのインデックス
   int<lower=0> Y_obs[N_obs]; // 欠測値を除いたデータベクトル
}

parameters {
   real r[N]; 
   real<lower=0> s_r; 
}

model {
   for (j in 2:N)
      r[j] ~ normal(r[j-1], s_r); // 1次のCAR model
   for (j in 1:N_obs) 
      Y_obs[j] ~ poisson_log(r[Idx_obs[j]]);
}

generated quantities { 
   real<lower=0> Y_mean[N];
   for (j in 1:N)
      Y_mean[j] = exp(r[j]);
}

R側のコードはこんな感じ。

# データを欠測させる
no <- c(6, 9, 12, 13, 26:30)
Idx_obs <- (1:50)[-no]
Y_obs <- Y[Idx_obs]
N <- length(Y)
N_obs <- length(Y_obs)

# サンプリング
dlist <- list(N = N, N_obs = N_obs, Idx_obs = Idx_obs, Y_obs = Y_obs)
fit2 <- stan("sec11_2.stan", data = dlist)

# 結果から中央値と信頼区間を算出
pred <- extract(fit2)
ci <- apply(pred$Y_mean, 2, quantile, probs = c(0.5, 0.025, 0.975))

# プロット
plot(Y, type = "n", ylim = c(0, 25))
polygon(c(1:N, N:1), 
        c(ci[2, ], rev(ci[3, ])),
        border = NA,
        col = "gray90")
points(Y)
points(no, Y[no], pch = 19)
lines(m, lty = 2)
lines(ci[1, ])
for(n in no){
  polygon(c(n-0.5, n+0.5, n+0.5, n-0.5),
          c(0, 0, 100, 100), 
          border = NA, 
          col = rgb(.3, .3, .3, .3))
}

するとこのようなものが出来る。
f:id:Rion778:20161008220158p:plain
空間相関を考慮しないバージョンはこんな感じ。

data {
   int<lower=1> N;     // 欠測値を含むデータ数
   int<lower=1> N_obs; // 欠測値以外のデータの数
   int<lower=1, upper=N> Idx_obs[N_obs]; // 欠測値以外のデータのインデックス
   int<lower=0> Y_obs[N_obs]; // 欠測値を除いたデータベクトル
}

parameters {
   real r[N]; 
   real<lower=0> s_r; 
   real mu;
}

model {
   for (j in 1:N)
      r[j] ~ normal(0, s_r);
   for (j in 1:N_obs) 
      Y_obs[j] ~ poisson_log(mu + r[Idx_obs[j]]);
}

generated quantities {
   real<lower=0> Y_mean[N];
   for (j in 1:N)
      Y_mean[j] = exp(mu + r[j]);
}

f:id:Rion778:20161008220315p:plain

参考

stanコードの記述にあたっては下記サイトが参考になったというかほとんどそのまんまです。

「データ解析のための統計モデリング入門」第9章メモつづき

前回。
rion778.hatenablog.com
これをrstanarmでやってみた。

準備

install.packages("rstanarm")

割と時間かかった。

library(rstan)
library(rstanarm)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
# データダウンロード
download.file(destfile = "ch9.Rdata", "http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/gibbs/d.RData")
load("ch9.Rdata")

実行

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

細かい設定は省略したけど1行で済むのは楽。あと早い。

プロット

rstanに入っている作図関数がそのまま使える。

plot(fit)
stan_dens(fit, separate_chains = TRUE) + geom_rug()
stan_trace(fit)

出力は特に違いがないので略。
stan_traceにinclude_warmup=TRUEを指定した時の挙動が若干おかしい気はする(thinが無視されている様子?)。

出力

デフォルトの出力はかなり要約されている。summaryで少し詳しい物が見られるが、$stanfitの中にrstanの時と同じ形式のやつが入っていたので、例えばextract()を使いたい時なんかはこれを使えば良い様だ。

fit                # デフォルト
summary(fit)       # 少し詳しく
fit$stanfit        # rstanと同じで

テキストの図を再現する

上述のように、$stanfitをextract()してやれば良い。変数名が何故かalphaとbetaだ。

post <- extract(fit$stanfit)
add.mean <- function(bb1, bb2, lty = 2, lwd = 1, ...){
  lines(d$x, exp(bb1 + bb2 * d$x), lty = lty, lwd = lwd, ...)
}
plot(d$x, d$y, type = "n", xlab = "x_i", ylab = "y_i")
for (i in 1:length(post$alpha)) {
  add.mean(post$alpha[i], post$beta[i], lty = 1, col = "#00000004")
}
points(d$x, d$y)
add.mean(median(post$alpha), median(post$beta), lty = 1, lwd = 2)

これも出力に差は無いので略。

geom_barにおけるNAの扱い

geom_barというかstat_countの方だと思うけど、例えばこんなデータがあったとして

data <- c(1, 2, 2, 3, 3, 3, NA, NA, NA, NA)

普通に棒グラフを描くとこうだ。

ggplot() + geom_bar(aes(data))

f:id:Rion778:20160825202906p:plain
ところがfactorにするとこうだ。

ggplot() + geom_bar(aes(factor(data)))

f:id:Rion778:20160825202941p:plain
どうも集計にnaを含むかどうかはfactorであるかどうかに依存しているらしい(正しい動作かどうかはよく分からない)。ちなみにgeom_barとstat_countにはna.rm引数が用意されているが、これは警告の有無を制御するだけでNAが除外されることには変わりがない。
これはLearnBayesパッケージに含まれるstudentdataのDrinkをプロットしようとして気付いたのだが、もしこのようなデータをstat_countに渡す必要が出てきたら、factorで無くするか、NAを除外する必要がある。逆にNAを含めたかったらfactorにしてやる。

install.packages("LearnBayes")
library(LearnBayes)
data(studentdata)
## NAを除外する方法
ggplot(subset(studentdata, !is.na(Drink)), aes(x = Drink)) + geom_bar()
## characterに変換する方法
ggplot(studentdata, aes(x = as.character(Drink))) + geom_bar()

「データ解析のための統計モデリング入門」第9章メモ

緑本9章の内容について、RStanでやってみる。
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

Stanコード

次の内容で作成し、sample.stanという名前で保存した。

data {
  int<lower=1> N;
  vector[N] x;
  int<lower=0> y[N];
  real mean_x;
}
parameters {
  real beta1;
  real beta2;
}
model {
  for(i in 1:N){
    y[i] ~ poisson(exp(beta1 + beta2 * (x[i] - mean_x)));
  }
}

コードはPyStanによるはじめてのマルコフ連鎖モンテカルロ法 - Qiitaを参考に少しだけいじった。やはり平均値は引いた方が少し早いようだ。Stanでは変数名にドットを使用できないので注意が必要。

Rコード

library(rstan)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

# データのダウンロード
download.file(destfile = "ch9.Rdata", "http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/gibbs/d.RData")
load("ch9.Rdata")

# データlistの準備
list.data <- list(
  x = d$x,
  y = d$y,
  N = nrow(d),
  mean_x = mean(d$x)
)

# stanで計算
fit <- stan(file = "sample.stan", data = list.data,
            iter = 1500, chains = 3, thin = 3, warmup = 100)

緑本メモ - 盆栽日記を参考にした。

結果の確認

> print(fit)
Inference for Stan model: sample.
3 chains, each with iter=1500; warmup=100; thin=3; 
post-warmup draws per chain=467, total post-warmup draws=1401.

        mean se_mean   sd   2.5%    25%    50%    75%  97.5%
beta1   1.97    0.00 0.08   1.81   1.92   1.98   2.03   2.13
beta2   0.08    0.00 0.07  -0.06   0.03   0.08   0.13   0.22
lp__  143.99    0.03 1.02 141.56 143.62 144.31 144.70 144.95
      n_eff Rhat
beta1  1287    1
beta2  1263    1
lp__   1315    1

Samples were drawn using NUTS(diag_e) at Mon Aug 15 00:45:10 2016.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
plot(fit)

f:id:Rion778:20160815014841p:plain

結果のプロット

トレースプロットはrstanパッケージのtraceplot()で描ける。バーンイン期間の表示の有無等も設定できるようだ。

traceplot(fit, inc_warmup = TRUE)

f:id:Rion778:20160815015315p:plain
また、サンプリング結果はextract()で数値として取り出せるので、これを用いることで事後分布のプロット等も描ける。

post <- extract(fit)

old <- par(mfrow = c(1, 2))
plot(density(post$beta1), main = "beta 1")
rug(post$beta1)
plot(density(post$beta2), main = "beta 2")
rug(post$beta2)
par(old)

f:id:Rion778:20160815020229p:plain
(※2016/08/17追記)
単に事後分布の密度推定をプロットしたいだけであればstan_dens()を使用するのが簡単だった。chain毎の推定結果を重ねて表示するオプションなんかもある。可視化のためのツールは一通り用意されている(?stan_plotを確認するとよい)。

stan_dens(fit1, separate_chains  = TRUE)

f:id:Rion778:20160817204956p:plain
(※追記ここまで)

テキストの図を再現する

テキストの図9.6については、生態学データ解析 - 本/データ解析のための統計モデリング入門のfig09_06.Rを少し修正すれば描くことができる。

# テキストの図を再現
add.mean <- function(bb1, bb2, lty = 2, lwd = 1, ...){
  lines(d$x, exp(bb1 + bb2 * (d$x - mean(d$x))), lty = lty, lwd = lwd, ...)
}

plot(d$x, d$y, type = "n", xlab = "x_i", ylab = "y_i")
for (i in 1:length(post$beta1)) {
  add.mean(post$beta1[i], post$beta2[i], lty = 1, col = "#00000004")
}
points(d$x, d$y)
add.mean(median(post$beta1), median(post$beta2), lty = 1, lwd = 2)

f:id:Rion778:20160815020806p:plain

WindowsにRStanをインストール

基本的にはRStan Getting Started (Japanese) · stan-dev/rstan Wiki · GitHubに従っていけた。

環境

Rtoolsのインストール

  • Building R for WindowsよりRtools34.exeをダウンロード。
  • セットアップ中、"Edit the system PATH"にチェック。

Makevarsファイルを作成するため、次のコードをRから実行。

dotR <- file.path(Sys.getenv("HOME"), ".R")
if (!file.exists(dotR)) dir.create(dotR)
M <- file.path(dotR, "Makevars")
if (!file.exists(M)) file.create(M)
cat("\nCXXFLAGS=-O3 -mtune=native -march=native -Wno-unused-variable -Wno-unused-function",
    file = M, sep = "\n", append = TRUE)

RStanのインストール

以下のコードをRから実行。

Sys.setenv(MAKEFLAGS = "-j4")
install.packages('rstan', repos = 'https://cloud.r-project.org/', dependencies = TRUE)

インストール完了後、Rを再起動。

C++が動くかどうかの確認

fx <- inline::cxxfunction( signature(x = "integer", y = "numeric" ) , '
  return ScalarReal( INTEGER(x)[0] * REAL(y)[0] ) ;
                           ')
fx( 2L, 5 )
# => 10

最初"cygheap base mismatch detected"が出たので、cygwin1.dllで検索してRtoolsに入っていたもの以外をリネームした(gnupackに入っていた)。
また、「もしg++のバージョンが4.9かそれ以降ならば...」に続くコードを実行してしまうと"Cannot export __gnu_lto_v1: symbol not defined"というようなエラーが出た。Makevarsファイルから該当箇所を削除したら直った(cf.Google グループ)

カテゴリカルデータ解析読書メモ(第10章)

決定木

データ全体を説明変数を用いて段階的にグループ分けする手法。

決定木を書く

テキストではmvpartパッケージを使用しているが、現在は公開が停止されているので代わりにpartykitパッケージを使用する。
rpart()の代わりに使う関数はctree()であるが、引数の入れ方はほとんど変わらない。結果はpartyオブジェクトになっているので、そのままプロットすれば決定木が表示される。

library(partykit)
library(vcd)
result <- ctree(survival ~ stage + operation + xray, weight = Freq,
                data = OvaryCancer)
result
# Model formula:
# survival ~ stage + operation + xray
#
# Fitted party:
# [1] root
# |   [2] stage in early: yes (n = 158, err = 19.6%)
# |   [3] stage in advanced: no (n = 141, err = 16.3%)
# 
# Number of inner nodes:    1
# Number of terminal nodes: 2
plot(result)

f:id:Rion778:20160807214654p:plain

量的変数を用いる

ctree()を使った結果はテキストと若干異なる。

# 量的変数を用いたグループ分け
data(iris)
result <- ctree(Species ~ ., data = iris)
result
# Model formula:
#   Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
# 
# Fitted party:
#   [1] root
# |   [2] Petal.Length <= 1.9: setosa (n = 50, err = 0.0%)
# |   [3] Petal.Length > 1.9
# |   |   [4] Petal.Width <= 1.7
# |   |   |   [5] Petal.Length <= 4.8: versicolor (n = 46, err = 2.2%)
# |   |   |   [6] Petal.Length > 4.8: versicolor (n = 8, err = 50.0%)
# |   |   [7] Petal.Width > 1.7: virginica (n = 46, err = 2.2%)
# 
# Number of inner nodes:    3
# Number of terminal nodes: 4
plot(result)

f:id:Rion778:20160807215145p:plain
また、predict()にtypeオプションの指定は必要ない。

iris$est <- predict(result)
xtabs( ~ Species + est, data = iris)
#             est
# Species      setosa versicolor virginica
#   setosa         50          0         0
#   versicolor      0         49         1
#   virginica       0          5        45

順序カテゴリカルデータ

パッケージrpartOrdinalも既に公開されていないが、これもctree()で解析できる。目的変数が順序カテゴリカルデータになっていれば特にオプションの指定は必要ない様だ。
データセットはパッケージRSADBEに同様の物が入っているのでこれを利用する。

library(RSADBE)
data("lowbwt")
lowbwt$Category <- factor(ifelse(lowbwt$BWT <= 2500, 3,
                                 ifelse(lowbwt$BWT <= 3000, 2,
                                        ifelse(lowbwt$BWT <= 3500, 1, 0))),
                          ordered = TRUE)
ord.rpart <- ctree(Category ~ AGE + LWT + RACE + SMOKE + PTL + HT + UI + FTV,
                   data = lowbwt)
ord.rpart
 
# Model formula:
#   Category ~ AGE + LWT + RACE + SMOKE + PTL + HT + UI + FTV
# 
# Fitted party:
#   [1] root
# |   [2] LWT <= 109: 3 (n = 42, err = 50.0%)
# |   [3] LWT > 109
# |   |   [4] UI <= 0
# |   |   |   [5] SMOKE <= 0
# |   |   |   |   [6] RACE <= 1: 0 (n = 39, err = 43.6%)
# |   |   |   |   [7] RACE > 1: 1 (n = 44, err = 63.6%)
# |   |   |   [8] SMOKE > 0: 3 (n = 47, err = 66.0%)
# |   |   [9] UI > 0: 3 (n = 17, err = 41.2%)
# 
# Number of inner nodes:    4
# Number of terminal nodes: 5

plot(ord.rpart)

f:id:Rion778:20160807215932p:plain