2chの人たちがした、RMP シミュレーション

捕鯨問題議論スレッド 36頭目(2013年ごろ…)で、Forty-Fourth Report of the International Whaling Commission のレポート p158-p166 の Fortran のコードとデータ・ファイルを、匿名の投稿者が再現して、 Catch.zip として、ファイルアップロードサイトに上げて(今は元ページはなくなっています…)、そしてそれを誰かが MSY の最小値を0固定だったものから変化させられる修正版 http://ideone.com/lCeftw (Catch_v2.f) にして上げました。それを使って行われた、シミュレーションの会話録です。当該プログラムをコンパイルし、計算結果を再現できました。

348 ：名無しさん＠お腹いっぱい。：2013/06/06(木) 12:15:09.90 ID:9MrPp+vO
RMPプログラムを入力したんで上げとくね
http://fast-uploader.com/file/6926043723658/
44th Annual report of the IWC
http://iwc.int/index.php?cID=1918&cType=document&download=1
p158～p166

Catch.exe　実行ファイル
Catch.f 　ソースリスト
CLC.PAR パラーメータファイル
CLC.DAT データファイル

Catch.exeを実行するとパラメータファイル、データファイルを読み込み、結果をCLC.OUTに書き込む。
パラメータファイル、データファイルは167ページのサンプルデータそのまま。

361 ：名無しさん＠お腹いっぱい。：2013/06/08(土) 21:25:58.61 ID:tqiHMkVs
別の数字を入れてみたいと思った時に、分からないのが資源量の次の数字。
167ページのInput data file CLC.DATでは
1988 27150.0 100.
となっているけど、最後の100って何？
上のExampleで One abundance estimate is available: 27150 in 1988.
The variance of the log estimate is 0.01 and hence information matrix consists of a single entry of value 1/0.01)
何の事やらさっぱり分からないので数学板で聞いてみた。

>511 名前：１３２人目の素数さん[sage] 投稿日：2013/06/05(水) 01:21:21.31
>>>505
>http://iwc.int/cache/downloads/b9t7r0k873wc4gc00wwso00kk/1992%20SC%20REP.pdf
>このページ66の表を参照
>N(CV) = 27150 (0.120)
>との数値があります。前後をそれほど読んではいませんが捕鯨関係だし同一起源のデータでしょう
>※CVは、変動係数（Coefficient of variation: 標準偏差を算術平均で割ったもの)

>weblio: 「対数正規分布」
>対数正規分布は(中略)、 生物群の個体数、(中略)の各現象に広く見られる統計分布であり…

>という事で、クジラ頭数に対して対数正規分布 LN(μ,σ^2) を仮定すると。
>CV = sqrt(V(x))/E(x) = sqrt( e^{σ^2} – 1 ) = 0.120
>σ^2 = ln( 0.120^2 +1) = 0.014 (≒0.01 はココから来ていると思われ)
>参考:Wikipedia 「対数正規分布」の項目「平均・分散」

>あと、CV=0.120 自体は何らかのモデリングによる推定っぽいです

実は数学はよくわからんのだけど、資源量NとCVがあって、資源量の次に1/ln(CV^2+1)を入れればよさそうだと分かった。
回答して下さった数学版の>511さん、ありがとう。

366 ：名無しさん＠お腹いっぱい。：2013/06/09(日) 20:09:42.28 ID:XqHqZIqb
2、3日まえに税金泥棒君がN+のスレで張ってくれてたレポートを利用させてもらうね
http://iwc.int/cache/downloads/6r8jq8llm4cgso0sc0k000w8c/2012%20SC%20REP.pdf

このp37のTable9の数字から捕獲枠を計算してみる。なんか数字が複数あるのでよく分からない。
なので適当にCPIIのトータルは720,054を使って、CPIIIのトータルは514,783を使う。
CVもinternalとwith AVの二種類あるので、両方使ってみる。

資源量 720,054、CV=0.08 → 157
資源量 720,054、CV=0.18 → 31

資源量 514,783、CV=0.09 → 124
資源量 514,783、CV=0.18 → 31

以上のデータをCLC.DATに入力して計算すると、結果は

720,054頭の場合　14,432頭(CV=0.08) 8,565頭(CV=0.18)
514,783頭の場合　 9,931頭(CV=0.09) 6,105頭(CV=0.18)

379 ：名無しさん＠お腹いっぱい。：2013/06/10(月) 20:31:57.05 ID:4V3XG0Ug
>>370
1993年のAllison Cの試算は年間最小1945頭、最大4490頭だったけど
なんか保守的な要素加えてたのかね

381 ：名無しさん＠お腹いっぱい。：2013/06/10(月) 21:26:39.41 ID:yC2zTTLb
>μとP(0)の組を指定された範囲内に満遍なく配置して(>>78)

