統計数理（2003）
第 51 巻 第 1 号 95–109
c 2003 統計数理研究所


特集「森林資源統計学」［原著論文］

森林資源構成表を用いた減反率の推定
藤掛 一郎

†

（受付 2002 年 7 月 29 日；改訂 2003 年 2 月 27 日）

要旨

我が国では森林が何齢級で伐採されるかを確率的に表現するために減反率という概念が示され，森林計画などにおいて利用されてきた．減反率の推定には 5 年間隔の相前後する森林資源
構成表を資料とすることが多い．しかし，この森林資源構成表から得られるデータの特性がこれまで十分に理解されてこなかった．本論文では，まず，森林資源構成表から減反率を推定する場合にデータに含まれる打ち切りや切断について論じた．そして，このデータ特性の把握に基づいて，減反率を最尤推定する方法を紹介し，また，既存の方法のうち，鈴木のモーメント法，Blandon の最尤法を再評価した．Blandon の最尤法と今回提案した方法は打ち切りや切断に対応しており，良好なシミュレーション結果を示したが，鈴木のモーメント法は問題を残した．また，これまで減反率の推定が期首の齢級構成に依存するのは好ましくないとされてきたが，打ち切りや切断への対応がなされておれば，そのこと自体は問題ではないと考えられた．
キーワード：減反率，最尤法，シミュレーション，伐採齢．
1.

はじめに

鈴木
（1972, 1979）はある地域，ある時点での伐採傾向を確率によって表現するために減反率という概念を提案した．減反率とは各齢級での伐採確率のことで，j 齢級の減反率は，ある林分の伐採齢が j 齢級である確率と定義される．その上で，鈴木はモーメント法によって，地域の森林資源構成表から，減反率を推定する方法を示した．このモーメント法による減反率の推定は今日まで地域森林計画における伐採量の予測に広く用いられている．しかし，この方法には問題点も指摘されており，いくつかの代替的な方法が示されてきた．
この報告では，森林資源構成表から減反率を推定する，もう一つの代替案を提示する．それは単純な最尤法を当てはめた推定法だが，森林資源構成表から得られるデータは切断や打ち切りのある特殊なデータであるため，この点に配慮して尤度を構成している．また，このデータ特性については，これまで十分な認識がなされてこなかった．そこで，これまで提案されてきた減反率推定法の中から，モーメント法と Blandon（1991）の方法という代表的な 2 つを採り上げ，データ特性への対応という観点から，改めてこれらの推定法を評価した．さらに，シミュレーションと実際の予測への応用を行い，今回提案する方法と鈴木のモーメント法，Blandonの最尤法を比較した．その結果，Blandon の最尤法と今回の方法は，標本の打ち切りや切断に対応しており，シミュレーションでも良好な結果を示したが，鈴木のモーメント法は打ち切り
や切断に対応できておらず，シミュレーション結果も芳しくないなど，いくつかの新たな知見を得ることができた．
†

宮崎大学 農学部：〒889–2192 宮崎市学園木花台西 1–1

96

統計数理

第 51 巻

第 1 号 2003

以下，2 節で森林資源構成表から得られるデータの特性について整理した後，3 節で新たな推定法を提案し，また，4 節で鈴木のモーメント法と Blandon の最尤法を紹介する．そして，
これらの推定法について，5 節でデータ特性への対応について論じ，6 節でシミュレーション，
7 節で予測への応用を試みる．
2.

森林資源構成表からのデータ

