植物版レッドリストの絶滅リスクの評価方法についての私の提案

松田裕之（東大海洋研）

　このMathematicaプログラムは“「植物レッドリスト」（種子植物・シダ植物）の調査・判定方法と判定結果の特徴”日本植物分類学会絶滅危惧植物問題専門第一委員会（文責：矢原徹一・加藤辰己）という文書を読んで、以前私が矢原氏に示したプログラムを改良したものである。実際の絶滅リスクの評価は矢原徹一氏と新潟大学の加藤辰巳氏の研究グループが議論して以前の私の提案を一部修正し、他の言語でプログラムを書き直して行った。したがって、下記のプログラムと完全に同一ではない。

　必要なファイルは入力用テキストファイル"redlist.dat"である。計算結果は画面と出力用テキストファイル"result.dat"に表示される。入力用テキストファイルは以下のように１行目が各欄の説明、２行目以降に各種の絶滅メッシュ数、減少率0.01，0.1，0.5，0.9，１のメッシュ数、個体数が1000、100、10，１以上のメッシュ数、不明メッシュ数の順にタブ区切りで１行１種で入力されている。

出力は種番号、推定個体数、存在メッシュ数（個体数不明メッシュを除く）、平均減少率、絶滅までの平均年数、10年後、20年後、・・・100年後の絶滅確率（1000回の試行のうちで絶滅した頻度）の順に計算している。減少率の推定方法はこのプログラムでは矢原・加藤案と異なり、最初の10年間の減少率を1000回の試行で平均して求めている。

　サクラソウ（仮にここでは種番号を１とした）の例では出力は

となる。

/*Extinction risk estimation program for Japanese Plant Redlist */

/* Programmed by Hiroyuki Matsuda, September 3, 1997 */

#include <stdio.h>

#include <math.h>

#include <limits.h>

#include <stdlib.h>

static unsigned long seed=1;

#define rnd() (seed *=69069UL, seed/(ULONG_MAX +1.0))

#define Ncrit 1.

#define Kmax 100 /* number of iterations */

#define Imax 1000 /* maximum number of patches */

#define Tmax 30 /* normally we assume Tmax=10 decades (=100yrs)*/

#define Smax 10 /* Number of species */

FILE *fp1, *fp2;

/***************************************************/

float u(int i) /*class of population decrease*/

{

float reply;

float uu[]={0.001, 0.01, 0.1, 0.5, 1., 1.};

if (i==5) reply = 1.;

else if (i==0) reply=uu[0];

else reply=uu[i-1]+(uu[i]-uu[i-1])*rnd();

return reply;

}

/***************************************************/

void initial(int ll[], float n[], int r[], float *np, int *t9, int *lp,

int *sp)

{

int i, j, p[6], pp[6], l[5];

int p0,p1,p2,p3,p4,p5, l1,l2,l3,l4,l5;

/******** frequency distribution of the rate of decrease*/

fscanf(fp1,"%d %d %d %d %d %d %d %d %d %d %d %d",

&sp, &p0,&p1,&p2,&p3,&p4,&p5,&l1,&l2,&l3,&l4,&l5);

p[0]=p0; p[1]=p1; p[2]=p2 ;p[3]=p3; p[4]=p4; p[5]=p5;

l[0]=0; l[1]=l1; l[2]=l2; l[3]=l3; l[4]=l4;

/* p[0]=0+1; p[1]=0; p[2]=0; p[3]=0; p[4]=1; p[5]=0;

l[0]=0; l[1]=0; l[2]=1; l[3]=0; l[4]=0;*/

/******** cumulative frequency */

pp[0]=p[0]; for(i=1; i<6; i++) {pp[i]=pp[i-1]+p[i];}

for(j=1;j<pp[0]+1;j++) r[j]=0;

for(i=0;i<5;i++){for(j=pp[i]+1;j<pp[i+1]+1;j++){r[j]=i+1;}}

/******** cumulative frequency of patch size desitribution */

ll[0]=0; for(i=1; i<5; i++) {ll[i]=ll[i-1]+l[i];}

/******** initial patch size of each class */

n[1]=sqrt(10.)*1000.; for(i=2;i<5;i++) {n[i]=n[i-1]/10.;}

/******** total population size at present */

*np=(float) l[1]*n[1] + (float) l[2]*n[2] +

(float) l[3]*n[3] +(float) l[4]*n[4];

/******** total number of patches */

*t9=p[0]+p[1]+p[2]+p[3]+p[4]+p[5];

*lp=l[1]+l[2]+l[3]+l[4];

}

