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