2017-07-20

Rによるハンバーガー統計学 #2

R Statistics

2.1 平均的ポテトを推定する

前回も説明したとおり、標準のvar()、sd()でそれぞれ不偏分散、標本標準偏差が計算される。

wakupote <- c(47, 51, 49, 50, 49, 46, 51, 48, 52, 49)
mean(wakupote); var(wakupote); sd(wakupote)

2.3 信頼区間/区間推定

テキスト通りにやるならこのように。

# t値を求める
qt(0.025, 9) # 自由度9の場合の95%信頼区間の下側t値 qはquantileから
qt(0.975, 9) # 同上側t値
qt(0.005, 9) # 自由度9の場合の99%信頼区間の下側t値
qt(0.995, 9) # 同上側t値

# 信頼区間を求める
se <- function(x){ sd(x) / sqrt(length(x)) } # 標本標準誤差を求める関数

# 95%信頼区間
mean(wakupote) + qt(0.025, 9) * se(wakupote) # 下側
mean(wakupote) + qt(0.975, 9) * se(wakupote) # 上側
# 99%信頼区間
mean(wakupote) + qt(0.005, 9) * se(wakupote) # 下側
mean(wakupote) + qt(0.995, 9) * se(wakupote) # 上側

もちろん、Rを使っているのであればt.test()に突っ込むだけという姿勢が正しい。必要であれば信頼区間だけを取り出すこともできる。

> t.test(wakupote, conf.level = 0.95)$conf.int
[1] 47.85957 50.54043
attr(,"conf.level")
[1] 0.95
> t.test(wakupote, conf.level = 0.99)$conf.int
[1] 47.27432 51.12568
attr(,"conf.level")
[1] 0.99

t値が計算できるのだから実用的な意味は特に無いがt分布表も作ってみよう。

t_table <- data.frame(
  df = (df <- c(1:30, 40, 60, 120, Inf)),
  p95 = qt(0.975, df),
  p99 = qt(0.995, df)
)
knitr::kable(t_table)

df	p95	p99
1	12.706205	63.656741
2	4.302653	9.924843
3	3.182446	5.840909
4	2.776445	4.604095
5	2.570582	4.032143
6	2.446912	3.707428
7	2.364624	3.499483
8	2.306004	3.355387
9	2.262157	3.249836
10	2.228139	3.169273
11	2.200985	3.105806
12	2.178813	3.054540
13	2.160369	3.012276
14	2.144787	2.976843
15	2.131449	2.946713
16	2.119905	2.920782
17	2.109816	2.898230
18	2.100922	2.878440
19	2.093024	2.860935
20	2.085963	2.845340
21	2.079614	2.831360
22	2.073873	2.818756
23	2.068658	2.807336
24	2.063899	2.796940
25	2.059539	2.787436
26	2.055529	2.778714
27	2.051830	2.770683
28	2.048407	2.763263
29	2.045230	2.756386
30	2.042273	2.749996
40	2.021075	2.704459
60	2.000298	2.660283
120	1.979930	2.617421
Inf	1.959964	2.575829

2.4 平均的チキンを推定する

🐔

chicken <- c(568, 530, 581, 554, 536, 518, 564, 552)
mean(chicken)
var(chicken)
sd(chicken)
se(chicken)
t.test(chicken, conf.level = 0.95)$conf.int
t.test(chicken, conf.level = 0.99)$conf.int

2.9 通過テスト

(1)
1. 母集団
2. 標本
3. 無作為抽出
4. 不明である

dについて実際どのような解答が想定されているか分からないが、問題文の情報量だけで度数分布の形状を推定することは不可能である。例えば、テストの総問題数が1問であった場合を考えてみれば想像できるだろう。実際にはどのような分布も想定すべきでないという姿勢で分析にかかるべきではないのか。

もっとうまいやり方があるような気がして仕方がない。

# 95%信頼区間
65 + qt(0.025, 500) * sqrt(60)/sqrt(500)
65 + qt(0.975, 500) * sqrt(60)/sqrt(500)
# 99%信頼区間
65 + qt(0.005, 500) * sqrt(60)/sqrt(500)
65 + qt(0.995, 500) * sqrt(60)/sqrt(500)

信頼区間の計算式においては確率pのみが変数であるので、逆に計算された信頼区間を変数として確率pが計算されるような関数を考えることができる。このとき、pの範囲が2.5%〜97.5%となるような区間が95%信頼区間、同様に0.5%〜99.5%となるような区間が99%信頼区間である。