/************************* BootStrap **********************/

void boot(int is[], float n[], int ll[], int r[], float np, int lp,

int t9, int sp, float *decm, float *d_cr, float *d_en, float *d_vu)

{

int i, j, k, t;

float patch[Imax]; /* patch size in patch i */

float ns,np1,dec;

/* is[t]: No. of sample in which the population is extinct

until the t-th decades */

*decm=0.; *d_cr=0.; *d_en=0.; *d_vu=0.;

for(t=0; t<= Tmax ; t++){is[t]=0;}

/* patch[i]: patch size of the i-th patch in the t-th decade*/

/* if the patch size is less than 1, this patch is recognized as extinct */

/* if the total size is 0, the population is recognized as extinct */

for(k=0; k< Kmax ; k++){

/* initial patch isze in each patch */

for(j=1;j<=4;j++){

for(i=ll[j-1]+1;i<=ll[j];i++){

patch[i]=n[j];

}

}

for(t=1;t<=Tmax;t++){

for(i=1;i<=lp; i++){

if(patch[i]>=1.)

patch[i]=u(r[(int)(rnd()*(float)t9)+1])*patch[i];

if(patch[i]<1.) patch[i]=0.;

}

ns=0.; for(i=1;i<=lp;i++){ns=ns+patch[i];}

if(ns<Ncrit) {is[t]=is[t]+1;}

if(t==1){

np1=0.;

for(i=1;i<=lp;i++) np1=np1+patch[i];

dec=1.-np1/np;

}

}

*decm=*decm+dec;

if(dec>0.2)*d_vu=*d_vu+1;

if(dec>0.5)*d_en=*d_en+1;

if(dec>0.8)*d_cr=*d_cr+1;

}

*decm=*decm/Kmax;

*d_vu=*d_vu/Kmax; *d_en=*d_en/Kmax; *d_cr=*d_cr/Kmax;

}

/********************************************************/

void result(int sp, int is[], float np, int lp,

float decm, float d_cr, float d_en, float d_vu)

{

int t;

float tt;

/* Expected extinction time (years) */

tt=0.;

for(t=1;t<=Tmax;t++){

tt=tt+(float)((is[t]-is[t-1])*t*10)/(float)(Kmax);}

tt=tt+(float)(Tmax*10.+5.)*(Kmax-is[Tmax])/(float)(Kmax);}

if(tt<=0.) tt=100.; else tt=tt-5.;

fprintf(fp2,"%5d, %7d, %5d, %7.3f, %10.3f, ",sp, (int)np, lp, decm, tt);

printf("%d, ",sp);

/* Cum. probability of extinction within 10 yrs,..., 100 yrs ***/

for(t=1;t<=Tmax;t++) {

printf("%.2f, ", (float)(is[t])/(float)(Kmax));

fprintf(fp2,"%.3f, ", (float)(is[t])/(float)(Kmax));

if(t/10*10==t) printf("\n");

}

printf("\n %.2f, %d, %.2f, %.2f, %.2f, %.2f, %.2f\n",

np, lp, decm, d_cr, d_en, d_vu, tt);

fprintf(fp2,"\n");

}

/********************************************************/

main()

