統計数理（2003）
第 51 巻 第 1 号 111–120
c 2003 統計数理研究所
特集「森林資源統計学」［研究ノート］

丸太価格に基づく減反率の推定
広嶋 卓也

†

（受付 2002 年 9 月 4 日；改訂 2003 年 2 月 5 日）

要旨

減反率とは，新植された林分がある齢級で伐採される確率を指し，我が国の伐採予測のために林野庁で用いられている代表的な手法である．従来，減反率を推定する際には，伐採齢の平均と分散の情報が必要とされていたが，そのような情報は一般に入手が困難である．そこで本論では伐採の予測を全国レベルに適用することを前提として，データの入手が容易な丸太価格に基づき減反率を推定する方法を検討した．その方法とは，径級別の丸太価格および立木本数データを利用して求めた丸太価格 林齢曲線に基づき，Yoshimoto 型減反率を推定するというものである．この方法の有効性を検証するため，関数型とパラメータの組み合わせに応じた数通りの Yoshimoto 型減反率と既存の伐採齢データから推定した従来型減反率との比較，さらには伐採予測への応用を試みた．結論として，本論の手法は，ただちに伐採予測に利用可能なものではないが，過去のデータに基づく精緻化を経ることにより利用可能となることが見込まれた．
キーワード：減反率，伐採予測，ポアソン過程，待ち時間，丸太価格．
1.

はじめに

1.1 減反率
減反率とは，新植された林分がある齢級で伐採される確率を指し，林野庁により我が国の伐採予測に用いられている代表的な手法である．減反率の発想の原点は，林分の寿命と α 粒子の飛跡の類似性に着目した鈴木（1959）に見ることが出来る．その後，鈴木（1972）は，ある林分が植栽されてから伐採されるまでの林木の成長過程をポアソン過程とし，伐採までの待ち時間の確率密度関数
Fk (t) = me−mt(1.1)(mt)k−1,(k − 1)!

を導いた．すなわち j 齢級の減反率 qj は

 j+1

(1.2)

qj =となる．また鈴木（1972）は，

∞

E(t) =0(1.4)

E(t2 ) =0†

∞

e−mt

0t2 Fk (t)dt =me−mt

jtFk (t)dt =

∞

 j+1

Fk (t)dt =j(1.3)m ∈ Ê > 0,

1m(k − 1)!

k∈Æ(mt)k−1dt(k − 1)!

Γ(k + 1)k(mt)kdt ==(k − 1)!
m(k − 1)!
m

∞0(mt)k+1 e−mt dt =k(k + 1)m2

東京大学大学院 農学生命科学研究科：〒113–8657 東京都文京区弥生 1–1–1; hiroshim@fr.a.u-tokyo.ac.jp

112

統計数理

第 51 巻

第 1 号 2003より，(1.5)

σ 2 (t) = E(t2 ) − E(t)2 =km2となるので，パラメータ m，k は平均伐採齢 E(t)，伐採齢分散 σ 2 (t) から求めることができるとした．すなわち鈴木は，伐採齢を確率変数と見なした場合の寿命分布がガンマ分布 Ga(k, m)に従うことを導き，観測された伐採齢の平均と分散から分布型を決定しようと考えたのである．
なお近年，鈴木（2002）は，林分はある林齢 b までは伐採されないことから密度関数を b だけ平行移動した
 Fk (t) = 0 ,

(1.6)

t≤b
k−1{m(t − b)} Fk (t) = me−m(t−b)(k − 1)!

,t>bなる分布の 3 つのパラメータを伐採齢の平均，分散，歪度から求めることを提案した．なお j
齢級の伐採率 cj は
qj
(1.7)
cj =j−11−qii=1と計算され，減反率 qj と区別される．
一方，Yoshimoto（1996a）は，非定常なポアソン過程に基づく減反率の確率密度関数(1.8)

Fk (t) = m(t)e−g(t)g(t)k−1(k − 1)!

