1.はじめに
1).発表の構成
1. はじめに（資源評価手法・問題の背景など）
2. 水産資源データ解析におけるNUOPTの利用（体長組成の年齢分解(成分数推定)・構造モ
デルとフルモデルの融合・生物情報の利用）
3. 水産資源評価におけるVMStudioの利用例（まぐろ類の漁獲効率(CPUE)に関する解析）
4. まとめ（NUOPTやVMStudio等に関する感想）
A. 付録（有限正規混合分布の理論的考察(続)）B. 補遺（サポートベクターマシーンのコンセプト）

水産資源データ解析におけるNUOPTおよびVMStudioの利用
庄野 宏（独)水産総合研究センター遠洋水産研究所）
Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

1

1.はじめに2).まぐろ・かつお類の水産資源解析

Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

2

1.はじめに
3).資源評価手法および取り扱う問題

a. 遠洋水産研究所（独)水研センター）の役割• 日本の領海外の漁業資源の水産資源解析→国際漁業委員会などで一括して管理する
b. 水産資源解析の2つの目的• 資源評価（海の中にいる魚の数を推定する）– プロダクションモデルやコホート解析などを使用• 資源管理（評価結果に基づき漁獲量を決定）– 再生産関係（親子関係）や海洋環境などを考慮
Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

3

Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

42.水産資源解析におけるNUOPT利用1).体長組成データの年齢分解問題

2.水産資源解析におけるNUOPT利用
2).過去の知見および生物情報の利用

a.キダイ体長組成(田中1956)

b. 年齢分解問題の条件設定

a. 過去の知見（〜1990年代）

b. 生物的な補助情報の利用

・有限混合正規分布を仮定して成分数およびパラメーター(平均・分散・混合比率)を推定• 成分数は目で見て判断！
• EMアルゴリズムの利用???
（成分数を固定後、各々の正規分布の平均・分散等の推定に際して効果を発揮!）• 耳石の年輪・日輪等の活用• 背鰭や尻鰭等の年齢形質
正規分布の平均に対して上記
VB成長曲線を仮定して母数（K, t0, L(∞)）を推定→構造モデル
L( )[1 exp{ K (t t0 )}]
t2500
2000
N um b er

1500
1000

• 成分数推定に際しては、カイ二乗検定が利用可能???

500
0
8 10 12 14 16 18 20 22 24 26 28 30 32 34 36
Length(cm)Hiroshi SHONO• 成分数推定に際しては、情報量規準は利用不可?????

・Catch-at-size→Catch-at-age
変換にて成分数推定が重要

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

v.s.フルモデル(母数制約無し)

水産資源分野の奇妙な知見
5

Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

62.水産資源解析におけるNUOPT利用3).成分数推定にχ2検定は使用不可

2.水産資源解析におけるNUOPT利用4).罰金付き最尤推定法(Leroux,1992)

例： X 1 ,L , X n (i.i.d.)〜N(x| 1 ,H0 :21n(1- )N(x| 2 ,)1 v.s. H1 : 012log n22) (025-21)( as n)?

l( | X ) (対数尤度関数)
i 1

L:

mlog

j

N(xi | j , 2j )j 1

ˆを
2 l( | X ) 2am,n とおき、成分数mの推定量m対数尤度の過剰な振る舞いゆえに不成立！
• 認定可能性の欠如(非識別性)

ˆ min{min L} の解として求める最小化問題 m• 帰無仮説の下でFisher情報行列が特異• 対数尤度が∞になる場合も存在(非有界性)

但し am,n は am,n 0, am 1,n am,n ,m0(n 0)nを満たす実数列 とする（注) はmに依存している）

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

Hiroshi SHONO

am,n7

Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

8

2.水産資源解析におけるNUOPT利用5).情報量規準による成分数の推定

2.水産資源解析におけるNUOPT利用
6).フルモデル(制約無しモデル)の特徴

罰金付き最尤法(Leroux,1992)により求められた成分数(m) は漸近的に過小推定でない → 一致性とは異なる
am ,n 3m 1→ AIC(3m 1) log n → BICa• 長所m ,n2正規分布以外mjh( x,j)(p.d.f.) でも適用可能!

j 1

Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

9– 柔軟なパラメーター推定が可能• 短所– モデルの正則性(有界性,識別性含む)の問題– 情報量規準による成分数推定の根拠の問題• Leroux (Annals of Statistics, 1992) の[学位]論文“罰則付き尤度関数(AICやBICを含む)による成分数の推定は、漸近的に過小推定にはならない”⇒ 過大推定の可能性は言及なし→片手落ちなのでは?
Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

10

2.水産資源解析におけるNUOPT利用
7).構造モデル(制約付モデル)の特徴

