FORMATH Vol. 10 (2011): 333–349

333

ベイズ法を用いた単木成長予測モデル

Bayesean Method for Predicting Forest Tree
Growth Curve
二宮嘉行・吉本敦

Ninomiya, Y. & Yoshimoto, A.
キーワード:

ベイズ統計, 成長モデル, 森林資源管理, 統計数理モデル

要約:

時系列的に変化する成長現象を予測するモデルの構築には, 観測データの定期的な収集が理想的であるが, なんらかの事情により不定期かつ十分な観測が行われない場合がある. 本稿では, そうした場合, 他で得られる十分な情報あるいは, なんらのモデルから得られる情報を事前情報として用いて, ベ
イズ統計により不十分な観測データしか得られないものに対する成長モデルの構築を行った.

Abstract:

Single-tree growth modeling is one of the major concerns in forest
growth analysis. If the measurement times for each tree are suﬃcient,
its modeling can be completed with ease by some nonlinear growth
function. If it is not the case, then modeling for each tree basedon only the limited data becomes almost impossible. In order toovercome this shortcoming, some values from the available standdensity management diagram can be incorporated into the growthmodeling using Baysean method. Those values from the diagram areused in the prior distribution, so that modeling for each tree can be
carried out from the posterior distribution.

Keywords:

Bayesean method, growth modeling, forest resource management,
mathematical modeling

Received Janunary 17, 2011; Accepted February 8, 2011

334

二宮嘉行・吉本敦

1. はじめに林分あるいは林木の成長を予測することは森林資源を管理する上で重要な課題の一つであり, 予測如何により計画の妥当性が左右される. これまで森林資源管理の分野では, 実務的に古くから規範的な収穫表の作成, あるいは安藤 (1966) によって提唱された林分密度管理図の作成が行われてきた. これらは表あるいは図面上で管理計画を調整するものである. それに対し, Konohira (1995) にあるように, コンピュータの普及に伴い, プログ
ラムシステムとしてのシステム収穫表が開発されている. その例には, 白
石 (1986) が開発した LYCS (Local Yield Table Creation System), 田中(1995) の “シルブの森”, 稲田 (2005) の “穂の国『創造』” 等があり, それ
ぞれが特有の方法論の展開によりシステムを構成している.
システムが開発される一方で, それを実務的に利用するための地域特有の係数推定が同時に行われている (Matsumoto, 1997, 嘉戸・田中, 1995).
すなわち, 個々の林分の成長は地域特有の様々な要因に左右されるため,それらの特定によるシステムの調整も必要不可欠となる. さらに, 同一林分内であっても, 個々の林木は局所的な地理条件の違いや林木間の競争などの影響を受け成長しているため, SPS (Arney, 1985), PROGNOSIS(Wykoﬀ et al., 1982), Vettenranta (1999), ORGANON (Hann, 2003),
SILVA (Pretzsch et al., 2002) などのように, 距離従属あるいは独立の単木レベルの成長モデルを構築する際には, それらの影響を考慮できるようにモデルの変数及びパラメータの調整も必要となっている.
林木の成長データは経時測定データであり時系列的に変化する成長現象を予測するモデルの構築には, 観測データの定期的な収集が理想的であるが, なんらかの事情により不定期かつ十分な観測が行われない場合がある.
その場合, 一般的には欠損値の取り扱いとして, 特別な対応をすることな
く, 従来の方法によりモデルの構築が行われる. しかしながら, その場合でも, 例えば同一林分内での他の事前情報が存在するような場合, 事前情報を使用したモデル構築も可能である. 本稿では, そうした場合を想定し, 他で得られる十分な情報あるいは, なんらのモデルから得られる情報を事前情報として用いて, ベイズ統計により不十分な観測データしか得られないも

335

ベイズ法を用いた単木成長予測モデルのに対する成長モデルの構築を行う.

2. 不十分な観測データからの推定
今, {Xt |t ∈ T } をある林木の観測データとする. ここで T は観測時点のセットとする. また, ここでは下記の Richards (Richards, 1958) の成長関数を使用する.

