R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), 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. 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 a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("Bhat-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('Bhat') > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "dfp" > > ### * dfp > > flush(stderr()); flush(stdout()) > > ### Name: dfp > ### Title: Function minimization with box-constraints > ### Aliases: dfp > ### Keywords: optimize methods > > ### ** Examples > > # generate some Poisson counts on the fly > dose <- c(rep(0,50),rep(1,50),rep(5,50),rep(10,50)) > data <- cbind(dose,rpois(200,20*(1+dose*.5*(1-dose*0.05)))) > > # neg. log-likelihood of Poisson model with 'linear-quadratic' mean: > lkh <- function (x) { + ds <- data[, 1] + y <- data[, 2] + g <- x[1] * (1 + ds * x[2] * (1 - x[3] * ds)) + return(sum(g - y * log(g))) + } > > # for example define > x <- list(label=c("a","b","c"),est=c(10.,10.,.01),low=c(0,0,0),upp=c(100,20,.1)) > > # call: > results <- dfp(x,f=lkh) Wed Jul 13 09:17:00 2005 starting at 1 fmin: 24607.17 10.00 10.00 0.01 reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size covariance matrix is not positive-definite or nearly singular dfp search intermediate output: iter: 0 fmin: -6156.999 nfcn: 23 start main loop: dfp loop: 1 reducing step size reducing step size reducing step size dfp search intermediate output: iter: 0 fmin: -25800.79 nfcn: 139 restart main loop: dfp loop: 2 reducing step size reducing step size reducing step size reducing step size dfp search intermediate output: iter: 0 fmin: -25805.55 nfcn: 183 restart main loop: DFP search completed with status 0 Optimization Result: iter: 14 fmin: -25834.62 nfcn: 338 label estimate low high 1 a 20.316295 19.230972 21.446605 2 b 0.49252516 0.43152055 0.56190643 3 c 0.050242832 0.045046771 0.055433653 > > > > cleanEx(); ..nameEx <- "dqstep" > > ### * dqstep > > flush(stderr()); flush(stdout()) > > ### Name: dqstep > ### Title: step size generator > ### Aliases: dqstep > ### Keywords: optimize iteration > > ### ** Examples > > ## Rosenbrock Banana function > fr <- function(x) { + x1 <- x[1] + x2 <- x[2] + 100 * (x2 - x1 * x1)^2 + (1 - x1)^2 + } > ## define > x <- list(label=c("a","b"),est=c(1,1),low=c(0,0),upp=c(100,100)) > dqstep(x,fr,sens=1) [1] 0.0504410 0.1010096 > > > > cleanEx(); ..nameEx <- "logit.hessian" > > ### * logit.hessian > > flush(stderr()); flush(stdout()) > > ### Name: logit.hessian > ### Title: Hessian (curvature matrix) > ### Aliases: logit.hessian > ### Keywords: array iteration methods optimize > > ### ** Examples > > ## Rosenbrock Banana function > fr <- function(x) { + x1 <- x[1] + x2 <- x[2] + 100 * (x2 - x1 * x1)^2 + (1 - x1)^2 + } > ## define > x <- list(label=c("a","b"),est=c(1,1),low=c(-100,-100),upp=c(100,100)) > logit.hessian(x,f=fr,del=dqstep(x,f=fr,sens=0.01)) reducing step size reducing step size reducing step size $df [1] -0.4985783196 0.0000999999 $ddf [,1] [,2] [1,] 2004611 -999800 [2,] -999800 499900 $nfcn [1] 9 $eigen [1] 2503510.783 1000.686 > ## shows the differences in curvature at the minimum of the Banana > ## function along principal axis (in a logit-transformed coordinate system) > > > > cleanEx(); ..nameEx <- "mcmc" > > ### * mcmc > > flush(stderr()); flush(stdout()) > > ### Name: mymcmc > ### Title: Adaptive Multivariate MCMC sampler > ### Aliases: mymcmc > ### Keywords: iteration methods optimize > > ### ** Examples > > # generate some Poisson counts on the fly > dose <- c(rep(0,50),rep(1,50),rep(5,50),rep(10,50)) > data <- cbind(dose,rpois(200,20*(1+dose*.5*(1-dose*0.05)))) > > # neg. log-likelihood of Poisson model with 'linear-quadratic' mean: > nlogf <- function (x) { + ds <- data[, 1] + y <- data[, 2] + g <- x[1] * (1 + ds * x[2] * (1 - x[3] * ds)) + return(sum(g - y * log(g))) + } > > # start MCMC near mle > x <- list(label=c("a","b","c"), est=c(20, 0.5, 0.05), low=c(0,0,0), upp=c(100,10,.1)) > # samples from posterior density (~exp(-nlogf))) with non-informative > # (random uniform) priors for "a", "b" and "c". > out <- mymcmc(x, nlogf, m1=2000, m2=2000, m3=10000, scl1=0.5, scl2=2, skip=10, plot=TRUE) Wed Jul 13 09:17:00 2005 trying a proposal COVM using the Hessian first pilot chain: using MLEs and log-transformed covariance n: 1000 acceptance: 0.461 -log lkh: -25832.822 20.1284 0.48921 0.0511068 n: 2000 acceptance: 0.452 -log lkh: -25830.097 19.6566 0.553955 0.0553785 second pilot run: eigen values: 0.1086997 0.0001322501 4.49005e-07 n: 3000 acceptance: 0.351 -log lkh: -25830.812 21.3629 0.447403 0.0464231 n: 4000 acceptance: 0.397 -log lkh: -25833.645 20.624 0.474851 0.0490967 n: 5000 acceptance: 0.414 -log lkh: -25833.518 20.5874 0.466866 0.0467717 n: 6000 acceptance: 0.363 -log lkh: -25834.061 20.2172 0.495489 0.0489941 n: 7000 acceptance: 0.382 -log lkh: -25827.983 20.2238 0.500651 0.0515664 n: 8000 acceptance: 0.4 -log lkh: -25831.509 20.773 0.442728 0.0456827 n: 9000 acceptance: 0.369 -log lkh: -25831.62 20.2164 0.462951 0.0451255 n: 10000 acceptance: 0.378 -log lkh: -25823.711 20.175 0.487742 0.0498267 n: 11000 acceptance: 0.375 -log lkh: -25829.975 19.7602 0.537283 0.0536095 n: 12000 acceptance: 0.385 -log lkh: -25829.933 20.5069 0.519423 0.0526417 > # start MCMC from some other point: e.g. try x$est <- c(16,.2,.02) > > > > cleanEx(); ..nameEx <- "newton" > > ### * newton > > flush(stderr()); flush(stdout()) > > ### Name: newton > ### Title: Function minimization with box-constraints > ### Aliases: newton > ### Keywords: optimize methods > > ### ** Examples > > # generate some Poisson counts on the fly > dose <- c(rep(0,100),rep(1,100),rep(5,100),rep(10,100)) > data <- cbind(dose,rpois(400,20*(1+dose*.5*(1-dose*0.05)))) > > # neg. log-likelihood of Poisson model with 'linear-quadratic' mean: > lkh <- function (x) { + ds <- data[, 1] + y <- data[, 2] + g <- x[1] * (1 + ds * x[2] * (1 - x[3] * ds)) + return(sum(g - y * log(g))) + } > > # for example define > x <- list(label=c("a","b","c"),est=c(10.,10.,.01),low=c(0,0,0),upp=c(100,20,.1)) > > # calls: > r <- dfp(x,f=lkh) Wed Jul 13 09:17:04 2005 starting at 1 fmin: 49929.56 10.00 10.00 0.01 reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size covariance matrix is not positive-definite or nearly singular covariance matrix is not positive-definite or nearly singular parameter c near upper boundary covariance matrix is not positive-definite or nearly singular parameter c near upper boundary parameter c near lower boundary parameter c near upper boundary parameter c near upper boundary covariance matrix is not positive-definite or nearly singular parameter c near upper boundary parameter c near lower boundary parameter c near upper boundary dfp search intermediate output: iter: 0 fmin: -33836.74 nfcn: 58 start main loop: DFP search completed with status 0 reducing step size Optimization Result: iter: 13 fmin: -51217.99 nfcn: 205 label estimate low high 1 a 20.127046 19.363575 20.912812 2 b 0.4943406 0.45047228 0.54236244 3 c 0.050167078 0.04648119 0.053851152 > x$est <- r$est > results <- newton(x,lkh) Wed Jul 13 09:17:04 2005 1 fmin: -51217.99 20.12704583 0.49434060 0.05016708 Bhat run: Wed Jul 13 09:17:04 2005 status: converged iter: 1 fmin: -51217.99 nfcn: 21 label estimate log deriv lower 95% upper 95% 1 a 20.127979 -0.006492 19.36447 20.913785 2 b 0.494256 -0.051552 0.450395 0.542269 3 c 0.05016 -0.022265 0.046475 0.053845 > > > > cleanEx(); ..nameEx <- "plkhci" > > ### * plkhci > > flush(stderr()); flush(stdout()) > > ### Name: plkhci > ### Title: Profile-likelihood based confidence intervals > ### Aliases: plkhci > ### Keywords: distribution multivariate > > ### ** Examples > > # generate some Poisson counts on the fly > dose <- c(rep(0,50),rep(1,50),rep(5,50),rep(10,50)) > data <- cbind(dose,rpois(200,20*(1+dose*.5*(1-dose*0.05)))) > > # neg. log-likelihood of Poisson model with 'linear-quadratic' mean: > nlogf <- function (x) { + ds <- data[, 1] + y <- data[, 2] + g <- x[1] * (1 + ds * x[2] * (1 - x[3] * ds)) + return(sum(g - y * log(g))) + } > > # for example define > x <- list(label=c("a","b","c"),est=c(10.,10.,.01),low=c(0,0,0),upp=c(100,20,.1)) > > # get MLEs using dfp: > r <- dfp(x,f=nlogf) Wed Jul 13 09:17:04 2005 starting at 1 fmin: 24607.17 10.00 10.00 0.01 reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size reducing step size covariance matrix is not positive-definite or nearly singular dfp search intermediate output: iter: 0 fmin: -6156.999 nfcn: 23 start main loop: dfp loop: 1 reducing step size reducing step size reducing step size dfp search intermediate output: iter: 0 fmin: -25800.79 nfcn: 139 restart main loop: dfp loop: 2 reducing step size reducing step size reducing step size reducing step size dfp search intermediate output: iter: 0 fmin: -25805.55 nfcn: 183 restart main loop: DFP search completed with status 0 Optimization Result: iter: 14 fmin: -25834.62 nfcn: 338 label estimate low high 1 a 20.316295 19.230972 21.446605 2 b 0.49252516 0.43152055 0.56190643 3 c 0.050242832 0.045046771 0.055433653 > x$est <- r$est > plkhci(x,nlogf,"a") neg. log. likelihood: -25834.62 will attempt to compute both bounds (+/- direction) trying lower bound ------------------------ starting at: 0.965577 initial guess: 20.8759 0.4669094 0.04943779 begin Newton-Raphson search for profile lkh conf. bounds: eps value for stop criterium: 0.001 nmax : 10 CONVERGENCE: 5 iterations chisquare value is: 3.841709 confidence bound of a is 21.44093 log derivatives: -0.01425335 -0.008702677 label estimate log deriv log curv 1 a 21.4409 -56.0714 5525.6 2 b 0.441749 -0.0142533 2823.65 3 c 0.0484816 -0.00870268 428.87 trying upper bound ------------------------ starting at: 0.9558562 initial guess: 19.76794 0.5195089 0.05104774 begin Newton-Raphson search for profile lkh conf. bounds: eps value for stop criterium: 0.001 nmax : 10 CONVERGENCE: 5 iterations chisquare value is: 3.841326 confidence bound of a is 19.22551 log derivatives: -0.01807002 -0.004989158 label estimate log deriv log curv 1 a 19.2255 55.7841 5772.52 2 b 0.547045 -0.01807 3132.6 3 c 0.0517257 -0.00498916 530.384 [1] 19.22551 21.44093 > plkhci(x,nlogf,"b") neg. log. likelihood: -25834.62 will attempt to compute both bounds (+/- direction) trying lower bound ------------------------ starting at: 0.9847576 initial guess: 19.87196 0.5260986 0.0521933 begin Newton-Raphson search for profile lkh conf. bounds: eps value for stop criterium: 0.001 nmax : 10 CONVERGENCE: 5 iterations chisquare value is: 3.841237 confidence bound of b is 0.5602529 log derivatives: 0.009452938 -0.04430035 label estimate log deriv log curv 1 a 19.4186 0.00945294 5779.12 2 b 0.560253 -29.8501 3150.3 3 c 0.0537823 -0.0443003 565.484 trying upper bound ------------------------ starting at: 0.9385299 initial guess: 20.76799 0.4610436 0.04829162 begin Newton-Raphson search for profile lkh conf. bounds: eps value for stop criterium: 0.001 nmax : 10 CONVERGENCE: 5 iterations chisquare value is: 3.841834 confidence bound of b is 0.4300321 log derivatives: 0.07238531 -0.0734946 label estimate log deriv log curv 1 a 21.2081 0.0723853 5525.9 2 b 0.430032 26.5941 2810.96 3 c 0.0458916 -0.0734946 390.857 [1] 0.4300321 0.5602529 > plkhci(x,nlogf,"c") neg. log. likelihood: -25834.62 will attempt to compute both bounds (+/- direction) trying lower bound ------------------------ starting at: 1.040933 initial guess: 20.14566 0.5174733 0.05284519 begin Newton-Raphson search for profile lkh conf. bounds: eps value for stop criterium: 0.001 nmax : 10 CONVERGENCE: 5 iterations chisquare value is: 3.841268 confidence bound of c is 0.05506524 log derivatives: 0.1309275 0.08963562 label estimate log deriv log curv 1 a 20.012 0.130927 5694.69 2 b 0.540247 0.0896356 3021.78 3 c 0.0550652 -21.4587 567.975 trying upper bound ------------------------ starting at: 0.8860591 initial guess: 20.48801 0.4687509 0.04763916 begin Newton-Raphson search for profile lkh conf. bounds: eps value for stop criterium: 0.001 nmax : 10 CONVERGENCE: 5 iterations chisquare value is: 3.841815 confidence bound of c is 0.04456647 log derivatives: 0.229306 0.1801300 label estimate log deriv log curv 1 a 20.6975 0.229306 5597.7 2 b 0.443058 0.18013 2922.08 3 c 0.0445665 15.0105 379.714 [1] 0.04456647 0.05506524 > # e.g. 90% confidence bounds for "c" > plkhci(x,nlogf,"c",prob=0.9) neg. log. likelihood: -25834.62 will attempt to compute both bounds (+/- direction) trying lower bound ------------------------ starting at: 0.7236085 initial guess: 20.17302 0.5133809 0.05242747 begin Newton-Raphson search for profile lkh conf. bounds: eps value for stop criterium: 0.001 nmax : 10 CONVERGENCE: 5 iterations chisquare value is: 2.705358 confidence bound of c is 0.05433885 log derivatives: 0.0970066 0.06557293 label estimate log deriv log curv 1 a 20.0574 0.0970066 5688.25 2 b 0.532647 0.0655729 3015.71 3 c 0.0543389 -17.6665 554.827 trying upper bound ------------------------ starting at: 0.632051 initial guess: 20.46033 0.4724961 0.04805727 begin Newton-Raphson search for profile lkh conf. bounds: eps value for stop criterium: 0.001 nmax : 10 CONVERGENCE: 5 iterations chisquare value is: 2.705799 confidence bound of c is 0.04554668 log derivatives: 0.1513282 0.1161526 label estimate log deriv log curv 1 a 20.6325 0.151328 5606.89 2 b 0.451033 0.116153 2932.11 3 c 0.0455467 12.9979 396.615 [1] 0.04554668 0.05433885 > > > > > ### *