2.水産資源解析におけるNUOPT利用
8).新モデル(Eguchi and Yoshioka,2001)• 長所

コンセプト：フルモデルと構造モデルの融合• 長所– モデルの安定性• パラメーターの次元が縮約されているため、識別性が保たれ、（対数）尤度関数がパラメーター空間で有界

• 短所

• 短所

e.g. VB曲線よりもゴンベルツ曲線が適当な場合 etc.

– 計算手順や背景となる理論が非常に複雑⇒成分数を固定すれば情報量規準により選択可能
数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)→罰則付き尤度関数(フルモデルの尤度と構造モデルの尤度の結合形)の有界性を数学的に証明！

– 成分数(山の数)が推定可能– 構造(VB型の成長曲線など)が推定可能– 成分数の推定が難しい– モデルの柔軟性（特定化に失敗する恐れも）Hiroshi SHONO– “対数尤度の過剰な振る舞い”の問題解決

11

Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

122.水産資源解析におけるNUOPT利用
9).新しいモデルの計算手順概略-(1)13

2.水産資源解析におけるNUOPT利用
10).新しいモデルの計算手順概略-(2)( ) の下で最尤法によりStep1: 構造モデルH:パラメーターの推定値 ( ˆ ) を求める

Step4: 最後に下式により成分数mの推定を行なう

Step2: 罰則付き尤度(構造及びフルモデルの融合)

・考えうる成分数を持つ複数のモデルを最初に仮定l ( )D ( J ) ( ( ˆ ), ))l ( )(1mˆ• クロス・バリデーション （全データから第k番目の観測値をx(k)(k=1,…,n)を抜いての計算実行） はデータ数nや成分数mが大きくなると実行不可能を最大にするパラメーターθの値を求める
Step3: クロスバリデーション(or近似手法)によりチューニングパラメーターλの値を求める
l( ) :D(J )↓近似的手法(ACV)やグリッドサーチの利用

フルモデルの対数尤度関数( , ) : 構造モデルとフルモデルのK-L情報量
数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

Hiroshi SHONO– 平均：VB成長曲線(t | inf, K , )e Kt )

einf (1推定値 L(∞)=48.7, K=0.13, t0= -1.08,
c=0.62, p=0.24, a=0.55, 成分数は４– 分散：平均のべき乗– 比率：減少率一定(t | c, p)c{ (t )} pRatio(t | a )e at /i l

14

2.水産資源解析におけるNUOPT利用
12).新しいモデルを含む計算機実験
X 1 ,L , X n(i.i.d )〜p.d . f .f(x)0.4N(x|10,5)

0.3N(x|14,6) 0.2N(x|19,7) 0.1N(x|25,8)(r1=0.41, r2=0.29, r3=0.19, r4=0.11)
0.08
0.175

真の成分数=4

キダイの例

0.15

0.06

(田中,1956)( L( )[1 exp{ K (t t0 )}]) 0.125
2

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

Hiroshi SHONO

2.水産資源解析におけるNUOPT利用11).モデルの仮定とNUOPTの適用事例• 構造モデルの仮定

min m CV ( , m ) or min m ACV ( , m )0.1

繰り返し数=100

0.04

0.075
0.05

e it

0.02

0.025

k i
5
Hiroshi SHONO

10152025

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

3035
15

2.水産資源解析におけるNUOPT利用13).計算機実験による結果の比較• 新しいモデル（サンプルサイズ=250）

10
Hiroshi SHONO

1520

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

2530162.水産資源解析におけるNUOPT利用14).豪州畜養ミナミマグロの体長組成a.2004年3-8月のage2-4の組成 b.2005年4-10月のage2-4の組成• フルモデル（サンプルサイズ=3000）
成分数 2AICBIC
Bayes*

3111490899459615*-Bayes型規準では混合比率に対するDilichlet事前分布を仮定
Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

17

Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

182.水産資源解析におけるNUOPT利用
15).文献（有限正規混合分布モデル）

3.資源評価におけるVMStudioの利用
1).CPUE(catch per unit effort)の定義• Chernoff, H. 1954: Ann. Math. Statist., 25, 573-578.
• Eguchi, S. and Yoshioka, K. 2001: The Institute ofMathematical Statistics., 809, 30pp.
• Henna, J. 1985: Ann. Inst. Math. Statist., 37, 235-240.
• Leroux, B. 1992: Ann. Statist., 20(3), 1350-1360.
• Shapiro, A. 1985: Biometrika, 72, 133-144.
• Shapiro, A. 1988; Inter. Statist. Rev., 56, 49-62.
• 庄野宏. 2006: 計量生物学, 27(1), 55-67.
• 田中昌一. 1956: 東海区水産研究所報告, 14, 1-4.

• CPUEとは?→漁獲効率を表す一般的な指標!

Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

