Individual Based Model

/*
Individual Based Model

RULE:

An event A that occurs with probability P is equivalent to

1) Generate a random number X that distributes uniformly between 0 and 1.

2) If X < P, we assume the event A has occured. Otherwise, the event A did not occur.

**************************************

Birth-death model
dn/dt=(a-b)n


We assume that the birth rate β is the probability that an
individual gives birth to an offspring and that the death
rate δ the probability that the individual dies within a unit
time. We also assume that within a short time interval Δt,
only one of the following three cases occurs mutually exclusively;
an individual

1) gives birth to an offspring,
2) dies,
or
3) neither gives birth nor dies.

This stochastic birth-death process could be implemented using the
algorithm with a constant time interval Δt small enough.

For all individuals repeat
1) Give birth to a new individual with probability βΔt.
2) Remove this individual with probability δΔt.

*/

#include 
#include 
#include 

#define ALPHA 0.03
#define BETA 0.02
#define INTV 10
#define DT 0.1  /*h 分割数*/

#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

double ran2();

unsigned long int idum =1024;

main()
{
  int pop_size, new_indiv, dead_indiv; /*population size*/
  int i,step; /*j;step number*/
  double prob1,prob2, ran; /*for probabilities*/
    
  prob1=ALPHA*DT;
  prob2=BETA*DT;
  
  pop_size=10; /*initial population size*/
  
  for(step=0;step<=1000;step++)
    {
      if(step % INTV==0)
	printf("%f %d\n",step*0.01, pop_size);
      new_indiv=0;
      dead_indiv=0;
      for(i=1;i<=pop_size;i++)
	{
	  ran = ran2(&idum);
	  if (ran < prob1)
	    new_indiv++;
	  else if (prob1 < ran && ran < prob1 + prob2)
	    dead_indiv++;
	  if (pop_size < 0) pop_size = 0; /*pop.size cannot be negative */
	}
      pop_size+=(new_indiv-dead_indiv);
    }
  return 0;
}

double ran2(long *idum)
{
  int j;
  long k;
  static long idum2=123456789;
  static long iy=0;
  static long iv[NTAB];
  float temp;

  if(*idum <= 0){
    if(-(*idum) < 1) *idum=1;
    else *idum = -(*idum);
    idum2=(*idum);
    for(j=NTAB+7; j>=0; j--){
      k=(*idum)/IQ1;
      *idum=IA1*(*idum-k*IQ1)-k*IR1;
      if(*idum < 0) *idum += IM1;
      if(j < NTAB) iv[j] = *idum;
    }
    iy=iv[0];
  }

  k=(*idum)/IQ1;
  *idum=IA1*(*idum-k*IQ1)-k*IR1;
  if(*idum < 0) *idum += IM1;
  k=idum2/IQ2;
  idum2=IA2*(idum2-k*IQ2)-k*IR2;
  if(idum2 < 0) idum2 += IM2;
  j=iy/NDIV;
  iy=iv[j]-idum2;
  iv[j]=*idum;
  if(iy<1) iy += IMM1;
  if((temp = AM*iy) > RNMX) return RNMX;
  else return temp;
}
出力結果
Back to C Language