森林資源構成表には，対象地域の森林資源の現況を示すために，齢級別の森林面積が表示される．森林資源構成表は様々な地域範囲を対象に，また様々なタイミングで作られるが，主要なものとして，民有林の地域森林計画区ごとに，地域森林計画の樹立に合わせ 5 年ごとに作られるものがある．
このような 5 年間隔の相前後する森林資源構成表から次のデータを得て，その 1 期間，5 年における減反率が推定される．まず，期首の森林資源構成表から，各齢級に属する森林面積 rjが得られる．次に期末の森林資源構成表から rj のうち伐採を受けずに，j + 1 齢級に進級した
面積 cj が分かる．最後に，両者の差をとって，期首に j 齢級にあった森林のうち 5 年間に伐採された面積 dj ≡ rj − cj が得られる．森林資源構成表では 15 齢級（近年では 21 齢級）
以上が一括表示されるため，期首に 13 齢級以下の森林について，これらの数値が得られる．つまり，
rj ，cj ，dj が j = 1, . . . , 13 について得られる．
減反率は，森林が植栽されてから伐採されるまでの期間を確率変数と見なし，その確率変数が従う確率分布を特定するものである．ところが，この植栽から伐採までの期間はかなりの長期に渡ることが多い．これに対して観測期間が短い場合には，データにおいて次のことが生じる可能性がある．一つは右からの打ち切りである．これは，観測対象が伐採される前に観測が終わる場合，言い換えると，観測終了時点でまだ伐採が起こっていない場合に生じる．また，もう一つは左からの切断である．これは，観測対象が観測期間の始めに既に成長の途中にある
場合，言い換えると，観測対象が植栽時点ではまだ観測の対象となっていない場合に生じる．
右からの打ち切りと左からの切断という概念について説明しよう．左から右へ時間が流れる
時間軸を考える．ある初めの出来事が起こってからある終わりの出来事が起こるまでの期間（例
えば，植栽から伐採までの期間）を観測する場合に，右からの打ち切りは，終わりの出来事が起こるのが，初めの出来事が起こってから一定時間経過後の，ある時点より右，つまりある時点より後であることは分かるが，その時点より右（後）のどこであるかは特定できないという観測が得られることを指す．この右からの打ち切りは，すべての観測対象について終わりの出来事が起こる前に，観測をやめてしまう場合などに生じる．また，左からの切断は，もしも，ある個体について，終わりの出来事が起こるのが，初めの出来事が起こってから一定時間経過後の，ある時点よりも左，つまりある時点より前であったなら，その個体は観測対象とはならな
い場合を指す．左からの切断は，観測対象が初めの出来事が起こった時点では観測の対象ではなく，一定時間経過して初めて観測対象となっており，もしその時間経過する前に終わりの出
来事が起こっていたならば，終わりの出来事の発生は確認できず，観測対象とはなり得なかったという場合に典型的に生じる．このような打ち切りや切断があるデータを使い，期間の従う確率分布について推定を行うには，打ち切りや切断を考慮して行うことが必要になる．この推定上の問題については，最尤法とモーメント法に則して，3 節と 4 節で述べる．
5 年間隔の相前後する森林資源構成表から rj ，dj ，cj を得る場合には，観測期間は 5 年，1
齢級分に過ぎない．しかも，期首に既に存在していた森林だけが観測対象となる．そのため，右からの打ち切りや左からの切断があまねく発生する．具体的には，森林資源構成表から得られるデータは図 1（A）のようなものと理解できる．図では個々の線は一群の森林についての観森林資源構成表を用いた減反率の推定

図 1.

97

切断と打ち切りについての概念図．注：○が伐採，×が打ち切りによる観測終了を，線の長さが観測期間を表す．

測結果を示す．○が伐採を，×が生存途中での観測打ち切りを，線の長さが観測期間を意味する．例えば，期首に 11 齢級にあった森林については，それが 12 齢級に移行する途中で，伐採されたか，あるいは生き延びて，12 齢級に進んだかが分かる．このような情報が期首に 1〜13
齢級にある森林について得られ，そしてそれが得られる情報のすべてである．打ち切りと切断については，まず，観測期間の 5 年間に伐採されず，進級した森林はすべて右から打ち切られたものと扱わねばならない．進級した森林については，伐採齢が進級した齢級以上であることは間違いないが，それ以上のどの齢級かは特定できない，との観測が得られたことになる．また，期首に 1〜13 齢級のデータを使用するが，これらは期首にすべて生存の途中にあるから，すべての観測対象は期首にある齢級で左から切断されたものとする必要がある．例えば，期首に 11 齢級にあって，伐採が観測された森林については，植栽されてから 11 齢級で伐採されるまでが観測されたわけではない．仮にもし伐採が 9 齢級で起こっていれば，この森林は観測対象とはならなかったはずである．この森林は 11 齢級になったところで初めて観測対象となり，12 齢級に進むことなく伐採されたことが観測されたに過ぎない．
このような切断や打ち切りのある標本は，植栽時点から伐採時点までがすべての観測対象について観測された図 1（C）のような標本とは違う．また，これまで森林資源構成表の観測打ち切りといえば，15 齢級以上の一括集計によるものばかりが注目され，得られる標本は図 1（B）に示したようなものと考えられてきた（田中（1979），Blandon（1991））
が，各齢級で期末に成立している森林はすべて打ち切りによって観測が終了したものと見なすべきである．
3.

最尤法による推定

減反率の推定は，森林の植栽から伐採までの期間，つまり伐採齢を確率変数 T で表すと，Tが Γ 分布に従うことを前提に，Γ 分布の 2 つのパラメータを推定することで行われてきた．こ
こでも Γ 分布の推定を前提とする．

98

統計数理

第 51 巻

第 1 号 2003

Γ 分布の確率密度関数を，パラメータ k，λ を用いて，次のように表す．
(3.1)

f (j; k, λ) =λ(λj)k−1 e−λjΓ(k)また，生存関数は森林がある齢級を超えて生き延びる確率を表すもので，次のように定義される．
F (j; k, λ) ≡ Pr(T ≥ j) =

(3.2)

∞

f (s; k, λ)dsj