[1]

iid

Xt = gt (Z) + εt , εt ∼ N (0, σ 2 ), gt (Z) = a(1 − e−b·t )c

なお, Z = (a, b, c)0 は gt 内の未知パラメータからなるベクトルであり, また誤差の分散 σ 2 も未知パラメータとする. まず上記観測データ{Xt |t ∈ T } を用いてパラメータ Z を推定する. 仮に観測データ数が少ない場合, すなわち T の要素数 |T | が小さいとき, 良好な推定を施すことは困難である. 特に |T | < 4 ではデータの数がパラメータ数より少なくなる
ため, 最尤推定法を使用することができないし, |T | < 3 では に対する最小二乗法を使用することもできない.
さて, 図 1 に 30 本の樹高データを用いて Richards 関数を当てはめた結果を示す. 仮にここで用いたデータを全データとすると推定されるパラメータの値は下記となる.

[2]

gt (Z) = 23.9(1 − e−0.0642·t )1.80

すなわち, Z = (23.9, 0.0642, 1.80)0 であり, これが 30 本全体のデータへの当てはめ結果である. これに対し, たとえ期間全体のデータを使用してもそれぞれの単体へ当てはめた結果は必ずしも同じにはならない. 図 2 に30 本全体から得られたパラメータを使用した曲線 (破線) と第 7 番目の林
木に当てはめた結果 (実線) を示す. この図からも明らかなように, 全体から得られた曲線と１つの単体から得られた曲線は必ずしも一致せず, 傾向
的にずれる場合も出てくる. この 1 つの単体から得られた曲線は, 全期間のデータを使用した結果であるため, 対象とする林木の成長曲線を推定す

336

敦









Height(m)

二宮嘉行・吉本









Age(year)

図 1. 全林木の樹高成長データへの成長関数当てはめ結果

るには十分である. では, データ数が少ない場合, たとえば 4 点あるいは 5
点のみの観測データしかない場合, そのデータに対する成長曲線の推定を行ってみる. 4 点では, 等間隔の 5, 10, 15, 20 年, 5 点では不規則な 6, 7,

8,17,19 年のデータをそれぞれ使用してみると, 図 3 となる. 図 3 の両図とも実線が得られた成長曲線で, 破線がその林木のデータ全てを用いた曲線,
そして点線が全ての林木データを使用した場合の曲線である. このように,データ数が少なくなったりすると, 推定される曲線が全データから得られる曲線から乖離していくことが分る. 特にサンプル外の部分においてその差は大きくなっている. すなわち, 少数のデータから得られる曲線に対する信頼性が疑わしくなる. そもそもデータ数が 4 より少なくなれば, 例え
ば最尤推定そのものができなくなる. そこで, ベイズ法により, その打開策
を考える.

3. ベイズ統計
今, 仮に図 4 のような 2 つの事象 A, B を考えると, A かつ B が起こる
確率 (同時確率) は下記となる.

ベイズ法を用いた単木成長予測モデル

図 2. 単木データへの当てはめ結果と全林木に対する成長曲線

図 3. データ観測が欠如した場合の少数データに対する成長曲線当てはめ結果[3]

P (A ∩ B) = P (B|A)P (A) = P (A|B)P (B)

337

338

二宮嘉行・吉本

敦

඲஦㇟㻌

஦㇟ 㻭㻌

஦㇟ 㻮㻌

図 4. ベイズの定理
すなわち, 事象 A が起こる確率 (周辺確率) P (A) と事象 A が起こったという条件下で事象 B が起こる確率 (条件付き確率) P (B|A) の積, あるいは事象 B が起こる確率 (周辺確率) P (B) と事象 B が起こったという条件下で事象 A が起こる確率 (条件付き確率) P (A|B) の積で表すことができる. このことから, 下記のベイズの定理が得られる （Bayes, 1763) .

[4]

P (A|B) =P (B|A)P (A)P (B)

