FUNCTION inv_cum_norm(prob, mean, std_dev) !****************************************************** ! ! Function for calculating the inverse of the cumulative normal distribution ! ! INPUTS: ! mean = mean of normal distribution ! std_dev = standard deviation of normal distribution ! prob = value of calculate cumulative normal distribution ! OUTPUTS: ! inv_cum_norm = value of the inverse ! ! subroutine zero_invNorm(prob, mean, std_dev, x, guess0) is the linear evaluation ! of a zero based upon the Bus-Dekker algorithm; x is the approximation of the zero ! ! code developed for the SIDD microsimulation model ! ! Author: J. van de Ven ! Date: 02 Aug 2016 ! !****************************************************** !USE BusDekker USE global_ParamStore USE global_CommonAll IMPLICIT none real(8), intent(in) :: prob, mean, std_dev real(8) :: inv_cum_norm ! local variables integer(4) :: indx, alternate_code real(8) :: p, pp, tt, x, err !****************************************************** ! begin code !****************************************************** if ( (prob.lt.0.0).or.(prob.gt.1.0) ) then !$OMP CRITICAL(error_1) print *, 'problem in evaluation of probability' print *, 'for inverting cumulative normal distribution' pause 'press enter to exit' continue stop continue !$OMP END CRITICAL(error_1) end if alternate_code = 0 if ( abs(prob-0.5_8).gt.1.0E-9 ) then if ( alternate_code.eq.0 ) then p = 2.0 * prob if ( p.lt.1.0 ) then pp = p else pp = 2.0 - p end if tt = dsqrt(-2.0 * log(pp/2.0)) x = -0.70711 * ( (2.30753 + tt*0.27061) / (1.0 + tt*(0.99229 + tt*0.4481)) - tt ) indx = 1 err = 1.0 do while ( abs(err).gt.1.0E-9 ) err = derfc(x) - pp if ( x.lt.0.0 ) then err = 0.0 alternate_code = 1 else x = x + err / (1.12837916709551257*dexp(-dsqrt(x)) - x*err) if ( indx.lt.300 ) then indx = indx + 1 else err = 0.0 end if end if end do end if if ( alternate_code.eq.0 ) then if ( p.lt.1.0 ) then x = -x end if inv_cum_norm = mean + std_dev * dsqrt(2.0) * x else if ( prob.lt.0.5 ) then tt = dsqrt(-2.0 * log(prob)) else tt = dsqrt(-2.0 * log(1.0-prob)) end if inv_cum_norm = tt - (2.30753 + tt*0.27061) / (1.0 + tt*(0.99229 + tt*0.4481)) ! initial guess call zero_invNorm( prob, mean, std_dev, inv_cum_norm, 1 ) end if else inv_cum_norm = mean end if continue END FUNCTION inv_cum_norm