1### SUMT (Sequential Unconstrained Maximization Technique) 2### borrowed from package 'clue' 3### 4### Adapted for linear constraints 5sumt <- function(fn, grad=NULL, hess=NULL, 6 start, 7 maxRoutine, constraints, 8 SUMTTol = sqrt(.Machine$double.eps), 9 # difference between estimates for successive outer iterations 10 SUMTPenaltyTol = sqrt(.Machine$double.eps), 11 # maximum allowed penalty 12 SUMTQ = 10, 13 SUMTRho0 = NULL, 14 printLevel=print.level, print.level=0, 15 SUMTMaxIter=100, 16 ...) { 17 ## constraints list w/components eqA and eqB. Maximization will 18 ## be performed wrt to the constraint 19 ## A %*% theta + B = 0 20 ## The user must ensure the matrices are in correct 21 ## form 22 ## maxSUMTiter how many SUMT iterations to perform max 23 ## 24 penalty <- function(theta) { 25 p <- A %*% theta + B 26 sum(p*p) 27 } 28 ## Penalty gradient and Hessian are used only if corresponding function 29 ## for the likelihood function is provided 30 gPenalty <- function(theta) { 31 2*(t(theta) %*% t(A) %*% A - t(B) %*% A) 32 } 33 hessPenalty <- function(theta) { 34 2*t(A) %*% A 35 } 36 37 ## strip possible arguments of maxRoutine and call the function thereafter 38 callWithoutMaxArgs <- function(theta, fName, ...) { 39 return( callWithoutArgs( theta, fName = fName, 40 args = names(formals(maxRoutine)), ... ) ) 41 } 42 43 SUMTMessage <- function(code) { 44 message <- switch(code, 45 "1" = "penalty close to zero", 46 "2" = "successive function values within tolerance limit", 47 "4" = "Outer iteration limit exceeded (increase SUMTMaxIter ?).", 48 paste("Code", code)) 49 return(message) 50 } 51 ## the penalized objective function 52 Phi <- function(theta, ...) { 53 llVal <- callWithoutMaxArgs( theta, "logLikFunc", fnOrig = fn, 54 gradOrig = grad, hessOrig = hess, sumObs = FALSE, ... ) 55 llVal <- llVal - rho * penalty( theta ) / length( llVal ) 56 g <- attributes( llVal )$gradient 57 if( !is.null( g ) ) { 58 if( is.matrix( g ) ) { 59 g <- g - matrix( 60 rep( rho * gPenalty( theta ) / nrow( g ), each = nrow( g ) ), 61 nrow = nrow( g ), ncol = ncol( g ) ) 62 } else { 63 g <- g - rho * gPenalty( theta ) 64 } 65 attributes( llVal )$gradient <- g 66 } 67 h <- attributes( llVal )$hessian 68 if( !is.null( h ) ) { 69 attributes( llVal )$hessian <- h - rho * hessPenalty( theta ) 70 } 71 return( llVal ) 72 } 73 ## gradient of the penalized objective function 74 if(!is.null(grad)) { 75 gradPhi<- function(theta, ...) { 76 g <- grad(theta, ...) 77 if(is.matrix(g)) { 78 g <- g - matrix( 79 rep( rho * gPenalty( theta ) / nrow( g ), each = nrow( g ) ), 80 nrow = nrow( g ), ncol = ncol( g ) ) 81 } else { 82 g <- g - rho * gPenalty( theta ) 83 } 84 return( g ) 85 } 86 } else { 87 gradPhi <- NULL 88 } 89 ## Hessian of the penalized objective function 90 if(!is.null(hess)) { 91 hessPhi <- function(theta, ...) { 92 return( hess(theta, ...) - rho*hessPenalty(theta) ) 93 } 94 } else { 95 hessPhi <- NULL 96 } 97 ## -------- SUMT Main code --------- 98 ## Note also that currently we do not check whether optimization was 99 ## "successful" ... 100 A <- constraints$eqA 101 B <- as.matrix(constraints$eqB) 102 ## Check if the matrices conform 103 if(ncol(A) != length(start)) { 104 stop("Equality constraint matrix A must have the same number\n", 105 "of columns as the parameter length ", 106 "(currently ", ncol(A), " and ", length(start), ")") 107 } 108 if(nrow(A) != nrow(B)) { 109 stop("Equality constraint matrix A must have the same number\n", 110 "of rows as the matrix B ", 111 "(currently ", nrow(A), " and ", nrow(B), ")") 112 } 113 ## Find a suitable inital value for rho if not specified 114 if(is.null(SUMTRho0)) { 115 rho <- 0 116 result <- maxRoutine(fn=Phi, grad=gradPhi, hess=hessPhi, 117 start=start, 118 printLevel=max(printLevel - 1, 0), 119 ...) 120 theta <- coef(result) 121 # Note: this may be a bad idea, 122 # if unconstrained function is unbounded 123 # from above. In that case rather specify SUMTRho0. 124 if(printLevel > 0) { 125 cat("SUMT initial: rho = ", rho, 126 ", function = ", 127 callWithoutMaxArgs( theta, "logLikFunc", 128 fnOrig = fn, gradOrig = grad, 129 hessOrig = hess, ... ), 130 ", penalty = ", penalty(theta), "\n") 131 cat("Estimate:") 132 print(theta) 133 } 134 ## Better upper/lower bounds for rho? 135 rho <- max( callWithoutMaxArgs( theta, "logLikFunc", fnOrig = fn, 136 gradOrig = grad, hessOrig = hess, ... ), 1e-3) / 137 max(penalty(start), 1e-3) 138 } 139 ## if rho specified, simply pick that and use previous initial values 140 else { 141 rho <- SUMTRho0 142 theta <- start 143 } 144 ## </TODO> 145 iter <- 1L 146 repeat { 147 thetaOld <- theta 148 result <- maxRoutine(fn=Phi, grad=gradPhi, hess=hessPhi, 149 start=thetaOld, 150 printLevel=max(printLevel - 1, 0), 151 ...) 152 theta <- coef(result) 153 if(printLevel > 0) { 154 cat("SUMT iteration ", iter, 155 ": rho = ", rho, ", function = ", callWithoutMaxArgs( theta, 156 "logLikFunc", fnOrig = fn, gradOrig = grad, hessOrig = hess, ... ), 157 ", penalty = ", penalty(theta), "\n", sep="") 158 cat("Estimate:") 159 print(theta) 160 } 161 if(max(abs(thetaOld - theta)) < SUMTTol) { 162 SUMTCode <- 2 163 break 164 } 165 if(penalty(theta) < SUMTPenaltyTol) { 166 SUMTCode <- 1 167 break 168 } 169 if(iter >= SUMTMaxIter) { 170 SUMTCode <- 4 171 break 172 } 173 iter <- iter + 1L 174 rho <- SUMTQ * rho 175 } 176 ## Now we replace the resulting gradient and Hessian with those, 177 ## calculated on the original function 178 llVal <- callWithoutMaxArgs( theta, "logLikFunc", fnOrig = fn, 179 gradOrig = grad, hessOrig = hess, sumObs = FALSE, ... ) 180 gradient <- attr( llVal, "gradient" ) 181 if( is.null( gradient ) ) { 182 gradient <- callWithoutMaxArgs( theta, "logLikGrad", fnOrig = fn, 183 gradOrig = grad, hessOrig = hess, sumObs = FALSE, ... ) 184 } 185 if( !is.null( dim( gradient ) ) ) { 186 if( nrow( gradient ) > 1 ) { 187 gradientObs <- gradient 188 } 189 gradient <- colSums( gradient ) 190 } else if( length( start ) == 1 && length( gradient ) > 1 ) { 191 gradientObs <- matrix( gradient, ncol = 1 ) 192 gradient <- sum( gradient ) 193 } 194 result$gradient <- gradient 195 names( result$gradient ) <- names( result$estimate ) 196 197 result$hessian <- callWithoutMaxArgs( theta, "logLikHess", fnOrig = fn, 198 gradOrig = grad, hessOrig = hess, ... ) 199 result$constraints <- list(type="SUMT", 200 barrier.value=penalty(theta), 201 code=SUMTCode, 202 message=SUMTMessage(SUMTCode), 203 outer.iterations=iter 204 ) 205 if( exists( "gradientObs" ) ) { 206 result$gradientObs <- gradientObs 207 colnames( result$gradientObs ) <- names( result$estimate ) 208 } 209 210 if( result$constraints$barrier.value > 0.001 ) { 211 warning( "problem in imposing equality constraints: the constraints", 212 " are not satisfied (barrier value = ", 213 result$constraints$barrier.value, "). Try setting 'SUMTTol' to 0" ) 214 } 215 return(result) 216} 217