3.1. ベイズ法「ベイズの定理」を使用してモデルのパラメータ推定を考える. ベイズの方法では, 観測データを与える X = (X1 , X2 , · · · , XT )0 が確率変数であるだけではなく, 成長曲線の未知パラメータ Z も確率変数であるとして考える. 同時に未知パラメータの分布 (事前分布) も考える. この分布については, なんらかの情報により事前に決めておく必要がある. つまり, 通常の統計分析の場合ではパラメータは定数として扱われるが, ベイズ統計ではパラメータは事前分布に従う確率変数として取り扱われる. その結果, 「ベイズの定理」を用いれば, 事前分布および全観測データに関する周辺確率分布をもとに, 全観測データを与えたもとでの未知パラメータの条件付き確率分布 (事後分布) が下記のように求められる.

[5]

fZ|X (z|x) =fX|Z (x|z)fZ (z)fX|Z (x|z)fZ (z)=RfX (x)
fX|Z (x|z)fZ (z)dz

ここで fZ|X (·) および fX|Z (·) はそれぞれ X を与えたもとでの Z の確率ベイズ法を用いた単木成長予測モデル

339(密度) 関数あるいは Z を与えたもとでの X の確率 (密度) 関数である. また, fZ (·) は Z に対する事前分布を意味する周辺確率 (密度) 関数である.
最初の等式が「ベイズの定理」によるものであり, 上述の「ベイズの定理」とは少し表現が異なるが, 事象 A と B をそれぞれ Z = z と X = x と大雑把にみなせば対応が理解できる. 通常, fZ|X (z|x) のモード, あるいは平均
値が Z に対するベイズ推定量となる. なお, fX|Z (·) は通常の尤度関数と同じ式であり, それゆえ, fZ (·) が与えられれば fZ|X (z|x) を得ることができる. すなわち, 事前分布と尤度関数により未知パラメータが推定される.
それでは, fZ (z) をいかに決定するかを考える. ここでは Z = (a, b, c)0の事前分布 fZ (z) を構築する方法として他のデータより推定する方法を用いて, かつ fZ (z) については正規分布 N(z0 , ΣZ ) の確率 (密度) 関数を仮定
する.

[6]

fZ (z) = (2π)−3/2 |ΣZ |−1 exp{−(z − z0 )0 Σ−1
Z (z − z0 )/2}

(z0 , ΣZ ) 及び σ 2 については他の林木から推定値を用いる. 具体的には, 林
木全部のデータを使用して推定したパラメータ値を z0 とし, 各林木からそれぞれ得られるパラメータ値を用いて ΣZ を算出する. なお, σ 2 は [1] で定義した成長関数の誤差項の分散である.

4. Z = (a, b, c)0 のベイズ推定量
推定には下記の通り z0 の周辺でのテイラー展開を使用し, まず非線形関数である Richards 関数を線形化する. そうすることにより, Z が正規分
布, ε も正規分布なので, X も正規分布になる.

[7]

Xt = gt (Z) + εt ≈ gt (z0 ) + At (z0 )(Z − z0 ) + εt for t ∈ T

ただし,

[8]

At (z0 ) = ∂gt (z)/∂z|z=z0

340

二宮嘉行・吉本敦

ここでベクトル表記を用いると下記となる.

[9]

X = g(Z) + ε ≈ g(z0 ) + A(z0 )(Z − z0 ) + ε

ただし, X = (X1 , . . . , XT )0 , g = (g1 , . . . , gT )0 , A = (A1 , . . . , AT )0 , ε =(ε1 , . . . , εT )0 である. 上記により Z のベイズ推定量が下記となる.

[10]

0
−1
z0 + (σ 2 Σ−1
A(z0 )0 (X − g(z0 ))Z + A(z0 ) A(z0 ))その導出は次の通りである.
まず Z を与えたもとでの X の条件付き平均は g(z0 ) + A(z0 )(Z − z0 )であり, また条件付き分散を ΣX|Z とすれば以下を得る. なお, この分散は通常推定で得られる誤差項の分散である (ΣX|Z = σ 2 I) .

[11]

