統計数理（2023）
第 71 巻 第 2 号 129–148©2023 統計数理研究所

特集「諸科学における統計数理モデリングの拡がりII」［研究ノート］

果樹豊凶の決定論的非線形予測——アンサンブル再構成と Lorenz 類推法による１年先収量予測——
酒井 憲司†（受付 2023 年 3 月 9 日；改訂 7 月 18 日；採択 7 月 20 日）

要旨決定論的カオスは軌道不安定性を有する低次元のダイナミクスから生成される．観測時系列からカオス性を判別することが非線形時系列解析（カオス時系列解析）の第 1 の目的である．その中でも，決定論的非線形予測はよく利用される手法の一つであるが，大きな時系列サイズが必要である．一方，作物の収量は年に 1 回しか計測できないため，取得できる収量時系列の大きさは非常に小さい．その反面，圃場や果樹園の作物個体数は大きいため，時系列集合のサイズは大きくなる．また，多くの果樹の繁殖様式が低次元の非線形ダイナミクスを有していると仮定できる．本稿では，果樹園から得られた収量時系列集合に対して行ったダイナミクスのアンサンブル再構成，決定論的非線形予測そして収量の 1 年先予測の実践例をウンシュウミカン
とピスタチオについて紹介する．
キーワード：カオス同期，カオス集団力学，決定論的非線形予測，隔年結果，アンサ
ンブル再構成，共通ノイズ同期．
1.

はじめに

果樹や森林樹木の果実および堅果（ドングリ）の生産量は個体ごとに大きく変動する．さらに，
これらの樹木個体群は種子・果実生産において豊凶同期現象を呈する．ある山全体のドングリが大豊作であると翌年は実をつけない．豊凶の周期は 2 年から数年，同期の空間的範囲は数km から数百 km まで樹種によってマスティングと呼ばれる多彩な同期現象を示す（Koenig and．柑橘類，ナッツ類，広葉果樹では豊凶が隔年交互に現れ（Sakai et al., 2008），
Knopes, 2005）
これを隔年結果と呼ぶ．樹木集団の同期メカニズムとしては，花粉による間接結合，共通ノイズそして生理融合による直接結合の 3 者が知られており，樹種と生育環境によって，これらがは物質収支モデル（Resource Budget Model，以下 RBM）によって複合する．Isagi et al.（1997）
個体の豊凶メカニズムを一次元写像のカオス振動子としてモデル化し，個体群の豊凶同期（マ
スティング）
を花粉結合による大域結合写像として定式化した．樹木豊凶同期は，個体間の相互作用の形態，RBM の制御パラメータおよび共通ノイズによって，同期強度・周期・空間分
布などに多彩な挙動を呈する．植物種を超えた多様な同期現象を考察するにあたって，RBMは強力な理論的基盤を与えてくれている（Satake and Iwasa, 2000）．さらに，RBM パラメータの同定についても様々なアプローチが精力的に試みられている（Rees et al., 2002; Crone et al.,2005; 佐竹, 2007; Abe et al., 2016）．また，農業や林業への応用を目指したカオス農学の立場†

東京農工大学 農学研究院：〒 183–8509 東京都府中市幸町 3–5–8

統計数理

130

第 71 巻 第 2 号 2023

から（酒井, 1997; Sakai, 2001），ウンシュウミカン（Sakai et al., 2008），ケヤキ（Sakai et al.,，ピスタチオ（Sakai et al., 2022），スギ・ヒノキ（Ishibashi et al., 2019）などの多様な同
2019）期のメカニズムが，カオス解析・同期解析および RBM に基づいた数理モデルによって解明さ．
れてきた（酒井, 2023）
果樹の豊凶データは一個体ごとに年 1 点しか得られないため，その時系列サイズは極めて小さい．しかし，幸いなことに，果樹園や林分内の個体数は数十から数万にも及ぶため，時系列
集合のサイズを十分大きくすることが可能である．また，RBM によって多くの果樹の繁殖様式が低次元の非線形ダイナミクスを有していると仮定できることから，時系列集合からのダイ
ナミクス再構成が可能となる．本稿では，最もシンプルな決定論的非線形予測法である Lorenz
類推法によって，数点という極小サイズの時系列に対して，統計的手法に拠らず，決定論的非線形予測を応用した 1 年先収量予測の実践例を示す．
2.

樹木豊凶の決定論的ダイナミクス

2.1 物質収支モデル（RBM）としてのテントマップ
Isagi et al.（1997）の RBM の模式図を図 1a に示した．樹体内には，なんらかの物質 S を貯蔵するプールがある．環境が一定であれば毎年 PS がプールに供給される．S + PS がプールの容量 LT を超えると，あふれた Cf によって開花数が定まる（図 1a 上段）．LT を超えなければ着花しない（図 1a 下段）．開花すると果実の結実のために Ca が消費される．この時，Ca /Cf– 2.3）
をコスト比 RC と呼ぶ．これを個体 i のサイズ M の個体群に拡張して
（2.1）
（
に定式化し
た．ここで，i は個体番号で i = 1, 2, . . . , M である．
LT − Cai (t)Si (t) + PS > LT
(2.1)
Si (t + 1) =Si (t) + PS
Si (t) + PS ≤ LT

但し，
(2.2)

Cai (t) = RC Cfi (t)(2.3)

Cfi (t) = Si (t) + PS − LT .

このダイナミクス Si (t)→Si (t + 1) は図 1b のリターンマップに示されるように 1 次元写像である．これは区分線形のテントマップで左側と右側の傾きは，それぞれ 1 と −RC である．RCを制御パラメータとする S の分岐図を図 1d に示した．一方，現実には S を観測することはで
きず，Ca を樹木個体の果実生産量として観測する．Ca のリターンマップと分岐図を図 1c と
図 1e に示した．Ca リターンマップから明らかなように，Ca (t) = 0 の場合には Ca (t + 1) は一．
意に定まらない（図 1c）
2.2 同期モデル
2.2.1 花粉結合モデル–（2.3）の RBM モデルでは個体間の相互作用がなく個体はそれぞれに独立に振る舞うた（2.1）
め同期は起こらない．自家不和合性の樹種は他の個体の花粉でないと受精できないため，開は花粉結合による平均場を（2.4）で構成し，（2.2）を（2.5）
花・結実に至らない．Isagi et al.（1997）
で置き換えることによって同期モデルとして定式化した．Y (t) は年ごとの集団全体の受粉成功率であり，β が大きいと受粉成功率は低くなり，同期は強くなる．

β
M

1
(2.4)
Y (t) =
Cfi (t)(t)(M − 1)P0
j=1,j=i

果樹豊凶の決定論的非線形予測131

図 1．物質収支モデル（RBM）の構造．a. 井鷺の物資収支モデル（RBM）の模式図，b. S のリ
，c. Ca のリターンマップ（RC = 1.4，
ターンマップ（RC = 1.4，P0 = 10，LT = 100）
P0 = 10，LT = 100），d. S の分岐図，e. Ca の分岐図．

(2.5)

