# tadpole predation exercise
library(bbmle)
library(emdbook)
data(ReedfrogFuncresp)
d <- ReedfrogFuncresp

plot(d$Killed ~ d$Initial)

d$freqkilled <- d$Killed/d$Initial
plot(d$freqkilled ~ d$Initial)

init <- c(5,10,15,20,30,50,75,100)
mfk <- sapply( init , function(z) mean(d$freqkilled[d$Initial==z]) )
plot(mfk ~ init , ylim=c(0,1) , xlab="Initial number" , ylab="Average freq killed" )

##############
### step 2
##############
# 3 intercept

x <- d$Initial
killed <- rbinom( length(x) , prob=0.5 , size=x )
plot( killed ~ x )


# 4

a <- 0.5
b <- -0.005
killed <- rbinom( length(x) , prob=a+b*x , size=x )

pkilled <- killed/x
plot( pkilled ~ x )

# 2

a <- 0.5
b <- -1/50
killed <- rbinom( length(x) , prob=a*exp(b*x) , size=x )
plot( killed/x ~ x )

# 1

a <- 0.75
b <- -.2
killed <- rbinom( length(x) , prob=a*x*exp(b*x) , size=x )
plot( d$Killed/x ~ x )
curve( a*x*exp(b*x) , from=0 , to=100 , add=TRUE , col="red" )


##############
### step 3
##############

# intercept model
sum( -dbinom(d$Killed,size=d$Initial,prob=0.25,log=TRUE) )

# linear model
sum( -dbinom(d$Killed,size=d$Initial,prob=a+b*x ,log=TRUE) )



#############
### step 4
#############

# intercept
x <- d$Initial
killed <- rbinom( length(x) , prob=0.15 , size=x )

m <- mle2( killed ~ dbinom( prob=a , size=x ) , start=list(a=0.5) )
coef(m)

rep.m <- function(a) {
    y <- rbinom( length(x) , prob=a , size=x )
    m <- mle2( y ~ dbinom( prob=ka , size=x ) , start=list(ka=0.5) , data=list(y=y) )
    coef(m)
}

a.reps <- replicate( 999 , rep.m(0.15) )

mean(a.reps)
sd(a.reps)
hist(a.reps,breaks=100)
lines( c(.15,.15) , c(0,10000) , lwd=2 , col="red" )


###########
## step 5
###########

# intercept
y <- d$Killed
x <- d$Initial

m0 <- mle2(  )

mricker <- mle2( y ~ dbinom( prob=a*x*exp(b*x) , size=x ) , start=list(a=0.005,b=0) )

plot( y/x ~ x )
k <- coef(mricker)
curve( k[1]*x*exp(k[2]*x) , from=min(d$Initial) , to=max(d$Initial) , col="red" , add=TRUE )