tただし， g(t) =m(τ )dτ0を導いた．これは，ポアソン過程で任意の時間 [0, t] における事象の発生回数（以下，カウント数）が微少区間 ∆t に増加する割合 m を m(t) に置き換え，それを時間積分した g(t) を物理量ととらえることにより減反率分布を求めるものである．たとえば Yoshimoto（1996b）は，g(t) に林分の現在価値を用いて木材価格，コスト，利率等の経済要因が変化した際の減反率の応答を分析した．ただし，Yoshimoto（1996a, 1996b）は（1.8）の確率密度関数では無限積分が必ずしも1 にならないことを指摘した．そこで Yoshimoto（2001a）は，増加率が Mitscherlich 型の成長関数に従う場合，非定常ポアソン過程のカウント数に関する確率の差分微分方程式を状態変数k について解くことによりこれは二項過程となることを示し，この場合，無限積分が必ず 1 になることを導いた．さらに Yoshimoto（2001b）は，成長関数として Logistic 式，Gompertz 式，Richards 式を用いた減反率を導いた．例として Richards 式に基づく減反率の確率密度関数を以下に示す．
(1.9)

Fk (t) = M f (t)e−g(t) M −1 Ck−1 (1 − e−g(t) )k−1 (e−g(t) )M −kabe−at (1 − e−at )b−1
1 − (1 − e−at )b
1g(t) = ln
1 − (1 − e−at )b

ただし， f (t) =a, b ∈ Ê > 0 ,M ∈Æ

これら Yoshimoto の減反率分布を推定するには観測値をもとに成長関数のパラメータを回帰分析等で推定し，さらにパラメータ k に適当な値を与える必要がある．

丸太価格に基づく減反率の推定

113

1.2 本論の目的
前述したように，（1.1）
を推定する際には，伐採齢の平均と分散の情報が必要となるが，そのような情報は減反率理論を適用する地域が広域になるほど入手が困難となる．そこで本論では，この理論を全国レベルに適用することを前提として，データの入手が容易な丸太価格に基づき（1.8），（1.9）を推定することを目的とし，減反率理論の適用事例を示す．
以下では，はじめに径級別丸太価格および径級別立木本数割合を利用して丸太価格 林齢の（，1.9）が前提とする 2 種類の成長関データセットを構築し，つぎにそれらデータに対して（1.8）
数を最小二乗法によりあてはめる．さらに 3 通りの k の値にもとづき（1.8）（，1.9）の減反率分布を推定し，分布（1.1）との比較を行う．比較の結果（1.1）に最も近い分布を示したものについて，最後に伐採面積を試算する．
2.

研究の方法

2.1 丸太価格 林齢曲線
一般に，丸太価格は末口径に応じて変化し，稲田（2002）によれば，20 cm 前後の適寸径でピークを持ちそれ以外では末口径の二乗に比例して増加するという．本論ではこのような特性を考慮して，丸太価格を径級別に扱うこととした．
さらに
（1.8）
（，1.9）を推定するには丸太価格の林齢に対する時系列変化，すなわち丸太価格 林齢曲線が必要となる．これまで丸太価格 末口径曲線に関する研究は上述の稲田（2002）の他に大洞 他（2001）などいくつか見られるが，径級別の丸太価格に着目して，丸太価格 林齢曲線を分析したものは見あたらない．
ある林分の価値を丸太価格で評価するには，林分の直径分布をもとに各径級に属する立木の本数割合を求め，それらを重みとする径級別丸太価格の加重平均を求めればよい．これを齢級別に求めることにより丸太価格の時系列変化が得られる．すなわち，第 j 齢級の丸太価格 Pj
は，(2.1)

Pj =s

si,j pi ,j = 1, 2, . . . , ri=1

ただし，pi ：第 i 径級の丸太価格，si,j ：第 i 径級，第 j 齢級の立木本数割合，r ：齢級数，s：径級数．丸太価格 林齢曲線は，Pj を観測値として適当な曲線を当てはめることにより求めることができる．曲線は，（1.8）に関しては任意であるが，上述の通り確率密度関数の無限積分が必ずしも 1 にならないという問題がある．これを避けるためには（1.8）
で曲線 g(t) が，g(0) = 0かつ limt→∞ g(t) = ∞ なる条件をみたせばよい．このような条件をみたす曲線として，n 次の多項式(2.2)y(t) =nai ti ,

an > 0i=1などが考えられる．
（1.9）に関しては，Richards 式(2.3)

y(t) = M (1 − e−at )b ,M, a, b ∈ Ê > 0に限定される．
2.2 パラメータ k の決定
定常なポアソン過程におけるカウント数は，その値が i − 1 から i へ増加する際 (i ≥ 2) の待ち時間が各々，独立に指数分布に従うステップ関数であり，減反率はカウント数が k となるま
2での待ち時間に関する伐採齢の確率分布である．ここで（1.3），（1.5）より k = E(t)と表される
σ 2 (t)