X|Z ∼ N (g(z0 ) + A(z0 )(Z − z0 ), ΣX|Z )ここで, Z ∼ N (z0 , ΣZ ) と仮定したので,

µ
[12]

¶
X
∼
Z
µµ
¶ µ
g(z0 )
ΣX|Z + A(z0 )ΣZ A(z0 )0N,z0ΣZ A(z0 )0

A(z0 )ΣZ
ΣZ¶¶を得る. なお, 前述したように A(z0 ) = ∂g(z)/∂z 0 |z=z0 である. その結果, X を与えたもとでの X の条件付き分布は下記となる. Z|X ∼N (µZ|X , ΣZ|X ) ただし, それぞれの項は以下の通りである.
ΣZ|X = ΣZ − ΣZ A(z0 )0[13]

× (ΣX|Z + A(z0 )ΣZ A(z0 )0 )−1 A(z0 )ΣZ
0 −1−1= (Σ−1
Z + A(z0 ) ΣX|Z A(z0 ))ベイズ法を用いた単木成長予測モデル

341µZ|X = z0 + ΣZ A(z0 )0× (ΣX|Z + A(z0 )ΣZ A(z0 )0 )−1 (X − g(z0 ))[14]= z0 + ΣZ A(z0 )0−1
0 −1
× (Σ−1
X|Z − ΣX|Z A(z0 )ΣZ|X A(z0 ) ΣX|Z )(X − g(z0 ))= z0 + ΣZ|X A(z0 )0 Σ−1
X|Z (X − g(z0 ))この結果, ΣX|Z = σ 2 I なので, [10] で与えた Z に対するベイズ推定量が算出できる. 前節において正規分布を事前分布として仮定したが, それは本節で導いたように事後分布も正規分布となり, ベイズ推定量が易しく求まるからである. このように, 事前分布と事後分布が同じ分布となるときの事前分布を共役事前分布という (例えば Carlin and Louis, 1996, Section

2.2 や古谷 2.1 節参照). 事前情報をより反映するような事前分布, あるいは事前情報がないときに用いられる Jeﬀreys の無情報事前分布 (Jeﬀreys,1946) やその改良版 (例えば Kass and Wasserman, 1996 参照) は必ずしも共役事前分布ではない. このような事前分布を用いるとき, ベイズ推定量は一般に数値的に求めなければならない. それがベイズ法の弱点とされてきたが, その数値解を与えるマルコフ連鎖モンテカルロ法が計算機の発展とともに開発され (例えば Carlin and Louis, 1996, Section 5.4 や古谷

3 章参照), それが克服されつつあるのが現状といえる. 本ケースでは過去の知見と計算の容易さから正規分布を用いたが, それとは異なる事前分布を用いて推定量や推定曲線を改良することは今後の課題の一つである.

4.1. 林木成長データへの適用
上記ベイズ法を用いて, 観測データが少ない場合の成長曲線の推定を行ってみる. まず Z の事前分布には全体のデータを使用して得られた推定値 z0 = (a0 , b0 , c0 )0 を用いる. 成長関数は Richards 曲線を使用するため,

342

二宮嘉行・吉本

Tree ID : 730

30

Tree ID : 7

25
10

15

Height(m)

20
15

5
0
20

30

40

50

0

10

20

30
Age(year)

Tree ID : 7

Tree ID : 7

40

50

40

50

30

30

Age(year)25

Single Data
All Data
Small Data
Bayes

10

15

Height(m)

20

20

25

Single Data
All Data
Small Data
Bayes

15

Height(m)

10
5
0

10

5
0

0

5

10

Height(m)Single Data
All Data
Small Data
Bayes

20

25

Single Data
All Data
Small Data
Bayes

0

0

10

20

30

40

50

0

Age(year)

10

20

30
Age(year)

図 5. 観測データの様々な欠損を想定した成長曲線の当てはめ結果

下記を得る.

µA(z0 ) =
[15]

敦



∂gt
∂a

∂gt
∂b

∂gt
∂c

¶¯
¯
¯
¯

z=z0