「信頼区間に母平均が含まれる確率が95%または99%である」という説明は信頼区間の説明としては正しくない(例えばRで楽しむ統計 (Wonderful R 1)に詳しい解説がある。上述の方法とは異なる信頼区間の構成方法も紹介されている。)。

2017-07-20

Rによるハンバーガー統計学 #1

R Statistics

(ハンバーガー統計学にようこそ！をRでやってみたものです。)

1.1 ポテトの長さの平均は

# ワクワクバーガーのポテトの長さデータ
waku <- c(3.5, 4.2, 4.9, 4.6, 2.8, 5.6, 4.2, 4.9, 4.4, 3.7,
          3.8, 4.0, 5.2, 3.9, 5.6, 5.3, 5.0, 4.7, 4.0, 3.1,
          5.8, 3.6, 6.0, 4.2, 5.7, 3.9, 4.7, 5.3, 5.5, 4.7,
          6.4, 3.8, 3.9, 4.2, 5.1, 5.1, 4.1, 3.6, 4.2, 5.0,
          4.2, 5.2, 5.3, 6.4, 4.4, 3.6, 3.7, 4.2, 4.8)

# 総和を求める
sum(waku)

# 平均を求める
sum(waku)/49

# 平均を関数で求める
mean(waku)

# モグモグバーガーのポテトの長さデータ
mogu <- c(4.5, 4.2, 3.9, 6.6, 0.8, 5.6, 3.2, 6.9, 4.4, 4.7,
          3.8, 3.0, 3.2, 4.9, 7.6, 3.3, 7.0, 3.7, 3.0, 4.1,
          5.8, 4.6, 4.0, 2.2, 7.7, 3.9, 6.7, 3.3, 7.5, 2.7,
          5.4, 5.8, 5.9, 3.2, 5.1, 3.1, 6.1, 4.6, 2.2, 4.0, 
          6.4, 5.2, 3.3, 6.4, 6.4, 2.6, 2.6, 5.2, 5.8)

# 平均値を求める
mean(mogu)

1.2 度数分布

実はhist()の戻り値の中に度数分布表は含まれていたりするのだけれど、このタイミングでそれもどうかなと思うので別解。

table(cut(waku, breaks = 0:8))
table(cut(mogu, breaks = 0:8))

cut()はnumericを指定の区間で区切ったfactorに変換する。これをtable()で集計すれば度数分布表が得られる。結果を見たほうが早いだろう。

> table(cut(waku, breaks = 0:8))

(0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] 
    0     0     1    14    19    13     2     0 
> table(cut(mogu, breaks = 0:8))

(0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] 
    1     0     7    13     8     9     8     3

1.4 分散と標準偏差

標準のvar()とsd()は不偏分散に基づいた標本分散と標本標準偏差を計算するので、テキストに沿った値が欲しければ関数を定義する必要がある。

my_var <- function(x){mean((x - mean(x))^2)}
my_sd  <- function(x){sqrt(my_var(x))}
# 組み込みのvar(), sd()は不偏分散に基づくので計算結果が異なる

# ワクワクバーガーの分散と標準偏差
my_var(waku)
my_sd(waku)

# モグモグバーガーの集計
my_var(mogu)
my_sd(mogu)

1.5 チキンで行こう

waku_c <- c(135, 142, 149, 146, 149, 144, 136, 138,
            156, 153, 150, 147, 136, 160, 142, 157)
mogu_c <- c(144, 143, 139, 166, 169, 144, 147, 138, 
            176, 133, 170, 137, 146, 140, 122, 177)
# 平均、分散、標準偏差の計算
mean(waku_c); my_var(waku_c); my_sd(waku_c)
mean(mogu_c); my_var(mogu_c); my_sd(mogu_c)

1.9 通過テスト

score <- data.frame(
  桜組 = c(78, 62, 81, 59, 72, 68, 75, 65, 80, 60, 78, 62, 70),
  桃組 = c(70, 72, 68, 75, 65, 71, 69, 76, 64, 80, 60, 73, 67),
  柳組 = c(57, 59, 55, 62, 52, 58, 56, 63, 51, 67, 47, 60, 54)
)

my_summary <- function(x){
  data.frame(
    平均     = round(mean(x), 2),
    分散     = round(my_var(x), 3),
    標準偏差 = round(my_sd(x), 3)
  )
}

apply(score, 2,  my_summary)

2017-07-06

