Project Euler Problem 11:15

R Project Euler

Problem 11: Largest product in a grid

問題自体は難しくないし計算時間もかからない。文字列の読み込みと添字操作の練習。

Problem 12: Highly divisible triangular number

gmpを使って無理矢理やっても2〜3秒だが、overviewにより効率的な方法が書いてある。約数の数同士の関係をちゃんと考えてやれば1秒以内に終わる。

Problem 13: Large sum

ただの和なのでgmpで普通に解ける。積とかなら工夫が居るかもしれない。

Problem 14: Longest Collatz sequence

試した範囲では普通にメモしながらrepeat回していくのが結局一番早かった。それでも10秒くらいかかる。

Problem 15: Lattice paths

ただの組み合わせ。

RGRからゴンペルツ曲線を推定する 3

R

ゴンペルツ曲線では相対成長速度と重量の関係はどうなってるのかと思って適当に回帰してみたら

r = a - b\ln w

が良く当てはまったので、この仮定(相対成長速度は重量の対数に比例して減少する)のもとに微分方程式

\frac{1}{w}\frac{dw}{dt} = a - b \ln w

を解いてみると、ちゃんとゴンペルツ関数が出て来る。

w_t = c_1c_2^{\exp (-bt)}

ここで

c_1 = \exp(a/b)

c_2 = \exp(\exp(bc)/b)

であり、cに時間に関するパラメータが分離されているので、始点の問題を棚上げできる。

あとは適当に。

gompertz <- function(t, c1, a, b){ # RGR vs Time
  c1 * exp(-(exp(a-b*t)/b))
}
t <- seq(0, 20, 0.5)
a <- 0.3
b <- 0.25
c1 <- 10
g <- gompertz(t, c1, a, b)
plot(t, g)

r_t <- function(w, t) sum(coef(lm(w ~ t + I(t^2)))[2:3] * c(1, 2*t[2])) / w[2]
r <- numeric(length(g))
r <- NA
for(i in 1:length(g)){
  r[i] <- r_t(g[c(i-1, i, i+1)], t[c(i-1, i, i+1)])
}

gompertz2 <- function(t, c1, a, b){ # RGR vs Weight
  exp(a/b - exp(b*c1 - b*t)/b)
}
result <- lm(r ~ log(g))
result_coef <- coef(result)
f <- function(c1) (sum(gompertz2(t, c1, result_coef[1], -result_coef[2]) - g))^2
c1 <- optimize(f, interval = c(-50, 100))$minimum
points(t, gompertz2(t, c1, result_coef[1], -result_coef[2]), col = "red", type = "l")

f:id:Rion778:20161023090414p:plain

RGRからロジスティック曲線を推定する

R

植物の重量に上限があったとして、相対成長速度が上限と今の重量の差に比例するとしてみよう。
{\displaystyle r \propto w_{max} - w }
これは適当な定数をもってきて書き換えると、要するに相対成長速度が重量の一次式で表現できるということになる。
{\displaystyle r = a - bw }
これを相対成長速度の定義に代入する。
{\displaystyle \frac{dw}{dt} = (a - bw)w }
そしてwについて解くとこうなる。
{\displaystyle w_t = \frac{a \exp (ac_1+a_t)}{b \exp (ac_1 + at) + 1} }
これは実はロジスティック方程式である(定数を整理してみると良い)。

aとbとc_1に適当な値を入れてプロットしてみよう。

a <- 0.4
b <- 0.2
c1 <- -15
t <- seq(0, 50, 1)

d <- logistic(t, c1, a, b)
plot(t, d)

f:id:Rion778:20161022203447p:plain
そして、重量と相対成長速度の関係は当然だが一次式になっている。

r_t <- function(w, t) sum(coef(lm(w ~ t + I(t^2)))[2:3] * c(1, 2*t[2])) / w[2]
r <- numeric(length(d))
r <- NA
for(i in 1:length(d)){
  r[i] <- r_t(d[c(i-1, i, i+1)], t[c(i-1, i, i+1)])
}
plot(d, r)

f:id:Rion778:20161022203615p:plain
点が不等間隔に並んでいるのも分かるだろう。重量の上限値でもって相対成長速度は0になるので、負の値になることはない。

ここからパラメータaとbを求めるのは簡単で、単に回帰すれば良い。またc_1はoptimizeを使えば求められる。ただし、c_1は曲線の立ち上がりを左右に移動させるだけなので、データのどこをt=0と定めるかによって左右される。本質的なのはaとbだ。

plot(t, d)
result <- lm(r ~ d)
result_coef <- coef(result)
result_coef # > 0.4074815  -0.2040751 
f <- function(c1) (sum(logistic(t, c1, result_coef[1], -result_coef[2]) - d))^2
c1 <- optimize(f, interval = c(-100, 100))$minimum
points(t, logistic(t, c1, result_coef[1], -result_coef[2]), col = "red", type = "l")

f:id:Rion778:20161022204522p:plain

ゴンペルツ曲線の例との違い

今回の例では相対成長速度と時間を切り離している。そのため、例えば測定データのセットが複数あった場合、相対成長速度vs重量というデータに変換してしまうことで複数の結果を束ねて解析することが容易になる。

