R version 3.6.0 (2019-04-26) -- "Planting of a Tree"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
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 an HTML browser interface to help.
Type 'q()' to quit R.
> ### BFGSR-related tests
>
> ## 1. Test maximization algorithm for convex regions
> ##
> ## Optimize quadratic form t(D) %*% W %*% D with p.d. weight matrix
> ## (ie unbounded problems).
> ## All solutions should go to large values with a message about successful convergence
> set.seed(0)
> options(digits=4)
> quadForm <- function(D) {
+ C <- seq(1, N)
+ return( - t(D - C) %*% W %*% ( D - C) )
+ }
> N <- 3
> # 3-dimensional case
> ## a) test quadratic function t(D) %*% D
> library(maxLik)
Loading required package: miscTools
Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
> W <- diag(N)
> D <- rep(1/N, N)
> res <- maxBFGSR(quadForm, start=D)
> all.equal(coef(res), 1:3, tolerance=1e-4)
[1] TRUE
> all.equal(gradient(res), rep(0,3), tolerance=1e-3)
[1] TRUE
> all.equal(nIter(res) < 100, TRUE)
[1] TRUE
> all.equal(returnCode(res) < 4, TRUE)
[1] TRUE
>
> ## Next, optimize hat function in non-concave region. Does not work well.
> hat <- function(param) {
+ ## Hat function. Hessian negative definite if sqrt(x^2 + y^2) < 0.5
+ x <- param[1]
+ y <- param[2]
+ exp(-(x-2)^2 - (y-2)^2)
+ }
>
> hatNC <- maxBFGSR(hat, start=c(1,1), tol=0, reltol=0)
> all.equal(coef(hatNC), rep(2,2), tolerance=1e-4)
[1] TRUE
> all.equal(gradient(hatNC), rep(0,2), tolerance=1e-3)
[1] TRUE
> all.equal(nIter(hatNC) < 100, TRUE)
[1] TRUE
> all.equal(returnCode(hatNC) < 4, TRUE)
[1] TRUE
>
> ## Test BFGSR with fixed parameters and equality constraints
> ## Optimize 3D hat with one parameter fixed (== 2D hat).
> ## Add an equality constraint on that
> hat3 <- function(param) {
+ ## Hat function. Hessian negative definite if sqrt((x-2)^2 + (y-2)^2) < 0.5
+ x <- param[1]
+ y <- param[2]
+ z <- param[3]
+ exp(-(x-2)^2-(y-2)^2-(z-2)^2)
+ }
> sv <- c(x=1,y=1,z=1)
> ## constraints: x + y + z = 8
> A <- matrix(c(1,1,1), 1, 3)
> B <- -8
> constraints <- list(eqA=A, eqB=B)
> hat3CF <- maxBFGSR(hat3, start=sv, constraints=constraints, fixed=3)
> all.equal(coef(hat3CF), c(x=3.5, y=3.5, z=1), tolerance=1e-4)
[1] TRUE
> all.equal(nIter(hat3CF) < 100, TRUE)
[1] TRUE
> all.equal(returnCode(hat3CF) < 4, TRUE)
[1] TRUE
> all.equal(sum(coef(hat3CF)), 8, tolerance=1e-4)
[1] TRUE
>
> proc.time()
user system elapsed
0.562 0.560 0.338
Generated by dwww version 1.15 on Thu May 23 22:39:23 CEST 2024.