1 2### Test battery for various optimization parameters for different optimizers. 3### 4### ... 5### 6library(maxLik) 7library(tinytest) 8 9tol <- .Machine$double.eps^(0.25) 10set.seed( 123 ) 11# generate a variable from normally distributed random numbers 12N <- 50 13x <- rnorm(N, 1, 2 ) 14 15## log likelihood function 16llf <- function( param ) { 17 mu <- param[ 1 ] 18 sigma <- param[ 2 ] 19 if(!(sigma > 0)) 20 return(NA) 21 # to avoid warnings in the output 22 N <- length( x ) 23 llValue <- -0.5 * N * log( 2 * pi ) - N * log( sigma ) - 24 0.5 * sum( ( x - mu )^2 / sigma^2 ) 25 return( llValue ) 26} 27 28# start values 29startVal <- c( mu = 0, sigma = 1 ) 30 31# 32expect_silent(ml <- maxLik( llf, start = startVal )) 33expect_equivalent(coef(ml), c(1.069, 1.833), tolerance=tol) 34## tol 35expect_silent(mlTol <- maxLik( llf, start = startVal, tol=1)) 36expect_equal(returnCode(mlTol), 2) 37 # tolerance limit 38expect_silent(mlTolC <- maxLik(llf, start=startVal, control=list(tol=1))) 39expect_equal(coef(mlTol), coef(mlTolC)) 40expect_equal(hessian(mlTol), hessian(mlTolC)) 41expect_equal(returnCode(mlTol), returnCode(mlTolC)) 42expect_silent(ml <- maxLik( llf, start = startVal, tol=-1)) 43 # negative tol switches tol off 44expect_silent(ml <- maxLik( llf, start = startVal, control=list(tol=-1))) 45expect_false(returnCode(ml) == 2) 46 # should not be w/in tolerance limit 47expect_error(ml <- maxLik( llf, start = startVal, tol=c(1,2)), 48 pattern="'tol' must be of length 1, not 2") 49expect_error(ml <- maxLik( llf, start = startVal, control=list(tol=c(1,2))), 50 pattern="'tol' must be of length 1, not 2") 51expect_error(ml <- maxLik( llf, start = startVal, tol=TRUE), 52 pattern="object of class \"logical\" is not valid for slot 'tol'") 53expect_error(ml <- maxLik( llf, start = startVal, control=list(tol=TRUE)), 54 pattern="object of class \"logical\" is not valid for slot 'tol'") 55 56## ----- reltol: play w/reltol, leave other tolerances at default value ----- 57expect_silent(mlRelTol <- maxLik( llf, start = startVal, reltol=1)) 58expect_equal(returnCode(mlRelTol), 8) 59mlRelTolC <- maxLik(llf, start=startVal, control=list(reltol=1)) 60expect_equal(coef(mlRelTol), coef(mlRelTolC)) 61expect_silent(ml0 <- maxLik( llf, start = startVal, reltol=0)) 62expect_true(nIter(ml0) > nIter(mlRelTol)) 63 # switching off reltol makes more iterations 64expect_silent(ml1 <- maxLik( llf, start = startVal, reltol=-1)) 65expect_equal(nIter(ml0), nIter(ml1)) 66expect_error(ml <- maxLik( llf, start = startVal, reltol=c(1,2)), 67 pattern="invalid class \"MaxControl\" object: 'reltol' must be of length 1, not 2") 68expect_error(ml <- maxLik( llf, start = startVal, control=list(reltol=c(1,2))), 69 pattern="invalid class \"MaxControl\" object: 'reltol' must be of length 1, not 2") 70expect_error(ml <- maxLik( llf, start = startVal, reltol=TRUE), 71 pattern="assignment of an object of class \"logical\" is not valid for slot 'reltol'") 72expect_error(ml <- maxLik( llf, start = startVal, control=list(reltol=TRUE)), 73 pattern="assignment of an object of class \"logical\" is not valid for slot 'reltol'") 74## gradtol 75expect_silent(mlGradtol <- maxLik( llf, start = startVal, gradtol=0.1)) 76expect_equal(returnCode(mlGradtol), 1) 77mlGradtolC <- maxLik(llf, start=startVal, control=list(gradtol=0.1)) 78expect_equal(coef(mlGradtol), coef(mlGradtolC)) 79expect_silent(ml <- maxLik( llf, start = startVal, gradtol=-1)) 80expect_true(nIter(ml) > nIter(mlGradtol)) 81 # switching off gradtol makes more iterations 82expect_error(ml <- maxLik( llf, start = startVal, gradtol=c(1,2)), 83 pattern="object: 'gradtol' must be of length 1, not 2") 84expect_error(ml <- maxLik( llf, start = startVal, control=list(gradtol=c(1,2))), 85 pattern="object: 'gradtol' must be of length 1, not 2") 86expect_error(ml <- maxLik( llf, start = startVal, gradtol=TRUE), 87 pattern="assignment of an object of class \"logical\" is not valid for slot 'gradtol' ") 88expect_error(ml <- maxLik( llf, start = startVal, control=list(gradtol=TRUE)), 89 pattern="assignment of an object of class \"logical\" is not valid for slot 'gradtol' ") 90## examples with steptol, lambdatol 91## qac 92expect_silent(mlMarq <- maxLik( llf, start = startVal, qac="marquardt")) 93expect_equal(maximType(mlMarq), 94 "Newton-Raphson maximisation with Marquardt (1963) Hessian correction") 95expect_silent(mlMarqC <- maxLik(llf, start=startVal, control=list(qac="marquardt"))) 96expect_equal(coef(mlMarq), coef(mlMarqC)) 97expect_error(ml <- maxLik( llf, start = startVal, qac=-1), 98 pattern = "assignment of an object of class \"numeric\" is not valid for slot 'qac'") 99 # qac should be "stephalving" or "marquardt" 100expect_error(ml <- maxLik( llf, start = startVal, qac=c("a", "b")), 101 pattern = "invalid class \"MaxControl\" object: 'qac' must be of length 1, not 2") 102expect_error(ml <- maxLik( llf, start = startVal, qac=TRUE), 103 pattern = "assignment of an object of class \"logical\" is not valid for slot 'qac'") 104mlMarqCl <- maxLik(llf, start = startVal, 105 control=list(qac="marquardt", lambda0=1000, lambdaStep=4)) 106expect_equal(coef(mlMarqCl), coef(mlMarq)) 107## NM: alpha, beta, gamma 108expect_silent(mlNMAlpha <- maxLik(llf, start=startVal, method="nm", beta=0.8)) 109expect_silent(mlNMAlphaC <- maxLik(llf, start=startVal, method="nm", control=list(beta=0.8))) 110expect_equal(coef(mlNMAlpha), coef(mlNMAlphaC)) 111 112## likelihood function with additional parameter 113llf1 <- function( param, sigma ) { 114 mu <- param 115 N <- length( x ) 116 ll <- -0.5*N*log( 2 * pi ) - N*log( sigma ) - 117 0.5*sum( ( x - mu )^2/sigma^2 ) 118 ll 119} 120 121## log-lik mixture 122logLikMix <- function(param) { 123 rho <- param[1] 124 if(rho < 0 || rho > 1) 125 return(NA) 126 mu1 <- param[2] 127 mu2 <- param[3] 128 ll <- log(rho*dnorm(x - mu1) + (1 - rho)*dnorm(x - mu2)) 129 ll 130} 131 132## loglik mixture with additional parameter 133logLikMixA <- function(param, rho) { 134 mu1 <- param[1] 135 mu2 <- param[2] 136 ll <- log(rho*dnorm(x - mu1) + (1 - rho)*dnorm(x - mu2)) 137 ll 138} 139 140## Test the following with all the main optimizers: 141pl2Patterns <- c(NR = "----- Initial parameters: -----\n.*-----Iteration 1 -----", 142 BFGS = "initial value.*final value", 143 BFGSR = "-------- Initial parameters: -------\n.*Iteration 1") 144for(method in c("NR", "BFGS", "BFGSR")) { 145 ## create data in loop, we need to mess with 'x' for constraints 146 N <- 100 147 x <- rnorm(N, 1, 2 ) 148 startVal <- c(1,2) 149 ## two parameters at the same time 150 ## iterlim, printLevel 151 expect_stdout(ml2 <- maxLik(llf, start=startVal, method=method, 152 iterlim=1, printLevel=2), 153 pattern = pl2Patterns[method]) 154 expect_stdout(ml2C <- maxLik(llf, start=startVal, method=method, 155 control=list(iterlim=1, printLevel=2)), 156 pattern = pl2Patterns[method]) 157 expect_equal(coef(ml2), coef(ml2C)) 158 ## what about additional parameters for the loglik function? 159 expect_silent(mlsM <- maxLik(llf1, start=0, method=method, tol=1, sigma=1)) 160 expect_silent(mlsCM <- maxLik(llf1, start=0, method=method, control=list(tol=1), sigma=1)) 161 expect_equal(coef(mlsM), coef(mlsCM)) 162 ## And what about unused parameters? 163 expect_error(maxLik(llf1, start=0, method=method, control=list(tol=1), 164 sigma=1, unusedPar=2), 165 pattern = "unused argument") 166 N <- 100 167 ## Does this work with constraints? 168 x <- c(rnorm(N, mean=-1), rnorm(N, mean=1)) 169 ## First test inequality constraints 170 ## Inequality constraints: x + y + z < 0.5 171 A <- matrix(c(-1, 0, 0, 172 0, -1, 0, 173 0, 0, 1), 3, 3, byrow=TRUE) 174 B <- rep(0.5, 3) 175 start <- c(0.4, 0, 0.9) 176 ## analytic gradient 177 if(!(method %in% c("NR", "BFGSR"))) { 178 expect_silent(mix <- maxLik(logLikMix, 179 start=start, method=method, 180 constraints=list(ineqA=A, ineqB=B))) 181 expect_silent(mixGT <- try(maxLik(logLikMix, 182 start=start, method=method, 183 constraints=list(ineqA=A, ineqB=B), 184 tol=1))) 185 expect_silent( 186 mixGTC <- try(maxLik(logLikMix, 187 start=start, method=method, 188 constraints=list(ineqA=A, ineqB=B), 189 control=list(tol=1))) 190 ) 191 ## 2d inequality constraints: x + y < 0.5 192 A2 <- matrix(c(-1, -1), 1, 2, byrow=TRUE) 193 B2 <- 0.5 194 start2 <- c(-0.5, 0.5) 195 expect_silent( 196 mixA <- maxLik(logLikMixA, 197 start=start2, method=method, 198 constraints=list(ineqA=A2, ineqB=B2), 199 tol=1, 200 rho=0.5) 201 ) 202 expect_silent( 203 mixAC <- maxLik(logLikMixA, 204 start=start2, method=method, 205 constraints=list(ineqA=A2, ineqB=B2), 206 control=list(tol=1), 207 rho=0.5) 208 ) 209 expect_equal(coef(mixA), coef(mixAC)) 210 expect_equal(hessian(mixA), hessian(mixAC)) 211 } 212} 213 214### Test adding both default and user-specified parameters through control list 215estimate <- function(control=NULL, ...) { 216 maxLik(llf, start=c(1,1), 217 control=c(list(iterlim=100), control), 218 ...) 219} 220expect_silent(m <- estimate(control=list(iterlim=1), fixed=2)) 221expect_stdout(show(maxControl(m)), 222 pattern = "iterlim = 1") 223 # iterlim should be 1 224expect_equal(coef(m)[2], 1) 225 # sigma should be 1.000 226## Does print.level overwrite 'printLevel'? 227expect_silent(m <- estimate(control=list(printLevel=2, print.level=1))) 228expect_stdout(show(maxControl(m)), 229 pattern = "printLevel = 1") 230 231## Does open parameters override everything? 232expect_silent(m <- estimate(control=list(printLevel=2, print.level=1), print.level=0)) 233expect_stdout(show(maxControl(m)), 234 pattern = "printLevel = 0") 235 236### does both printLevel, print.level work for condiNumber? 237expect_silent(condiNumber(hessian(m), print.level=0)) 238expect_silent(condiNumber(hessian(m), printLevel=0)) 239expect_silent(condiNumber(hessian(m), printLevel=0, print.level=1)) 240