なんでn-1なんですか

R Statistics

不偏分散はなぜnではなくn-1で割るのか、という問題。

以前似たようなことやったけどきちんと証明していなかったのと今ならどう書くか試してみたかったので。

証明をする

確率変数 $X$ の実現値を $X_i$ 、平均を $\bar{X}$ とする。

まず、 $\sum\left(X_i - \bar{X}\right)^2$ の期待値を求めよう。

任意の数 $\mu$ を足して引く。

$E\left[\sum\left((X_i - \mu) - (\bar{X} - \mu)\right)^2\right]$

ここで、

$\sum(x_i - \bar{x})^2 = \sum x_i^2 - n\bar{x}^2$

であることに注意すれば、

$\sum E \left[ (X_i - \mu)^2 \right] - nE \left[ (\bar{X} - \mu)^2 \right]$

$\sum E \left[ (X_i - \mu)^2 \right] - nE \left[ \left( \frac{1}{n} \sum (X_i - \mu) \right)^2 \right]$

ここで任意の数 $\mu$ が母平均に等しかったとすれば、

$\sum E \left[ (X_i - \mu)^2 \right] = \sigma^2$

であるので、

$n\sigma^2 - \sigma^2 = (n-1)\sigma^2$

従って、 $\sigma^2$ の不偏推定量 $\hat{\sigma^2}$ は

$\hat{\sigma^2} = \frac{1}{n-1}\sum\left(X_i - \bar{X}\right)^2$

となる。

証明をしない

標準正規分布からのサンプリングに対し、nで割った場合とn-1で割った場合の分散の分布をシミュレーションしてプロットしてみよう。

library(dplyr)
library(ggplot2)

# (n-1)/nを乗ずることでnで割った場合の分散を求める
varp <- function(x){
  var(x) * (length(x) - 1)/length(x)
}

# n = 2:20において各1000回ずつvar(rnorm(n))及びvarp(rnorm(n))を求める
d <- data.frame()
N <- 1000
for(n in 2:20){
  d <- rbind(d,
             data.frame(n = n,
                        type = "var",
                        var = replicate(N, var(rnorm(n)))))
  d <- rbind(d,
             data.frame(n = n+0.3, # プロット時にちょっと右にずらすため調整
                        type = "varp",
                        var = replicate(N, varp(rnorm(n)))))
}

ggplot(d, aes(x = n, y = var, color = type)) + 
  geom_point(alpha = 0.1) +
  stat_summary(fun.y = mean, geom = "line") + # 平均をプロット
  theme_classic() +
  geom_hline(yintercept = 1, lty = 2, col = "gray", alpha = 0.3) # 真の分散=1を示す線

f:id:Rion778:20170706233312p:plain

このようにnで割った分散の期待値は真の分散(=1)より少し小さくなる傾向がある。一方でn-1で割った場合の期待値は真の分散に一致する。また、nで割った場合でもnが大きくなれば真の分散に近づく傾向があることが分かる。実際、nで割った場合でもn→∞で真の分散に一致するので、nで割った分散は不偏推定量ではないが一致推定量ではある。

2017-06-06

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

Math

ユーザーの平均継続期間は「解約率=一定」と仮定すれば「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_t$ 、 $dt$ 時間経過後のユーザー数を $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) - t$ 、 $dx_t = x_{t+dt} - x_t$ であるため、 $dx_t$ は $x_t$ が時間とともに減少するという仮定のもとでは負の値であり、これにマイナスを付けたものが解約ユーザー数となる

2017-05-29

Project Euler Problem 22

問

Problem 22 - PukiWiki

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

解答

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

2017-05-29

Project Euler Problem 21

R Project Euler

問題の概要

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以上の場合の処理方法が分からないが、とりあえず今回の範囲ではそのような友愛数のペアは出現しない。

2017-05-08

スペインの施設園芸

Agriculture

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

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

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

このエリアの施設面積だけで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年もたてば状況も変わっている可能性が高い。また時間をみつけて新しい情報を確認したいところだ。

2017-05-02

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

R ggplot2 Physics

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

プランクの法則

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

プランクの法則 - 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でやればパッと作れる。最初forとrbindで適当にやったら大変時間がかかった。

プロット

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

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

2017-05-01

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

Agriculture R

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

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

有機質資材中の窒素は基本的には微生物の分解を受けることで植物が利用しやすい形態となる。有機質資材はそれぞれ特有の割合で炭素と窒素を含んでいて、これを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)")