複数の葉や果実の生長、肥大を測定した場合に、どこをt=0として測定を開始するのかというのは難しい問題である。似たような大きさのものを選んだとしても、特にサイズが小さいうちは指数関数的に生育が進むので、初期の誤差が生育曲線の位相の大きなずれとなって現れてくる。この問題を棚上げしたまま共通のパラメータを推定できるのは大きなメリットだろう。

ちなみにゴンペルツ曲線とロジスティック曲線はそもそも形が違う(ゴンペルツ曲線は対称ではない)ので、その点は注意が必要である。

RGRからゴンペルツ曲線を推定する 2

R

前回はRGRを重量の対数の変化速度として求めたが、そもそもRGRの定義は

{\displaystyle r = \frac{1}{w}\frac{dw}{dt} }

であったので、複数のデータポイントがあれば何らかの方法で平滑化して微分すれば、より精度良くRGRが推定できるだろう。例えば前後2点を含めた3点をラグランジュ補間して得た二次関数を微分するとか。

Rでラグランジュ補間をするパッケージは多分あるが面倒なので二次関数を当てはめてしまう。

r_t <- function(w, t) sum(coef(lm(w ~ t + I(t^2)))[2:3] * c(1, 2*t[2])) / w[2]

r <- numeric(length(g))
r <- NA
for(i in 1:length(g)){
  r[i] <- r_t(g[c(i-1, i, i+1)], t[c(i-1, i, i+1)])
}

そしてあとは前回のヤツを多少いじって実行する。

result <- glm(r ~ t, family = gaussian(link="log"), start = c(1, 0))
result_coef <- coef(result)
plot(t, g)
f <- function(c1) (sum(gompertz(t, c1, result_coef[1], -result_coef[2]) - g))^2
c1 <- optimize(f, interval = c(0, 100))$minimum
points(t, gompertz(t, c1, result_coef[1], -result_coef[2]), col = "red", type = "l")

すると少しはマシな推定値が出る。
f:id:Rion778:20161022170713p:plain

ただ、繰り返しになるがゴンペルツ曲線を求めるなら最初からゴンペルツ曲線をデータに当てはめたほうが良い。例えばgに適当な誤差を与えて推定を繰り返してみると、マズさが分かるだろう。

特定のタイミングにおけるRGRの推定値がほしければ、ゴンペルツ関数を微分したものに推定したパラメータを渡せば良い。微分した結果を関数として使うにはderiv()を使って引数にfunc=TRUEを指定する。

dg <- deriv(~ c1 * exp(-(exp(a-b*t)/b)), "t", func = TRUE)
attr(dg(t), "gradient")
r <- attr(dg(t), "gradient")[,1] / gompertz(t, c1, a, b)
plot(r)

f:id:Rion778:20161022172949p:plain

RGRからゴンペルツ曲線を推定する

R

植物の重量wについて、瞬間的な成長速度は植物の重量に比例すると仮定する。

\frac{dw}{dt} = rw

このときrを相対成長速度(RGR: Relative Growth Rate)と呼ぶ。もしrが変化しないのであれば、植物は指数関数的に増大を続けてしまい、地球は程なくして滅ぶだろう。

w_t = w_0 e^{rt}

しかし、植物は無限に大きくなるわけではないので、一般的には相対成長速度は時間とともに減少する。ここに例えば適当な定数aとbを用意して、rは次のように減少すると仮定してみよう。

r= e^{a-bt}

先の式に代入すれば次のようになる。

\frac{dw}{dt} = we^{a-bt}

そして面倒なのでWolfram Alphaにまかせて微分方程式を問いてもらうとこうなる。

w_t = c_1 e^{-\frac{e^{a-b t}}{b}}

さらに定数項を少し整理して、c_2をこのように定義する。

c_2 = e^{-\frac{e^a}{b}}

そして得られた式をゴンペルツ関数(Gompertz function)と呼ぶ。

w_t = c_1c_2^{e^{-bt}}

瞬間的なRGRの推定

最初の仮定と若干違ってしまうのだが、短い時間\Delta tで植物の体重を測定したためにrが変化しなかったとしよう。すると、最初の式を\Delta tの間で積分することで次の式が得られる。

r = \frac{\Delta \ln w}{\Delta t}

実際の測定で推定できそうな雰囲気がしてきただろう。これで植物体重の変化からRGRの変化が推定できる。そして、RGRの変化についてはr=e^{a-bt}と仮定したので、定数aとbは、対数リンク関数を設定したglmで推定できる。

RGRからゴンペルツ曲線を推定する

まず、aとbが残ったバージョンでゴンペルツ関数を定義し、適当なパラメータでゴンペルツ曲線を描いてみよう。

gompertz <- function(t, c1, a, b){
  c1 * exp(-(exp(a-b*t)/b))
}

t <- seq(1, 20, 0.5)
a <- 0.3
b <- 0.25
c1 <- 10
plot(t, gompertz(t,c1, a, b), type = "l")

f:id:Rion778:20161020222308p:plain