Cai (t) = RC Cfi (t)Y (t)–（2.4）の P0 は秋にプールに供給される物質の内的ポテンシャルであり，（2.1）（2.3）式の RBM
では供給量は毎年一定としているため，（2.1）
および（2.3）式において PS = P0 である．
花粉の到達距離内での同期は大域結合写像（Kaneko, 1990）として，それより広い範囲では
局所結合写像としてモデル化されることになる．β が大きくなり完全同相同期になると，接線分岐と境界衝突分岐を交互に繰り返す周期増加分岐となる（Sakai et al., 2019）．周期増加分岐
では制御パラメータである RC の増加に従って卓越する周期が 1 つずつ増加する．これに対して，周期倍分岐では周期は偶数毎に増加する．隔年結果やマスティングでは奇数年の周期も現れるため，RBM が周期増加分岐のダイナミクスを内包していることの重要性を強調しておきたい．なお，不連続点を有しない単峰型の 1 次元写像とした場合には（Ye and Sakai, 2016;，周期増加分岐ではなく周期倍分岐となる．
Esmaeili et al., 2021）
2.2.2 共通ノイズモデル
ピスタチオなどの雌雄異株の樹種では，雌株に対して十分な花粉を提供できるように果樹園
内に雄株を配置する．そのため，花粉結合による同期メカニズムを仮定できないので，ピスタ
チオの同期は共通ノイズ同期であると考えてきた（Lyles et al., 2009; Rosenstock et al., 2010;
Esmaeili et al., 2022; Sakai et al., 2022）．微気象変動などの外乱が果樹園内の各個体に同一に

132

統計数理

第 71 巻 第 2 号 2023

作用すると仮定し（2.6）
によって共通ノイズ PS (t) を導入する．PS (t) は t 年毎に変動するが，常にすべての個体に同一に作用する．
–（2.3）の PS を（2.6）
式で置き換えて RBM の共通ノイズ同期モデルを定式化した（Sakai（2.1）．
et al., 2022）(2.6)

PS (t) = P0 {1 + eC σ(t)}ここで，σ(t) は平均値 0，標準偏差 1 の正規乱数であり，eC が共通ノイズの強度である．
共通ノイズが印加された場合のリターンマップの幾何学的形状は図 1b のテント写像と変わらず，t 年毎の共通ノイズの大きさ eC σ(t) に応じて左右に平行移動する（Sakai et al., 2022）．
2.2.3 拡散結合モデル
ピスタチオ果樹園における収量の空間分布に 1/3 べき乗則をともなう顕著な長距離相関の存，統計物理の Ising モデルとのアナロジーから，近在が明らかになり（Nobel et al., 2015, 2018）
接する個体間における相互作用の存在が示唆された．そこで筆者らは，根系ネットワークによる個体間の物質交換（Klein et al., 2016; Lev-Yadun, 2011; Simard et al., 1997; Gaion et al.,を仮定して，（2.7）
によって，距離 γ = 11 m にある近傍個体間の拡散結合による結合写像
2018）として定式化した（Sakai et al., 2022）．
格子（Kaneko, 1985）(2.7)

SAi (t) = Si (t) +1
N E(i)Mεj,i [Sj (t) − Si (t)]i=jd(i,j)<γ

ただし，εj,i を樹木 i, j 間の拡散結合，d(i, j) を樹木 i, j 間の距離，N E(i) を樹木 i の γ 近傍の樹木の数，SAi (t) を着花前の i 個体内の貯蔵物質としている．（2.3）を（2.8）で置き換え（2.9）が得られる．
SAi (t) + PS (t) − LT
SAi (t) + PS (t) > LT
(2.8)
Cfi (t) =0
SAi (t) + PS (t) ≤ LT
(2.9)

Si (t + 1) = SAi (t) + PS (t) − Cfi (t) − Cai (t)ここで示した拡散結合を有する RBM のリターンマップは，結合に起因する揺らぎがテント–（2.9）を追加すること
マップ上に現れるが，骨格構造は変わらない（Sakai et al., 2022）．（2.7）
で，カオス振動子 RBM の結合写像格子（Kaneko, 1985）の共通ノイズ同期として定式化した．
–（2.9）の結合形式では同相同期となるが，拡散結合がおこる生活史上のタイミング
なお，
（2.7）
が異なれば逆相同期にもなる（Prasad et al., 2017）．
2.2.4 RBM 制御パラメータの空間勾配
RBM の制御パラメータは RC であるが，ピスタチオにおいては受粉成功率 Z の空間分布を考慮する必要があった（Sakai et al., 2022）．果樹園の西端からの距離 l における受粉成功率
Z(l) と，新たな制御パラメータとしての作物係数（Cropping coeﬃcient）
mC を以下のように定める．
(2.10)

mC = RC Z(l)（2.2）式を（2.11）
式に置き換え，RBM 制御パラメータ mC を定式化した．
(2.11)

なお，

Cai (t) = mC Cfi (t)果樹豊凶の決定論的非線形予測Z(l) =(2.12)

1
l
+α
RC
WE

α=1−

(2.13)

133



1
RC

（2.12）中の WE は果樹園の東西距離（777 m）
である．空間勾配を α で表し，西端では mC = 1，
（2.12）と（2.13）を定めた．
東端では mC = RC となるように
3.

決定論的非線形予測

3.1 決定論的非線形時系列解析とその農業技術上の意義
決定論的カオス（以下，カオス）は低次元の決定論的ダイナミクスによって生成される一見複雑で不規則な変動である．カオスの軌道不安定性（初期値への鋭敏な依存性）
によって誤差が指数関数的に拡大するため，初期値への鋭敏な依存性を示し，往々にして実験における再現性が困難となる．変動の主要部分がカオスであれば，短期予測や制御の可能性があるので，観測時系列からカオスか否かを判別することが重要になる．その際，カオスの属性である決定性と軌
道不安定性を利用するのであるが，カオスは十分条件が定義できないため，決定論的非線形予測，Lyapunov 指数，相関次元，推移誤差，Wayland テスト，False Nearest Neighbor 法，サロゲート法などからカオスの必要条件を多面的に調べ，現象に関する専門的知見と照合しながら，カオス性の判定をおこなう（池口 他, 2000）．これを非線形時系列解析もしくはカオス時と呼んでいる．農業現場は不規則で複雑な変動や再現性が困難な系列解析（池口・合原, 1997）
現象に囲まれている．にもかかわらず，伝承・蓄積されてきた農業技術や農業者の経験と技によってこれらをうまく取り扱いながら作物生産を行ってきた．変動現象にカオスが関与していれば，予測手法や制御手法の導入可能性が生まれ，農業技術への寄与が期待できる．これが，農業におけるカオス理論の意義でありカオス農学の目的でもある．
本稿では果樹収量の 1 年先予測手法の実践を紹介するが，まず 3 章では 1 本の観測時系列からのダイナミクスの再構成法と Lorenz 類推法による決定論的非線形予測手法を紹介した上で，4 章でそれを果樹収量の将来予測手法に拡張する．
3.2 時間遅れ埋込によるダイナミクスの再構成
農業現場で見られる複雑な変動現象の微分方程式や差分方程式が不明であることは一般である．しかし，「時間遅れ埋込法」が発明されたことで，1 つの観測時系列のみを用いて多次元のダイナミクスを再構成することが可能となり，カオス理論の現場応用に扉が開かれた．Takensの埋め込み定理は，得られた観測時系列 x(t) からその時間遅れを成分とした m 次元ベクトルx(t) = {x(t), x(t + τ ), . . . , x(t + (m − 1)τ ), t = 1, 2, . . . , N } によって描いた相軌道が，元のダイナミクスとトポロジカルに等価なダイナミクスを再構成していることを保証している（池口・．実際の解析においては，時間遅れ τ は卓越する周期の 1/4 や自己相関が最初にゼ
合原, 1977）
ロになる時刻などが経験的に用いられている．埋め込みの次元 m をダイナミクスの次元より大きく設定しないと誤った最近傍ベクトルを選択するため予測精度が大きく低下する．
3.3 Lorenz 類推法による決定論的非線形予測{x(t), 1 ≤ t ≤ N } が与えられたとき，m 次元の相空間に予測器と予測対象をそれぞれ（3.1）と（3.2）で構成する．
(3.1)

