DrawFromDistri.f90

Go to the documentation of this file.

!==============================================================================
!
!  Name   :  DrawFromDistri
!
!  Purpose:  draws random value from  binomial, poisson or normal distribution
!
!  Remarks:  contains the SUBROUTINEs
!
!            DrawFromDistri (<SeedsDispFromThisCellAndSpec<SeedsDispersedFromThisCell<
!                              LocalDynOneTimeStep<SpatialDynOneTimeStep<TimeLoop<TreeMig)
!
!            DrawFromPoissonDist }   (<GetVariate<GetBioclim<GetEnvFactors<
!            DrawFromBinomDist   }      SpatialDynOneTimeStep<TimeLoop<TreeMig AND/OR
!            DrawFromNormalDist  }      DrawFromDistri<SeedsDispFromThisCellAndSpec<
!                                       SeedsDispersedFromThisCell<LocalDynOneTimeStep<
!                                       SpatialDynOneTimeStep<TimeLoop<TreeMig)
!
!==============================================================================
! design          : H. Lischke, N. Zimmermann
! author(s)       : H. Lischke, N. Zimmermann
! implementation  : H. Lischke, N. Zimmermann
! cleaner         : T.J. Loeffler
! copyright       : (c) 1999, 03 by H. Lischke
!==============================================================================
 
!===================================================================== 
!===============================================================
SUBROUTINE drawfromdistri(n, p, ix)
 
    ! <CALL modules for TYPE definitions and +- fixed variables>---------------------------------------------------------------
    IMPLICIT NONE
 
    ! <Passed variables>--------------------------------------------------------------------------------
    INTEGER, INTENT(in)  :: n      ! parameter
    REAL, INTENT(in)  :: p      ! parameter
    INTEGER, INTENT(out) :: ix     ! drawn value
 
    ! <Local variables>---------------------------------------------------------------------------------
    REAL :: randomvalue, & ! [0,1] random value
            startdist, mu, sigma, rx
 
    ! <Here the SUBROUTINE starts>**********************************************************************
    if (p .ge. 1.0) then
        ix = n
        return
    end if
    if ((p  <=  0.0) .or. (n  <=  0)) then
        ix = 0.0
        return
    end if
 
    call random_number(randomvalue)
    startdist = (1.0 - p)**n
    if ((startdist == 0.0) .or. (startdist == 1.0) .or. (n .gt. 50)) then
        ! for n->infinity, the binomial distribution is approximated by a poisson or normal distribution 
        mu = p*n 
        if (mu  <  20) then ! for small p, the Binomial Distribution is approximated by a  Poisson distribution
            call drawfrompoissondist(mu, ix, randomvalue)
        else  ! otherwise by a normal distribution
            sigma = sqrt(mu)
            call drawfromnormaldist(mu, sigma, rx, randomvalue)
            ix = rx
        end if
    else
        call drawfrombinomdist(n, p, ix, randomvalue)
    end if
end SUBROUTINE drawfromdistri
 
! -------------------------------------------------------------------------------------------------
!===================================================================== 
!===============================================================
SUBROUTINE drawfrompoissondist(mu, x, randomvalue)
 
    ! <CALL modules for TYPE definitions and +- fixed variables>---------------------------------------------------------------
    IMPLICIT NONE
 
    ! <Passed variables>--------------------------------------------------------------------------------
    REAL, INTENT(in)  :: mu, &   ! parameter                     [=..]
                         randomvalue     ! [0,1] random value            [=..]
    INTEGER, INTENT(out) :: x               
 
    ! <Local variables>---------------------------------------------------------------------------------
    double precision :: distrival, distrisum, ri, fak, mud
    INTEGER          :: i
 
    ! <Here the SUBROUTINE starts>**********************************************************************
    x = 0
    mud = mu
    distrival = exp(-mud)
    distrisum = distrival
 
    do i = 1, 10000
        ri = i
        if ((distrisum .ge. randomvalue) .or. (distrisum .ge. 1.0)) exit
        fak = mu/ri
        distrival = distrival*fak
        distrisum = distrisum + distrival
    end do
 
    x = ri - 1
end SUBROUTINE drawfrompoissondist
! ****<Here the SUBROUTINE ends>*******************************************************************
 
!===================================================================== 
!===============================================================
SUBROUTINE drawfrombinomdist(n, p, x, randomvalue)
 
    ! <CALL modules for TYPE definitions and +- fixed variables>---------------------------------------------------------------
    IMPLICIT NONE
 
    ! <Passed variables>--------------------------------------------------------------------------------
    INTEGER, INTENT(in)  :: n              ! number of species
    REAL, INTENT(in)  :: p, &  ! probability                      [=..]
                         randomvalue    ! [0,1] random value               [=..]
    INTEGER, INTENT(out) :: x               
 
    ! <Local variables>---------------------------------------------------------------------------------
    INTEGER          :: i
    REAL             :: pTq
    double precision :: distrival, distrisum, rn, ri, fak, fak1, fak2
 
    ! <Here the SUBROUTINE starts>**********************************************************************
    x = 0
    distrival = (1.0 - p)**n
    distrisum = distrival
    rn = n
    ptq = p/(1.0 - p)
    do i = 1, n
        if (distrisum .ge. randomvalue) exit
        ri = i
        fak1 = (rn - ri + 1.0)   
        fak2 = fak1/ri  ! this gives a product, from which the binomial coefficient n!/(i!(n-i)!) can be calculated, with (n-i+1)! = n! / (n-i)!
        fak = ptq * fak2 ! by multiplying the product over all pTq, i.e. pTq**i with the first distrival, i.e. (1.0 - p)**n, gives the second part of the binomial distribution density
        distrival = distrival*fak
        distrisum = distrisum + distrival
        x = i  ! this is the sampled value
    end do
