FRA-SA-2024-ABCWG02-04

令和 6(2024) 年度 再生産関係の推定・管理基準値計算・将来予測シミュレーションに関する技術ノート
資源評価高度化作業部会

2024 年 6 月 21 日

1 令和 5 年度版からの変更・追加点• 編集上のミス等の微小な修正

2 基本的な式と記号の説明
本資料で使われる主要な数式記号の定義を表 1 に，ホッケー・スティック (HS)[1]，ベバートン・ホルト(BH)[2]，リッカー (RI)[4] 型の再生産関係式，個体群動態の基本的な式を以下に示した．
ホッケー・スティック型再生産関係 (HS)(
R(SBt |a, b) =aSBt−Sminab(SBt−Smin < b)(SBt−Smin ≥ b)(1)

ただし，親魚量と加入量が線形関係の場合には過去最大親魚量を HS の折れ点と仮定する．一方，観察された親魚量の範囲で加入が一定となる場合には，HS の折れ点を過去最小親魚量とする．結果として HS の折
れ点 (b) が過去最小親魚量または過去最大親魚量となった場合，b は仮定値をあてていることになるが，bを過去最小親魚量から過去最大親魚量の範囲で探索した結果としての推定値と考え，AIC などの情報量規準の計算時はどのような場合でも a, b については 2 パラメータ分が推定されているものとした．
ベバートン・ホルト型再生産関係 (BH)

R(SBt |a, b) =aSBt−Smin
1 + bSBt−Smin(2)リッカー型再生産関係 (RI)

R(SBt |a, b) = aSBt−Smin exp (−bSBt−Smin )1(3)

その他の再生産関係式
その他の再生産関係式についても使用についての十分な科学的根拠と合意があれば使用できる．たとえば、Shepherd (式 4)[6]、Cushing (式 5)[6] などの再生産関係式がある．

R(SBt |a, b, γ) =aSBt−Smin
1 + b(SBt−Smin )γ(4)bR(SBt |a, b) = aSBt−S
min

(5)

また、Hockey-Stick 型の再生産関係は折れ点の部分で微分不可能であるため、パラメータの最尤推定が難しい場合がある。それを解決するために、連続的な Hockey-Stick 型の再生産関係として Mesnil (式 6)[7]も提案されている。この場合のパラメータ γ の推定は難しいため、適当な値を固定して与えるのが一般的である。

R(SBt |a, b, γ) = 0.5a(SBt−Smin +pp
b2 + 0.25γ 2 − (SBt−Smin − b)2 + 0.25γ 2(6)個体群動態の式
個体群動態の決定論的な基本式は以下のような年齢構造モデルで表す．


R(SBt |a, b)
Ns,t = Ns−1,t−1 exp(−Ms−1,t−1 − Fs−1,t−1 )
Ns−1,t−1 exp(−Ms−1,t−1 − Fs−1,t−1 ) + Ns,t−1 exp(−Ms,t−1 − Fs,t−1 )s = SminSmin < s < Smaxs = Smax(7)

また、プラスグループを考慮しない場合は(
Ns,t =

R(SBt |a, b)
Ns−1,t−1 exp(−Ms−1,t−1 − Fs−1,t−1 )s = SminSmin < s ≤ Smax

(8)

となる。

3 再生産関係式のパラメータ推定
再生産関係式のパラメータ推定においては、まず、残差の自己相関を考慮せずにパラメータ a，b を推定する (3.1 節)．このときは最小二乗法または最小絶対値法を使用する．その後，残差の自己相関も考慮するかどうかを検討する (3.2 節)。

3.1 パラメータ a，b の推定
観測値と再生産関係式からの予測値との残差 εt (T1 ≤ t ≤ T2 ) を
obs
εt = log(NSobs) − log R(SBt−S|a, b)
min ,tmin(9)と定義し，最小二乗法の場合には
T2
X
t=T1

2

ε2t

(10)

最小絶対値法の場合には
T2
X

|εt |

(11)

t=T1

を最小化することで，パラメータ a，b を推定する．最小絶対値法の方が，外れ値の影響を受けにくく，頑健な推定値が得られやすい．AIC 等の情報量規準で使用される対数尤度は，最小二乗法の場合は正規分布，
最小絶対値法の場合はラプラス分布を用いて計算される．どちらの最適化法を用いた場合でも，残差の標準偏差 σ は

v
u
u
σ=t

2
X
1
ε2t
T2 − T1 + 1

T

(12)

t=T1

と定義する．

3.2 残差の自己相関の有無を判断する推定された再生産パラメータ a，b のもとで，加入の予測値に対する観測値の残差 εt に自己相関があるかどうかを検討する．自己相関を導入したほうがよいかどうかは，ρ が統計的に有意かどうかや AIC 等を用いて判断する．また，簡便には R[10] で実装されている関数 acf で出力される図で示される ρ の 95% 信頼
区間等からも判断できる．
自己相関を推定する方法としては，再生産関係式のパラメータ a，b を推定してから，残差についての自己回帰モデル (式 13) で自己相関係数 ρ を推定する二段階の手法（二段階推定）が簡便である。

εt = ρεt−1 + ξt

(13)ξt ∼ Normal(0, (1 − ρ2 )σ 2 ))(14)

ここで ξt は

v "
u
u1
1
σ=t
ε2T1 +n1 − ρ2T2X#(εt − ρεt−1 )2(15)

t=T1 +1

二段階推定は安定的なパラメータ推定が得られる場合が多い．
一方、a，b，ρ を含む正確尤度 (exact likelihood, 式 9, 13, 14, 15) を最小化することで，a，b，ρ を同時に推定する方法（同時推定）もある。同時推定によるパラメータ推定は，最尤推定値を得ることが難しい，
推定値にバイアスが生じる，推定値が不安定になるなどの問題点が指摘されている [9]．

ρ の推定手法として 2 段階推定を使うか、同時推定を使うかは、パラメータの安定性や局所解の存在、
ブートストラップ等、一連のモデル診断を通じて判断する。ただし、マアジ太平洋系群を模したシミュレー
ションの結果では、自己相関係数が非常に大きい場合、二段階推定は真の自己相関係数を過小推定する傾向が見られた ([11])。そのため、自己相関が非常に大きいと思われ、それが加入変動に大きく影響を与える可能性がある場合、同時推定を用いたほうが加入変動の不確実性をより正確に捉えられると考えられる。ただし、その場合でも、推定パラメータの安定性を十分確認したり、MSE などを用いて誤りに対する影響を評価する必要はある。
自己相関の影響を除いたときの加入の残差の分散は (1 − ρ2 )σ 2 で，自己相関でも説明しきれなかったランダムな加入変動の大きさを示す．平成 31 年度研究機関会議で SD として示されている数値はに対応し、それ以降の研究期間会議で示された SD は式 12 に従う。

3p(1 − ρ2 )σ 2

図 1 HS に解の不安定性がある場合の例。最尤推定値と思われる値周辺にランダムに初期値を設定してパラメータ推定したときに推定された HS のパラメータ b(x 軸) とそのときの対数尤度（y 軸、左）、
試行した中で得られた最大の対数尤度と個々の推定で得られた対数尤度との差の絶対値の対数（log10）をとったもの（y 軸、右)。右図の y 軸は差の log10 をとったものなので、-3 は最大対数尤度との差が0.001 であることを示す。図の黒点と赤点は、異なる乱数のシードからランダムに初期値を変えてそれ
ぞれ 500 回繰り返し計算を行った場合の結果。左図で示すように尤度が高い b=5 万トン付近では尤度の値がほとんど変わらないパラメータの範囲が存在する。右図で、尤度の差の対数をとると 10−6 以下の差で尤度の違いが見られるが、尤度の違いの傾向はスムーズではなく、高い初期値依存性を示唆する。

3.3 HS を用いるときの注意点
連続的に変化しない HS 再生産関係においては a, b のパラメータ推定の際、初期値の値によって異なるパラメータが推定され、また、その目的関数の値がほとんど同じとなるような状況が見られることがある（図 1）。特に HS で最小絶対値法を用いる場合に生じやすい。このような状況の場合、推定パラメータが初期値に依存し、単独の推定では頑健で再現性のあるパラメータ推定値を得られない問題点がある。このような状況が生じた場合には、図 1 のように初期値を変えた解析を実施し、推定されたパラメータの中央値に最も近いパラメータを得る初期値を利用して推定したパラメータ推定結果を用いることにする。

4 管理基準値の推定
管理基準値は，推定された再生産関係パラメータと確率的な加入変動等のもとで，漁獲係数を一定にした場合の確率的な平衡状態 (時系列的に平均値が一定) となるときの漁獲量や資源量の平均値から推定する．
通常は，近年の代表的な年齢別漁獲係数 Fcurrent,s に対して，選択率を保ったまま一律に漁獲圧を削減または増加させるような漁獲係数一定のシナリオを想定する．この場合の将来の漁獲係数は Fcurrent,s に特定の乗数 α を乗じた αFcurrent,s となる．ここで 1 − α は，例えば α < 1 の場合，目的とする管理基準値を得るために削減すべき漁獲努力量の削減率と解釈される．Fcurrent,s は，通常，資源評価最終年を含む直近 l 年4図2

将来予測における誤差分布の仮定の違い．(左) 対数正規分布の誤差を仮定した場合．(右) 残差のリサンプリングによる誤差を仮定した場合

分の年齢別漁獲係数の平均値が用いられ，以下のように計算される．
T3
X

Fcurrent,s =t=T3 −l+1

Fs,t
l

(16)

確率的シミュレーションにおける加入尾数は，再生産関係式から計算される加入尾数にランダムな誤差

exp(εkt ) を乗じたものとして計算される．εkt は，対数正規分布の誤差を仮定する場合には(Normal(0, σ 2 )
ρ=0
εkt ∼
k
2 2
Normal(ρεt−1 , (1 − ρ )σ ) ρ > 0(17)となる．また，残差リサンプリングによるシミュレーションを実施する場合にはεkt ∼ (random draw from εt=T1 ,…,T2 )(18)となる（図 2）．
管理基準値計算のために実施する確率的シミュレーションにおける個体群動態式は以下となる．

kkR(SBt |a, b) exp(εt + κ)kkNs,t
= Ns−1,t−1exp(−Ms−1,t−1 − αFcurrent,s−1 ) kk
Ns−1,t−1 exp(−Ms−1,t−1 − αFcurrent,s−1 ) + Ns,t−1exp(−Ms,t−1 − αFcurrent,s )s = SminSmin < s < Smaxs = Smax(19)

プラスグループを考慮しない場合には、

(
k
Ns,t
=

R(SBtk |a, b) exp(εkt + κ)kNs−1,t−1exp(−Ms−1,t−1 − αFcurrent,s )s = SminSmin < s ≤ Smax

(20)となる。ここでは t ≥ T3 で，t = T3 のとき Ntk = Ntobs ，Ftk = Ftobs とする．κ は漁獲なしのときの親魚
量 (初期資源量)SB0 の決定論的な値 (式 7 または 8) と，確率的予測値 (式 19 または 20) の平均を一致させるためのバイアス補正項で，対数正規分布を仮定する場合には
κ = −0.5σ 2 ，残差リサンプリングをおこな∑う場合には − log(

T2
t=T1 exp(εt )n) となる．

管理基準値は，α を様々に変え，十分平衡状態に達したと判断される T4 年における状態から推定する．
経験的には，世代時間 G の 20 の年数 T4 = T3 + 20G であれば十分平衡状態に達するようだが，系群によって平衡状態に十分達したと判断される年数を用いる．

4.1 MSY 管理基準値
MSY 管理基準値は，K 回シミュレーションを実施した場合の T4 年における漁獲量 CTk4 の平均値

∑K

(

k
k=1 CT4

K) が最大になるような α を探索的に求め，推定する．CTk4 は，資源量推定に Baranov の漁獲方5

程式を使っている場合，

CTk4 =SX
max
s=SminαFcurrent,skexp(−αFcurrent,s − Ms,T4 )Ns,Tv
4 s,T4αFcurrent,s + Ms,T4(21)

Pope の近似式を使っている場合
CTk4 =

SX
max(1 − exp(−αFcurrent,s )) exp(−s=Smin

Ms,T4
k
)Ns,T
v
4 s,T42(22)kで計算される．このときの平均漁獲量を MSY, αFcurrent,s を Fmsy,s (Fmsy ) と定義する．また，Fs,t=PK
Fmsy,s で K 回シミュレーションを行ったときの t = T4 における親魚資源量の平均 ( k=1 SBTk4 /K) がPK
SBmsy ，MSY/ k=1 BTk4 K −1 を MSY を与える漁獲率 (Umsy ) とする．K は乱数の違いによって推定値が大きく変わらないくらい十分に大きい数を用いる（自己相関なしの場合は最低 2000 回，ありの場合は 1

万回以上が望ましい）．

Fmsy,s は年齢別の漁獲係数のベクトルとして与えられ，特定の選択率を前提とする．そのため，神戸プ
obs
ロットなどのために，過去に見られた年齢別漁獲係数 Ft,sが Fmsy,s に比べてどの程度大きかった・小さかったか (F/Fmsy ) を計算するためには，MSY 計算で前提とされた選択率と過去の選択率の違いを考慮すobs
る必要がある．具体的には Fmsy,s のときの SPR に一致するような xFt,sを求め，1/x = F/Fmsy と定義する．

4.2 その他の管理基準値
その他，平衡状態において特定の平均漁獲量 (たとえば，限界管理基準値の標準値となる MSY の 60%や，禁漁水準を与える MSY の 10%、あるいは特定の年齢群の漁獲量の最大化など) や平均親魚量 (たとえば SB0 の 30% など) を達成するような管理基準値を計算する場合も同様に，t = T4 において平均漁獲量または親魚量が目的とする値と一致するような α を探索的に求め，得られた漁獲係数のもとで平衡状態に達したときの平均親魚量を資源量の管理基準値として用いる．

4.3 再生産関係に基づかない生物学的管理基準値
再生産関係の不確実性が高い場合には %SPR や YPR など、再生産関係の仮定を用いなくても計算できる生物学的管理基準値から管理基準値の提案をすることもできる（1B ルール）．%SPR と YPR は加入

1 尾あたりが生涯に期待される個体数 (Ns , numbers per recruit, NPR) から、加入 1 尾あたりの親魚量(Spawning per recruit, SPR) や漁獲量 (yield per recruit, YPR) を計算することで定義される。NPR(Ns )は漁獲係数 Fs 、自然死亡係数 Ms のもとで、

(
Ns =

1
Ns−1 exp(−Ms−1 − Fs−1 )として計算される。このときの SPR(Fs ) は

s = SminSmin < s < ∞

(23)

P∞

s=smin Ns ws ms となる。また、漁獲がない場合の SPR
SPR(Fs )（SPR(0)）に対する割合（%SPR(Fs ) = SPR(0) × 100) が x% となるような Fs を Fx%SP R と呼ぶ。同様に、Fs のもとでの YPR(Fs ) は、Baranov の漁獲方程式を用いる場合には∞X
s=SminFsexp(−Fs − Ms )Ns
Fs + Ms

6(24)

Pope の近似式を使っている場合には∞X(1 − exp(−Fs )) exp(−s=Smin

Ms
)Ns
2

(25)

から計算される。Fmax は YPR(Fs ) を最大にするような Fs として定義される。F0.1 は YPR(Fs ) 関数の傾きが原点における YPR(Fs ) の傾きの 1/10 となるときの Fs として定義される。これらの漁獲係数をもとにした管理基準値が、1B ルールで Fmsy の代替値として提案される。また、上の SPR, YPR の定義式
における Fs , ws , Ms はすべて時間の次元 t が省略されているが、Fs , ws , Ms が時間によって変化する場合には、適宜適切なパラメータを用いることとする。
また、1B ルールのもとでの将来予測では、再生産関係が不明なため、親が減ると子が減るような親子関係が仮定されず、将来の加入尾数が親魚量によらず一定になるような加入で将来予測が実施される。ただし、加入の確率的な変動を表現するため、加入の確率分布としては過去に観測された加入尾数のリサンプリ
ングや対数正規分布などが仮定される。将来予測における加入の平均尾数が N̄Smin で、Fx%SPR が Fmsy の代替値であるとき、SBmsy の代替値を SPR(Fx%SP R ) × N̄Smin と定義する。また、漁獲がない場合の親魚資源量 SB0 の代替値は SPR(0) × R̄ となる。
限界管理基準値の代替値は SBmin （過去最小親魚量）または SB0 の代替値の 10% 程度の値が用いられる。ここで、SBmin = mint=1,...,T3 SBt である。

5 ABC 算定や目標達成確率計算のための将来予測
漁獲管理規則のもとでの中長期的な将来予測でも，加入誤差については，基本的に式 17，18 などから計k算された εkt を用いて計算する．一方，将来の漁獲係数は漁獲管理規則に従って毎年異なる値 Fs,tを用いる．したがって，短期予測のための将来予測の個体群動態式は

kkR(SBt |a, b) exp(εt + κ)
k
k
k
Ns,t
= Ns−1,t−1exp(−Ms−1,t−1 − Fs−1,t−1
)

 k
k
k
k
Ns−1,t−1 exp(−Ms−1,t−1 − Fs−1,t−1 ) + Ns,t−1exp(−Ms,t−1 − Fs,t−1)

s = SminSmin < s < Smaxs = Smax(26)となる．また、プラスグループを考慮しない場合は

(
k
Ns,t
=

R(SBtk |a, b) exp(εkt + κ)kk
Ns−1,t−1exp(−Ms−1,t−1 − Fs−1,t−1)

s = SminSmin < s ≤ Smax

(27)kは，通常，管理の開始を T3 + 2 年目からと想定し，となる．Fs,t

(
k
Fs,t
=

Fcurrent,sβγFmsy,s = αβγFcurrent,st = T3 + 1
t > T3 + 1(28)kと仮定する．t = T3 + 1 年目の Fs,tについては必ずしも Fcurrent,s を用いる必要はなく，T3 + 1 年目の漁獲状況を最もよく代表していると判断される漁獲係数をあてる．t > T3 + 1 では，漁獲管理規則から計算される β ，γ を Fmsy,s に乗じる．漁獲管理規則における β は 1 以下の安全係数で一定値をとり，γ は図 3に示された漁獲管理規則に従って以下のように計算される．



1
γ=

SBlimit < SBtk
SBtk −SBban

SB−SBban limit0

SBban < SBtk ≤ SBlimit
SBtk ≤ SBban
7(29)図 3 1 系資源の漁獲制御ルールの模式図

ここで SBlimit は限界管理基準値，SBban は禁漁水準である．

ABC は通常 T3 + 2 年目の漁獲量の平均値として算出される．資源量推定に Baranov の漁獲方程式を使っている場合にはCTk3 +2 =

SX
max

k
Fs,T
3 +2

s=SminkFs,T+ Ms,T3 +2
3 +2

k
k
exp(−Fs,T− Ms,T3 +2 )Ns,Tv
3 +2
3 +2 s,T3 +2(30)

Pope の近似式を使っている場合にはCTk3 +2 =

SX
max

k
(1 − exp(−Fs,T)) exp(−3 +2

s=SminMs,T3 +2 k
)Ns,T3 +2 vs,T3 +2
2

(31)

PK

k
CT
3 +2
が ABCT3 +2 として提案される．これは、式 26、式k=1
K
27 から計算される将来の個体群動態とは乖離しているが、２年分の将来予測における漁獲量 CTk3 +2 の平均より計算された漁獲量 CTk3 +2 の平均値

値を ABC とする手順については Okamura et al (2020)[12] による MSE で頑健性が確認されている。また市野川ほか (2022)[13] では、ABC の計算が将来予測と実際で異なることへの影響を調べ、将来予測における平均値などの平均的な挙動にはほとんど影響を与えないものの、ABC の誤差を考慮した場合には将来予測の予測区間がより大きくなることを示している。

6 その他
6.1 再生産関係の誤りによる ABC の推定誤差を考慮した簡易 MSE
ABC は式 30、 式 31 により、資源評価最終年 T3 年の 2 年後 T3 + 2 年の漁獲量の平均値として与えられる。この値は、（資源評価の誤差がないと仮定すると）既知である T3 年までの資源状態から、T3 + 1 年と T3 + 2 年分の将来予測を実施した結果を利用している。したがって ABC を計算するさいには、T3 + 1
年と T3 + 2 年の加入尾数・年齢別体重・自然死亡係数・成熟率、また、T3 + 1 年の年齢別漁獲係数、さらに、HCR を用いて漁獲量を計算する T3 + 2 年には HCR のパラメータとして SBlimit ・SBban ・Fmsy,s が必要となる。
資源評価の誤差がない理想的な条件においても、再生産関係を一つに決定するための十分な情報が得られない場合がある。また、選択した再生産関係とは異なる再生産関係がもし真であった場合、管理基準値
や ABC のバイアスを引き起こす恐れがある。そこで、再生産関係の選択の誤りから生じる管理基準値やABC のバイアスが引き起こすリスクを評価するため、再生産関係の選択の誤りがあるような状況を想定したシミュレーションを実施する。
資源量推定や管理誤差の不確実性を考慮したシミュレーションは一般に管理戦略評価 (Management

Strategy Evaluation, MSE) と呼ばれているが [3]、本シミュレーションでは不確実性を考慮する範囲が狭8いため、あくまで簡易的な MSE である。本 MSE で正しいと仮定する数値と考慮する不確実性を以下にまとめた。ここでは、将来予測の任意の年 T6 年における漁獲量を決定する場合である。

• 正しく知っているとするもの：T6 − 1 年と T6 年の年齢別体重・自然死亡係数・成熟率、T6 − 1 年の年齢別漁獲係数• 間違った再生産関係を用いた場合に間違う値：T6 − 1 年と T6 年の加入尾数 (NSmin ,T6 −1 , NSmin ,T6 )、SBlimit・SBban (つまり γ)・Fmsy,s 、ひいては T6 − 1 年 ABC(ABCT6 −1 )。真の値に対して、将来予′′′測による推定値を NS′ min ,T6 −1 , NS′ min ,T6 ）、SBlimit
・SBban・Fmsy,s・γ ′ ・ABCT′ 6 とする。

シミュレーションでは上記の方法によって推定された ABC(ABCT′ 6 ) が T6 年の TAC として与えられ、ABCT′ 6 どおりの漁獲が行われるとする。つまり、T6 年の漁獲量はkk
C ′ T6 = ABC ′ T6 =SX
max

F ′ s,T6kkexp(−F ′ s,T6 − Ms,T6 )Ns,Tv
6 s,T6k

F ′ ks,T6 + Ms,T6
s=Smin(32)

Pope の近似式を使っている場合にはkk
C ′ T6 = ABC ′ T6 =SXmax

(1 − exp(−F ′ s,T6 )) exp(−ks=Smin

Ms,T6
k
)Ns,T
v
6 s,T62(33)となる。ここで、F ′ s,T6 = cβγ ′ Fmsy,s で、c は、式 32 または式 33 の右辺と左辺が一致するように探索的kkを F ′ s,T6 で置き換えたものを用いる。
に求めた定数である。真の個体群動態は式 26 または 27 の Fs,T
6k

ここで、簡易 MSE の中で ABC 計算のさいに推定された加入尾数や ABC の相対誤差 (Relative error)を評価するため、以下のような相対誤差の指標を計算する。加入尾数の推定誤差はRErec =

(RT6 − RT′ 6 )
RT 6

(34)

(ABCT6 − ABCT′ 6 )
ABCT6

(35)

とした。また、ABC の相対誤差は

REABC−true =と(ABCT6 ,pseudo − ABCT′ 6 )REABC−pseudo =
ABCT6

(36)

の 2 通りを計算した。ABCT6 は全てのパラメータを真の値で計算したときの ABC、ABCT6 ,pseudo はABC 計算のさいに T6 − 1 年と T6 年の加入尾数 (NSmin ,T6 −1 , NSmin ,T6 ) だけ真の値を用いて計算した
ABC（管理基準値のみ間違っている）である。簡易 MSE の適用例と結果の解釈については市野川 (2020,
FRA-SA2020-BRP01-07) で詳述している。

6.2 過去にさかのぼって管理をおこなった場合のシミュレーション
過去に実際に見られた加入のパターンのもとで、もし MSY 管理基準値に基づいた HCR に従って漁獲を調整していたら過去の資源量はどのように変化していたかをシミュレーションする。本シミュレーション
では、過去の任意の年 (T5 年) の年齢別資源尾数 Ns,T5 から資源評価最終年 T3 まで、式 26 または式 27 に従った前進計算をおこなう。そのさいの生物パラメータ (Ms,t 、ms,t 、ws,t 、vs,t ) は過去に見られた実際の値 (VPA への入力値) を用いる。T5 ≤ t ≤ T3 の期間の漁獲係数 Fs,t は、式 28 における t > T3 + 1 の場合

9

の条件の場合の式を用いる。これはつまり、T5 年から HCR に従って漁獲係数を調整した場合を想定している。

T5 + 1 年からは、過去の再生産関係から推定された再生産関数とパラメータ (R(SBTk |a, b)) を期待値とした加入が発生すると考える。毎年の加入の残差は、過去に見られた残差をそのまま使った場合 (シナリオ 1)に加えて、残差を確率分布で発生させる場合 (シナリオ 2) も計算した。シナリオ１もシナリオ 2 も、HCRを用いて漁獲係数を調整したことにより、親魚量が過去に実現されたものとは異なってくるため、過去の加
入尾数の期待値も変化する。シナリオ 1 では、実際の加入尾数の期待値からの確率的なずれを、（環境の良し悪しには漁獲のコントロールが影響しないと考えて）過去に見られた残差の値をそのまま用いる。シナリ
オ 2 では、対数残差が正規分布に従うと考えて確率的に値を発生させるシミュレーションを繰り返す (t 年

(T5 ≤ t ≤ T3 ) の時点で、将来の生物パラメータを既知として将来予測を実施することに相当する）。シナ
リオ 1 の結果は、シナリオ 2 のうちの一つの実現値にあたる。

• 1. 3 節で推定された再生産関係のパラメータ a,b,rho に従って過去の T 年 (T5 ≤ t ≤ T3 ) の加入尾数の期待値 R(SBT |a, b) を計算し、その予測値に過去に実現された加入尾数の残差 εT を指数変換したものを乗じる (NSmin ,t = R(SBt |a, b) exp(εT ))。

• 2.

T 年 の 加 入 尾 数 の 期 待 値 シ ナ リ オ 1 と 同 じ 方 法 で 計 算 す る が 、加 入 の 残 差 は εkT ∼Normal(ρεkT −1 , (1 − ρ2 )σ 2 ) とする。

図4

マイワシ対馬暖流系群で過去にさかのぼった場合のシミュレーションを実施した例。1990 年(start from1990)、2000 年 (start from2000)、2005 年 (start from2005) から、通常レジームの管理
基準値を使って管理を開始した 3 通りの結果の比較。細い実線は過去に見られた残差をそのまま用いたもの（シナリオ 1）、塗りの範囲は過去の残差を確率分布からランダムに発生させたシミュレーションの 90 ％信頼区間、太い実線はシミュレーションの結果の平均値。親魚量が禁漁水準を割っている 2000
年、2005 年では、管理開始当初の漁獲量はゼロになったのち、親魚量が回復したのちに漁獲が再開している。

10

また、マイワシ対馬暖流系群のように、年によって再生産関係のパラメータが異なると仮定されている場
合、本シミュレーションで仮定される過去の再生産関係は、パラメータ推定のさいに仮定した同じ年に再生産関係が切り替わるとしている。図 4 は対馬暖流系群で過去にさかのぼったシミュレーションを実施した
ときの例を示す。

6.3 将来予測において漁獲量一定方策を取る場合
将来予測においてある一定期間、F でなく漁獲量を一定と仮定して漁獲を行うシミュレーションが必要になる場合がある [17]。たとえば将来予測年の t 年目に一定漁獲量 GCt (given catch) で漁獲するとする。そkの場合、式 30, 31 に変わって、以下の式を満たすような係数 x を探索的に求め、xFs,tを Tt 年の漁獲係数

として個体群動態式の中で用いる。
SX
max

GCt =

s=Smink
xFs,tkk
exp(−xFs,t− Ms,t )Ns,t
vs,tk
xFs,t + Ms,t

(37)

Pope の近似式を使っている場合には
GCt =

SX
max

k
(1 − exp(−xFs,t)) exp(−s=SminMs,t k)Ns,t vs,t2(38)

また、資源量が少なくなりすぎて目的とする漁獲量 GCt まで漁獲できない場合には、式 37, 38 においてk
= exp(100) とした場合の漁獲量 GCt,F =100 の 99% まで漁獲をおこなうような x を探索的に求めるこFs,tととする。ここの割合は暫定的に 99% としているが、目的や漁業の実態に応じて適当な値を設定することも可能である。

6.4 TAC の持ち越しを許容する場合
実際の漁業現場では、与えられた ABC を目一杯まで漁獲することはなく、毎年、何割かの ABC が未消化のままとなる状況がよく見られる。そこで、未消化分の ABC を翌年のに付け替えるような ABC の運用の仕方が提案されることがある。ここでは、そのような ABC の持ち越しがどの程度まで可能かを検討する
ためのシミュレーションを以下のようにして実行する。（ここでは暗黙に ABC=TAC として考えるものとする。
）
シミュレーション内で想定されている毎年の ABC(ABCkt ) は、取り残しを考えない場合 Ctk と一致する。
通常の将来予測シミュレーションの場合、式 30 または 31、簡易 MSE の場合、式 32 または 33 から計算される。しかし、取り残しがあってそれが持ち越される場合、毎年の ABC は本来の ABC(ABCkt ) に加えて、持ち越し分が足された ABC+kt となる。その上で、t 年の取り残し率 (reserve rate) を RRt と定義すると、

t 年に取り残す漁獲量は RRt × ABC+kt となる。その場合、t 年の実際の漁獲量は (1 − RRt ) × ABC+kt とkなる。そこで、t 年の漁獲圧 xFs,tは 6.3 節における GCt を (1 − RRt ) × ABC+kt とすることで、式 37 または 38 から求められる。
上の段落から計算された t 年の取り残し量は RRt × ABC+kt であるが、翌年に持ち越せる ABC は次のルールから計算される持ち越し量のうち、小さいほうの値とする。

• ABC の持ち越し分は、t 年にそれが消化されなかったとしても、その分は翌年には持ち越せない。
この条件から計算される ABC の持ち越し分は ABCkt − GCt = ABCkt − (1 − Reservet ) × ABC+ktとなる。ただし、マイナス値になる場合には 0 に置き換えられる。

11

図 5 ABC の持ち越しを許容する場合の結果の図• 持ち越せる漁獲量の上限はもともとの ABCkt に対して一定率 (carry rate)（CRt ）までである。この条件から計算される ABC の持ち越し分は ABCkt × CRt である。
以上のルールに従って将来予測をおこなった場合の例が図 5 である（取り残し率、持ち越し率をともに0.3 とした場合）。2021 年の管理開始年には繰越分が存在しないため、もとの ABC（赤）ともとの ABC ＋
繰越分（緑）は一致し、実漁獲（青）が 3 割小さくなっている。2 年目以降は、もとの ABC ＋繰越分がABC を常に上回っているが、実漁獲は ABC ＋繰越分の 3 割を取り残しているため、実質的にもとの ABCよりも小さい漁獲が実現され、それが翌年に繰り越されている。同様の繰り越し設定がスケトウダラ太平洋
系群・日本海系群で検討されてことがあり、実際の管理シナリオを設定する時の根拠となった [15][16].

6.5 将来予測において漁獲量の変動を一定の幅以内に抑える場合
将来予測において、基本的には漁獲管理規則に従うものの、漁獲量が急激に変動する場合にはその変動に上限（CVlimupper ; 0 ≤ CVlimupper < 1）や下限（CVlimlower ; 0 ≤ CVlimlower < 1）をつける必要がある
場合がある [13]。漁獲量の変動に上限をつける場合には、式 30 または 31 から計算される CTk3 +2 と CTk3 +1
年の漁獲量を比較し、(′CTk3 +2 =CTk3 +2
CTk3 +1 × (1 + CVlimupper )12

CTk3 +2 ≤ CTk3 +1 × (1 + CVlimupper )CTk3 +2 > CTk3 +1 × (1 + CVlimupper )(39)とする。下限をつける場合には(′CTk3 +2 =CTk3 +2
CTk3 +1 × (1 − CVlimlower )CTk3 +2 ≥ CTk3 +1 × (1 − CVlimlower )CTk3 +2 < CTk3 +1 × (1 − CVlimlower )(40)′とする。下限と上限をつける場合には上式の両方を適用する。将来予測では T3 + 2 年の漁獲量が CTk3 +2 とkなるような xFs,tを探索的に求め、これを T3 + 2 年の漁獲係数として個体群動態式の中で用いる。

参考文献[1] Clark CW, Charles AT, Beddington JR, Mangel M (1985) Optimal capacity decisions in a developing ﬁshery. Marine Resource Economics 2:25–53[2] Beverton RJH, Holt SJ (1957) On the Dynamics of Exploited Fish Populations. Her Majesty s
Stationary Oﬃce, London.
[3] Punt et al (2016) Management strategy evaluation: Best practives. Fish and Fisheries. 17:303–334[4] Ricker WE (1954) Stock and recruitment. Journal of Fisheries Research Board of Canada 11:559–
623[5] Shepherd JG (1982) A versatile new stock-recruitment relationship for ﬁsheries, and the construction of sustainable yield curves. Journal du Conseil International pourl’Exploration de la Mer,
40:67–75[6] Cushing DH (1971) The dependence of recruitment on parent stock in diﬀerent groups of ﬁshes.
Journal du Conseil International pourl’Exploration de la Mer. 33: 340-362[7] Mesnil B, Rochet MJ (2010) A continuous hockey stick stock–recruit model for estimating MSY
reference points. ICES Journal of Marine Science, 67: 1780—1784[8] Thorson JT, Jensen OP, Zipkinc EF (2014) How variable is recruitment for exploited marine ﬁshes?
Can J Fish Aquat Sci 71:973–983[9] Johnson KF, Council E, Thorson JT, Brooks E, Methot RD, Punt AE (2016) Can autocorrelated
recruitment be estimated using integrated assessment models and how does it aﬀect population
forecasts? Fish Res 183:222–232[10] R Core Team (2018) R: A language and environment for statistical computing. R Foundation for
Statistical Computing, Vienna, Austria. https://www.R-project.org/[11] FRA-SA2020-BRP01-06. 市野川桃子, 西嶋翔太, 岡村 寛 (2020) シミュレーションを用いた自己相関係数の同時推定手法の推定バイアス評価（非公開資料）[12] Okamura H, Ichinokawa M, Hilborn R (2020) Evaluating a harvest control rule to improve thesustainability of Japanese ﬁsheries. bioRxiv. DOI: 10.1101/2020.07.16.207282[13] 市野川 桃子・西嶋 翔太・ 向 草世香・ 黒田 啓行・ 大下 誠二 (2022) 改正漁業法下での様々な代替漁獲管理規則の検討：マイワシ 2 系群を例に. 日本水産学会誌. 88: 239-255.

[14] FRA-SA2020-BRP01-07. 市野川桃子 (2020) 簡易的 MSE を用いた複数の管理基準値の頑健性の比較・HCR の検討（非公開資料）[15] FRA-SA2020-BRP03-3 (2020) スケトウダラ日本海北部系群における資源管理方針に関する検討会（第１回）からの依頼への対応.https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/HW_suketou_n_20201014.pdf

13

[16] FRA-SA2020-BRP03-4 (2020) スケトウダラ太平洋系群における資源管理方針に関する検討会（第１回）からの依頼への対応.https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/HW_

suketou_p_20201014.pdf[17] FRA-SA2020-BRP03-5 (2020) スケトウダラ太平洋系群における資源管理方針に関する検討会（第２回）からの依頼への対応.https://www.fra.affrc.go.jp/shigen_hyoka/SCmeeting/2019-1/HW_

suketou-pa_20201211.pdf

14

表1

記号

t

本資料で用いられる数式の記号の定義

説明
時間 (time) の次元を表す添え字．通常は年なので，ここでは t 年といった表記を使用する．資源評価開始年を t = 1, 資源評価最終年を t = T3 とする．
将来予測開始年は t = T3 + 1，一定の漁獲圧で漁獲して平衡状態に達したときの年を T4 とする．また，再生産関係の推定のためにデータを用いる最初の加入年を T1 ，最後の加入年を T2 とする．通常は，1 ≤ T1 < T2 ≤ T3 ≪ T4 ．

s

成長段階 (stage) を表す添え字．通常は年齢なので，ここでは s 歳といった表記を使用する．s = Smin , ..., Smax で，Smin は加入年齢，Smax は Smax
歳以上をひとまとめにしたグループ (プラスグループ).

Ns,t

s 歳 t 年の資源尾数．たとえば加入尾数は Nsmin ,t となる．

Fs,ts 歳 t 年における漁獲死亡係数

Ms,ts 歳 t 年における自然死亡係数ms,t

s 歳 t 年における成熟率

ws,t

親魚量を計算するときの s 歳 t 年の個体あたりの体重vs,t

漁獲量を計算するときの s 歳 t 年の個体あたりの体重

Bt

PSmax
t 年の親魚資源量．SBt = s=Smin
Ns,t ws,t ms,t ．
PSmax
t 年の総資源量．Bt = s=Smin Ns,t ws,t ．

k

確率的な将来予測シミュレーションにおける各試行に対する添字．将来予

SBt

測の期間で各パラメータの右上につけて各試行ごとの値を示す．例えば，将来予測における加入尾数の場合は NSkmin ,t , 親魚量の場合は SBtk のように使用する．k = 1, …, K.

obs

資源評価モデルにより得られた過去の資源評価値に対する添字．再生産関係を推定する際には観測値 (observation) として取り扱うため，obs と表記した．例えば過去の加入尾数は Nsobsとなる．
min ,tG

P∞

世代時間 (年)．

s=0 (s exp(−

Ps Ms,t )ms,t )/(exp(−

Ps Ms,t )ms,t ) と定

義する．

n

再生産関係推定に用いるデータの数．n = T2 − T1 + 1．

εtt 年の再生産関係からの加入尾数の予測値と観測値の対数残差σ

再生産関係からの加入尾数の予測値と観測値の対数残差における標準偏差ρ

加入の対数値の残差の自己相関係数

κ

将来予測における加入尾数のバイアス補正のための係数a, b

再生産関係で推定されるパラメータ

R(SBt |a, b)

パラメータ a，b と親魚資源量 SBt を持つ再生産関係式，またはそこから計算される加入尾数の期待値

Reservet

TAC の取り残し率15