ユーザーの平均継続期間が解約率の逆数であることを等比数列を使わずに理解する

ユーザーの平均継続期間は「解約率=一定」と仮定すれば「1/解約率」で示せる。

これについてはユーザの平均継続期間が「1/解約率」で求められることの数学的証明 - it's an endless world.等比数列を用いたわかりやすい証明があるが、「解約率=一定」とした場合の平均継続期間は「故障率=一定」とした場合のMTBFと計算上は同じなので、MTBFの計算に関する考え方を流用すると等比数列を使わずに示すこともできる。

MTBF(平均故障間隔)

MTBF(Mean Time Between Failures = 平均故障間隔)はあるシステムが故障するまでの平均動作時間で、システムの稼働時間/故障回数として求める。

例えば100個の装置を100時間観察し、その間に2つ装置が故障した場合、MTBF


\mathrm{MTBF} = \frac{100 \times 100}{2} = 5,000\ (\mathrm{hours})

のように延べ稼働時間から求める(cf. 平均故障間隔 - Wikipedia)。要するに、100個の装置100時間稼働を、1個の装置10000時間稼働と同一視する。

ユーザーの平均継続期間

現在のユーザー数をx_tdt時間経過後のユーザー数を x_{t+dt}とすれば、その間の解約ユーザー数はx_t - x_{t+dt} = -dx_tである*1。このとき、平均継続期間Uは、次式で求められる。

U = \frac{x_t dt}{-dx_t}

要するに、ユーザー1人についてx_t dtという期間中に-dx_t回解約があるとみなして、平均契約期間を計算する。サービスを解約してもなぜかすぐに再開する1人のユーザーを仮定するのだ。これは上述のMTBF計算の場合と同様の考え方である。

次に、単位時間中のユーザーの離脱率が一定値rだとすると、ユーザー数は現在のユーザー数x_tに離脱率rを乗じた速度で減少する。

すなわち、ユーザー数の変化速度\frac{dx_t}{dt}は次式で表される。

\frac{dx_t}{dt} = -rx_t

これを先程の式と比較すれば、

U = \frac{1}{r}

であり、平均継続期間は離脱率の逆数であることが示せる。

*1:dt = (t+dt) - tdx_t = x_{t+dt} - x_tであるため、dx_tx_tが時間とともに減少するという仮定のもとでは負の値であり、これにマイナスを付けたものが解約ユーザー数となる

Project Euler Problem 22

Problem 22 - PukiWiki

  • 名前のリストが入ったテキストファイルを読み込み、ソートする。
  • アルファベットを1から26までの整数に対応付けてその和を計算し、さらに名前の順位を乗ずることで「名前のスコア」を算出する。
  • 例えばCOLINなら3 + 15 + 12 + 9 + 14 = 53で、順位は938なので、53*938=49714がCOLINのスコアになる。
  • 全ての名前のスコアの和を求めよ。

解答

データの読み込みと基本的な文字列操作と演算ができれば特に問題はない…。

続きを読む

Project Euler Problem 21

問題の概要

Problem 21 - PukiWiki

  • d(n)をnの「真の約数」の和とする。
  • 「真の約数」はnの約数のうちnを除外したものと定義する。
  • d(a) == b かつ d(a) == bであるようなaとbをAmicable numbers(友愛数)と呼ぶ。
  • 完全数(すなわちa == b)は除外する。
  • 10000以下の友愛数の総和を求めよ。

解答

最初完全数を除外する条件を見逃していてハマった。

ちなみに問題文からはaが10000以下でbが10001以上の場合の処理方法が分からないが、とりあえず今回の範囲ではそのような友愛数のペアは出現しない。

続きを読む

スペインの施設園芸

スペインにはハウスがたくさんある

昔作った資料をもとに話す機会があったので、スペインの施設園芸について手持ちの文献を再度調べていた。

スペインの南部のアルメリア県にはハウスが集中している地域がいくつかあって、特にアルメリア市の西部に位置するアドラ市とエルエヒド市を中心とするエリアは明確にハウスが密集している。衛星写真に切り替えて確認してみてほしい。

このエリアの施設面積だけで2万haあると言われており、日本全国の施設面積の1/3程度に相当する。

スペインの施設はブドウの棚から発展した

この地域でこのように施設園芸が普及するまでの経緯が面白い。もともとこの地域はブドウ(ワイン用ではなく生食用だったそうだ)の栽培が盛んで、棚仕立てでブドウを栽培していたのだが、ある時生産者がこの棚の上にビニールを載せることでその下で色々な野菜を作れる、ということに気付いたそうだ。

そこから施設園芸が広まったので、ハウスの形状も独特であった。今もこの地域で最も普及しているハウスはRaspa y Amagado(ラスパ・イ・アマガド)と呼ばれるもので、ベースはビニールの載ったブドウの棚なのだが、雨水を排水するために屋根がゆるくギザギザと波打っている。Raspaは魚の骨、Amagadoは屈曲を意味し、天井部の高い部分がRaspaで谷になっている部分がAmagadoなのだそうだ。

この話を「施設と園芸」の2012年の夏号(No.158)で読んで、一度見てみたいものだ(写真は載っていたがモノクロで分かりにくかった)と思っていたのだが、昨日改めて航空写真などを確認していてそう言えばストリートビューというものがあったと気付き、適当な場所を確認してみたら記事の説明そのもののRaspa y Amagadoを確認することができた。天井部の高い部分はまさに背骨という感じだ。

上に示した場所の場合は最新の写真が2011年7月のものだが、ハウスの中は空っぽだ。この地域では夏は生産物単価が下がるのと、ハウス内部が高温になって作業が大変という理由から基本的に9月から翌年5月までの作型で栽培をするそうだ。同じ場所の2009年1月の写真を遡って確認できるので見てみると、ハウスの中にトマトが植わっているのが分かる。

夏の間は前作の片付けや次作の準備を進めつつビーチで優雅に過ごすそうだ。

ハウスはなんだか安っぽく見える(実際建造費は安いのかもしれない)が、開口部はかなり目合の細かいネットで完全に覆われているのが伺える。出入り口は基本的に二重扉となっており、外部からの害虫の侵入を防止している。中には送風装置を備えて体に付着した害虫を払い落とすような仕組みを備えた施設もあるそうだ。日本の先進地域並の装備である。

加えてスペインの農業で特徴的なのは天敵の利用が一般的であるという点だ。天敵利用は高度な技術を必要とするが、ヨーロッパでは有機農業への需要が高いこともあり、天敵利用技術が発達しており、天敵などの生物農薬を生産販売するメーカーも複数ある。販売する苗にあらかじめ天敵を定着しておくということも行われているそうだが、これは生産者の天敵利用に対する深い理解と技術力がなければ到底不可能なことだ。天敵利用に関しては明らかに日本の先を行っている。

かつては気候に恵まれ安く野菜が生産できるという点だけが利点で、オランダのような施設園芸の先進国は天敵利用技術に先んじていたため、多少値段が高くても有意性を維持できていた。というようなことが古い教科書には書いてある。

そのスペインに天敵利用が急速に広まり、HACCPやグローバルGAPを取得する生産者も続々出てきて施設園芸の先進国も苦しい立場に…というのも少し昔の話で、今ではモロッコやイスラエルが農業新興国として安くて高品質な野菜を輸出するようになっており、今度はスペインの立場も危うくなっているそうだ。

などということが書いてあった記事に関心していたのももう5年も前の話だ。5年もたてば状況も変わっている可能性が高い。また時間をみつけて新しい情報を確認したいところだ。

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

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

プランクの法則

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で解いた。