新たに提案する推定法は最尤法の応用である．最尤法で用いる尤度関数は標本に関する同時確率であり，個々の森林の伐採が独立であるとすると，個々の森林の観測された結果に対する確率をすべての森林について掛け合わせたものとなる．
ここで標本を構成する個々の観測対象は何かが問題になる．同じ施業を施されるひとまとまりの林分を一つの観測対象とすることが考えられるが，森林資源構成表のように集計された資料を用いる場合には，それは不可能である．代わりに，例えば 1 ha の森林といったように単位面積の森林を独立に伐採を受ける 1 つの観測対象とみなす方法が採られてきている．後に紹介する鈴木のモーメント法，Blandon の最尤法でもそうである．ここで提案する最尤法の応用に際してもその方法を踏襲し，1ha ずつの森林が独立に取り扱われた結果，森林資源構成表にまとめられた結果が生じたものと仮定する．そして，rj ，dj ，cj は ha 単位の整数で表示されるものとする．
次に，個々の森林について観測された事象は何かが問題となるが，切断や打ち切りのあるデー
タに対して最尤法を使う場合には，ここで切断や打ち切りを考慮に入れて，観測された事象を特定し，それに基づいて，尤度を構成せねばならない．まず，観測期間中に伐採された観測対象については，T = j が観測されたことになるが，観測が打ち切られた観測対象については，
その森林の伐採齢 T が期末にある齢級かそれ以上であることが観測された事象である．よって，尤度の構成においては Pr(T ≥ j + 1) を考えないといけない．また，切断のある観測対象に関しては，その観測対象の尤度への貢献を，期首にある齢級まで生き延びた上での，という条件付き確率の形で表現する．Pr(· | T ≥ j) である．既に述べたとおり，すべての森林について切断があるから，どの観測対象についても尤度への貢献は条件付き確率の形をとる．
また，伐採のあった森林については T = j が観測されたとしたが，T が Γ 分布に従うとし，連続変数だと仮定していることから，むしろ，伐採が観測された森林については，j ≤ T < j +1 が観測されたとする方が整合的である．つまり，個々の森林の伐採齢 T は連続変数だが，j ≤ T < j + 1
であった場合に，森林資源構成表から得られる齢級単位のデータでは，j 齢級で伐採されたことになり，dj に含められるとしておく．
以上のことから，標本に対する尤度関数を次のように表現する．
(3.3)

L(k, λ) =1≤j≤13=1≤j≤13Pr(j ≤ T < j + 1 | T ≥ j)dj Pr(T ≥ j + 1 | T ≥ j)cj

F (j; k, λ) − F (j + 1; k, λ)F (j; k, λ)dj 

F (j + 1; k, λ)F (j; k, λ)cj

このように，打ち切りや切断のあるデータを取り扱う場合，尤度を，伐採齢がある齢級以上の確率として表現したり，条件付き確率を使って表すことが必要となるのである．なお，既述の通り，森林資源構成表で 15 齢級以上が一括表示されていることを前提にしているので，尤度の構成に用いるのは期首に 13 齢級以下の標本である．

99

森林資源構成表を用いた減反率の推定（3.3）の最後の式では，生存関数を使って尤度関数を表現しており，尤度はパラメータ k，λの関数となっている．dj ，cj に観測値を入れ，尤度を最大化する k，λ を最尤推定量とする．なお，T は連続変数だが，齢級単位でしか観測されない場合のように，観測が区間単位でしか行われない場合に，ここで行っているように最尤法を適用する方法については，Lawless（（1982）,（1997）, Section 3.5）は打ち切りや切断が
pp. 259–261）に紹介がある．Klein and Moeschberger（ある場合や，区間単位の観測がなされた場合の尤度の作り方を簡潔にまとめている．また，この種の最尤推定量の性質については，Lawless（（1982）, Section 1.4, Appendix E）, Kalbﬂeisch（1980）, Section 3.4）に詳しい．
and Prentice（4.

既存の 2 つの推定法

4.1 鈴木のモーメント法
鈴木（1972，1979）は一般にモーメント法といわれる推定法を当てはめ，減反率を推定する方法を示した．これは次のような方法である．
伐採齢 T が Γ 分布に従う場合，その平均と分散は次式のように，2 つのパラメータで表さ
れる．
(4.1)

Var(T ) = k/λ2

E(T ) = k/λ ,これは母平均と母分散である．もし，標本から伐採齢の標本平均と標本分散が得られるならば，これらの式の母平均と母分散を標本平均，標本分散で置き換えて，パラメータ k，λ について解くことで，k，λ の推定値を得ることができる．
得られる標本が図 1（C）のような打ち切りも切断も含まないもので，そこから標本平均，標本分散が計算されるならば，この方法を用いて，パラメータの推定が問題なく行える．ところが，5 年間隔の相前後する森林資源構成表から得られるデータは図 1（A）
についてのものであについてのデータからモーメント法によって減反率を推定する場合，（ A）のうち○印，る．（A）
つまり伐採で終わっている観測対象面積 dj についての平均，分散によって，（4.1）の母平均，母分散が置き換えられ，パラメータが推定される．つまり，次式を k，λ について解くことで，k，λ を推定する．


(4.2)



1≤j≤131≤j≤13

jdj

dj = k/λ1≤j≤13


2

j dj


1≤j≤13
dj −
1≤j≤13


jdj



2
dj

= k/λ21≤j≤13

