読者です 読者をやめる 読者になる 読者になる

Project Euler Problem 16:20

Problem 16: Power digit sum

gmpで解いた。

Problem 17: Number letter counts

どの数のときにどれだけの文字数がどこに追加されるかを考えて、1000個のベクトルにまとめて加算していった。
最初綴りを間違えていて解けなかった。恥ずかしい。

Problem 18: Maximum path sum I

文字列の読み込みは直接書くよりも一旦ファイルに保存してread.table()とかでやるほうが楽なようだ。
あとは動的計画法ボトムアップ方式で解く。

Problem 19: Counting Sundays

月初のみチェックすれば良いので、各月の日数を加算していけばよい。

Problem 20: Factorial digit sum

gmpで解いた。

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}
}
これは実はロジスティック方程式である(定数を整理してみると良い)。

abc_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になるので、負の値になることはない。

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

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}

しかし、植物は無限に大きくなるわけではないので、一般的には相対成長速度は時間とともに減少する。ここに例えば適当な定数abを用意して、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}と仮定したので、定数abは、対数リンク関数を設定したglmで推定できる。

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

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

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でabを推定する。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である。いい加減にして欲しい。