Issue in Windows/M1 with type = "hh" after lamW update
petersonR opened this issue · 8 comments
petersonR commented
Hello, first off, thank you again for this package!
As you know, bestNormalize depends on lamW via LambertW. After the recent update to lamW, I began seeing errors in bestNormalize CRAN tests on Windows and M1 Macs for the Lambert function in bestNormalize. (see check results here). It appears that the Lambert transformation used by bestNormalize (see code here) is not inverting correctly, but this only occurs only on Windows and M1 macs.
In order to keep bestNormalize on CRAN I'll need to take these tests out, but I thought I would give you a heads up. I dug a bit into this and it looks like the type = "hh" doesn't seem to be working on Windows (see possible convergence issue below):
data(iris)
train <- iris$Petal.Width
summary(train)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.100 0.300 1.300 1.199 1.800 2.500
lambert_obj <- LambertW::Gaussianize(train, type = "hh")
summary(lambert_obj)
Y1.X
Min. :1.898
1st Qu.:1.898
Median :1.898
Mean :1.898
3rd Qu.:1.898
Max. :1.899
attributes(lambert_obj)
$dim
[1] 150 1
$dimnames
$dimnames[[1]]
NULL
$dimnames[[2]]
[1] "Y1.X"
$`Gaussianized:mu`
[1] 1.898343
$`Gaussianized:sigma`
[1] 6.4e-05
$`Gaussianized:delta_l`
[1] 7.545524e+31
$`Gaussianized:delta_r`
[1] 2.558503
$`Gaussianized:alpha_l`
[1] 1
$`Gaussianized:alpha_r`
[1] 1aadler commented
Hi. I'm the author and maintainer of lamW.
The only calculation change made for version 2.2.0 was to the principal/positive branch for -6.4e-3 <= x <= 6.4e-3, but abs(x) > 1e-16, where instead of using Fritsch's Pade estimate I implemented a degree-six polynomial minimax estimate. When abs(x) <= 1e-16 it still returns (x). So there is a slight change for those values.
However, when I recompiled on my (Windows) machine using Fritsch, or even using the original Halley step, I still get the same issue with LambertW::Gaussianize(train, type = "hh") (see below), which makes me think the issue isn't with lamW. What else could it be?
Using Fritsch (as per versions 2.0.0–2.1.2)
double lambertW0_CS(double x) {
if (x == R_PosInf) {
return(R_PosInf);
} else if (x < -M_1_E) {
return(R_NaN);
} else if (std::abs(x + M_1_E) <= EPS) {
return(-1.0);
} else if (std::abs(x) <= 1e-16) {
/* This close to 0 the W_0 branch is best estimated by its Taylor/Pade
expansion whose first term is the value x and remaining terms are below
machine double precision. See
https://math.stackexchange.com/questions/1700919
*/
return(x);
} else {
double w;
if (std::abs(x) <= 6.4e-3) {
// double p = std::sqrt(2.0 * (M_E * x + 1.0));
// double Numer = (0.2787037037037037 * p + 0.311111111111111) * p - 1.0;
// double Denom = (0.0768518518518518 * p + 0.688888888888889) * p + 1.0;
// return(HalleyIter(x, Numer / Denom));
// // Minimax Approximation calculated using R package minimaxApprox 0.1.0
// return((((((-1.0805085529250425e1 * x + 5.2100070265741278) * x -
// 2.6666665063383532) * x + 1.4999999657268301) * x -
// 1.0000000000016802) * x + 1.0000000000001752) * x +
// 2.6020852139652106e-18);
/* Use equation (5) in Fritsch */
w = ((1.33333333333333333 * x + 1.0) * x) /
((0.83333333333333333 * x + 2.33333333333333333) * x + 1.0);
} else if (x <= M_E) {
/* Use expansion in Corliss 4.22 to create (2, 2) Pade approximant.
* Equation with a few extra terms is:
* -1 + p - 1/3p^2 + 11/72p^3 - 43/540p^4 + 689453/8398080p^4 - O(p^5)
* This is just used to estimate a good starting point for the Fritsch
* iteration process itself.
*/
double p = std::sqrt(2.0 * (M_E * x + 1.0));
double Numer = (0.2787037037037037 * p + 0.311111111111111) * p - 1.0;
double Denom = (0.0768518518518518 * p + 0.688888888888889) * p + 1.0;
w = Numer / Denom;
} else {
/* Use first five terms of Corliss et al. 4.19 */
w = std::log(x);
double L_2 = std::log(w);
double L_3 = L_2 / w;
double L_3_sq = L_3 * L_3;
w += -L_2 + L_3 + 0.5 * L_3_sq - L_3 / w + L_3 / (w * w) - 1.5 * L_3_sq /
w + L_3_sq * L_3 / 3.0;
}
return(FritschIter(x, w));
}
}library(lamW)
train <- iris$Petal.Width
summary(train)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.100 0.300 1.300 1.199 1.800 2.500
lambert_obj <- LambertW::Gaussianize(train, type = "hh")
summary(lambert_obj)
Y1.X
Min. :1.898
1st Qu.:1.898
Median :1.898
Mean :1.898
3rd Qu.:1.898
Max. :1.899
> attributes(lambert_obj)
$dim
[1] 150 1
$dimnames
$dimnames[[1]]
NULL
$dimnames[[2]]
[1] "Y1.X"
$`Gaussianized:mu`
[1] 1.898343
$`Gaussianized:sigma`
[1] 6.4e-05
$`Gaussianized:delta_l`
[1] 2.71722e+15
$`Gaussianized:delta_r`
[1] 2.558455
$`Gaussianized:alpha_l`
[1] 1
$`Gaussianized:alpha_r`
[1] 1Using HalleyStep (as per versions prior to 2.0.0)
double lambertW0_CS(double x) {
if (x == R_PosInf) {
return(R_PosInf);
} else if (x < -M_1_E) {
return(R_NaN);
} else if (std::abs(x + M_1_E) <= EPS) {
return(-1.0);
} else if (std::abs(x) <= 1e-16) {
/* This close to 0 the W_0 branch is best estimated by its Taylor/Pade
expansion whose first term is the value x and remaining terms are below
machine double precision. See
https://math.stackexchange.com/questions/1700919
*/
return(x);
} else {
double w;
if (std::abs(x) <= 6.4e-3) {
double p = std::sqrt(2.0 * (M_E * x + 1.0));
double Numer = (0.2787037037037037 * p + 0.311111111111111) * p - 1.0;
double Denom = (0.0768518518518518 * p + 0.688888888888889) * p + 1.0;
return(HalleyIter(x, Numer / Denom));
// // Minimax Approximation calculated using R package minimaxApprox 0.1.0
// return((((((-1.0805085529250425e1 * x + 5.2100070265741278) * x -
// 2.6666665063383532) * x + 1.4999999657268301) * x -
// 1.0000000000016802) * x + 1.0000000000001752) * x +
// 2.6020852139652106e-18);
// /* Use equation (5) in Fritsch */
// w = ((1.33333333333333333 * x + 1.0) * x) /
// ((0.83333333333333333 * x + 2.33333333333333333) * x + 1.0);
} else if (x <= M_E) {
/* Use expansion in Corliss 4.22 to create (2, 2) Pade approximant.
* Equation with a few extra terms is:
* -1 + p - 1/3p^2 + 11/72p^3 - 43/540p^4 + 689453/8398080p^4 - O(p^5)
* This is just used to estimate a good starting point for the Fritsch
* iteration process itself.
*/
double p = std::sqrt(2.0 * (M_E * x + 1.0));
double Numer = (0.2787037037037037 * p + 0.311111111111111) * p - 1.0;
double Denom = (0.0768518518518518 * p + 0.688888888888889) * p + 1.0;
w = Numer / Denom;
} else {
/* Use first five terms of Corliss et al. 4.19 */
w = std::log(x);
double L_2 = std::log(w);
double L_3 = L_2 / w;
double L_3_sq = L_3 * L_3;
w += -L_2 + L_3 + 0.5 * L_3_sq - L_3 / w + L_3 / (w * w) - 1.5 * L_3_sq /
w + L_3_sq * L_3 / 3.0;
}
return(FritschIter(x, w));
}
}library(lamW)
train <- iris$Petal.Width
summary(train)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.100 0.300 1.300 1.199 1.800 2.500
lambert_obj <- LambertW::Gaussianize(train, type = "hh")
summary(lambert_obj)
Y1.X
Min. :1.898
1st Qu.:1.898
Median :1.898
Mean :1.898
3rd Qu.:1.898
Max. :1.899
> attributes(lambert_obj)
$dim
[1] 150 1
$dimnames
$dimnames[[1]]
NULL
$dimnames[[2]]
[1] "Y1.X"
$`Gaussianized:mu`
[1] 1.898343
$`Gaussianized:sigma`
[1] 6.4e-05
$`Gaussianized:delta_l`
[1] 2.643938e+15
$`Gaussianized:delta_r`
[1] 2.558455
$`Gaussianized:alpha_l`
[1] 1
$`Gaussianized:alpha_r`
[1] 1Session Info
sessionInfo()
R Under development (unstable) (2023-08-09 r84924 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8 LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C LC_TIME=English_United States.utf8
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lamW_2.2.0
loaded via a namespace (and not attached):
[1] vctrs_0.6.3 cli_3.6.1 rlang_1.1.1 stringi_1.7.12 generics_0.1.3 RcppParallel_5.1.7 glue_1.6.2
[8] colorspace_2.1-0 plyr_1.8.8 scales_1.2.1 fansi_1.0.4 grid_4.4.0 munsell_0.5.0 tibble_3.2.1
[15] MASS_7.3-60.1 LambertW_0.6.7-1 reshape2_1.4.4 lifecycle_1.0.3 stringr_1.5.0 compiler_4.4.0 dplyr_1.1.2
[22] RColorBrewer_1.1-3 Rcpp_1.0.11 pkgconfig_2.0.3 rstudioapi_0.15.0 R6_2.5.1 tidyselect_1.2.0 utf8_1.2.3
[29] pillar_1.9.0 magrittr_2.0.3 tools_4.4.0 gtable_0.3.3 ggplot2_3.4.2 aadler commented
Also, I run the following tests via tinytest every time I make a change, and all three methods (Halley, Fritsch Pade, and polynomial minimax) pass.
# Copyright (c) 2015, Avraham Adler All rights reserved
# SPDX-License-Identifier: BSD-2-Clause
tol <- sqrt(.Machine$double.eps)
# Test that functions return proper values
principleBranchAnswers <- runif(5000, min = -1, max = 703.22703310477016)
principleBranchTests <- principleBranchAnswers * exp(principleBranchAnswers)
secondaryBranchAnswers <- runif(5000, min = -714.96865723796657, max = -1)
secondaryBranchTests <- secondaryBranchAnswers * exp(secondaryBranchAnswers)
# Test that function works properly in general
expect_equal(lambertW0(principleBranchTests), principleBranchAnswers,
tolerance = tol)
expect_equal(lambertWm1(secondaryBranchTests), secondaryBranchAnswers,
tolerance = tol)
# Test that function works properly for larger numbers
expect_equal(lambertW0(1000) * exp(lambertW0(1000)), 1000, tolerance = tol)
# Test that function behaves properly near 0
V0 <- seq(-2e-2, 2e-2, 2e-6)
V0E <- V0 * exp(V0)
expect_equal(lambertW0(V0E), V0, tolerance = tol)
# Test that W0 behaves properly VERY close to 0
expect_identical(lambertW0(1e-275), 1e-275)
expect_identical(lambertW0(7e-48), 7e-48)
expect_identical(lambertW0(-3.81e-71), -3.81e-71)
# Test that function behaves properly near -1/e
expect_identical(lambertW0(-1 / exp(1)), -1)
expect_identical(lambertWm1(-1 / exp(1)), -1)
# Test that function behaves properly near its asymptotes
L <- seq(1e-6 - exp(-1), -0.25, 3e-6)
V0 <- lambertW0(L)
vm1 <- lambertWm1(L)
expect_equal(V0 * exp(V0), L, tolerance = tol)
expect_equal(vm1 * exp(vm1), L, tolerance = tol)
vm1 <- seq(-714, -714.96865, -3e-5)
vm1E <- vm1 * exp(vm1)
expect_equal(lambertWm1(vm1E), vm1, tolerance = tol)
# Test that function behaves properly at its asymptotes
expect_identical(lambertW0(Inf), Inf)
expect_identical(lambertWm1(0), -Inf)
# Test that NaNs are returned for values outside domain
expect_true(is.nan(lambertW0(-Inf)))
expect_true(is.nan(lambertW0(-1)))
expect_true(is.nan(lambertW0(c(1, -1)))[[2]])
expect_true(is.nan(lambertWm1(-Inf)))
expect_true(is.nan(lambertWm1(Inf)))
expect_true(is.nan(lambertWm1(-0.5))) # x < -M_1_E
expect_true(is.nan(lambertWm1(1.2))) # x > 0
# Test that integers are converted to reals for principle branch
expect_identical(lambertW0(c(-1L, 0L, 1L, 2L, 3L, 4L)),
lambertW0(c(-1, 0, 1, 2, 3, 4)))
aadler commented
@gmgeorg and @petersonR, unfortunately, I have even more bad news. I went to the CRAN archives and installed lamW 2.1.2 from source on my machine---for which there were no errors for either of your packages, I believe---then reinstalled LambertW from source, and I still get the same problem. See below for details. I don't think is a lamW issue. I tried to trace how Gaussianize works, but this is not my area of expertise. Please see below for some cursory findings.
I am very happy to try and help debug the issue, but unless I've missed something, I think the timing of the error and the lamW update is more serendipitous than anything else. Thanks.
Cursory debugging
Problem case
Now a little investigation shows this only occurs when type = "hh" and 'method = "IGMM"`. This is borne out by the following:
> LambertW::IGMM(train, "hh")
Call: LambertW::IGMM(y = train, type = "hh")
Estimation method: IGMM
Input distribution:
mean-variance Lambert W x F type ('h' same tails; 'hh' different tails; 's' skewed): hh
Parameter estimates:
mu_x sigma_x delta_l delta_r
1.898343e+00 6.437686e-05 2.712248e+15 2.558455e+00
Obtained after 100 iterations.When type = "hh", there are a number of calls affected in IGMM.
mLambertW
if (any(type == c("h", "hh"))) {
lamW.U.sd <- mLambertW(theta = list(beta = c(0, 1), delta = mean.default(tau.init[grepl("delta",
names(tau.init))])), distname = "normal")$sdSome set up:
y <- train
type <- "hh"
skewness.x = 0
kurtosis.x = 3
tau.init = LambertW::get_initial_tau(y, type)
robust = FALSE
tol = .Machine$double.eps^0.25
location.family = TRUE
not.negative = NULL
max.iter = 100
delta.lower = -1; delta.upper = 3
LambertW::mLambertW(theta = list(beta = c(0, 1), delta = mean.default(tau.init[grepl("delta", names(tau.init))])), distname = "normal")$sd
[1] 1I don't know if 1 is improper here, but this is the result.
lp_norm
This enters the equation here
else if (type == "hh") {
tau.trace <- rbind(0, tau.init[c("mu_x", "sigma_x", "delta_l",
"delta_r")])
colnames(tau.trace) <- c("mu_x", "sigma_x", "delta_l",
"delta_r")
while (lp_norm(tau.trace[kk + 1, ] - tau.trace[kk, ],
p = 2) > tol && kk < max.iter) {
zz <- normalize_by_tau(y, tau.trace[kk + 1, ])
DEL <- delta_GMM(zz, delta.init = tau.trace[kk +
1, c("delta_l", "delta_r")], kurtosis.x = kurtosis.x,
tol = tol, not.negative = not.negative, type = "hh",
lower = delta.lower, upper = delta.upper)
delta.hat <- DEL$delta
uu <- W_2delta(zz, delta.hat)
xx <- normalize_by_tau(uu, tau.trace[kk + 1, ], inverse = TRUE)
tau.trace <- rbind(tau.trace, c(mean(xx), sd(xx),
delta.hat))
if (!location.family) {
tau.trace[nrow(tau.trace), "mu_x"] <- 0
}
total.iter <- total.iter + DEL$iterations
kk <- kk + 1
}
se <- c(1, sqrt(1/2), 1, 1)/sqrt(out$n)
zz.tmp <- normalize_by_tau(y, tau.trace[nrow(tau.trace),
])
delta.hat <- delta_GMM(zz.tmp, delta.init = tail(tau.trace[,
c("delta_l", "delta_r")], 1), kurtosis.x = kurtosis.x,
skewness.x = skewness.x, tol = tol, not.negative = not.negative,
type = "hh")$delta
tau.trace <- rbind(tau.trace, c(tail(tau.trace[, c("mu_x",
"sigma_x")], 1), delta.hat))
}I'm not sure how to test this one.
delta_GMM
Also appears in the above. Looks ok?
> LambertW::delta_GMM(train, "hh")
$iterations
[1] 8
$delta
[1] 0.00000000 -0.02733716Reinstalling old lamW and failure remains
> install.packages("C:/R/lamW_2.1.2.tar.gz", repos = NULL)
* installing *source* package 'lamW' ...
** package 'lamW' successfully unpacked and MD5 sums checked
** using staged installation
** libs
using C++ compiler: 'g++.exe (GCC) 12.2.0'
g++ -std=gnu++17 -I"C:/R/RCurrent/R-devel/include" -DNDEBUG -I'C:/R/RCurrent/R-devel/library/Rcpp/include' -I'C:/R/RCurrent/R-devel/library/RcppParallel/include' -I"C:/rtools43/x86_64-w64-mingw32.static.posix/include" -DRCPP_PARALLEL_USE_TBB=1 -O2 -march=native -pipe -mno-rtm -c RcppExports.cpp -o RcppExports.o
g++ -std=gnu++17 -I"C:/R/RCurrent/R-devel/include" -DNDEBUG -I'C:/R/RCurrent/R-devel/library/Rcpp/include' -I'C:/R/RCurrent/R-devel/library/RcppParallel/include' -I"C:/rtools43/x86_64-w64-mingw32.static.posix/include" -DRCPP_PARALLEL_USE_TBB=1 -O2 -march=native -pipe -mno-rtm -c lambertW.cpp -o lambertW.o
g++ -std=gnu++17 -shared -O2 -march=native -pipe -mno-rtm -s -static-libgcc -o lamW.dll tmp.def RcppExports.o lambertW.o -LC:/R/RCurrent/R-devel/library/RcppParallel/lib/x64 -ltbb -ltbbmalloc -LC:/rtools43/x86_64-w64-mingw32.static.posix/lib/x64 -LC:/rtools43/x86_64-w64-mingw32.static.posix/lib -LC:/R/RCurrent/R-devel/bin/x64 -lR
installing to C:/R/RCurrent/R-devel/library/00LOCK-lamW/00new/lamW/libs/x64
** R
** inst
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
** building package indices
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (lamW)
> packageVersion("lamW")
[1] ‘2.1.2’
> install.packages("LambertW", type = "source")
trying URL 'https://cran.rstudio.com/src/contrib/LambertW_0.6.7-1.tar.gz'
Content type 'application/x-gzip' length 439401 bytes (429 KB)
downloaded 429 KB
* installing *source* package 'LambertW' ...
** package 'LambertW' successfully unpacked and MD5 sums checked
** using staged installation
** libs
using C++ compiler: 'g++.exe (GCC) 12.2.0'
g++ -std=gnu++17 -I"C:/R/RCurrent/R-devel/include" -DNDEBUG -I'C:/R/RCurrent/R-devel/library/Rcpp/include' -I'C:/R/RCurrent/R-devel/library/lamW/include' -I"C:/rtools43/x86_64-w64-mingw32.static.posix/include" -O2 -march=native -pipe -mno-rtm -c RcppExports.cpp -o RcppExports.o
g++ -std=gnu++17 -I"C:/R/RCurrent/R-devel/include" -DNDEBUG -I'C:/R/RCurrent/R-devel/library/Rcpp/include' -I'C:/R/RCurrent/R-devel/library/lamW/include' -I"C:/rtools43/x86_64-w64-mingw32.static.posix/include" -O2 -march=native -pipe -mno-rtm -c W_Cpp.cpp -o W_Cpp.o
g++ -std=gnu++17 -I"C:/R/RCurrent/R-devel/include" -DNDEBUG -I'C:/R/RCurrent/R-devel/library/Rcpp/include' -I'C:/R/RCurrent/R-devel/library/lamW/include' -I"C:/rtools43/x86_64-w64-mingw32.static.posix/include" -O2 -march=native -pipe -mno-rtm -c W_delta_Cpp.cpp -o W_delta_Cpp.o
g++ -std=gnu++17 -I"C:/R/RCurrent/R-devel/include" -DNDEBUG -I'C:/R/RCurrent/R-devel/library/Rcpp/include' -I'C:/R/RCurrent/R-devel/library/lamW/include' -I"C:/rtools43/x86_64-w64-mingw32.static.posix/include" -O2 -march=native -pipe -mno-rtm -c W_gamma_Cpp.cpp -o W_gamma_Cpp.o
g++ -std=gnu++17 -I"C:/R/RCurrent/R-devel/include" -DNDEBUG -I'C:/R/RCurrent/R-devel/library/Rcpp/include' -I'C:/R/RCurrent/R-devel/library/lamW/include' -I"C:/rtools43/x86_64-w64-mingw32.static.posix/include" -O2 -march=native -pipe -mno-rtm -c kurtosis.cpp -o kurtosis.o
g++ -std=gnu++17 -I"C:/R/RCurrent/R-devel/include" -DNDEBUG -I'C:/R/RCurrent/R-devel/library/Rcpp/include' -I'C:/R/RCurrent/R-devel/library/lamW/include' -I"C:/rtools43/x86_64-w64-mingw32.static.posix/include" -O2 -march=native -pipe -mno-rtm -c lp_norm_Cpp.cpp -o lp_norm_Cpp.o
g++ -std=gnu++17 -I"C:/R/RCurrent/R-devel/include" -DNDEBUG -I'C:/R/RCurrent/R-devel/library/Rcpp/include' -I'C:/R/RCurrent/R-devel/library/lamW/include' -I"C:/rtools43/x86_64-w64-mingw32.static.posix/include" -O2 -march=native -pipe -mno-rtm -c normalize_by_tau_Cpp.cpp -o normalize_by_tau_Cpp.o
g++ -std=gnu++17 -I"C:/R/RCurrent/R-devel/include" -DNDEBUG -I'C:/R/RCurrent/R-devel/library/Rcpp/include' -I'C:/R/RCurrent/R-devel/library/lamW/include' -I"C:/rtools43/x86_64-w64-mingw32.static.posix/include" -O2 -march=native -pipe -mno-rtm -c skewness.cpp -o skewness.o
g++ -std=gnu++17 -shared -O2 -march=native -pipe -mno-rtm -s -static-libgcc -o LambertW.dll tmp.def RcppExports.o W_Cpp.o W_delta_Cpp.o W_gamma_Cpp.o kurtosis.o lp_norm_Cpp.o normalize_by_tau_Cpp.o skewness.o -LC:/rtools43/x86_64-w64-mingw32.static.posix/lib/x64 -LC:/rtools43/x86_64-w64-mingw32.static.posix/lib -LC:/R/RCurrent/R-devel/bin/x64 -lR
installing to C:/R/RCurrent/R-devel/library/00LOCK-LambertW/00new/LambertW/libs/x64
** R
** data
** inst
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
** building package indices
** installing vignettes
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (LambertW)
The downloaded source packages are in
‘D:\Rtemp\Rtmp0ACopw\downloaded_packages’
> packageVersion("lambertW")
[1] ‘0.6.7.1’
> data(iris)
>
> train <- iris$Petal.Width
> summary(train)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.100 0.300 1.300 1.199 1.800 2.500
> lambert_obj <- LambertW::Gaussianize(train, type = "hh")
> summary(lambert_obj)
Y1.X
Min. :1.898
1st Qu.:1.898
Median :1.898
Mean :1.898
3rd Qu.:1.898
Max. :1.899
> attributes(lambert_obj)
$dim
[1] 150 1
$dimnames
$dimnames[[1]]
NULL
$dimnames[[2]]
[1] "Y1.X"
$`Gaussianized:mu`
[1] 1.898343
$`Gaussianized:sigma`
[1] 6.4e-05
$`Gaussianized:delta_l`
[1] 2.712248e+15
$`Gaussianized:delta_r`
[1] 2.558455
$`Gaussianized:alpha_l`
[1] 1
$`Gaussianized:alpha_r`
[1] 1aadler commented
Any findings?
gmgeorg commented
Thanks for looking into this in so much detail.
I could replicate this issue on my machine w/ existing LambertW package (w/o lamW update). So this does seem like a coincidence. Not sure why the error just happens now in bestNormalize package and not previously.
I root-caused this to the starting values of the optimization from IGMM. IGMM() diverges to delta_l -> infinity
LambertW::IGMM(train, type="hh")
Call: LambertW::IGMM(y = train, type = "hh")
Estimation method: IGMM
Input distribution:
mean-variance Lambert W x F type ('h' same tails; 'hh' different tails; 's' skewed): hh
Parameter estimates:
mu_x sigma_x delta_l delta_r
1.898343e+00 6.437686e-05 2.973560e+15 2.558455e+00
Obtained after 100 iterations.
The issue is that IGMM() (and related functions) don't properly pass the delta.lower and delta.upper limits (which default to -1 and 3 respectively). This means the MLE gets starting values from IGMM() that are already degenerate and can't get out of this degenerate solution by likelihood optimization.
I verified that if I pass starting values that overwrite the degenerate (delta_l, delta_r) values from IGMM, the MLE gives non-trivial/non-degenerate results
tt = LambertW::get_initial_theta(train, distname="normal", type="hh", method="IGMM")
tt$delta = c("delta_r" = 0.001, "delta_l" = 0.001)
mod <- LambertW::MLE_LambertW(train, distname="normal", type="hh", theta.init = tt)
Call: LambertW::MLE_LambertW(y = train, distname = "normal", type = "hh",
theta.init = tt)
Estimation method: MLE
Input distribution: normal
mean-variance Lambert W x F type ('h' same tails; 'hh' different tails; 's' skewed): hh
Parameter estimates:
mu sigma delta_r delta_l
1.38384 0.66008 0.99165 -0.12424
Fix is up #5 ; need some more testing/review (and then push up to CRAN)
gmgeorg commented
@petersonR @aadler I merged PR and submitted v0.6.8 to CRAN. Should be updated in next couple of days.
gmgeorg commented
Closing this out as v0.6.8 should fix this issue now. @petersonR feel free to re-open if your tests still run into this issue after the upgrade
petersonR commented
Great work - thank you both for looking so hard into this and addressing!