プランクの法則を綺麗にプロットしたい

技術的に難しい部分とか特に無いのでこちらで。

プランクの法則

Wikipediaで確認しよう。ちなみに英語版のほうが詳しいぞ。

あまりきちんと説明するとボロがでるのでいい加減に説明すると、暖かいものほど強い光を出して、暖かいものほど光の波長のピークが短波長側にずれるというやつだ。

よくあるプロット

こういうのよく見かけるだろう。

f:id:Rion778:20170502233943p:plain

この手のプロットだと英語版Wikipedia(Planck’s law - Wikipedia)にあるやつが割と好きだ。太陽の表面温度くらいの放射スペクトルがちょうど可視光になってるのがよく分かる。

ほしいプロット

温度も連続変数のはずだ。だとすれば、なんかこう、グラデーションっぽい感じでプロットしてほしい。

やっていく

パッケージ

今回はggplot2だけです。

library(ggplot2)

計算

k = 1.38e-23 # Boltzmann constant
h = 6.62e-34 # Planck constant
c = 3e8      # speed of light

p_rad <- function(lambda, t){
  2*h*c^2/(lambda^5) * 1/(exp((h*c)/(lambda*k*t)) - 1)
}

lim = 5000
t = lim - (1:(lim^(1/1.3)))^1.3
lambda = 10^-9 * 1:2000
param = expand.grid(t=t, lambda=lambda)
u = cbind(param, rad = p_rad(param$lambda, param$t))

温度と波長について少しずつ値を変えながら色々な組合せでの放射強度を計算していくわけだけれども、温度については先程のプロットからも分かるように、高温側では少し値が変わっただけで放射強度も大きく変わり、低温側ではそれほど変わらないという関係がある。

従って、高温側の隙間を小さく、低温側の隙間を大きくするような感じで調整してやると多分ポイントのムラが減る。

また、「色々な組み合わせ」はexpand.gridでやればパッと作れる。最初forrbindで適当にやったら大変時間がかかった。

プロット

jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

ggplot(u, aes(x = lambda, y = rad, color = t)) +
  geom_point(size = 0.5) + 
  scale_color_gradientn(name = "Temperature(K)", colors = jet.colors(7)) +
  theme_classic() +
  labs(x = "Wavelength (m)", y = "Intensity (W/sr/m^2)") +
  theme(legend.position = "right")

jetというのは少し前までMATLABのデフォルトだったカラーマップで、色々なところでよく見かける。作り方はUsing jet colormap in R and ggplot2 – Hi!!を参考にした。

そしてgeom_pointでプロットする。とても、とても時間がかかる。

f:id:Rion778:20170502235451p:plain

ほかに上手いやり方があるような気がして仕方がないが、思いつかなかった。

有機質資材からの窒素放出量を予測したい

有機質資材からの窒素放出

ほ場に有機質資材、例えば牛ふん堆肥であったり稲わらであったりというものを施用した場合、これらに含まれている窒素というのは必ずしもすぐには使いやすい形にならない。

有機質資材中の窒素は基本的には微生物の分解を受けることで植物が利用しやすい形態となる。有機質資材はそれぞれ特有の割合で炭素と窒素を含んでいて、これをC/N比と言う。そして、有機物を分解する微生物も炭素と窒素を体に含んでいて、微生物のC/N比はだいたい20くらいである。従って、C/N比が20より小さいときは、微生物が有機物を体内に取り込むと、窒素のほうが多いので窒素が土壌中に放出される。アバウトな説明だがこれが有機物の肥料効果が現れる流れである。従って、C/N比が低い有機物ほど肥料効果が高く、早い。

では20より大きいと炭素が放出されるかというとそうではなくて、微生物は土壌中にもともとあった窒素を利用して増殖をはじめる。従って土壌中の窒素分が吸われてしまい、負の肥料効果とでも言うべきものが現れる。C/N比が大きな有機物を投入することで土壌中の窒素がなくなってしまう現象を窒素飢餓と呼ぶ。

しかし、微生物に取り込まれた窒素も、微生物が死ねば土壌中に放出される。炭素は呼吸に伴って炭酸ガスとして減少していくので、基本的には有機物中の窒素はいつかは全て放出される(ガス化する場合もあるが)。

したがって、毎年有機物を入れ続ければ、土壌中に放出される窒素の量は年々増えていき、放出量は投入量に漸近していく。しかし、平衡状態に達するには通常長い年月がかかる。何年か連用したときの窒素放出量を知るすべはないだろうか。

土壌中における有機物の分解モデル

土壌中の有機物の分解モデルはいくつか考えられているものがあって、有機物が半減期の異なるいくつかの種類のコンパートメントの複合したものであるという前提でモデル化したもので比較的良好な結果がえられているものが多い。

前置きはこの程度にしておいて

今回はhttp://agriknowledge.affrc.go.jp/RN/2010322673のモデルを使う*1

Rでやっていこう。

まず、paramというデータテーブルが次のように用意してあるものとする。