X (i) = {x(i), x(i + τ ), . . . , x(i + (m − 1)τ ), 1 ≤ i ≤ K}

統計数理

134

第 71 巻 第 2 号 2023

X (j) = {x(j), x(j + τ ), . . . , x(j + (m − 1)τ ), K ≤ j ≤ N − p − (m − 1)τ }(3.2)

K = N/2 − (m − 1)τ であれば時系列の前半分と後半分で予測器と予測対象になる．ここで，（3.1）から探索する．そうすると，予測（3.2）の予測対象 X (j) の最近傍ベクトル X (k) を予測器
対象 X (j) の p ステップ先の予測ベクトルは X̂ (j + p) = X (k + p) として算出される．また，
j + p における実ベクトルは X (j + p) である．予測ベクトルと実ベクトルの第 m 成分はそれぞれ，{x̂(j +p+(m−1)τ )} および {x(j +p+(m−1)τ )} であり，これらを K ≤ j ≤ N −p−(m−1)τ
に対して求めると，予測区間 [K + p + (m − 1)τ, N ] における予測値が得られたことになる．予測ステップ p の関数として x と x̂ の相関係数 r(p) や正規化平均自乗誤差（Relative Root Mean
RRMSE（p）などを求める．
Square of Error）

N −p−(m−1)τ
{x(j + p + (m − 1)τ ) − x̂(j + p + (m − 1)τ )}2
j=K
(3.3)
RRMSE (p) =N −p−(m−1)τ
{x(j + p + (m − 1)τ ) − x}2
j=K

但し，
(3.4)

x=

1
N − p − (m − 1)τ − K + 1

N −p−(m−1)τ



x(j + p + (m − 1)τ )j=K

ここで，RRMSE (p) = 0 は完全な予測を示し，RRMSE (p) = 1 は平均値以外全く予測できな（カオス），カラードノイズ，
かったことを意味する．sin 波，sin 波 + ホワイトノイズ，Lorenz63（図 2b）で示した．
ホワイトノイズの決定論的非線形予測結果を相関係数（図 2a）と RRMSE (p)
そこで短期予測可能性と長期予測不能性を示し，決定性と軌道不安定性を確認したとしてカオスであることを主張するのである．例えば，sin の時系列では長期予測も完全に可能であり，逆に，一様乱数では短期予測も不能であるので当然ながらどちらもカオスではない．Lorenz63 とは，（3.5）式の Lorenz 方程式のパラメータを，σ = 10, r = 28, b = 8/3 として求めた時系列であり，典型的なカオスとして知られている（Lorenz, 1963）．
ẋ = σ(y − x)(3.5)

ẏ = −xz + rx − yż = xy − bz

Lorenz63 では短期予測は可能で，長期予測は不能である．これによってカオスの必要条件が確認できる．しかし，一様乱数を積分した有色ノイズでも短期予測可能性と長期予測不能性が示されており，このままではカオスと誤判定してしまう．これを回避するために，推移誤差，サ
ロゲート法，正規化決定論的非線形予測等（Sakai et al., 2019）の様々な手法が提案されている（池口 他, 2000）．
決定論的非線形予測では対象とする変動が定常であり，ダイナミクスのパラメータも変化しないと仮定して，決定性を検出することが目的である．そのため，予測器と予測対象の定め方（3.1）を予測対象
は多様に可能である．よって，前半と後半を逆にして，（3.2）を予測器としてとしても差し支えない．しかし，本稿では果実収量の 1 年先予測への応用を意図しているので，現在およびできるだけ近い過去のデータセットのみから未来の収量を予測しなければならない．次節ではそのことを踏まえて予測器と予測対象の設定を工夫する．
4.

果樹収量予測のための決定論的非線形予測

4.1 ダイナミクスのアンサンブル再構成
1 年生の作物（イネ，ムギなど）は収穫後には耕起されて更地になるので，年次間の収量変動

果樹豊凶の決定論的非線形予測135

図 2．決定論的非線形予測．a. 相関係数 r(p) による評価，b. RRMSE (p) による評価．

にダイナミクスの存在を仮定できないので本稿の考察対象としない．一方，果樹作物は図 1 のRBM で示したように，現在までの収量と翌年の収量にダイナミクスの存在が仮定できる．更に，果樹園の各個体はクローンとして同一のダイナミクスを有すると考えてよいので，果樹園の個体数と同じサイズの時系列集合を用いてダイナミクスの再構成（アンサンブル再構成）
が以下のように可能となる．
個体数 M の果樹園で N 年間の収量を観測した場合，時系列 xi = {xi (t), t = 1, 2, . . . , N } を要素とする集合 A = {xi , i = 1, 2, . . . , M } が得られる．3.1 節で述べた時間遅れを τ = 1 とする．
(4.1)

Xi (t, m) = {xi (t), xi (t + 1), . . . , xi (t + (m − 1)), t = 1, 2, . . . , N − m + 1}

（4.1）で示した M 個の相軌道によってダイナミクスを m 次元相空間内に (N − m + 1)×M のベクトル点で再構成する．これをダイナミクスのアンサンブル再構成と呼んでいる（Sakai et．
al., 2008）
4.2 果樹収量予測のための予測器構成法
アンサンブル再構成されたダイナミクスに対して Lorenz 類推法による決定論的非線形予測（年）であり，果樹収量の将来予測を行うため，予測器 L は予測を行を試みる．時間遅れ τ は 1
う時点での現在年とその過去のみから構成する．
予測を行う対象年 tN において p 年先予測を施して予測年 tN + p の収量を予測する．m 次元のダイナミクスであるとして，対象年における対象ベクトル Xj (tN , m) を m 年間の収量を成分として（4.2）のように構成する．
(4.2)

Xj (tN , m) = {xj (tN − m + 1), xj (tN − m + 2), . . . , xj (tN )}なお，（4.1）においてはベクトルの時間表示 t を第 1 成分の時間表示と一致させていたが，（4.2）
以降ではこれを第 m 成分の時間表示 t に一致させている．この表記の方が予測年 tN + p と対
象年 tN ，予測値および実値の関係が分かりやすくなるためである．
（4.3）で構成する．
次に，予測器 L(tN , p, m) を設定する．対象 tN における各個体の相軌道を(4.3)

xi (tN ) = {xi (1), . . . , xi (tN ), i = 1, 2, . . . , M }これを要素として，予測器 L(tN , p, m) を（4.4）
で設定する．
(4.4)