19

CPUE=

Catch(漁獲量:漁獲尾数や漁獲重量で表わす)
Effort(努力量:延縄船‑針数,巻網船‑操業日数等)• 人間の側から見ると→漁獲効率・漁獲能率• 生物の側から見ると→資源密度（ミクロ的）・相対資源量（マクロ的）(年トレンドの把握)• CPUE標準化とは?→資源の年変動に対応する部分を取り出す作業（加工していないCPUEから資源密度以外の様々な要因を取り除く）• 標準化CPUEは資源評価モデルに利用される
Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

20

3.資源評価におけるVMStudioの利用
2).統計モデルを利用したCPUE解析

3.資源評価におけるVMStudioの利用
3).マイニング手法による解析の動機• LogCPUE-Normalモデル(共分散分析)
log(CPUEijk)=Intercept+Yeari+Areaj+Monthk+(Year*Area)ij+errorijk, errorijk〜N(0,σ2)
Year, Area, Month：順序のないカテゴリカル変数• Catch-Poisson(or NegativeBinomial)モデル(GLM)
E[Catchijk]=Effort*exp(Intercpet+Yeari+Areaj+M
onthk+Year*Areaij),Catchijk~Po(λ) or NB(a,b)• CPUE年トレンド抽出法（要因分析法）• 統計的学習理論(データマイニング)への関心• 線形モデル(統計モデル)から非線形モデルへ• 教師付き学習(応答変数が存在)でのイメージ

CPUE yHiroshi SHONO

exp IC+YEAR y

AREA (YEAR*AREA) y L

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

21

3.資源評価におけるVMStudioの利用
4).Case Study:計算手順およびデータ
1. 出力変数と入力変数を以下のように設定– 出力変数：SBT-CPUE(catch at age4+/1000hooks)– 入力変数：年(1969-2007),月(4-9), 統計海区(4-9:経度の情報含む), 緯度(30°S-50°S: 5°刻み)

2. データを学習用(80%)と検証用(20%)に分割
3. 学習用データで3手法によるルールを作成– サポートベクター回帰– ニューラルネットワーク（誤差逆伝搬法:教師付）– 樹形モデル（CART algorithm, If…thenで2分割）

4. 上記ルールにより検証用データで性能評価
Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

23– 回帰問題(neural network>サポートベクター回帰)– 判別問題(サポートベクターマシン>neural network)• 実データにおけるvalidationを通じた性能比較• 予測値に基づく簡便な要因分析法(CPUE年傾向抽出法)の提案->ブラックボックスからの脱却
Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

22

3.資源評価におけるVMStudioの利用
5).Case Study: 計算に関する条件設定
サポートベクター ニューラルネット
回帰ワーク
カーネルの種類 Active Function(Gaussian Kernel) (Sigmoid function)
上記パラメータ Objective Function(Sigma=1.0)(SS: sum of square)
Slack変数# of hidden layer(C=1.0)(5)
誤差項の大きさ # of repeat=1000(Error=1.0)
Decay=0
Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

樹形モデル
Algorithm(CART)
分岐方法(Info Gain)
複雑度係数
(1)
節点の不純度(0.01)
24

3.資源評価におけるVMStudioの利用
6).Case Study:モデルの性能評価-(1)3.資源評価におけるVMStudioの利用
7).Case Study:モデルの性能評価-(2)• 観測値と（それに対応する）予測値との差異• 観測値と予測値の差異をPearson相関係数と順位相関係数(Kendoll-τ, Spearman-ρ)で判Neural network
Tree regression model
SBT 断 Support vector machine– MAE= Σ|CPUEpre(i)-CPUEobs(i)|/N （絶対誤差）– MSE=Σ(CPUEpre(i)-CPUEobs(i))2/N （平方誤差）
SBT
data
Training
Testing
Total

data Pearson KendollτSpearmanρ Pearson Kendollτ SpearmanρPearson KendollτSpearmanρ
Training 0.910 0.498
0.676 0.937 0.564
0.724 0.565 0.362
0.490
Testing 0.829* 0.474* 0.653* 0.821 0.470
0.637 0.61 0.374
0.516
Total
0.894 0.494
0.672 0.913 0.544
0.706 0.574 0.364
0.494

Suppot vector machine
Neural network
Tree regression model
MSE
MAE
MSE
MAE
MSE
MAE
103
8.07
67.5
5.81
373
12.7
175*
9.48
188
9.27*
339
12.1
118
8.35
91.76.5036612.6注）*印(黄色塗りつぶし)のマイニング手法の性能が一番良い（値が小さいほど性能が良い）
Hiroshi SHONO

数理システムユーザーコンファレンス2010 （東京：六本木ヒルズ） 2010/11/19 (Friday)