「μを範囲内に満遍なく配置する」の意味について。
CLC.PARの”MAX MSY % (PYMAX)” がμ、繁殖値の最大値、0.05が指定されている。
CLC.PARの”NO. OF MSY STEPS (PNYSTE)” が満遍なく配置するための数字、200が指定されている。

つまり5％を200で割った数字0.025％刻みで0から5％まで200段階で数字を満遍なく配置するという意味。
(実際は少し違い、下限の数字が0.025の半分の0.0125～最大5-0.0125の4.9875まで0.025％刻み)

0.000125
0.000375
0.000625
0.000875
・
・
0.049375
0.049625
0.049875

次に計算したいのは調査捕鯨で繁殖値がもっと絞られたら捕獲枠がどう変化するか。
ところが残念なことに繁殖値の最大値はパラメーターファイルで入力することになっているが、最低値はプログラムの中で0に固定されている。
従って繁殖値の範囲をもっと絞るためにはプログラムの改造が必要になる。

388 ：名無しさん＠お腹いっぱい。：2013/06/11(火) 00:08:50.05 ID:Yxl6keb1
>>381
修正してみた。これでどう？
http://ideone.com/lCeftw

パラメータはこんな感じで設定
PROBABILITY LEVEL (PPROB) 0.4102
MAX MSY % (PYMAX) 0.05
MIN MSY % (PYMIN) 0.01
NO. OF MSY STEPS (PNYSTP) 200.0
KSTEP (PKSTEP) 0.05
DEPLETION STEP (PDSTEP) 0.01
BIAS MIN (PBMIN) 0.0
BIAS MAX (PBMAX) 1.6667
NO. OF BIAS STEPS (PNBSTP) 100.0
SCALE FACTOR (PSCALE) 4.0
PHASEOUT PERIOD (PHASET) 8.0
PHASEOUT PROPORTION (PHASEP) 0.2
ASSESSMENT CYCLE (PCYCLE) 5.0
INTERNAL PROTECTION LEVEL 0.54
CATCH CONTROL SLOPE (PSLOPE) 3.0

391 ：名無しさん＠お腹いっぱい。：2013/06/11(火) 14:07:03.82 ID:swwZXUbm
>>388
サンクス、サンクス。
どうせだから結果も書いてちょうだい。
調査捕鯨で少なくとも繁殖力が1%よりは大きいと分かったら、2%より、3%よりは大きいと分かったら。それぞれの場合結果がどうなるのか。
または少なくても4％以下だと分かったら、3%以下だと分かったら。
または1％から4％の間、2%から3%の間とか、いろんなパターンでお願い。

392 ：名無しさん＠お腹いっぱい。：2013/06/11(火) 14:32:20.09 ID:c069KU5Q
MAX MSY % (PYMAX) 0.04
MIN MSY % (PYMIN) 0.00
Year: 1994 Catch limit: 401
Year: 1995 Catch limit: 401
Year: 1996 Catch limit: 401
Year: 1997 Catch limit: 321
Year: 1998 Catch limit: 241

MAX MSY % (PYMAX) 0.04
MIN MSY % (PYMIN) 0.01
Year: 1994 Catch limit: 581
Year: 1995 Catch limit: 581
Year: 1996 Catch limit: 581
Year: 1997 Catch limit: 465
Year: 1998 Catch limit: 349

MAX MSY % (PYMAX) 0.04
MIN MSY % (PYMIN) 0.02
Year: 1994 Catch limit: 741
Year: 1995 Catch limit: 741
Year: 1996 Catch limit: 741
Year: 1997 Catch limit: 593
Year: 1998 Catch limit: 445

393 ：名無しさん＠お腹いっぱい。：2013/06/11(火) 14:33:26.71 ID:c069KU5Q
MAX MSY % (PYMAX) 0.05
MIN MSY % (PYMIN) 0.00
Year: 1994 Catch limit: 493
Year: 1995 Catch limit: 493
Year: 1996 Catch limit: 493
Year: 1997 Catch limit: 394
Year: 1998 Catch limit: 296

MAX MSY % (PYMAX) 0.05
MIN MSY % (PYMIN) 0.01
Year: 1994 Catch limit: 678
Year: 1995 Catch limit: 678
Year: 1996 Catch limit: 678
Year: 1997 Catch limit: 542
Year: 1998 Catch limit: 407

MAX MSY % (PYMAX) 0.05
MIN MSY % (PYMIN) 0.02
Year: 1994 Catch limit: 847
Year: 1995 Catch limit: 847
Year: 1996 Catch limit: 847
Year: 1997 Catch limit: 677
Year: 1998 Catch limit: 508