しかし，（ A）
において各齢級で観測された伐採面積 dj と，（C）の打ち切りも切断もない標本における齢級別伐採面積とは異なるものである．打ち切りや切断があるデータから，伐採のあった森林だけについて単純に伐採齢の平均や分散を求めても，それらは，伐採齢の母平均や母分散を直接推定する資料にはならないはずである．
例えば，ある標本では，左からの切断はないが，すべての森林が 10 齢級までで観測終了となっていたとしよう．この時，打ち切られた観測対象の面積を無視して，10 齢級までに伐採された面積だけから計算される伐採齢の平均は必ず 10 齢級以下になる．10 齢級を超えて生き延びる森林がいくら多くても，そうなる．これは，打ち切られた観測対象の伐採齢は 11 齢級以上であるのに，その情報を無視して平均を計算するからである．そして，このように計算される平均の値は，打ち切り齢級の設定に依存し，それゆえ，打ち切りのない場合の平均伐採齢とは異なるものである．

100

統計数理

第 51 巻

第 1 号 2003

また，右からの打ち切りはないが，左からの切断があり，すべての森林が 6 齢級になって初めて観測対象となる場合，その観測結果から得られる伐採齢の平均は必ず 6 齢級以上となる．
実際には 6 齢級未満で伐採される森林がいくら多くても，そうなる．観測対象に含められた森
林は 6 齢級未満で伐採されようがないという情報を無視して，平均を計算するからである．この場合にも，切断があるにもかかわらず，それを無視して計算された平均伐採齢は，切断がない場合に得られる平均伐採齢とは違うものであり，モーメント法の推定において母平均を置き換えてよいものではない．
このように，打ち切りや切断のあるデータと，そうでないデータとでは，異なる情報を含んでおり，異なる取り扱いが必要である．森林資源表からデータを得る場合，各齢級で伐採と打ち切り，切断があるため，分かりにくいが，図 1（A）の dj と（C）における齢級別伐採面積を同等に扱うべきではない．モーメント法は，打ち切りや切断のない（C）の状態のデータから，平均や分散などの標本モーメントを得ることを前提としている．これまで減反率の推定で行われてきたように，打ち切りや切断のある（ A）の状態で得られたデータから dj だけをモーメント法
による推定に用いることには無理があるといわざるを得ない．
4.2 Blandon の最尤法
Blandon が提案した方法は最尤法を用いる方法である．ただし，今回提案した方法とは異なり，そこで用いる尤度関数は，切断がなく，打ち切りは 15 齢級の一括表示によってだけ生じる
ような観測がなされた場合に対応するものである．つまり，図 1（B）のような観測を前提とした尤度の構成となっている．そして，Blandon の最尤法ではまず，dj ，cj (J = 1, . . . , 13) から，d∗j (j = 1, . . . , 13)，c∗ を導く．この d∗j (j = 1, . . . , 13)，c∗ は，現実の伐採傾向とある意味で同一の伐採傾向のもとにある森林を対象に図 1（B）のような観測を行った場合に得られる，13 齢級までの各齢級の伐採面積と，14 齢級における打ち切り面積になっている．この d∗j ，c∗ を得た後，打ち切りを考慮して，尤度関数を次のように構成する．
(4.3)

L∗ (k, λ) =∗Pr(T = j)dj Pr(T ≥ 14)c∗1≤j≤13

その上で，Blandon（1991）はこの右辺の確率表現を確率密度関数を用いて，(4.4)

Pr(T = j) = f (j; k, λ)Pr(T ≥ 14) = 1 −

(4.5)



f (j; k, λ)1≤j≤13と表している．しかし，本論では，もともと連続変数である T について齢級という区間単位でしか観測されないという考えに従って，Blandon の最尤法における尤度を次のように表現することにする．つまり，先に行ったように，d∗j ，c∗ についてそれぞれ，(4.6)

Pr(j ≤ T < j + 1) = F (j; k, λ) − F (j + 1; k, λ)(4.7)

Pr(T ≥ 14) = F (14; k, λ)が尤度となるとし，標本全体の尤度を，(4.8)

L∗ (k, λ) =∗[F (j; k, λ) − F (j + 1; k, λ)]dj F (14; k, λ)c1≤j≤13と表すことにする．これを最大化する k，λ が Blandon の推定値となる．

∗

森林資源構成表を用いた減反率の推定

5.

101

データ特性への対応