次に重量を対数変換し、その変化量としてRGRの近似値を得てみよう。

g <- gompertz(t, c1, a, b)
r <- diff(log(g))/0.5
plot(t[-1], [r)

f:id:Rion778:20161020224446p:plain]

本当の植物の葉っぱなどを実際に測定をしてみる(重量測定は破壊が伴うので、例えば葉面積で代替する)と、RGRは実際にこのような減衰をすることが多い。

あとはglmでaとbを推定する。c_1についてはoptimize()で頑張って求める。

newt <- t[1:length(r)]+0.5
result <- glm(r ~ newt, family = gaussian(link="log"), start = c(1, 0))
result_coef <- coef(result)
plot(t, g)
f <- function(c1) (sum(gompertz(newt, c1, result_coef[1], -result_coef[2]) - g[-1]))^2
c1 <- optimize(f, interval = c(0, 100))$minimum
c1 <- 10
points(newt, gompertz(newt, c1, result_coef[1], -result_coef[2]), col = "red", type = "l")

するとこうだ。
f:id:Rion778:20161020225004p:plain
測定間隔の間はRGRが変化しないと仮定しているので、少しずれる。測定間隔が広いほどずれは大きくなる。

重量データに直接当てはめる

ちなみに、もし本当にゴンペルツ曲線のパラメータがほしいのであれば、重量の測定データに直接当てはめるべきで、このような回り道をすべきではない。
例えばoptim()を使うとこうなる。

f <- function(c) (sum((gompertz(t, c[1], c[2], c[3]) - g)^2))
res <- optim(c(10, 0.5, 1), f)$par
res # > [1] 9.9993114 0.3007776 0.2500798
plot(t, g)
points(t, gompertz(t, res[1], res[2], res[3]), col = "red", type = "l")

f:id:Rion778:20161020233012p:plain

Excelで表形式に入力されてしまったデータを整形するにはINDEX関数が便利

Excel

f:id:Rion778:20161017191321p:plain
こういう「お前は何がしたかったんだ」みたいなファイルが引き継がれることありますよね？ない？そいつはよかったな！！！！！

クソみたいなデータを整形する

普通データはこう入力するだろう。こうでなくては何も始められない。
f:id:Rion778:20161017191525p:plain
だが不幸なことに現実はこうではないのだ。何かを始めるためには現実を修正しなければならない。

INDEX関数

データの整形にはINDEX関数を使う。

INDEX関数はこのように引数を指定する。

=INDEX(配列, 行番号, 列番号)

すると、配列中の所定の行、列の位置にある値を返す。

Excelに配列などあったかという話だが、ある。例えば、{}の中にコンマ","もしくはセミコロン";"で区切って要素を指定することで直接値を指定して配列を作成できる。コンマで区切ると列が、セミコロンで区切ると行が増える。例えば、次の式は6を返す。

=INDEX({1, 2, 3; 4, 5, 6}, 2, 3)

セル内の数式はデバッグしにくいのでこのようなデータ入力はすべきではないが、もしやるなら改行しておくと分かりやすいだろう。

=INDEX({1,2,3;
        4,5,6},
        2, 3)

ちなみにセル内改行はWindowsだとAlt+Enterだが、MacではCommand+Option+Enterである。Windows版とMac版ではこういう細かい違いがたくさんあって大変不快である。

そして、セルの範囲は配列としても振る舞うので、とりあえずこの場では配列=セル範囲という認識でよい。

作業列を追加する

あとは元の表に対応する行番号と列番号を生成するだけだが、これは作業列を追加すると便利である。作業列などを追加するのはなんだかダサい気がするが、INDEX関数の中でごちゃごちゃやると何かミスが有っても気づきにくいので、デバッグのしやすさを考慮すると作業列(あるいは作業行)は積極的に使っていったほうが良い。

作業列を作成するコツは最初の1ループ分(今回の例では表の1行目分)は手で入力し、それに続く要素は最初のループを参照するように指定することだ。
f:id:Rion778:20161017193915p:plain
そして、必要であれば参照に修正を加える。今回の例では行番号が3ずつ増えるよう修正している。オートフィルは予期せぬ動作をする例が多いので、あまり当てにしてはならない。また、行の追加・削除をする場合には、全体をコピーして値貼り付けして値に変換しておいた方が良い(Excelが空気を読んで調整してくれる場合も多いが、信用してはならない)。

作業列が完成したらINDEX関数を組んでいく。配列にするセル範囲を絶対参照にするのを忘れてはならない。完成したら下に伸ばしていくだけだ。作業列がちゃんと完成していれば、アクティブセルの右下をダブルクリックするだけで適切な範囲までオートフィルしてくれる。ドラッグによるオートフィルは時間がかかる上にミスしやすいので、避けられる場面では避けよう。
f:id:Rion778:20161017195121p:plain
なお、絶対参照にするショートカットはWindowsではF4だが、MacではCommand+tである。いい加減にして欲しい。

Project Euler Problem 1:10

R Project Euler

久しぶりに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章メモつづき

R

前回。
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の扱い

R ggplot2

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()