* Given different daily mortality rates (predation/accident etc.)

* Given a constant failure rate and a defined number of eliminations then the expected survival rate follows the negative binomial distribution.

* Negative binomial results are as follows:

* mean = pr/(1-p)

* r is number of failures:

* r = 1

* mortality rate = .01

* if p=.99 ->

di "E= " .99/.01

* mortality rate = .015

* if p=.985 ->

di "E= " .985/.015

* mortality rate = .02

* if p=.98 ->

di "E= " .98/.02

* mortality rate = .025

* if p=.975 ->

di "E= " .975/.025

* mortality rate = .03

* if p=.97 ->

di "E= " .97/.03

* mortality rate = .035

* if p=.965 ->

di "E= " .965/.035

* Let's see if the simulation gives us the expected results

clear

set obs 12000

gen mort_rate = mod(_n,6)/200+.01

* This will split the data into 6 groups of 200 with each group having mortality rates from .01 to .035

tab mort_rate

* An indicator that the mosquito is alive

gen living=1

* An indicator that the mosquito survived

gen survived=0

* Sample 500 days, see how many mosquitos survive to live for each day

forv i = 1/500 {

* Draw a single temporary binomial draw to see if this day the mosquitos die

qui gen died = rbinomial(1,mort_rate) if living==1

* If the mosquitos died then change their "living" status

qui replace living=0 if living==1 & died == 1

* Replace survived to indicate the current day if they are still alive

qui replace survived=`i' if living==1

* Displays user feedback telling what round the days currently are on and what percentage of mosquitoes live.

qui sum living

di "Round `i' :" r(mean)

* Drop the temporary variable "died"

drop died

}

hist survived, by(mort_rate)

bysort mort_rate: sum survived

## No comments:

## Post a Comment