114

第 51 巻統計数理

第 1 号 2003

ことから，
（1.1）の k は無次元量であることがわかる．鈴木（2002）は，この k に物理的な意味はなく，その値を伐採行為と関連づけて観測することはできないとしている．
∞∞
一方，（1.8），（1.9）では，E(t) = 0 tFk (t)dt や E(t2 ) = 0 t2 Fk (t)dt ，さらにはモーメント
 ∞ xt母関数 MX (x) = 0 e Fk (t)dt の一般型を導くことはできないため，k は伐採齢の平均，分散，歪度などから求められない．Yoshimoto（1996b, 2001b）は，（1.8）（，1.9）の k を，森林所有者の伐採行動の意思決定に基づき（1.8）の g(t) や（2.3）の成長関数に関連して決定される，何かしらの物理量であるとしている．そこで本論もこの考えを踏襲し，k は，丸太価格曲線に関連して，径級別の価格に一定の重みをかけた加重平均として決定することとした．すなわち
(2.4)

k=s

wi pii=1

ただし，wi ：第 i 径級の重み．wi の決定方法については以下の 3 つのケースを想定した．
ケース 1. wi は，各径級に属する立木の現存量に比例する．すなわち1gi,j ,r j=1r

(2.5)

wi =

i = 1, 2, . . . , s

ただし gi,j ：第 i 径級における第 j 齢級の立木本数割合．
ケース 2. wi は，各径級に属する立木の伐採量に比例する．すなわち

r
j=1 fj gi,jwi = r

(2.6)

j=1 fjただし，fj ：第 j 齢級の伐採率．ここでは伐採量は立木本数割合と伐採率の積に比例するものと仮定している．
ケース 3. wi は，各径級に属する丸太の市場出荷額に比例する．
この場合の重みは市場公表値をそのまま用いることとした．
3.

結果と考察

3.1 丸太価格 林齢曲線使用したデータは，農林水産省統計情報部が発行する「統計情報
木材価格」のスギ径級別全国価格と，全国人工林の齢級別直径分布データである．丸太価格に関しては 2002 年 5 月時
点で m3 あたり，スギ小丸太：8〜14 cm が 10,800 円，スギ中丸太 I：14〜24 cm が 14,060 円，ス
ギ中丸太 II：24〜30 cm が 15,820 円，スギ大丸太：30 cm〜が 19,940 円となっており，立木本数割合に関してもこの 4 つの径級別に求めた．
（2.1）にもとづく計算結果を表 1 に示す．表 1 最右列の丸太価格は齢級とともに増加し 15,000 円前後で頭打ちとなるシグモイド型の変化を示した．このことから丸太価格 林齢曲線は（2.2）（，2.3）の曲線でうまく近似できる．曲線のあてはめにあたっては，（2.2），（2.3）とも原点を通る曲線のため，観測値全体から初期値を引き，単位を 100 円／ m3 とするなどの換算をした．最小二乗法によるあてはめの結果，（2.2）
では，3 次多項式の a1 は 3.062，a2 は −0.068，a3 は 0.0005 となり，
（2.3）では M は 44.794，a は 0.113，bは 1.519 となった．結果を図 1 に示す．なお図 1 の縦軸の値は元の価格に再換算した．
3.2 k の値
3 つのケースの径級別重みと k の値を求めた結果を表 2 に示す．なおケース 2 の伐採率は，
森林基本計画研究会（1997）の林産物需給の長期見通しにおけるすう勢見通しシナリオに示され丸太価格に基づく減反率の推定

115

表 1. 径級別素材割合と丸太評価価格．

図 1.

丸太価格 林齢曲線．

た 2002 年時の全国人工民有林，伐採齢平均と分散に基づいて求めた．ケース 3 の重みは，日本銀行公表の国内卸売物価指数算定時の重みをそのまま用いた．ケース 1 および 2 の k は，丸太価格 林齢曲線で前提とされる資源の現存量に基づいて求めたものであるのに対し，ケース

116

統計数理
表 2.

第 1 号 2003

各ケースにおける k の値．

図 2.

表 3.

第 51 巻減反率分布の比較．

各モデルの平均と分散．