資材名 a c f
乾燥牛ふん 0.10 0.58 0.32
牛ふん堆肥 0.04 0.15 0.81
牛ふんオガクズ堆肥 0.39 0.11 0.85
豚ぷん堆肥 0.02 0.39 0.59
豚ぷんオガクズ堆肥 0.02 0.34 0.64
バーク堆肥 0.14 -0.08 0.96
未熟稲ワラ堆肥 -0.06 0.26 0.80
中熟稲ワラ堆肥 0.03 0.18 0.79
完熟稲ワラ堆肥 0.15 0.11 0.74
稲ワラ -0.40 0.95 0.45
小麦ワラ -1.83 1.63 1.20
モミガラ 0.19 -0.11 0.92
水稲 -0.35 0.63 0.72
オガクズ -0.03 -2.77 3.80
製紙カス -1.15 0.49 1.66
発酵製紙カス 0.43 -0.37 0.94
余剰汚泥 0.65 0.28 0.07

パッケージ

こんな感じで。data.tableは気分。

library(dplyr)
library(data.table)
library(ggplot2)
library(ggrepel)

残存率、蓄積率、放出率を求める

毎年の有機物投入量を1として、t年後の残存率、累積率、放出率は、

## functions ---------------
residual_rate <- function(a, c, f, t){ # 残存率
  a * 0.01^t + c * 0.63^t + f * 0.955^t
}

accum_rate <- function(a, c, f, t){ # 累積率
  a * 0.01 * (1-0.01^t)/(1-0.01) + 
    c * 0.63 * (1-0.63^t)/(1-0.63) +
    f * 0.955 * (1-0.955^t)/(1-0.955)
}

release_rate <- function(a, c, f, t){ #放出率
  1 - (a * 0.01^t + c * 0.63^t + f * 0.955^t)
}

そしてtをベクトルで与えられるように関数を整備して、プロットに備える。

sim_fun_gen <- function(t, name, fun){
  tmp_p <- param[資材名==name]
  a <- tmp_p$a
  c <- tmp_p$c
  f <- tmp_p$f
  sapply(t, function(t) fun(a, c, f, t))
}

resid_rate_sim <- function(t, name){
  sim_fun_gen(t, name, residual_rate)
}
accum_rate_sim <- function(t, name){
  sim_fun_gen(t, name, accum_rate)
}
rel_rate_sim <- function(t, name){
  sim_fun_gen(t, name, release_rate)
}

窒素放出量を予測する

実際に一番知りたいのは窒素の放出量だ。連用何年目に何kg窒素が出てくるのかということだ。

放出量は、放出率に実際に資材が持っている窒素の量を乗ずれば出てくる。

t = 0:50
shizai <- c("豚ぷん堆肥", "豚ぷんオガクズ堆肥", "牛ふん堆肥", "牛ふんオガクズ堆肥",
            "稲ワラ", "中熟稲ワラ堆肥", "小麦ワラ", "オガクズ")
n_rate <- c(2.7, 1.4, 1.1, 0.8, 0.6, 0.5, 0.4, 0.1)
names(n_rate) <- shizai

result <- data.table(資材名=character(0), 連用年数=numeric(0), 窒素放出量=numeric(0))

for(i in 1:length(n_rate)){
  result <-
    rbind(result,
          data.table(資材名 = names(n_rate)[i],
                     連用年数 = t,
                     窒素放出量 = rel_rate_sim(t, names(n_rate)[i])* n_rate[i] * 10)
          )
}

プロットする

result %>%
  ggplot(aes(x = 連用年数, y = 窒素放出量, color = 資材名)) +
  geom_line() + 
  theme_bw(base_family = "Osaka") +
  geom_text_repel(
    data = subset(result, 連用年数==max(連用年数)),
    aes(label = 資材名),
    family = "Osaka",
    nudge_x = 10,
    size = 3,
    segment.size = 0.25,
    segment.alpha = 0.5
  ) +
  lims(x = c(min(result$連用年数), max(result$連用年数)*1.3)) + 
  theme(legend.position = "none") + 
  geom_hline(aes(yintercept = 0), lty = 2, cex = 0.1) +
  labs(title = "有機質資材を1t/10aずつ連用した水田での無機体窒素放出経過",
       x = "連用年数(年)",
       y = "無機体窒素放出量(kg/10a)")

f:id:Rion778:20170501021522p:plain

どの資材も平衡状態とみなせるような状態になるまでに50年以上はかかる。

そして、C/N比の低い資材は最初窒素を吸収する性質があることもわかる。オガクズなどは30年近く連用してやっと窒素が出てくる。このような資材を使う場合は窒素を補ってやらないと窒素飢餓が起きる危険性がある

一方、豚ぷん堆肥などはもともと窒素含量が多いというのもあるが、5年も連用しないうちに窒素の放出量は10kg/10aを超え、無視できない量になっている。このような資材の場合は投入量を制限するか、連用を控えるなどの対策をしなければ窒素過剰*2となる危険性があることが分かるだろう。

*1:実は元の論文はあまり確認してなくて、パラメータ等は堆肥・有機質肥料の基礎知識を参考にした

*2:カリウムやリン酸など他の成分が過剰となる危険性もあり、特に畜ふん堆肥は成分をよく把握して使用することが重要。

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

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 = 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からロジスティック曲線を推定する

植物の重量に上限があったとして、相対成長速度が上限と今の重量の差に比例するとしてみよう。
{\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

前回は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からゴンペルツ曲線を推定する

植物の重量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関数が便利

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

久しぶりに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まで指定できるようだ(あまり大きいと終わらなくなるが)