{

int i, is[Tmax+1], sp, lp, t9, r[Imax], ll[5];

float np, n[Imax], decm, d_cr, d_en, d_vu;

fp1=fopen("red.dat","r");

fp2=fopen("result.dat","w");

fprintf(fp2,"Critical total population size is %d\n", (int)Ncrit);

for(i=1;i<=Smax;i++){

initial(ll, n, r, &np, &t9, &lp, &sp);

boot(is, n, ll, r, np, lp, t9, sp, &decm, &d_cr, &d_en, &d_vu);

result(sp, is, np, lp, decm, d_cr, d_en, d_vu);

}

printf("Critical total population size is %d\n", (int)Ncrit);

}

なお、Mathematicaプログラムを以下に示す。こちらは減少率と絶滅待ち時間は計算しない。

(*class of population decrease*)

uu[0]=0.001;uu[1]=0.01;uu[2]=0.1;uu[3]=0.5;uu[4]=1;uu[5]=1;

u[i_]:=If[i==0,uu[0],uu[i-1]+(uu[i]-uu[i-1])*Random[]];

fi="a:\redlist.txt";fo="a:result.txt";

k9=1000; (* the number of iteration *)

t8=10; (* time limit is 10 decades *)

Read[fi,{Word,Number,Number,Number,Number,

Number,Number,Number,Number,Number,Word}];

Do[ str=Read[fi,{Word,Number,Number,Number,Number,Number,

Number,Number,Number,Number,Number,Number}];

l[0]=0;

Do[ p[i]=str[[i+2]],{i,0,5}];

Do[ l[i]=str[[i+7]],{i,1,4}];

(* cumulative frequency *)

Do[pp[i]=Sum[p[j],{j,0,i}],{i,0,5}];

(* the number of total sample of decreasing rate *)

t9=Sum[p[j],{j,0,5}];

(* artificial time series of decreasing rate *)

Do[r[t]=0,{t,1,pp[0]}];

Do[Do[r[t]=i,{t,pp[i-1]+1,pp[i]}],{i,1,5}];

(* cumulative frequency of patch size desitribution *)

ll[0]=0;Do[ll[i]=ll[i-1]+l[i],{i,1,4}];

(* patch size of each class *)

n[4]=N[Sqrt[10]];Do[n[i]=n[i+1]*10,{i,3,1,-1}];

(* total population size at present *)

np=Sum[l[i]*n[i],{i,1,4}];

(* total number of patches *)

lp=Sum[l[i],{i,1,4}];

(* BootStrap *)

(* is[t]: No. of sample in which the population is extinc until the t-th decades *)

Do[is[t]=0,{t,0,t8}];

(* patch[i]: patch size of the i-th patch in the t-th decade*)

(* if the patch size is less than 1, this patch is recognized as extinct *)

(* if the total size is 0, the population is recognized as extinct *)

Do[ Do[ Do[ patch[i]=n[j],{i,ll[j-1]+1,ll[j]}],{j,1,4}];

Do[ Do[ pti=u[r[Round[Random[]*t9+1]]]*patch[i];

patch[i]=If[pti<1,0,pti],

{i,1,lp}];

n1[k,t]=Sum[patch[i],{i,1,lp}];

is[t]=is[t]+If[n1[k,t]>1,0,1],

{t,1,t8}],

{k,1,k9}];

(* Expected extinction time (years) *)

tt=N[Sum[(is[t]-is[t-1])*t,{t,1,t8}]*10/k9-5];

(* Cumulative probability of extinction within 10 yrs, 20 yrs, ..., 100 yrs *)

pt=Table[N[is[t]/k9],{t,1,t8}];

(* Mean rate of population decrease in the deterministic model *)

s=1-(1/np)^(1/tt);

Write[fo,{sp,np,lp,tt,s,pt,Table[p[i],{i,1,5}],

Table[l[i],{i,1,4}]}];

Print[sp," ",np," ",lp," ",pt],{sp,1,10}];

Close[fi];Close[fo];