0
(1 − e)=  a0 b0 c0 e−b0 t (1 − e−b0 t )c0 −1 
a0 log(1 − e−b0 t )(1 − e−b0 t )c0−b0 t c0343

ベイズ法を用いた単木成長予測モデル

30

Tree ID : 7

15
0

5

10

Height(m)2025Single Data
All Data
Small Data
Bayes

0

10

20

30

40

50

Age(year)

図 6. 観測データが１つしかない場合の成長曲線当てはめ結果

なお,

 P30
[16]

Σz = 

200i=1 (ai − a0 ) /30
P30
2

(b
−
b
)
/30
0
0
0
i=1 i
P30
2
0
0
i=1 (ci − c0 ) /30

なお, zi = (ai , bi , ci )0 は第 i 番目の林木の全期間データを使用して推定したパラメータである. 前述したように全体のデータに対し推定されたパラメータ値は z0 = (23.9, 0.0642, 1.80)0 である.
図 5 に第 7 番目の林木について観測データが 4 つの場合 (5, 10, 15, 20年), 3 つの場合 (5, 15, 20 年), 2 つの場合 5, 20 年), 1 つの場合 20 年)を想定して推定した結果を示す. 推定したパラメータは, 全期間のデータを使用したもの (Single Data), 上記の個数に限られたもの (Small Data),
及びベイズ法で算出したもの (Bayes) である. なお, 図中の All Data は事
前分布に使用された平均値を示す. パラメータ推定に当たっては, 関数の形をシグモイドにするために, べき乗項以外のパラメータを指数関数にて変換して推定した.

344

二宮嘉行・吉本敦

30

Tree ID : 73

15
0

5

10

Height(m)2025Single Data
All Data
Small Data
Bayes

0

10

20

30

40

50

Age(year)

図 7. 成長に変異がある場合の１観測データに対する成長曲線当てはめ結果

図 5 からも分かるように, データ数が少ない場合でも, 事前分布の影響を受け, 比較的もっともらしい曲線になっていることが観察できる. 他の林木についても下記のように同じような結果であった.
次に１つのデータだけが得られる場合を想定し, 5 年目から 1 年間隔で20 年までの 1 つのデータに対し, 同様にベイズ推定量を求めてみた (図 6).
図 6 からも分かるように, 1 つのデータのみであるが, ある程度の推定結果を得ることができている. 図 7 は第 72 番の林木への適用である. 成長の途中において, 急激な増加が観察されているため, 2 層のデータになっているが, このようなデータに対してもある程度, 修復ができていることが分かる.

5. おわりに
新たな問題に対しモニタリングを始める場合, データ不足は常に研究者の頭を悩ましている. 事前に様々な情報があるにも関わらず, 入手できる

ベイズ法を用いた単木成長予測モデル

345

データに依存してしまえば, 十分な成果を上げることができない. そう言った場合, ここで示したようなベイズ法を用いればある程度の修正が可能になる. 上記の例においても, たった一つのデータであっても, 事前情報を用いることにより, ある程度もっともらしい成長関数を推定することができる. なお, ここでは計算を単純にするために, 非線形の成長関数を線形化し, 対応した.
今後の課題としては, 林分に複数の成長パターンがある (例えば Yanagi-

hara et al., 2008) ときに, その情報を活かすようなベイズ法を構築することが挙げられる. つまり, 各林木に対してどの成長パターンの群に属しているかを表すベイズ推定量を与え, またその群の成長パターンを事前情報として用いようということである. このとき, パラメータの事前分布として考えるのは混合 (正規) 分布モデルであり, 他の林木からこの分布を推定する必要がある. この際に問題となるのは, 混合分布モデルが局所錐モデ
ル (Dacunha-Castelle and Gassiat, 1997) と呼ばれるものの典型例であり, その非正則性ゆえに通常のモデル選択理論の適用が妥当でないことである (例えば二宮ら, 2007 参照). そこで, 罰則付尤度を用いた方法 (Chenet al., 2001) を拡張して妥当なモデル選択理論を導く, ということがこの課題の際に必要となる.
引用文献
安藤