L(tN , p, m) = {xi (tN ), i = 1, 2, . . . , M }予測器 L(tN , p, m) において，対象ベクトル Xj (tN , m) の最近傍ベクトル Xk (tN,j,k , m) を（4.5）

統計数理

136

第 71 巻 第 2 号 2023

として探索する．
(4.5)

Xk (tN,j,k , m) = {xk (tN,j,k − m + 1), xk (tN,j,k − m + 2), . . . , xk (tN,j,k )}但し，m ≤ tN,j,k ≤ tN − p．p 年先予測ベクトル X̂j (tN + p, m) は次式となる．
(4.6) X̂j (tN + p, m) = Xk (tN,j,k + p, m)= {xk (tN,j,k + p − m + 1), xk (tN,j,k + p − m + 2), . . . , xk (tN,j,k + p)}（4.7）によって予測値 x̂j (tN + p) は，Xk (tN,j,k + p, m) の第 m 成分 xk (tN,j,k + p) となる．
(4.7)

x̂j (tN + p) = xk (tN,j,k + p)

なお，個体 j の tN + p 年の収量の実値は xj (tN + p) である．
この操作を M 個の個体に実施し，予測年 tN + p において M 個の実値 {xj (tN + p), j =1, 2, . . . , M )} と予測値 {x̂j (tN + p), j = 1, 2, . . . , M )} とから相関係数 r(tN + p, p) と正規化平均
自乗誤差 RRMSE (tN + p, p) を計算する．さらにこれを，m + p ≤ tN ≤ N − p における対象年
tN に対して実施し求めた予測値と実値によって相関係数 r(p) と RRMSE (p) を求め p 年先予測の精度を評価する．
4.3 RBM モデルを用いた共通ノイズの影響評価
決定論非線形予測では時系列に定常性を前提としている．しかし，実データにおいては厳密な定常性は期待できない．よって，実データ解析の立場からは，予測精度が目的に応じて満たされているか否かで妥当性を判断することになる．
果樹個体群の同期現象としての隔年結果は花粉結合と共通ノイズの 2 者が大きく，根系ネッ
トワークによる拡散結合もこれに加わる．共通ノイズの印加によって定常性が失われるし，大域結合写像や結合写像格子も厳密には大規模な自由度を持つため，低次元の決定論的カオスであるとする仮定を逸脱している．しかし，この仮定が実用上許容できるかどうかは，実データへの適用結果で判断すればよい．
5 章の実践例に先立って，共通ノイズの影響を RBM 時系列集合のモデルによって見積もっておく．RBM は図 1b に示すように貯蔵物質 S の 1 次元写像であるので，m = 1 として決定論的非線形予測は可能である．しかし，これまでのところ S を測定することは困難で間接的な推定が試みられているが（Rees et al., 2002; Crone et al., 2005），直接かつ簡便に観測できるのは個体の果実収量 Ca である．図 1c に示すように，Ca のリターンマップでは，Ca (t) = 0 において Ca (t + 1) が一意に定まらない．よって，少なくとも m = 2 以上としなければならない．
制御パラメータ RC = 1.9 について決定論的非線形予測を実施した．ここでは 100 個体について初期値から 180 回計算を繰り返し，181〜210 までの計算値でサイズ N = 30 の時系列集合を得た．初期値を [LT − P0 RC , LT ] の一様乱数で与えている．それぞれの個体が全くばらばらに不規則に振る舞うことに留意してほしい．わかりやすいように 100 個体の内の 4 個体を図 3aに示した．決定論的非線形予測は図 3b のように高い短期予測可能性と長期予測不能性が明瞭（p = 1）における予測値
（中空）
と実値
（中実）
を t = 3 − 30 において比である．図 3c に短期予測較した．ほとんどの予測値は実値と重なっている．図 3d に実値 vs. 予測値プロットを示した．
r(1) = 1.0 および RRMSE (1) = 0.150 とほとんど完ぺきな予測が得られている．このように，
100 個の RBM のカオス振動子は個体毎に独立に振る舞うが，全個体のダイナミクスは同一で，ノイズもないため極めて高い短期予測が実現できている．図 3e に長期予測（p = 20）における
予測値
（中空）
と実値
（中実）
を t = 21 − 30 において比較した．ほとんど予測値と実値は重ならない．図 3f に実値 vs. 予測値プロットを示した．r(20) = 0.0717 および RRMSE (1) = 1.35 と

果樹豊凶の決定論的非線形予測137

図 3．共通ノイズが無い場合の RBM 集合の決定論的非線形予測．a. RC = 1.9 における 100
個体のうち 4 個体の収量の時系列を表示（t = 1 − 30）．b. 決定論的非線形予測の精度（p = 1）における実値（中実）およを相関係数 r(p) と RRMSE (p) で表示．c. 短期予測び予測値（中空）時系列（t = 3 − 30）．100 個体中 4 個体を例示し，t = 31 では予測値（p = 1）における実値対予測値プロット（t = 3 − 30）．
（中空）のみプロット．d. 短期予測
e. 短期予測（p = 1）における実値（中実）および予測値（中空）時系列（t = 21 − 30）．100
個体中 4 個体を例示し，t = 31 − 50 では予測値（中空）のみプロット．f. 長期予測
（t = 21 − 30）．RBM のすべての個体に対し（p = 20）における実値対予測値プロットて P0 = 10，LT = 100，eC = 0 とした．

カオスの軌道不安定性によって長期予測不能となっている．なお，ここでは共通ノイズが無いにおける共通ノイズの強度は eC = 0 として計算した．
場合であるので，（2.6）
5 章で紹介するウンシュウミカンやピスタチオの実データでは，共通ノイズが外生的ノイズとして印加されているため定常性が確保できているとは言えない．そこで，RBM 集合の短期予測に対する共通ノイズの影響を調べておく．制御パラメータは RC = 1.4，（2.6）における共通ノイズの強度は eC = 0.2 として N = 30 の時系列 M = 100 個によって時系列集合を構成した．予測器は予測対象の 1 年以上過去のデータしか用いない．共通ノイズは年ごとに変動するダイナミカルノイズなので予測精度の低下は当然であるが，その程度が実用上の関心事となる．図 4a のように収量予測では相関係数 r = 0.700，RRMSE = 0.826 であり，図 4b の予測と実値の時系列においても，共通ノイズによって収量の短期予測精度が低下している．一方，
図 4c に示す収量差分の相関係数は r(1) = 0.921，RRMSE (1) = 0.438 と高い予測精度を示し,
図 4d の収量差分時系列の実値と予測値比較からも 1 年先予測もおよそ成功している．このよ（eC = 0; 図 3d）に
うに，共通ノイズ eC = 0.2 を印加した場合の収量については共通ノイズ無し
比較して予測精度は低下したが，RRMSE (1) = 0.826 と 1 よりも小さい値に収まっている．さらに，差分収量の短期予測精度は r(1) = 0.912，RRMSE (1) = 0.438 と十分高い．以上から，この程度の共通ノイズの場合には，果樹栽培の実用上許容できる予測が期待できる．本節での 1
年先予測手法を実証するために，ウンシュウミカンとピスタチオの実データに対して 5 章で実践例を解説する．

