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.


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

Immigration-Emigration model
dn/dt=(a-b)

We also assume that the change of the population size within time
interval Δt is at most ±1. This means that only one of the
following three cases occurs mutually exclusively during Δt;

1) an individual joins,
2) an individual quits,
3) no change, with probability αΔt, βΔt, 1 - αΔt - βΔt,

respectively.

This would be satisfied if we assume the interval Δt is so small that
the population size changes ±1.


According to the above rule,
1) we draw a uniformly distributed random number X,
2) If X falls between 0 and αΔt the population size is incremented
by one,
or
3) If X falls between αΔt and αΔt + βΔt, it is decremented by one,
and
4) Otherwise, the population size does not change.

*/

#include 
#include 
#include 

#define ALPHA 0.3
#define BETA 0.2
#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; /*population size*/
  double step; /*j;step number*/
  double prob1,prob2, rand_uniform; /*for probabilities*/

  prob1=ALPHA*DT;
  prob2=BETA*DT;
  
  pop_size=10; /*initial population size*/
  
  for(step=0.0;step<1000.0;step++)
    {
	printf("%f %d\n",step*0.1, pop_size);
      
      rand_uniform = ran2(&idum);
      if (rand_uniform < prob1)
	pop_size++;
      else if (prob1 < rand_uniform && rand_uniform < prob1 + prob2)
	pop_size--;
      if (pop_size < 0) pop_size = 0; /*pop.size cannot be negative */
    }
  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