貴 (1966) 林分密度管理図とその使い方, 126p, 農林出版, 東京.

Arney, J.D. (1985) A modeling strategy for the growth projection ofmanaged stands, Can. J. Forest Res. 15: 511–518.
Bayes, T.R. (1763) An essay towards solving a problem in the doctrineof chances, Philosophical Transactions of the Royal Society, 53:
370–418. Reprinted in Biometrika, 45: 243–315, 1958.
Carlin, B.P. and Louis, T.A. (1996) Bayes and Empirical Bayes Methods
for Data Analysis, Chapman & Hall.
Chen, H., Chen, J. and Kalbﬂeisch, J.D. (2001) A modiﬁed likelihood
ratio test for homogeneity ﬁnite mixture models, J. Roy. Statist.

346

二宮嘉行・吉本敦

Soc., Series B, 63: 19–29.
Dacunha-Castelle, D. and Gassiat, E. (1997) Testing in locally conic
models and application to mixture models, ESAIM Probability andStatistics 1: 285–317.
古谷知之 (2008) ベイズ統計データ分析, 朝倉書店, 東京.

Hann, D.W. (2003) ORGANON User’s manual Edition 7.0. Dept. ofForest Resources, Oregon State Univ., Oregon, USA.
稲田充男 (2005) 穂の国「創造」<http://www.vector.co.jp/soft/win95/

business/se360304.html>
Jeﬀreys, H. (1946) An invariant form for the prior probability in estimation problems, Proceedings of the Royal Society of London, Ser.
A. 186: 453–461.
Kass, R.E. and Wasserman, L. (1996) The selection of prior distributions by formal rules, J. Amer. Statist. Assoc., 91: 1343–1370.
嘉戸昭夫・田中和博 (1995) タテヤマスギの樹幹解析資料から推定したシステム収穫表の成長パラメータ, 日林論 105: 185–188.

Konohira, Y. (1995) Deﬁnition of the system yield table, J. Forest
Plann. 1: 63–67
Matsumoto, M. (1997) Construction of yield tables for sugi (Cryptomeria japonica) in Kumamoto district using LYCS, J. Forest Plann.
3: 55–62
二宮嘉行・柳原宏和・吉本敦 (2007) 非正則な統計モデルに基づく林分成
長分析, 森林資源管理と数理モデル, 6: 43–56.

Pretzsch, H., Biber, P. and Dursky, J. (2002) The single tree-based
stand simulator SILVA: construction, application and evaluation,
Forest Ecol. Manag. 162: 3–21
Richards, F.J. (1958) A ﬂexible growth function to empirical use, J.
Exp. Bot. 10: 290–300.
白石則彦 (1986) 同齢単純林の成長予測に関する研究, 東大演報 75:

199–256.

ベイズ法を用いた単木成長予測モデル347

田中和博 (1995) 林分表と樹高曲線から将来の林分表と樹高曲線を予測するシステム (シルブの森), 科研報告書「システム収穫表プログラム」(No.0456019), pp22–32, 文部省.
Vettenranta, J. (1999) Distance-dependent models for predicting thedevelopment of mixed coniferous forests in Finland, Silva Fennica
33: 51–72.
Wykoﬀ, W.R., Crookston, N.D. and Stage, A.R. (1982) User’s guideto the stand prognosis model, USDA For. Serv. Gen. Tech. Report
INT-133. 112p.
Yanagihara, H., Ninomiya, Y. and Yoshimoto, A. (2008) Analysis ofgrouped growth patterns in even-aged sugi forest stand within theframework of mixture model, FORMATH vol. 7: 39-60.

348
二宮

嘉行・

二宮嘉行・吉本

敦

原

敦

宏和・吉本

付録
ベイズ推定のための R プログラム例
付録
・FullDataResN2.txt:
ベイズ推定のための 観測データ
プログラム例・