田中（1979, 1981）や Blandon（1991）は鈴木のモーメント法がデータ特性への対応という点で，好ましからざる性質を持っていることを示している．また，Blandon の最尤法に対しては，Blandon 自身が鈴木のモーメント法と比べて，その性質を論じている．しかしながら，これまでの議論は，森林資源構成表を用いた観測が図（A）のような打ち切りと切断を受けたものであることを正確には認識していなかったと思われる．そこで，データ特性への対応という観点から，鈴木のモーメント法，Blandon の最尤法，今回の方法を比べ，それぞれの推定法の性質について整理する．
森林資源構成表から得られるデータ特性への対応としては，次の 4 点がポイントである．
（a） 15 齢級以上面積の一括表示（b）全齢級での打ち切り（c）すべての観測対象の切断（d） rj 間の相対的な大きさへの依存
これまで田中や Blandon らは，森林資源構成表を使う場合に注意すべき点として，2 点を指摘してきた．1 つは，観測の打ち切りに関する（a）の 15 齢級以上の一括表示の問題である．また，
もう 1 つは，（d）の推定結果が rj 間の相対的な大きさに依存するが，rj 間の相対的な大きさは観測期間以前の過去の植伐の歴史によって決まるもので，これが推定結果に影響を与えるのは好ましくないこと，であった．鈴木のモーメント法を使った場合に，観測期間中の伐採傾向がある意味で変わらなくとも，それ以前の植伐の影響で生じる齢級構成の偏りなどが，rj 間の相対的な大きさを通して，この 5 年間における減反率の推定に影響を与える点が指摘された．しかし，図 1 の整理に基づくと，これ以外に（b）と（c）のデータ特性への対応も大事であるといえる．つまり，打ち切りはすべての齢級で生じていること，また，すべての観測対象は切断を受けていることである．
表 1 は鈴木のモーメント法，Blandon の最尤法と今回の方法が
（a）
〜
（ d）
に対応しているかど
うかをまとめたものである．まず，（a）〜（c）について見ていこう．鈴木のモーメント法はこれ
ら 3 つのデータ特性を考慮することなく作られている．森林資源構成表から得られるデータには図 1（A）のように切断と打ち切りがある．しかし，鈴木のモーメント法は，
（a）
〜
（c）
の切断と
打ち切りのいずれも存在しないかのようにデータを扱っている．具体的には，各齢級で打ち切られた cj 分の観測を推定に利用することなく，また，図 1（A）
で○で表された各齢級での伐採された分の面積である dj を図 1（C）の切断も打ち切りもない状態での各齢級での伐採面積であるかのように扱って，推定に用いている．
表 1.

データ特性への対応．

102

統計数理

第 51 巻

第 1 号 2003

これに対し，Blandon の最尤法と今回の方法は（a）〜（c）に対応している．Blandon の最尤法
が（a）に対応していることは尤度関数の作り方から分かる．また，（ b）
各齢級での打ち切りと（c）
全観測対象の切断については，そのような特性を持つデータから，それを持たない仮想データを導くことで対応している．つまり，切断と打ち切りを考慮に入れ，観測期間内に観測された情報だけを使って，仮想データが導かれる．また，今回の方法は，期首に 13 齢級以下の cj , dj
しか使わないことから，（a）に対応している．また，
（ b）
と（c）
には尤度の構成において対応していることは先に説明した通りである．
最後に，（d）については Blandon の最尤法だけが対応している．Blandon の最尤法は，初めに仮想データ d∗j ，c∗ を導くとした．ただし，尤度関数（4.8）の対数を取ったものが d∗j ，c∗ に関して線形であることから推察されるように，推定において意味があるのは 14 個の d∗j ，c∗ の相対的な大きさだけである．そこで，dj と cj の比率という 13 個の情報から，d∗j ，c∗ の相対比について 13 個の情報が導かれる．このように，パラメータの推定に用いるのは，13 個の dj と cj の比率だけである．その意味で，Blandon の最尤法では，rj 間の相対的な大きさとは関係なく，推定値が得られる．
反対に，鈴木のモーメント法や今回の方法では，各齢級における dj と cj の比率が同じであるという意味で，観測期間における伐採傾向が一定でも，推定に用いる標本によって rj 間の相対的な大きさが違えば，推定結果は変わる．鈴木のモーメント法では，パラメータの推定は，dj と cj の比率ではなく，13 個の dj に依存する．rj 間の相対的な大きさが違えば，13 個の djの相対的な大きさが変わるので，推定はこのことに影響を受ける．また，今回の方法では，それが用いる尤度関数（3.3）から推察されるように，パラメータの推定は 26 個の dj と cj に依存する．従って，今回の方法でもやはり 13 個の dj と cj の比率が一定でも rj 間の相対的な大きさによって，推定は影響を受ける．
しかし，推定値が rj 間の相対的な大きさに依存するのは，考え方によってはむしろ自然なことである．今回の方法が基づく最尤法では，一つ一つの観測対象に同じ重みを与えて，それらの同時確率分布を考える．このように実際に観測された観測対象の 1 つ 1 つに同じ重みを与える立場からすれば，推定が全齢級での観測を集約して表す 26 個の dj と cj に依存するのは自然だが，そのうち，13 個の dj と cj の比率だけに依存しなければならない必然性はない．もともと（ d）
は，鈴木のモーメント法が持つ問題点として指摘され，Blandon はこれを克服する推定法として 13 個の dj と cj の比率だけに依存する方法を示したのだが，（b），（c）の打ち切りと切断への対応が施された推定量であれば，（ d）の依存が生じていても問題はないのではないだ
ろうか．この問題については，シミュレーションの結果を得た後で，再び論じることにする．
6.

シミュレーション

今回提案した最尤法による推定法に加え，鈴木のモーメント法，Blandon の最尤法の推定法としての性質を確かめるために，シミュレーションを行った．初めに，シミュレーションの設定を紹介する．まず，Γ 分布に従う乱数を発生させた．Γ 分布のパラメータは k = 4.8，λ = 0.4とした．これは後に出てくる愛知県ヒノキ人工林についての 1966 年と 1971 年の森林資源構成表から推定される 66–71 年の k，λ の推定値に近いものである．この場合，伐採齢の平均と分散は 12 齢級と 30 齢級 2 となる．また，乱数の発生では，発生させる乱数の数が 1,000 個の場
合，10,000 個の場合，50,000 個の場合の 3 つの場合を試みた．こうして発生した乱数の小数点以下を切り下げて整数に置き換えた．つまり，j ≤ T < j + 1 であれば，T = j と置き換えた．
これは，森林資源構成表から得られる伐採齢は整数であることに対応するための処理であり，（3.3），（4.8）の尤度関数の構成は，このようなデータの発生方法に対応している．

