55# multivariate extenstion. Works for lm and lmer
66# 2017-12-25: bug fix with multivariace bcnPower
77# 2019-03-07: bug fix in estimateTransform.bcnPowerlmer, thanks to wouter@zoology.ubc.ca
8- # 2019-11-14,15: change class(x) == "y" to inherits(x, "y") and likewise for !=
8+ # 2025-05-15: The argument gamma.min was not used, not it is
9+ # 2025-05-18: change class() == to inherits()
910
1011bcnPower <- function (U , lambda , jacobian.adjusted = FALSE , gamma ) {
1112 if (is.matrix(U )){
@@ -65,7 +66,8 @@ bcn.sv <- function(X, Y, weights, itmax=100, conv=.0001, verbose=FALSE,
6566 gamma.1d <- function (Y , weights , lambda , gamma , xqr ){
6667 fn1 <- function (gam ) bcnPowerllik(NULL , Y , weights , lambda = lambda ,
6768 gamma = gam , xqr = xqr )$ llik
68- f <- optimize(f = fn1 , interval = c(0.01 , max(Y )), maximum = TRUE )
69+ # f <- optimize(f=fn1, interval=c(0.01, max(Y)), maximum=TRUE) bug, ignores gamma.min
70+ f <- optimize(f = fn1 , interval = c(gamma.min , max(Y )), maximum = TRUE )
6971 list (lambda = lambda , gamma = f $ maximum , llik = f $ objective )
7072 }
7173 # get qr decomposition
@@ -84,14 +86,15 @@ bcn.sv <- function(X, Y, weights, itmax=100, conv=.0001, verbose=FALSE,
8486 last.value <- res
8587 res <- lambda.1d(Y , weights , res $ lambda , res $ gamma , xqr )
8688 res <- gamma.1d(Y , weights , res $ lambda , res $ gamma , xqr )
87- if (res $ gamma < 1.5 * gamma.min ){
89+ if (res $ gamma < 1.01 * gamma.min ){
90+ # if(res$gamma < 1.5 * gamma.min){
8891 gamma.ok <- FALSE
8992 res <- lambda.1d(Y , weights , res $ lambda , gamma.min , xqr )
9093 }
9194 crit <- (res $ llik - last.value $ llik )/ abs(res $ llik )
9295 if (verbose )
9396 print(data.frame (Iter = i , gamma = res $ gamma ,
94- lambda = res $ gamma , llik = res $ llik , crit = crit ))
97+ lambda = res $ lambda , llik = res $ llik , crit = crit ))
9598 }
9699 if (i == itmax & conv > crit )
97100 warning(paste(" No convergence in" itmax , " iterations, criterion =" crit , collapse = " "
@@ -163,10 +166,15 @@ estimateTransform.bcnPower <- function(X, Y, weights,
163166 w <- if (is.null(weights )) 1 else sqrt(weights )
164167 xqr <- qr(w * as.matrix(X ))
165168# if d = 1 call bcn.sv and return, else call bcn.sv to get starting values.
166- if (d == 1 ) bcn.sv(X , Y , weights , start = FALSE ) else {
169+ if (d == 1 ) bcn.sv(X , Y , weights = weights ,
170+ itmax = itmax , conv = conv , verbose = verbose ,
171+ start = FALSE , gamma.min = gamma.min ) else {
167172# The rest of this code is for the multivariate case
168173# get starting values for gamma
169- sv <- apply(Y , 2 , function (y ) unlist(bcn.sv(X , y , weights , start = TRUE )))
174+ sv <- apply(Y , 2 , function (y ) unlist(
175+ bcn.sv(X , y , weights = weights ,
176+ itmax = itmax , conv = conv , verbose = verbose ,
177+ start = TRUE , gamma.min = gamma.min )))
170178 res <- as.list(as.data.frame(t(sv ))) # output to a list
171179# gamma.estimated converted to numeric, so fixup
172180 res $ gamma.estimated <- ifelse(res $ gamma.estimated == 1 , TRUE , FALSE )
@@ -197,7 +205,7 @@ estimateTransform.bcnPower <- function(X, Y, weights,
197205 res $ gamma [! gamma.ok ] <- gamma.min
198206 if (all(gamma.ok )){
199207 hess <- try(optimHess(c(res $ lambda , res $ gamma ), fn4 ))
200- res $ invHess <- if (inherits(hess , " try-error" NA else solve(- hess )
208+ res $ invHess <- if (inherits(hess , " try-error" NA else solve(- hess )
201209 } else {
202210 fn4a <- function (lam ) fn4(c(lam , res $ gamma ))
203211 hess <- try(optimHess(res $ lambda , fn4a )) # hessian for lambda only
0 commit comments