end SUBROUTINE drawfrombinomdist
!****<Here the SUBROUTINE ends>********************************************************************
 
!===================================================================== 
!===============================================================
SUBROUTINE drawfromnormaldist(mu, s, rx, randomvalue)
 
    ! <CALL modules for TYPE definitions and +- fixed variables>---------------------------------------------------------------
    IMPLICIT NONE
 
    ! <Passed variables>--------------------------------------------------------------------------------
    REAL, INTENT(in)  :: mu, &  ! mean               [=..]
                         s, &  ! standard deviation [=..]
                         randomvalue    ! [0,1] random value [=..]
    REAL, INTENT(out) :: rx             ! REAL output value  [=..]
 
    ! <Local variables>---------------------------------------------------------------------------------
    INTEGER          :: i, index, istart
    REAL             :: xtransf, &
                        v1, f1, &
                        v2, f2, &
                        nd(0:201)     !                                       [=NormDistFuncArr(index)]
 
    ! <Here the SUBROUTINE starts>*********************************************************************
    ! <First, initialize the Standard Normal Distribution (mu=0, sigma=1)  array nd(..), in the range of -5 to +5 >------------
    ! it has 202 (0:201) elements, element 100 corresponds to the median and mean, namely 0.  0 and 201 correspond to -5, 5 respectively.          
    nd = (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0000106885, 0.0000133457, &
           0.0000166238, 0.0000206575, 0.0000256088, 0.0000316712, 0.0000390756, 0.0000480963, 0.0000590589, &
           0.000072348, 0.0000884173, 0.0001078, 0.00013112, 0.000159109, 0.000192616, 0.000232629, &
           0.000280293, 0.000336929, 0.000404058, 0.000483424, 0.000577025, 0.000687138, 0.000816352, &
           0.000967603, 0.00114421, 0.0013499, 0.00158887, 0.00186581, 0.00218596, 0.00255513, 0.00297976, &
           0.00346697, 0.00402459, 0.00466119, 0.00538615, 0.00620967, 0.00714281, 0.00819754, 0.00938671, &
           0.0107241, 0.0122245, 0.0139034, 0.0157776, 0.0178644, 0.0201822, 0.0227501, 0.0255881, 0.0287166, &
           0.0321568, 0.0359303, 0.0400592, 0.0445655, 0.0494715, 0.0547993, 0.0605708, 0.0668072, 0.0735293, &
           0.0807567, 0.088508, 0.0968005, 0.10565, 0.11507, 0.125072, 0.135666, 0.146859, 0.158655, 0.171056, &
           0.18406, 0.197663, 0.211855, 0.226627, 0.241964, 0.257846, 0.274253, 0.29116, 0.308538, 0.326355, &
           0.344578, 0.363169, 0.382089, 0.401294, 0.42074, 0.440382, 0.460172, 0.480061, 0.5, 0.519939, &
           0.539828, 0.559618, 0.57926, 0.598706, 0.617911, 0.636831, 0.655422, 0.673645, 0.691462, 0.70884, &
           0.725747, 0.742154, 0.758036, 0.773373, 0.788145, 0.802337, 0.81594, 0.828944, 0.841345, 0.853141, &
           0.864334, 0.874928, 0.88493, 0.89435, 0.9032, 0.911492, 0.919243, 0.926471, 0.933193, 0.939429, &
           0.945201, 0.950529, 0.955435, 0.959941, 0.96407, 0.967843, 0.971283, 0.974412, 0.97725, 0.979818, &
           0.982136, 0.984222, 0.986097, 0.987776, 0.989276, 0.990613, 0.991802, 0.992857, 0.99379, 0.994614, &
           0.995339, 0.995975, 0.996533, 0.99702, 0.997445, 0.997814, 0.998134, 0.998411, 0.99865, 0.998856, &
           0.999032, 0.999184, 0.999313, 0.999423, 0.999517, 0.999596, 0.999663, 0.99972, 0.999767, 0.999807, &
           0.999841, 0.999869, 0.999892, 0.999912, 0.999928, 0.999941, 0.999952, 0.999961, 0.999968, 0.999974, &
           0.999979, 0.999983, 0.999987, 0.999989, 0.999991, 0.999993, 0.999995, 0.999996, 0.999997, 0.999997, &
           0.999998, 0.999998, 0.999999, 0.999999, 0.999999, 0.999999, 1.0, 1.0, 1.0, 1.0/)
    rx = 0
    ! if random value == 0 then a value at the very left tail of the distribution is chosen
    if (randomvalue == 0.0) then
        rx = -5.0 * s + mu
        return
    end if
    !   0.5 is the median of the distribution. therefore for smaller random values the first half, for larger random values the second half of the array can be searched
    if (randomvalue .gt. 0.5) then
        istart = 100
    else
        istart = 0
    end if
    !  search in the table where the random value belongs to
    do i = istart, istart + 101
        if (nd(i) .ge. randomvalue) then
            index = i
            exit
        end if
    end do
    !  interpolate 
    f1 = nd(index - 1)
    f2 = nd(index)
    ! transform the indices from 1 to 200 to (0,1) normal distribution. 
    v1 = (index - 1)/20.0 - 5.0
    v2 = (index)    /20.0 - 5.0
    xtransf = v1 + (v2 - v1)*(randomvalue - f1)/(f1 - f2)
    ! transform the (0,1) normal distribution to the (mu,sigma) normal distribution
    rx = xtransf*s + mu
 
end SUBROUTINE drawfromnormaldist
!done