統計数理

138

第 71 巻 第 2 号 2023

図 4．共通ノイズが印加された RBM 集合の決定論的非線形予測（p = 1 の短期予測）．a. 短期（Ca ）対予測値（Ĉa ）プロット．b. 短期予測
（p = 1）に予測（p = 1）における収量の実値．100 個体中の 9 個体
おける収量の実値（中実）および予測値（中空）時系列（t = 3 − 30）（中空）のみプロット．c. 短期予測
（p = 1）における収量差分を例示．t = 31 では予測値（ΔĈa ）プロット．d. 短期予測
（p = 1）における収量差分の実値の実値（ΔCa ）対予測値（中実）および予測値時系列（中空）．100 個体中の 9 個体を例示．t = 31 では予測値（中空）のみプロット．RBM のすべての個体に対して RC = 1.4，P0 = 10，LT = 100，eC = 0.2．

5.

果樹の決定論的収量予測

5.1 ウンシュウミカン
ウンシュウミカンの観測データには個体ごとの変異や外乱，さらに比較的弱い共通ノイズ同期が想定される．ウンシュウミカンは単為結実のため，受粉する必要がなく花粉結合同期を起こすことはないが，共通ノイズと考えられる弱い同相同期を示す．
（個体毎の収穫果実数）に対し
図 5 のようにウンシュウミカンの 48 個体 7 年間の収量データて 1 年先予測を以下のように実施する．なお，1998 から 2004 年はそれぞれ t = 1 から t = 7に対応させている．1 年先予測であるから予測ステップ p = 1 で次元 m = 2 である．個体 j
において予測の対象年 tN = 3 として，予測年 tN +1 = 4 の収量の予測値 x̂j (4) を求めよう．
予測の対象ベクトル xj (3) = {xj (2), xj (3)} の最近傍ベクトル xk (2) = {xk (1), xk (2)} を予測器
L(3, 1, 2) から得たとする．
（4.5）
から tN,j,k = 2 である．
（4.7）から x̂j (4) = xk (3) となるので，
個体 j の予測年 4 の収量の予測値は x̂j (4) = xk (3) となる．即ち，2000 年（t = 3）より過去の収（t = 4）の収量予測を行ったことになる．図 5a で示したように
量データのみから 1 年先 2001 年

果樹豊凶の決定論的非線形予測139

図 5．48 個体のウンシュウミカン個体群の決定性検出 - 収量時系列集合に対する Lorenz 類推（p = 1）における収量実値対予測法の適用において最近傍 1 点のみ使用．a. 短期予測
（p = 1）における収量実値対予測値プロッ
値プロット（2001–2004）．b. 年毎の短期予測（中実）および 1 年先予測値 x̂(t)（中空）の時系列比較．横ト．c. 個体毎の収量実値 x(t)軸の t = 1 から 7 はそれぞれ 1998 年から 2004 年に対応する．

2001 から 2004 年（t = 4 − 7）まで通年の予測精度は r(1) = 0.665，RRMSE (1) = 0.814 である．
（中実）と
図 5b は 2001 年から 2004 年まで年ごとの予測精度を示した．図 5c に個体毎に実収量
1 年先予測収量（中空）を比較した．なお, 予測器 L(tN , 1, 2) が用意している探索可能ベクトル
点数は (tN − m)×M なので，2001 年から 2004 年までそれぞれ M = 48 から 4M = 192 点である．2001 年のみ RRMSE (1) = 1.10 > 1 であり収量の変動予測は失敗している．これは, 探索
可能ベクトル点数が M = 48 と最も小さいことが要因と考えられる．なお，2005 年（t = 8）の収量は測定されていないので予測値のみを中空で示してある．

140

統計数理

第 71 巻 第 2 号 2023

図 6．48 個体のウンシュウミカン個体群の決定性検出 - 収量差分時系列集合に対する Lorenz
類推法の適用において最近傍 1 点のみ使用．a. 短期予測
（p = 1）における収量差分の（p = 1）における収量差分の実値対予測値プロット（2001–2004）．b. 年毎の短期予測
実値対予測値プロット．c. 個体毎の収量差分実値 Δx(t)（ 中実）および 1 年先予測値
Δx̂(t)（ 中空）の時系列を比較．横軸の 1 から 6 はそれぞれ [2003–2001] [2001–2002][2002–2003] [2003–2004] の期間に対応する．

また，実データの決定論的非線形予測においてカオスの誤判別があった場合，差分値の予測精度が非常に低くなることが指摘されている（池口 他, 2000）．誤判別の有無を確認するために，収量差分の予測精度を図 6 に示した．ここで，収量差分 Δx(t + 1) = x(t + 1) − x(t) に対して収量差分予測値は Δx̂(t + 1) = x̂(t + 1) − x(t) で定義する．図 6a のように通年で相関
係数 r(1) = 0.883 および RRMSE (1) = 0.478 と高い予測精度が得られた．また [2000–2001]，

果樹豊凶の決定論的非線形予測141

図 7．48 個体のウンシュウミカン個体群の決定性検出 - 収量時系列集合に対する Lorenz 類推法の適用において最近傍 10 点を使用し精度向上を図った．a. 収量の予測値対実値プ
ロット，b. 収量差分の実値対予測値プロット．

[2001–2002]，[2002–2003] および [2003–2004] においてもそれぞれ予測精度は高い．
図 5 および図 6 で示したように，最近傍 1 点のみでもある程度の予測精度が得られており，
図 4 で予め調べておいた RBM 共通ノイズモデルにおける短期予測結果と合致している．このことは，ウンシュウミカンの隔年結果現象は決定論的ダイナミクス（内生要因）
が第 1 義的であって，環境要因（外生要因）
は 2 次的要因であることを示唆している．
Lorenz 類推法の改良として複数の近傍ベクトルに対して，予測対象ベクトルとの距離に応じた重み平均を用いたり，ヤコビ行列の推定などが試みられている．ここでは 10 個の近傍ベク
トルを用いて，単純に算術平均して予測ベクトルを推定し精度向上を試みた．収量と収量差分についての予測結果を図 7a および図 7b に示した．収量については，相関係数 r(1) = 0.719 お
よび RRMSE (1) = 0.698 と図 5 に比較して予測精度は向上している．また，収量差分について
も，r(1) = 0.913 および RRMSE (1) = 0.410 と図 6 に比較して精度が向上している．このように，1 年先予測など現場応用を目指す場合には，近傍点数や近傍の設定など試行して経験的に予測精度の向上を図ることになる．
5.2 ピスタチオ
ピスタチオは換金作物としてカリフォルニア州での栽培が盛んである．強い同相同期を呈し，州全体の収穫量からも 2 年周期の隔年結果が顕著に観察できる（Rosenstock et al., 2010,．ピスタチオは雌雄異株であり，雌株に十分な花粉が供給されるよう
2011; Noble et al., 2018）に雄株が一定間隔で配置されている．当該の果樹園は東西 777 m, 南北 416 m で 9,526 本の個体が 5.2 m × 6.4 m で配置されている．雄株は 26 m 間隔で配置され，すべての雌株は雄株から14 m 以内に存在し十分な花粉が供給されている．そのため，花粉結合による同期ではなく，これまで共通ノイズ同期としてモデル化してきた．ところが，収量の空間分布に強い距離相関の存在が確認された．Ising モデルによるアナロジーから，当該果樹園において近接する個体間に相互作用，例えば，根系ネットワークによる物質交換の存在が示唆された．これを受けて，2.2.3–
項で示したように，RBM に
（2.7）
（2.9）
式を追加することで，RBM の結合写像格子（Kaneko,–
における共通ノイズ同期モデルを提案した（Sakai et al., 2022）．また，（2.10）（2.13）式に
1985）
よって RBM の制御パラメータの空間勾配 α も含んだ 4 つのパラメータ (RC , eC , ε, α) 推定は，同相同期強度（Prasad et al., 2017; Sakai et al., 2019）を評価関数としたパラメータスタディに