森林資源構成表を用いた減反率の推定

図 2.

103

シミュレーションにおける植栽面積の仮定．注：横軸は期首から n 期前，縦軸は各期の植栽面積の比率を表す．

次に，この整数化した乱数を各齢級に割り振った．各齢級への割り振りの割合について 4 つのパターンを用意した．図 2 はその 4 パターンを示したものである．1 つ目のパターンは 13 の齢級に同じ割合で割り振るもので，これは，観測期間の前 13 期間に渡って同じ植栽面積が続いた場合を表す．このパターンを法正状態と呼ぶ．これに対して，残る 3 つのパターンは過去
13 期間の植栽面積に変動がある場合を表すものである．非法正状態（1）としたパターンは 7，8
期間前に植栽のピークがある場合，非法正状態（2）は植栽面積が減り続けている場合，非法正状態（3）は植栽面積が増え続けている場合を表している．このように，いくつかのパターンを用意したのは，標本が同じ Γ 分布に従って発生してくるが，期首における rj の相対的な大きさが様々な場合に，それぞれの推定法がどのような結果を出すかを見るためである．なお，現在の日本の森林の齢級構成に近いのは非法正状態（1）のパターンである．
このように仮定された割合に従って，整数化した乱数を各齢級に割り振った後，整数化した乱数の値が割り振られた齢級よりも小さかったならば，その観測対象は観測が始まる前に既に伐採されたものとして，推定に用いる標本からは除外した．よって，発生させた乱数の数，
1,000 個，10,000 個，50,000 個はシミュレーションにおいて実際に k，λ の推定に使う標本サイ
ズとは異なる．このシミュレーションでは繰り返しのたびに，k，λ の推定に使う標本サイズは異なる．一例として，50,000 個の乱数を発生させ，上述の処理を行った結果，k，λ の推定に使う標本に含まれる乱数の割合と，標本の齢級別相対頻度（発生させた乱数の総数 50,000 に対する相対頻度）
は図 3 のようになった．
最後に，標本の中で，整数化した乱数の値が割り振られた齢級と同じならば，その観測対象
は観測期間中に伐採されたものとし，整数化した乱数の値が割り振られた齢級よりも大きけれ
ば，その観測対象は打ち切りによって観測が終了するものとして，前者は dj に，後者は cj にカウントした．こうして，dj ，cj を 13 齢級分について得て，それぞれの推定法をあてはめ，k，λ の推定値を得た．
以上の作業を各乱数発生数，過去の植栽面積についての各パターンに対して 1,000 回ずつ繰

104

統計数理

図 3.

第 51 巻

第 1 号 2003

シミュレーションにおける期首の齢級構成の例．注：横軸は齢級，縦軸は発生させた乱数の総数に対する相対頻度を表す．各図中の％は標本として残った乱数の割合．
表 2.

シミュレーションの結果（パラメータ k の推定）．

森林資源構成表を用いた減反率の推定

105

表 3. シミュレーションの結果（パラメータ λ の推定）．

り返した．その結果を表 2 と表 3 にまとめた．表 2 は k の推定結果，表 3 は λ の推定結果についてまとめたもので，Bias，MSE はそれぞれ 1,000 回の平均絶対バイアスと平均二乗誤差を示している．
すべての設定において鈴木のモーメント法は大きくバイアスを含んだ推定値を生んでおり，
他の 2 つの方法と比べ，格段に成績が悪い．今回の方法と Blandon の最尤法は，法正状態で乱数発生数が 50,000 個の場合に，ほぼ同じか，Blandon の最尤法がわずかによい結果を残している以外は，すべての設定で，今回の方法が若干優れた結果を残している．
今回の方法と Blandon の最尤法は，いずれも乱数発生数，従って標本サイズが大きくなると
ともに，推定の精度が上がり，両者の差も縮まる．地域森林計画区やそれ以上の範囲では，スギやヒノキといった樹種については 1 万 ha や 5 万 ha，あるいはそれ以上の面積が存在するのが普通であり，そのような状況では，1 ha ごとの観測対象に対して独立に伐採齢が決まるという仮定が成立するならば，今回の方法と Blandon の最尤法はほぼ同様の推定値を与えるものと考えてよいだろう．
ところで，鈴木のモーメント法を用いた場合，過去の植栽面積についての設定によって推定結果が大きく変わっている．鈴木のモーメント法に対する，（ d）
rj の相対的な大きさによって推定値が影響を受けてしまうという批判は，このような事態を指したものであった．しかし，今回の方法は（d）の問題があるにもかかわらず，（d）の問題のない Blandon の最尤法と同程度かが必ずしも問題ではないことを示している．鈴木の若干優れた結果を残している．これは
（ d）

106

統計数理

第 51 巻

第 1 号 2003