394 ：名無しさん＠お腹いっぱい。：2013/06/11(火) 14:35:11.28 ID:c069KU5Q
MAX MSY % (PYMAX) 0.10
MIN MSY % (PYMIN) 0.00
Year: 1994 Catch limit: 901
Year: 1995 Catch limit: 901
Year: 1996 Catch limit: 901
Year: 1997 Catch limit: 721
Year: 1998 Catch limit: 541

MAX MSY % (PYMAX) 0.10
MIN MSY % (PYMIN) 0.01
Year: 1994 Catch limit: 1111
Year: 1995 Catch limit: 1111
Year: 1996 Catch limit: 1111
Year: 1997 Catch limit: 889
Year: 1998 Catch limit: 667

MAX MSY % (PYMAX) 0.10
MIN MSY % (PYMIN) 0.02
Year: 1994 Catch limit: 1312
Year: 1995 Catch limit: 1312
Year: 1996 Catch limit: 1312
Year: 1997 Catch limit: 1050
Year: 1998 Catch limit: 787

395 ：名無しさん＠お腹いっぱい。：2013/06/11(火) 17:26:04.66 ID:swwZXUbm
>>392-394
ありがとさんです。

>事実として繁殖力の　「　常　識　値　」　を入力した場合と、同じく繁殖力の　「　実　測　値　」　を入力した場合では、

>　入　力　諸　元　の　値　に　差　が　あ　る　んだから、

>そこから　導　か　れ　る　結　果　（　捕　獲　頭　数　）　にも　差　が　出　る　のは

>　　　　ア　　　　タ　　　　リ　　　　マ　　　　エ　　　　♪

2、3年前にトリパン君がこうレスしているのを見て、実際どれくらい差が出るのかと思ったのが始まり。
最初は自分でプログラムを組めないかと思ったけど、尤度計算のあたりが全く分からず挫折。
最近になってやっとソースプログラムに辿りついて、今日ついに目標達成。(中身はまだよく分からんままだが)
なんか肩の荷が下りた感じ。

396 ：名無しさん＠お腹いっぱい。：2013/06/11(火) 19:53:47.72 ID:Yxl6keb1
もう少しパターンを。範囲の広さ狭さでどの程度影響があるか試してみた
MAX MSY % (PYMAX) 0.05
MIN MSY % (PYMIN) 0.04
Year: 1994 Catch limit: 1140
Year: 1995 Catch limit: 1140
Year: 1996 Catch limit: 1140
Year: 1997 Catch limit: 912
Year: 1998 Catch limit: 684

MAX MSY % (PYMAX) 0.06
MIN MSY % (PYMIN) 0.03
Year: 1994 Catch limit: 1109
Year: 1995 Catch limit: 1109
Year: 1996 Catch limit: 1109
Year: 1997 Catch limit: 887
Year: 1998 Catch limit: 665

MAX MSY % (PYMAX) 0.07
MIN MSY % (PYMIN) 0.02
Year: 1994 Catch limit: 1044
Year: 1995 Catch limit: 1044
Year: 1996 Catch limit: 1044
Year: 1997 Catch limit: 835
Year: 1998 Catch limit: 626

397 ：名無しさん＠お腹いっぱい。：2013/06/11(火) 19:54:22.26 ID:Yxl6keb1
MAX MSY % (PYMAX) 0.10
MIN MSY % (PYMIN) 0.09
Year: 1994 Catch limit: 2421
Year: 1995 Catch limit: 2421
Year: 1996 Catch limit: 2421
Year: 1997 Catch limit: 1937
Year: 1998 Catch limit: 1453

MAX MSY % (PYMAX) 0.11
MIN MSY % (PYMIN) 0.08
Year: 1994 Catch limit: 2401
Year: 1995 Catch limit: 2401
Year: 1996 Catch limit: 2401
Year: 1997 Catch limit: 1921
Year: 1998 Catch limit: 1441

MAX MSY % (PYMAX) 0.12
MIN MSY % (PYMIN) 0.07
Year: 1994 Catch limit: 2361
Year: 1995 Catch limit: 2361
Year: 1996 Catch limit: 2361
Year: 1997 Catch limit: 1889
Year: 1998 Catch limit: 1417

398 ：名無しさん＠お腹いっぱい。：2013/06/11(火) 20:08:54.45 ID:Yxl6keb1
それから、CLC.DATの この部分、
abundance data:
No. of sightings estimate 1
(I4,F10.0,10F8.0)
1988 27150. 100.
の sightings estimate 件数とその観測値も結果に影響が出るのではないかな。
プログラム上は、ここに複数の年度とその観測データを置けるようになっています。

ただ、どういうデータを用意すれば良いかが分からない。