25注）*印(黄色塗りつぶし)のマイニング手法の性能が一番良い（値が大きいほど性能が良い）
Hiroshi SHONO

3.資源評価におけるVMStudioの利用
8).予測値ベースの年トレンド抽出法• 要因分析が非常に難しい→ブラックボックス• 簡便な手法を提案（→質的説明変数に着目）
preCPUEyear

NmNl1
1 Na 1
preCPUEyear, month, area, latitude
Nm month 1 Na area 1 Nl latitude 1

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

26

4.まとめ
1).NUOPTに関する個人的な印象• 大規模問題に有利な内点法を使用している
ためか、初期値依存度が低く、性能が良い！
– c.f. 準ニュートン法(収束は早いが初期依存性大)– c.f. シンプレックス法(収束が遅くて時間がかかる)• 数式に似た記述法(文法)が非常に分かり易く、制約条件や解の制限等の追加記述も容易！
• パラメータの区間推定機能を装備して欲しい（Hessian, Profile-Likelihood, Bootstrap etc.）
Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

27

Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

28

4.まとめ2).VMStudioに関する個人的な印象

4.まとめ3).統計ソフトに関する個人的な印象• 装備しているデータマイニング手法が幅広く、各々の手法に関する機能やオプションも多彩

Softwar Advantage
Disadvantage
eS-Plus/R カスタマイズが得意・ インタプリタゆえに、最先端の機能を装備 実行速度が遅い*SAS
大規模な定型的処理 ライセンス料が高額
SPSS
ユーザーフレンドリー CUIが使いにくくモデなGUI機能を有する ル選択などが難しい*-性能向上(Win7/64bit/12GBメモリ)による解決– e.g. SVMやSVM-Rに関する多彩なカーネル関数– e.g. Naïve-Bayes・Boosting等最先端機能を装備– e.g. 複数手法の結果をまとめるモデル統合機能• GUIがユーザーフレンドリーで操作性が高い• NUOPTの最適化アルゴリズムおよびTrial andErrorによる停止規則などを取り入れてほしい
Hiroshi SHONO

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

29

S-Plus ・Validation/Rに無い機能も(多変量GARCH等)
R
Hiroshi SHONO

・最先端の機能が多く存在するがバグも多い
数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)30

A.付録（有限混合分布の理論的背景）
1).カイ二乗バー分布(Shapiro,1985)

END

特別な場合に 2log n ( n:尤度比検定統計量)の漸近分布がカイ二乗分布の重み付き平均の分布(カイ二乗バー分布)になる

Thank you for your kind attention!

例： X 1 ,L , X n (i.i.d.)〜 p.d . f .
【謝辞:数理システム株式会社の皆様方】
特に田辺様 、佐藤様、橋本様、徐様、中園様、田澤様、小木様、福田様、友廣様

0j1,成分数 3つ
数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

Hiroshi SHONO

31A.付録（有限混合分布の理論的背景）
2). 罰金付き最小距離法(Chen,1993)
mの推定量mˆ min[min{Fn ( x), F (x| )) cnm

d ( Fn ( x), F (x| )) op (cn ),*cno(1), cnmj 1logj}]

Hiroshi SHONOmj 1j

12 logmj 1jg j ( x)分布が既知の場合
162012211322

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

32

B.補遺1).サポートベクターマシンのイメージ
サポートベクターマシンのコンセプト(計算手順)

0, cn : 実数列

F (x| ) : 分布関数,. Fn ( x) : 経験分布関数 (例:cn 1/ n)c.f. Minimum Distance (Henna,1985)mˆ

2min m N | minm
( n)n2( n)
Hiroshi SHONO

,( n)
n2

F ( x | ) Fn ( x) dFn ( x )• 超平面（1本の直線）で分離したい→How?
1. カーネル関数 e.g. K(x, y)=exp(-|x-y|2/σ2) etc.
2. 高次元空間へ移す i.e. 地球から月に移動( n)exp{ 2 2 ( n)}n20( n),n 1

数理システムユーザーコンファレンス2010（東京：六本木ヒルズ） 2010/11/19 (Friday)

33

B.補遺2).サポートベクターマシンのイメージ(続)Hiroshi SHONO

3. マージン最大化を用いて一意に分離可能（道路の幅を広くする）• 高次元空間での操作を元の空間からカーネ
ル関数を通じて行う（月での操作を地球上からリモートで行う）• カーネルトリック（検証データに対する予測性能が非常に良い）
数理システムユーザーコンファレンス2010 （東京：六本木ヒルズ） 2010/11/19 (Friday)

35

Hiroshi SHONO

数理システムユーザーコンファレンス2010 （東京：六本木ヒルズ） 2010/11/19 (Friday)

34