3 の k はそのようなメカニズムを無視して求めたものである．
ケース 1 および 2 の k は図 1 の収束値 15,000 円以下のため（1.9）の確率分布を求めることができるが，平均市場取引価格にあたるケース 3 の k は収束値を越えてしまうため確率分布を求めることができない．
3.3 減反率分布
以下の分析ではケース 3 を除外し，（1.8）のケース 1，（1.8）のケース 2，（1.9）のケース 1，（1.9）のケース 2 をそれぞれ丸太価格モデル A，B，C，D と呼ぶ．さらに比較のため，ケース 2 で用いた伐採齢データに基づきパラメータの値を決定した（1.1）を取り上げ，これを伐採齢モデルと呼ぶ．図 2 に丸太価格モデルおよび伐採齢モデルの減反率分布を，表 3 にこれら分布の平均と分散をそれぞれ示す．
関数型の違いを比較すると，（1.8）に従う丸太価格モデル A，B では，平均値が伐採齢モデルに近くなると分散も伐採齢モデルに近くなったが，平均値が大きくなると分散が極端に大きくなり伐採齢モデルと乖離した．
（1.9）に従う丸太価格モデル C，D では，平均値が伐採齢モデルに近くなると分散が極端に小さくなり伐採齢モデルと乖離した．すなわち，正確に平均伐採齢を再現する k の値を決定できるのであれば（1.8）を選択するべきである．
ケースの違いを比較すると，ケース 1 に従う丸太価格モデル A，C では，平均値が伐採齢モデルに近くなったのに対し，ケース 2 に従う丸太価格モデル B，D では過大となった．k を決定する際に，ケース 2 は伐採率という形で伐採齢の情報を利用するのに対しケース 1 は利用し丸太価格に基づく減反率の推定
表 4.

117

丸太価格モデルの伐採齢モデルに対する適合度．ていないことを考えれば，ケース 1 の方が実践的である．すなわちケース 1 は，事前に伐採齢の平均と分散といった，分布を規定する情報が無くとも，前提とする成長関数に関する情報のみで正確な分布を再現できることを示唆している．ここで再度ケース 1 における k の算定式
（2.5）を見てみると，これは丸太価格 林齢曲線の観測値の平均値，すなわち全国人工民有林の平均丸太価格を求めているに他ならない．ポアソン過程のカウント数 k がそのような平均値と一致するときに伐採が実現されれば，そこから導かれる減反率分布は，現実の伐採齢に基づく減反率分布とよく対応するものと推測される．
総じて，（1.8）の k をケース 1 の手法で決定した丸太価格モデル A が伐採齢モデルに最も近い分布となった．
3.4 伐採予測への応用
最後に丸太価格モデル A の伐採予測への適用性を調べるため，伐採齢モデルと種々の比較を行った．表 4 は，丸太価格モデル A の伐採齢モデルに対する適合度を検定するために，丸太価格モデル A の減反率分布に従う 1000 個の乱数を 10 回発生させた際の平均度数と，伐採齢モデルの減反率分布を理論確率とした期待度数を示したものである．適合度検定の結果，自由度22 で χ2 値が 167 となり，両者には 1％水準で有意な差が認められた．この差は主としてモデルの前提となる伐採齢データの精度に起因すると考えられ，丸太価格と伐採齢に関して観測年のそろった過去のデータを利用することにより，減少すると考えられる．

118

統計数理

第 51 巻

第 1 号 2003

図 3.

伐採率の比較．

図 4.

伐採面積の比較．

つぎに（1.7）に従い伐採率を計算した結果を図 3 に示す．ただし伐採率の計算にあたっては，丸太価格モデル A については 33 齢級以降で，伐採齢モデルについては 26 齢級以降でそれぞれ減反率が 0.001 以下となったので，それらの齢級ですべて伐採されるものとした．丸太価格モデル A の伐採率は，17 齢級までは増加傾向で伐採齢モデルとほぼ同じ値となったが，それ以降は齢級の増加とともにいったん減少傾向となり，依然，増加傾向の伐採齢モデルと乖離した．この差は，図 2 の分布の違いに起因するものである．現実には 20 齢級前後の高齢級の林分は何らかの理由で伐採されず伐採率が低くなるため，丸太価格モデル A の方が現実の伐採性向に合うと考えられる．
これら伐採率に基づき全国人工民有林の伐採面積を試算した結果を図 4 に示す．全国約 800万 ha の人工民有林に対し，丸太価格モデル A に基づく総伐採面積は約 64 万 ha，伐採率モデルは約 83 万 ha となった．現実には 20 齢級以上の高齢林がほとんど存在しないため，図 3 に見られる伐採率の差は伐採面積へはさほど反映されなかった．
以上の比較から鑑みて，丸太価格モデル A は，ただちに伐採予測に利用可能なものではないが，過去のデータに基づく精緻化を経ることにより十分利用可能なものとなると言える．すな
わち丸太価格に基づき減反率の推定，ひいては伐採予測が可能となる．
4.