142

統計数理

第 71 巻 第 2 号 2023

図 8．ピスタチオ実データとモデル A およびモデル B に対する 1 年先予測．実データ
収量の時系列集合とモデル A およびモデル B の時系列集合に対する Lorenz 類推（2.7）–（2.11）式で与えられ，パラ法による 1 年先収量予測を行った．モデル A は
メータ (RC , eC , ε, α) = (1.4, 0.2, 0.05, 0.2857) では空間勾配 α を含んでいる．モデ（2.7）–（2.9）式で与えられ，パラメータ (RC , eC , ε, α) = (1.4, 0.2, 0.05, 0) は
ルBは
空間勾配 α を含まない．a. 実データ収量の実値対予測値プロット，r(1) = 0.143，
RRMSE (1) = 1.15．b. 実データ収量差分の実値対予測値プロット，r(1) = 0.747，
RRMSE (1) = 0.666．c. モデル A 収量の実値対予測値プロット，r(1) = 0.377，
RRMSE (1) = 1.41．d. モデル A 収量差分の実値対予測値プロット，r(1) = 0.802，
RRMSE (1) = 0.728．e. モデル B 収量の実値対予測値プロット，r(1) = 0.760，
RRMSE (1) = 0.693．f. モデル B 収量差分の実値対予測値プロット，r(1) = 0.939，
RRMSE (1) = 0.365．実データは（Sakai et al., 2022）の Appendix で Excel ファイ
ルとして公開しており，単位は kg/樹木個体である．M = 9,526 本であるが欠損値を含んでいる．モデル A およびモデル B においては，すべての M = 9,526 の個体に対して P0 = 10，LT = 100 として数値実験を行い，モデルパラメータ同定における評価
関数は収量ではなく，同期強度の空間分布によって行った．

よった．空間勾配を考慮したモデル A：(RC , eC , ε, α) = (1.4, 0.2, 0.05, 0.2857) と考慮しないモ
デル B：(RC , eC , ε, α) = (1.4, 0.2, 0.05, 0) による数値実験によってそれぞれの時系列集合を得た．実測時系列集合とモデル A およびモデル B の時系列集合に対して，Lorenz 類推法による
1 年先収量予測を実施した（図 8）．
当該果樹園においては M = 9,526 個体の収量データ（kg/個体）を 6 年間取得したが，栽培管理上排除した個体などがあり，6 年間欠損なく観測できたのは 3,127 個体であったので，観測（本）である．一方，モデル A とモデル B の数
時系列集合のサイズは，N = 6（年），M = 3,127，M = 9,526（本）である．
（4.1）式で得られ値実験における時系列集合のサイズは，N = 10（年）た再構成ダイナミクスは，(N − m + 1)×M 個のベクトル点から成っている．
実測収量の短期予測精度は r(1) = 0.143 および RRMSE (1) = 1.15 であり予測はできていな

果樹豊凶の決定論的非線形予測

143

い（図 8a）．一方，実測収量差分の短期予測精度は r(1) = 0.747 および RRMSE (1) = 0.666 であ，RRMSE (1) が 1.0 よりもある程度小さいため，年次間の増減については予測が出来り（図 8b）ている．
実データ（収量）に対する同期強度を評価関数としたパラメータ推定の結果，空間勾配 α を組み込んだモデル A：(RC , eC , ε, α) = (1.4, 0.2, 0.05, 0.2857) が適当とされたが，モデル A の収量の予測精度は r(1) = 0.377 および RRMSE (1) = 1.41 となり予測不能であった（図 8c）．一方，収量差分に関しては，r(1) = 0.802 および RRMSE (1) = 0.728 とある程度の予測精度を示した（図 8d）．これらの結果は図 8a–b の実データの場合と同様であった．
一般的な空間勾配が無い果樹園を想定しているモデル B：(RC , eC , ε, α) = (1.4, 0.2, 0.05, 0) においては，収量の短期予測精度は r(1) = 0.760 および RRMSE (1) = 0.693 と，ある程度の予測．収量差分の短期予測では r(1) = 0.939 および RRMSE (1) = 0.365 とが出来たといえる（図 8e）
高い予測精度を示した（図 8f）．これらは，図 7 で示した 48 個体のウンシュウミカンの短期予測と同等の精度である．
以上の結果から，制御パラメータの空間勾配 α の存在は，当該果樹園もしくは当該地域に特異的な属性と考えられる．図 8e–f で示したモデル B の予測結果は，一般的なピスタチオ果
樹園では短期予測が可能であることを示唆しており，今後，現場応用を通じての検証が期待される．
ウンシュウミカンやピスタチオ等の樹木作物では，周期 2 年の隔年結果が栽培上問題となる．一方，気象データなど外生的な変数を用いず，過去の収量データのみから 1 年先の差分値を予測できていることは，隔年結果が内生的なダイナミクスに因る変動が外生要因に優越していることを示している．また，栽培現場では，もちろん高い精度での収量予測が望ましいが，個体毎に翌年の増減が収量差分として高精度で予測できる提案手法は，栽培技術・栽培計画・資材調達など営農現場において有益である．
5.3 アンサンブル再構成法による収量予測技術の課題
バタフライ効果で知られる Lorenz はポアンカレと並ぶカオス理論の生みの親と言えよう（Lorenz, 1963）．Lorenz 類推法は最初の決定論的非線形予測法として紹介されているが，彼が気象データの予測に用いたところ期待した予測精度が得られなかったとされる（池口 他,
2000）．しかし，彼の提案した手法は，本稿で示したようにカンキツ類やナッツ類などの豊凶（p = 1）とは最短時間単位であるが，農作業栽培
予測には極めて有効である．また，1 ステップの観点からは p = 1 年は十分に長期であるので営農の観点からは将来予測が出来ている．果樹
収量時系列集合からのダイナミクスのアンサンブル再構成が有効であれは，決定論的非線形予測だけではなく，決定論的非線形時系列解析の諸手法の活用が期待できるが，個体レベルおよびそれよりマクロまたミクロな諸現象においての課題を下記に示す．
5.3.1 個体レベル
多数の集合時系列を収集・蓄積すれば，集合サイズ（個体数）や時系列サイズ（観測年数）の拡大によって予測器の性能向上を図ることができる．気象データなど外生要因をモデルに組み込む工夫をすれば，予測性能を向上させることができる．ダイナミクスが再構成できれば，本稿で紹介した決定論的非線形予測の他に，相関次元，Lyapunov 指数，推移誤差等の諸量を求めることができる，そしてこれらが，果樹園栽培の新たな評価変数になるかもしれない．また，
再構成ダイナミクスへのカオス制御（Sakai and Noguchi, 2009）は，熟練農業者による吟味を経て，摘花，摘果，剪定，接ぎ木，肥培管理などの栽培技術に寄与できるかもしれない．