モーメント法では，打ち切りや切断が存在しないかのようにして，標本平均，標本分散を計算するが，そうすると，得られる標本平均，標本分散が，大きく期首の齢級構成に依存することになる．しかし，打ち切りや切断を考慮して，観測期間中に得られる情報にだけ基づく推定量を作るならば，期首の齢級構成に依存すること自体は問題のないことではないかと考えられる．
なお，（d）の問題に関して，田中（1981）はこの問題があるために，鈴木のモーメント法は毎期の植栽面積が等しい場合にだけ利用できるとしている．しかし，ここでのシミュレーション
結果から，鈴木のモーメント法では，法正状態の場合であっても，推定がうまくいく保証はな
いことが分かる．
また，Blandon は，この鈴木のモーメント法では推定値が期首の齢級構成によって大きく変わるという問題に対処するため，（d）の問題を排除し，13 個の dj と cj の比率にだけ依存する推定法を考えたのだった．しかし，上で述べたように，今回の方法が良好なシミュレーション結果を残していることから見ても，（ d）の問題を排除しなければならないわけではないと考えられる．Blandon の最尤法に関しては，（ d）の問題を排除しようとして，13 個の dj と cj の比率にだけ依存する推定法を考えた結果，
（ b）
（
，c）
の打ち切りと切断に対応したことが評価されるべきではなかろうか．
ただし，今回の方法も Blandon の最尤法も，非法正状態（1），非法正状態（2）の場合に，法正状態，非法正状態（3）の場合よりも推定の精度が落ちている．非法正状態（1），非法正状態（2）の場合には若齢級の標本サイズが小さいことが問題ではないかと推測されたが，今回の結果だ
けからは確かなことは分からなかった．
7.

予測への応用例

最後に，鈴木のモーメント法と Blandon の最尤法，今回の方法を，伐採面積の予測に応用した結果を紹介する．田中（1981）と Blandon（1991）で用いられた 1966，1971，1976 年の愛知県のヒノキ人工林の森林資源構成表を資料とした．まず，それぞれの方法を使って 66–71 年の減反率を表す Γ 分布を推定した．次に，この推定結果と 1971 年の齢級構成とから 71–76 年の伐
採面積を予測した．具体的には，66–71 年の減反率を表す Γ 分布のパラメータ推定値が k̂，λ̂であり，71 年の j 齢級の森林面積を xj (j = 3, . . . , 13) とすると，71–76 年の j 齢級の伐採面
積を，
(7.1)

[F (j; k̂, λ̂) − F (j + 1; k̂, λ̂)] × xj

によって予測した．なお，1，2 齢級では，通常伐採は起こらないはずであるにもかかわらず，
資料の不備のために相前後する森林資源構成表からは大きな伐採面積が出てくるなど，この齢級での実現値を正確に得ることが難しい．そこで，伐採面積の予測は 3 齢級から 13 齢級までについてだけ行い，71 年，76 年の森林資源構成表から得られる伐採面積の実現値と比べた．表 4はその結果をまとめたものである．なお，表 4 の鈴木のモーメント法と Blandon の最尤法を適用した場合の予測結果は，Blandon（1991）が同様の作業をして報告した結果とは一致しない．
これは，最尤法を使う際の尤度の構成の仕方が違うことや，予測を行う際に（7.1）の代わりにBlandon が用いた方法が違うからではないかと想像される．Blandon（1991）
は（7.1）に相当する予測の方法を明示していない．しかし，いずれにせよ，表 4 の結果と Blandon（1991）が示した
結果はそれほど大きく変わるものではなく，いずれの結果を用いても，以下で述べる結論に変わりはない．
シミュレーションの結果から推察されるように，Blandon の最尤法と今回の方法は同じようなパラメータ推定値を与えており，そのため，予測結果も類似している．これに対し，鈴木のモーメント法だけは，パラメータ推定値，予測結果とも他の 2 つとはかけ離れている．鈴木の森林資源構成表を用いた減反率の推定

107

表 4. 齢級別伐採面積予測結果の比較．

モーメント法は平均，分散ともかなり小さな推定値を与えている．つまり，伐採が比較的低い
齢級で集中して起こる状態を推定している．このようになる理由としては，鈴木のモーメント法が他の 2 つの方法と違って，各齢級で打ち切りによって観測を終了した観測対象の分，すな
わち，cj を使わないことが大きいのではないかと推察される．
しかし，鈴木のモーメント法ほどではないにせよ，Blandon の最尤法，今回の方法も過大な伐採面積予測を与えた．これは，予測に際して，観測期間中に成立したのと同じ減反率が将来の予測期間に対しても成り立つことを仮定しているが，これが成立していないことに一因があると考えられる．実際，日本ではこれまでいろいろな方法で推定された減反率の平均や分散が時系列的に変化していることは田中（1980），松下（1992），Blandon（1994）によって確かめられている．高度成長期以降，日本では，安価な外材の輸入が拡大する中で，実質木材価格が著しく下落するなど，林業経営を取り巻く環境が大きく変化してきた．その中で，伐採，再造林を先送りする森林所有者が増えてきている．従って，過去に成立した減反率が将来にも成立すると仮定して，それをそのまま将来の予測に用いれば，将来の伐採面積を過大に予測する可能性が高い．従って，この応用例で Blandon の最尤法や今回の方法を用いた場合にも，なお埋めきれなかった予測と現実との乖離を埋めていくためには，林業経営を取り巻く経済環境の変化に対応して，減反率を調整する仕組みをモデル内部に含んだ，予測モデルが必要だと考えられる．

