R version 2.4.1 (2006-12-18)
Copyright (C) 2006 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> invisible(options(echo = TRUE))
> 
> library(bbmle)
> 
> funcresp <-
+ structure(list(Initial = as.integer(c(5, 5, 10, 10, 15, 15, 20, 
+ 20, 30, 30, 50, 50, 75, 75, 100, 100)), Killed = as.integer(c(1, 
+ 2, 5, 6, 10, 9, 7, 10, 11, 15, 5, 21, 32, 18, 25, 35))), .Names = c("Initial", 
+ "Killed"), class = "data.frame", row.names = c("1", "2", "3", 
+ "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", 
+ "16"))
> 
> attach(funcresp)
> 
> binomNLL2 = function(p) {
+   a = p[1]
+   h = p[2]
+   ##  cat(a,h,"\n")
+   p = a/(1+a*h*N)
+   -sum(dbinom(k,prob=p,size=N,log=TRUE))
+ }
> 
> N=0; k=0
> parnames(binomNLL2) = c("a","h")
> m2a = mle2(binomNLL2,start=c(a=0.5,h=0.0125),
+   data=list(N=Initial,k=Killed))
> p1a = profile(m2a); p1a
There were 50 or more warnings (use warnings() to see the first 50)
Likelihood profile:

$a
            z   par.vals.a   par.vals.h
1  -4.1986714  0.295147615 -0.002889674
2  -3.1509952  0.341378731  0.002607598
3  -2.2323085  0.387609846  0.007024647
4  -1.4128689  0.433840962  0.010701993
5  -0.6718644  0.480072077  0.013862508
6   0.0000000  0.526303193  0.016643616
7   0.6315243  0.572534308  0.019110825
8   1.2143907  0.618765424  0.021394466
9   1.7613448  0.664996539  0.023488432
10  2.2783229  0.711227654  0.025466822
11  2.7703242  0.757458770  0.027346030
12  3.2417836  0.803689885  0.029159198
13  3.6967744  0.849921001  0.030931996

$h
            z    par.vals.a    par.vals.h
1  -3.7508090  0.3274414241 -0.0023561362
2  -3.1636217  0.3548018002  0.0008104892
3  -2.5551771  0.3848461640  0.0039771145
4  -1.9288107  0.4174290358  0.0071437399
5  -1.2890406  0.4522236670  0.0103103652
6  -0.6413233  0.4887395904  0.0134769906
7   0.0000000  0.5263031926  0.0166436159
8   0.6539107  0.5645080828  0.0198102413
9   1.2903768  0.6025673699  0.0229768666
10  1.9132078  0.6400950914  0.0261434920
11  2.5191552  0.6767631680  0.0293101173
12  3.1060150  0.7123612724  0.0324767427
13  3.6724861  0.7467818281  0.0356433680

> c2a = confint(p1a); c2a
        2.5 %     97.5 %
a 0.402494943 0.68249515
h 0.006987242 0.02638530
> 
> binomNLL2b = function(p,N,k) {
+   a = p[1]
+   h = p[2]
+   ##  cat(a,h,"\n")
+   p = a/(1+a*h*N)
+   -sum(dbinom(k,prob=p,size=N,log=TRUE))
+ }
> parnames(binomNLL2b) = c("a","h")
> m2b = mle2(binomNLL2,start=c(a=0.5,h=0.0125),
+   data=list(N=Initial,k=Killed))
> c2b = confint(m2b)
There were 50 or more warnings (use warnings() to see the first 50)
> 
> N=Initial; k=Killed
> m2c = mle2(binomNLL2,start=c(a=0.5,h=0.0125))
> c2c = confint(m2c); c2c
There were 50 or more warnings (use warnings() to see the first 50)
        2.5 %     97.5 %
a 0.402494943 0.68249515
h 0.006987242 0.02638530
> 
> detach(funcresp)
> 
> 
> proc.time()
[1] 9.040 0.108 9.217 0.000 0.000
>