観測データ#-----------------------------------------expGC.Chapman.Richards <- function(a,b,c,t){ exp(a)*( 1 - exp(-exp(b)*t) )^c }F <- read.table("FullDataResN2.txt")
names(F) <- c("ID", "Age", "DBH", "Height", "Z")ID <c(7,12,15,31,37,52,59,62,72,73,76:78,83,86,87,93,96,100,102,104,107,111,116,120,126,127,131,133,135)
tp <- c(23.89,0.064,1.8) # true values of parameters for heightm <- c(5,10,15,20) # years used to make curvees <- matrix(0,length(ID),3); xv <- NULL
j <- 1x <- F$Age[F$ID==ID[j]]y <- F$Height[F$ID==ID[j]]da <- data.frame(x,y)
RSS.f <- function(b) { temp <- y-expGC.Chapman.Richards(b[1], b[2], b[3], x);t(temp)%*%temp }opt <- optim(c(log(tp[1]),log(tp[2]),log(tp[3])), RSS.f, control=list(maxit=10000))
beta <- opt$pary.hat <- c(1:length(x))
y.hat[1:length(x)] <- expGC.Chapman.Richards(beta[1], beta[2], beta[3], x)
conv<- opt$convergence
if(conv != 0) stop
es[j,] <- c(exp(beta[1:2]),beta[3])err <- y - y.hat
xv <- c(xv,err)}
xv <- mean(xv^2)
va <- mean((es[,1]-tp[1])^2)
vb <- mean((es[,2]-tp[2])^2)
vc <- mean((es[,3]-tp[3])^2)### plotfor(j in 1:1){
xx <- F$Age[F$ID==ID[j]]yy <- F$Height[F$ID==ID[j]]t <- 0:50### est: MLEs of parameters obtained from restricted measurementsaa <- xx==m[1]for(ii in 1:length(m)) aa <- aa | xx==m[ii]x <- xx[aa]; y <- yy[aa] # data for restricted time points
da <- data.frame(x,y)
est <- c(0,0,0)
RSS.f <- function(b) { temp <- y-expGC.Chapman.Richards(b[1], b[2], b[3], x);t(temp)%*%temp }opt <- optim(c(log(es[j,1]),log(es[j,2]),log(es[j,3])), RSS.f, control=list(maxit=10000))
beta <- opt$pary.hat <- c(1:length(x))
y.hat[1:length(x)] <- expGC.Chapman.Richards(beta[1], beta[2], beta[3], x)
conv<- opt$convergence
if(conv != 0) stop
est <- c(exp(beta.6[1:2]),beta.6[3])ベイズ法を用いた単木成長予測モデル
ベイズ法を用いた単木成長予測モデル349### bp: Bayes estimators obtained from restricted measumentsaa <- matrix(c((1-exp(-tp[2]*x))^tp[3],tp[1]*tp[2]*tp[3]*(1-exp(-tp[2]*x))^(tp[3]-1)*exp(-tp[2]*x),tp[1]*(1-exp(-tp[2]*x))^tp[3]*log(1-exp(-tp[2]*x))),length(x),3)
bp <- solve(t(aa)%*%aa/xv+diag(c(1/va,1/vb,1/vc)))%*%(t(aa)%*%(y-tp[1]*(1-exp(-tp[2]*x))^tp[3])/xv)+tp
title <- paste("Tree ID :",ID[j])plot(x,y,lwd=8,xlim=c(0,50),ylim=c(0,30),xlab="Age(year)",ylab="Height(m)",type="p",main=title)lines(t,es[j,1]*(1-exp(-es[j,2]*t))^es[j,3],lwd=3,lty=1,col=1)
lines(t,tp[1]*(1-exp(-tp[2]*t))^tp[3],lty=2,lwd=3,col=4)
lines(t,est[1]*(1-exp(-est[2]*t))^est[3],lty=3,lwd=3,col=3)
lines(t,bp[1]*(1-exp(-bp[2]*t))^bp[3],lty=4,lwd=3,col=2)
legend(25,12, c("Single Data","All Data","Small Data","Bayes"),lwd=3,lty=c(1:4),col=c(1,4,3,2))}