108

8.

統計数理

第 51 巻

第 1 号 2003

おわりに5 年間隔の相前後する森林資源構成表から減反率を推定する場合に，得られる標本の打ち切りと切断について明らかにするとともに，それに基づいて，新たな推定法を示し，また，既存の方法の中から鈴木のモーメント法と Blandon の最尤法を採り上げて，その性質を再評価した．
今回の方法と Blandon の最尤法は打ち切りと切断に対応しており，シミュレーションでも良好な結果を示したが，鈴木のモーメント法はこれらのデータ特性に対応しておらず，シミュレー
ションでもよい結果を残せなかった．また，これまで鈴木のモーメント法に関して，推定値が期首の齢級構成，rj の相対的な大きさに依存することを問題視する見方があり，Blandon の最尤法では，13 個の dj と cj の比率だけを利用する推定法が採られたが，標本の打ち切りや切断に対応しておれば，推定が rj の相対的な大きさに依存すること自体は問題ではないと考えられた．
鈴木のモーメント法は，今日でも，都道府県が地域森林計画を樹立する際の伐採量予測などに使われている．しかし，今回の結果はそのやり方が再検討されるべきであることを示している．ただし，伐採予測に関しては，たとえ Blandon の最尤法や今回の方法を用いたとしても，
7 節で論じたように，これらの方法も減反率の時系列変化に対応できないことから，今日のよ
うに林業経営を取り巻く環境が変化する中では，精度の高い予測を期待することはできない．
減反率の変化を扱えるような予測法の開発が待たれるが，このことについて一言付け加える
ならば，Blandon の最尤法よりも今回の方法のほうが共変量を導入するなど，モデルの拡張に対応しやすいのではないかと想像される．共変量を導入する方法については，例えば Lawless（（1982）, Chapter 6）に，パラメータを 3 つ含み，そのうちの 1 つが特定の値を取る時に Γ 分布となる，一般 Γ 分布を使い，共変量を導入する方法が示されている．本稿で紹介した方法をより実用的なものとしていくためには，このようなモデルの拡張を利用して，精度の高い予測ができないか，検討していくことが課題である．

参 考 文 献Blandon, P.（1991）. Gentan probability and censored sample theory（I）, Journal of Japanese Forestry
Society, 73, 187–196.
Blandon, P.（1994）. Movements in the parameters of the Gentan probability distribution, Journal ofJapanese Forestry Society, 76, 1–9.
. The Statistical Analysis of Failure Time Data, Wiley,Kalbfleisch, J. D. and Prentice, R. L.（1980）
New York.
. Survival Analysis, Springer, New York.
Klein, J. P. and Moeschberger, M. L.（1997）
Lawless, J. F.（1982）. Statistical Models and Methods for Lifetime Data, Wiley, New York.
松下幸司（1992）. 伐採齢級の平均と標準偏差の時系列的変化に関する研究
北海道におけるカラマツを事例に，鹿児島大学演習林報告, 20,191–211.
鈴木太七（1972）. 林業における確率過程論の応用（I），日本林学会誌，54, 234–243.
鈴木太七（1979）. 『森林経理学』，朝倉書店，東京．
田中和博（1979）. 時間打ち切り標本からの減反率の推定について, 日本林学会大会発表論文集，90,
123–124.
田中和博（1980）. 愛知県における最近の伐採傾向について，28 回日本林学会中部支部大会講演集，
189–191.
田中和博（1981）. C 方式による減反率計算の一例，林業統計研究会誌，6, 19–24.

Proceedings of the Institute of Statistical Mathematics Vol. 51, No. 1, 95–109 (2003)

109

Estimation of Gentan Probability Based on the Forest Resource Table
Ichiro Fujikake(Faculty of Agriculture, Miyazaki University)

The concept of Gentan Probability has been introduced and used in forest planning inJapan to describe the age class in which a forest stand is cut. Gentan Probability is usuallyestimated on the basis of data from the forest resource table. However, the characteristicsof the data obtained from such a source have not been explored thoroughly. This paper first
discusses the censoring and truncation in such data in estimating the Gentan Probability.
Based on the discussion, a method of maximum likelihood is developed in estimating the
Gentan Probability, and the properties of the procedure developed by Suzuki based on
the method of moments and the method introduced by Blandon are reconsidered. Both
Blandon’s method and the newly developed method are considered to be suitable fordata with the censoring and truncation and performed well in a simulation. However, theapplication of the method of moments appeared to have problems and performed worse.
Although it has been said that the estimation of the Gentan Probability shouldn’t dependon the age structure at the beginning of the observation period, it is concluded that thedependence itself is not a problem if censoring and truncation are adequately taken into
account.

Key words: Cutting age, Gentan Probability, maximum likelihood method, simulation.