統計数理

144

第 71 巻 第 2 号 2023

5.3.2 マクロレベル
個体レベルの RBM ダイナミクスはカオス振動子でもある．カオス同期（Pecora and Carroll,
によって上位の階層にもダイナミクスが創発しうる．現に，
1990, 2015; Rosenblum et al., 1996）
ウンシュウミカンでは，個体→果樹園→集荷組織→市町村→府県→全国という少なくとも 6 つの過程において隔年結果が存在しているのであるが，個体から集荷組織レベルまでは RBM のもつ内生的要因と外生的要因（気象）という自然要素に拠っていると考えてよい．一方，農林業
センサスからも明らかなように，全国のウンシュウミカン収穫量は明確な 2 年周期の変動を示「成り年」と「不成り年」
が交互にし，年次間の変動量は 10% から 20% 程度で，ほとんど確実に繰り返される．そして, 当該業界ではこれがあたかも自然の法則のように受け止められている．
しかし，個体の平均収量から類推して約 1 億本のウンシュウミカンが栽培されているが，日本列島の西半分全域で共通のノイズ印加が毎年実現しているとは到底考えられない（Sviridova
and Sakai, 2016）．見込み予測行動による創発という行動経済学的仮説についての検証を試みているところであるが，そのためには，6 つの過程のどこかで機能しているスイッチング行動を実データから探索しなければならない．いずれにせよ，全国収穫量の 2 年周期の隔年結果は事実であり，このメカニズムの解明は大変興味深い研究課題である．
5.3.3 ミクロレベル
個体レベルより下位の階層において RBM を引き起こすメカニズムが当然に存在しているはずである．開花遺伝子の発現（Shalom et al., 2012），花芽形成ホルモン等の寄与が無ければRBM は振動子として作動できない．よって，この機能を担う物質とそのメカニズムの解明は，それとともに，2.2.3 項の拡散結合モデル
大変興味深いテーマである（Satake and Kelly, 2021）
として紹介したように，根系ネットワークによる物質交換が明らかに樹木豊凶に影響を与えている（Sakai et al., 2022）．また，牧草地は 1 年生の牧草からなっているが，収穫後に耕起し翌春に播種というのではなく，収穫後に残された地下茎ネットワーク（ルートマット）によって再生される．そのため，ダイナミクスの存在を仮定することは的外れではい（Tilman and Wedin,
．
1991）
6.

カオスと農業

農業機械の研究者にとっては，播種・耕うん・肥培管理また病害虫・野生動物防除など通じて密度依存性は身近な概念である．密度依存性をモデル化した Logistic 写像（May, 1976）の周期倍分岐のカオスも農業現場での経験と平仄があう．トラクタ非線形動力学の要素過程である Bouncing Ball Dynamics は Logistic 写像と全く等価な分岐構造（周期倍分岐）を持っている（Watanabe and Sakai, 2019）．本稿で紹介した隔年結果の数理は昇圧型 DC-DC 変換器の数理と完全に等価で（Banerjee et al., 2000; 酒井, 2023），周期増加分岐のカオスルートを示す．これらは単なる偶然の一致なのであろうか．低次元の非線形ダイナミクスは農業のいたるところ．農業は必然的に人間が見て触れるに転がっている（酒井, 1993, 1995, 1997, 2023; Sakai, 2001）
量の変化を扱う．そこでは，再現性と予測および類推可能性を備えたものでなければ農業技術
としては成立していない．科学革命後の産業技術基盤となった線形理論建設の舞台裏で，様々な非線形現象への関心も着実に育まれていたのであろうし, ようやくポアンカレが 3 体問題として可視化してくれたともいえる（高橋 他, 1998）．その後さらに 1 世紀を要したが．カオス，
非線形力学，同期理論，グラフ，ネットワーク，粒子系のモデリングと解析，そしてリモート
センシングなどによるオブザーバの収集能力，さらに手元のパソコンの計算パワーによって実証可能な技術基盤を農業技術として構築できる道具立てを近代科学はようやく揃えてくれたよ
うに思う．本稿で提案したアンサンブル再構成法は果樹栽培に限らず，カオス集団力学の決定

果樹豊凶の決定論的非線形予測

145

論的時系列解析法としての一般性を有するように定式化できている．農業とカオスは親和性が高いのである．
謝辞

文科省科学研究費 20K21347 の補助を受けた．

参

考

文

献

Abe, T., Tachiki, Y., Kon, H., Nagasaka, A., Onodera, K., Minamino, K., Han, Q. and Satake, A.
(2016). Parameterisation and validation of a resource budget model for masting using spatiotemporal flowering data of individual trees, Ecology Letters, 19, 1129–1139.
Banerjee, S., Karthik, M. S., Yuan, G. and Yorke, J. A. (2000). Bifurcations in one-dimensional piecewise smooth maps-theory and applications in switching circuits, IEEE Transactions on Circuits
and Systems I, 47, 389–394.
Crone, E. E., Polansky, L. and Lesica, P. (2005). Empirical model of pollen limitation, resource acquisition, and mast seeding by a bee-pollinated wildflower, American Naturalist, 166, 396–408.
Esmaeili, S., Hastings, A., Abbott, K., Machta, J. and Nareddy, V. R. (2021). Density dependent
resource budget model for alternate bearing, Journal of Theoretical Biology, 509, 110498.
Esmaeili, S., Hastings, A., Abbott, K. C., Machta, J. and Nareddy, V. R. (2022). Noise-induced versus
intrinsic oscillation in ecological systems, Ecology Letters, 25, 814–827.
Gaion, L. A. and Carvalho, R. F. (2018). Long-distance signaling: What grafting has revealed?, Journalof Plant Growth Regulation, 37, 694–704.
池口徹, 合原一幸 (1997). 力学系の埋め込み定理と時系列データからのアトラクタ再構成（〈特集〉力学系, 応用数理, 7(4), 260–270.
理論 - 応用数理における新しい展開）
池口徹, 小室元政, 山田泰司 (2000).『カオス時系列解析の基礎と応用』（合原一幸 編）, 産業図書, 東京.
Isagi, Y., Sugimura, K., Sumida, A. and Ito, H. (1997). How does masting happen and synchronize?,
Journal of Theoretical Biology, 187, 231–239.
Ishibashi, A. and Sakai, K. (2019). Dispersal of allergenic pollen from Cryptomeria japonica andChamaecyparis obtusa: Characteristic annual fluctuation patterns caused by intermittent phasesynchronisations, Scientific Reports, 9, 11479.
Kaneko, K. (1985). Spatiotemporal intermittency in coupled map lattices, Progress of Theoretical
Physics, 74(5), 1033–1044.
Kaneko, K. (1990). Globally coupled chaos violates the law of large numbers but not the central-limit
theorem, Physical Review Letters, 65, 1391–1394.
Klein, T., Siegwolf, R. T. W. and Körner, C. (2016). Belowground carbon trade among tall trees in atemperate forest, Science, 352(6283), 342–344.
Koenig, W. D. and Knopes, J. M. H. (2005). The mystery of masting in trees, American Scientist, 93,
340–347.
Lev-Yadun, S. (2011). Why should trees have natural root grafts?, Tree Physiology, 31, 575–578.
Lorenz, E. N. (1963). Deterministic nonperiodic flow, Journal of Atmospheric Sciences, 20, 130–141.
Lyles, D., Rosenstock, T. S., Hastings, A. and Brown, P. H. (2009). The role of large environmental noise in masting: General model and example from pistachio trees, Journal of Theoretical
Biology, 259(4), 701–713.
May, R. (1976). Simple mathematical models with very complicated dynamics, Nature, 261, 459–467.
Noble, A. E., Machta, J. and Hastings, A. (2015). Emergent long-range synchronization of oscillating