おわりに

本論では Yoshimoto（1996a, 1996b, 2001b）の導いた減反率理論を用いて，従来のように伐採

丸太価格に基づく減反率の推定

119

齢の平均・分散を用いることなく丸太価格から全国レベルの減反率分布を求める方法を検討した．具体的には丸太価格 林齢曲線という新しい概念に基づく 2 通りの成長曲線を推定し，さらに減反率のパラメータについて 3 通りの推定方法を提案し，それらの組み合わせの中から現実の伐採性向に近い減反率分布を再現できるものを見出した．本論はある種の事例研究ではあるが，Yoshimoto の減反率を実践的な問題へ適用する際の礎になるものと考える．
今後の課題としては，伐採齢の平均・分散と径級別丸太価格に関して，観測時点のそろった過去の実績値を用いて，本論の手法で同様に減反率分布を再現できるか検証すること，その上で径級別丸太価格の予測モデルと組み合わせて，丸太価格の変動から減反率をダイナミックに予測するシステムを構築することなどが挙げられる．

参 考 文 献稲田充男
（2002）. 素材価格曲線の検討，『森林資源管理と数理モデル
21 世紀ニューミレニアムに向けて（吉本敦，松村直人，近藤洋史 編）』，第 1 巻，59–72，森林計画学会出版局，東京．
大洞智宏，渡邉仁志，横井秀一（2001）. 岐阜県東農地域のヒノキ林における長伐期施業導入の条件
木材市場価格調査からの考察，日林学術講，112，p. 153．
森林基本計画研究会 編（1997）. 『21 世紀を展望した森林・林業の長期ビジョン
持続可能な森林経営の推進』，地球社，東京．
鈴木太七（1959）. 遷移確率行列による収穫予定，日林関東支論，10，36–38．
鈴木太七（1972）. 林業における確率過程論の応用（I），日林誌，54，234–243．
鈴木太七（2002）. 減反率の推定について，『森林資源管理と数理モデル
21 世紀ニューミレニアムに向けて』（吉本敦，松村直人，近藤洋史 編），第 1 巻，1–27，森林計画学会出版局，東京．
Yoshimoto, A.（1996a）. A new stochastic model for harvesting behavior with application to nonstationary forest growth and supply, Canadian J. For. Res., 26, 1967–1972.
Yoshimoto, A.（1996b）. Economic analysis of harvesting behavior using the modified Gentan probability theory, J. For. Res., 1, 67–72.
. Gentan probability analysis with a state-dependent discrete forest growing
Yoshimoto, A.（2001a）model, J. For. Res., 6, 101–110.
Yoshimoto, A.（2001b）. Application of the logistic, Gompertz, and Richards growth functions toGentan probability analysis, J. For. Res., 6, 265–272.

120

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

Estimation of Gentan Probability Based on Price of Logs
Takuya Hiroshima(Graduate School of Agricultural and Life Science, The University of Tokyo)

Gentan probability is defined as the probability that a newly planted forest stand will
be harvested after a certain period. This probability is derived from the Poisson process
and the waiting time theory, and has been applied to yield prediction in Japan. For
estimating the probability, the mean and variance of cutting age are required, although
it is hard to obtain such a data in practice. In this paper, the analysis of Yoshimoto’s
Gentan probability is conducted using the price function of logs. First, two types of timedependent log price curves are estimated by the OLS method for the growth functionof Yoshimoto’s Gentan probability. Next, parameters of the Gentan probability density
function are determined on the basis of three hypothetical cases. Finally, six cases of theGentan probability distribution are derived and are compared with the actual distribution
derived from the observed mean and variance of the cutting age. The result is that one
case of the distributions of Yoshimoto’s Gentan probability is almost consistent with theactual one. This implies that the Gentan probability based on the price function of logs
may be useful for yield prediction.

Key words: Gentan probability, log price, Poisson process, waiting time, yield estimation.