146

統計数理

第 71 巻 第 2 号 2023

ecological populations without external forcing described by Ising universality, Nature Communications, 6, 6664.
Noble, A. E., Rosenstock, T. S., Brown, P. B., Machta, J. and Hastings, A. (2018). Spatial patterns oftree yield explained by endogenous forces through a correspondence between the Ising modeland ecology, The Proceedings of the National Academy of Sciences, 115(8), 1825–1830.
Pecora, L. M. and Carroll, T. L. (1990). Synchronisation in chaotic systems, Physical Review Letters,
64, 821–824.
Pecora, L. M. and Carroll, T. L. (2015). Synchronisation of chaotic systems, Chaos, 25, 097611.
Prasad, A., Sakai, K. and Hoshino, Y. (2017). Direct coupling: a possible strategy to control fruit
production in alternate bearing, Scientific Reports, 7, 39890.
Rees, M., Kelly, D. and Bjørnstad, O. N. (2002). Snow tussocks, chaos, and the evolution of mast
seeding, American Naturalist, 160, 44–59.
Rosenblum, M. G., Pikovsky, A. S. and Kurths, J. (1996). Phase synchronization of chaotic oscillators,
Physical Review Letters, 76(11), 1804–1807.
Rosenstock, T. S., Rosa, U. A., Plant, R. E. and Brown, P. H. (2010). A reevaluation of alternate
bearing in pistachio, Scientia Horticulturae, 124, 149–152.
Rosenstock, T. S., Hastings, A., Koenig, W. D., Lyles, D. J. and Brown, P. H. (2011). Testing Moran’stheorem in an agroecosystem, Oikos, 120, 1434–1440.

酒井憲司 (1993). 非線形農学の萌芽, 数理科学, 31(9), 43–48.
酒井憲司 (1995). 耕地のカオス, 数理科学, 33(3), 67–72.
酒井憲司 (1997).『カオス農学入門』, 朝倉書店, 東京.
Sakai, K. (2001). Nonlinear Dynamics and Chaos in Agriculture Systems, Elsevier Science, Amsterdam.
酒井憲司 (2023). カオス同期に基づく自然への理解と農業応用, 電気学会誌, 143(1), 29–32.
Sakai, K. and Noguchi, Y. (2009). Controlling chaos (OGY) implemented on a reconstructed ecologicaltwo-dimensional map, Chaos, Solitons and Fractals, 41(1), 630–641.
Sakai, K., Noguchi, Y. and Asada, S. (2008). Detecting chaos in a citrus orchard: Reconstruction of
nonlinear dynamics from very short ecological time series, Chaos, Solitons and Fractals, 38(5),1274–1282.
Sakai, K., Hoshino, Y., Prasad, A., Fukamachi, A. and Ishibashi, A. (2019). Period-3 dominant phasesynchronisation of Zelkova serrata: Border-collision bifurcation observed in a plant population,
Scientific Reports, 9, 15568.
Sakai, K., Brown, P., Rosenstock, T., Upadhyaya, S. and Hastings, A. (2022). Spatial phase synchronisation of pistachio alternate bearing: Common-noise-induced synchronisation of coupledchaotic oscillators, Chaos, Solitons and Fractals, 165(2), 112764.
佐竹暁子 (2007). 理論と実証分析の相互フィードバック：植物の繁殖同調モデルを例に, 日本生態学会
誌, 57, 200–207.
Satake, A. and Iwasa, Y. (2000). Pollen coupling of forest trees: Forming synchronized and periodic
reproduction out of chaos, Journal of Theoretical Biology, 203, 63–84.
Satake, A. and Kelly, D. (2021). Studying the genetic basis of masting, Philosophical Transactions of
the Royal Society B, 376(1839), 20210116.
Shalom, L., Samuels, S., Zur, N., Shlizerman, L., Zemach, H., Weissberg, M., Ophir, R., Blumwald,
E. and Sadka, A. (2012). Alternate bearing in citrus: Changes in the expression of flowering
control genes and in global gene expression in ON- versus OFF-crop trees, PLoS One, 7(10),e46930.
Simard, S. W., Pery, D. A., Jones, M. D., Myrold, D. D., Durall, D. M. and Molina, R. (1997). Net
transfer of carbon between ectomycorrhizal tree species in the field, Nature, 388, 579–582.

果樹豊凶の決定論的非線形予測

147

Sviridova, N. and Sakai, K. (2016). Common noise induced synchronization on collective dynamics of
citrus production, Journal of the Japanese Society of Agricultural Machinery and Food Engineers, 78(3), 221–226.
高橋陽一郎, 山口昌哉, 辻下徹 (1998).『数学 — 複雑系の科学と現代思想』, 青土社, 東京.
Tilman, D. and Wedin, D. (1991). Oscillations and chaos in the dynamics of a perennial grass, Nature,353, 653–655.
Watanabe, M. and Sakai, K. (2019). Impact dynamics model for a nonlinear bouncing tractor during
inclined passage, Biosystems Engineering, 182, 84–94.
Ye, X. and Sakai, K. (2016). A new modified resource budget model for nonlinear dynamics in citrus
production, Chaos Solitons and Fractals, 87, 51–60.

148

Proceedings of the Institute of Statistical Mathematics Vol. 71, No. 2, 129–148 (2023)

Deterministic Nonlinear Prediction of Tree Crops—Ensemble Reconstruction of Discrete Dynamics and a Lorenz Method for
One-year Forward Prediction of Individual Yield—
Kenshi Sakai
Institute of Agricultural Science, Tokyo University of Agriculture and Technology

Deterministic chaos is generated from low-dimensional dynamics with orbital instability. Discriminating chaotic behavior from observed time series is the main objectiveof nonlinear time series analysis. Deterministic nonlinear prediction is a powerful tool
for this purpose but requires large time series sizes. Because crop yields can only bemeasured once per year, the observed yield time series size is very small. However, thenumber of crop individuals in a field or orchard is large; thus, the size of the time series
sets is large. In addition, the reproductive patterns of many fruit trees can be assumed tohave low-dimensional nonlinear dynamics. In this paper, we apply deterministic nonlinear
prediction to yield time series sets obtained from orchards and present practical examples
of one-year-ahead forecasting of yields for Citrus unshiu and Pistachio (Pistacia vera).

Key words: Chaotic synchronization, deterministic nonlinear prediction, alternate bearing, ensemble reconstruction of dynamics.