From e588dd26b7356ad7d029941ab7cbe30ebe31ff4f Mon Sep 17 00:00:00 2001 From: "Daniel P. Palomar" Date: Sat, 9 Apr 2022 10:11:45 +0800 Subject: [PATCH] Making besselK more robust... --- R/fit_mvst.R | 84 +++++++++++++++---------- tests/testthat/fitted_mvst_check.RData | Bin 712 -> 708 bytes 2 files changed, 50 insertions(+), 34 deletions(-) diff --git a/R/fit_mvst.R b/R/fit_mvst.R index f2e3eb6..0c8597d 100644 --- a/R/fit_mvst.R +++ b/R/fit_mvst.R @@ -271,7 +271,7 @@ dmvst <- function(X, nu = 3, gamma = 1, mu = 0, scatter = 1) { kappa <- sqrt(nu + rowSums(Xc * (Xc %*% scatter_inv))) first_term <-Xc %*% solve(scatter) %*% gamma - (N/2) * log(2*pi) - 0.5 * sum(log(eigen((scatter))$values) ) second_term <- log(2) + (nu/2) * log(nu/2) - lgamma(nu/2) - third_term <- -((nu+N)/2) * log(kappa/delta) + log_besselK(delta*kappa, -((nu+N)/2)) + third_term <- -((nu+N)/2) * log(kappa/delta) + log_besselK(delta*kappa, -((nu + N)/2)) return(first_term + second_term + third_term) } } else # Gaussian case @@ -300,20 +300,21 @@ Estep_mvst <- function(X, nu, gamma_unscaled, mu, scatter_unscaled, alpha) { # kappa <- sqrt(nu + stats::mahalanobis(x = X, center = mu, cov = scatter_inv, inverted = TRUE)) if (delta == 0) { - E_tau <- ((kappa**2)/2)*lambda - E_invtau <- (1/((kappa**2)/2)) * (1/lambda) - E_logtau <- digamma(lambda) - log((kappa**2)/4) + E_tau <- (nu + N)/(kappa**2) + E_invtau <- (kappa**2)/(nu + N -2) + E_logtau <- digamma(lambda) - log((kappa**2)/2) } else { - tmp <- besselK_ratio(delta * kappa, lmd = lambda) + tmp <- besselK_ratio(delta * kappa, lambda) E_tau <- (delta / kappa) * tmp #E_invtau <- (kappa / delta) * 1/besselK_ratio(delta * kappa, lmd = lambda - 1) E_invtau <- (kappa / delta) * (tmp - (2*lambda)/delta/kappa) # this saves computing bessel functions again if (lambda < 150) { - dev_cal <- function(val) numDeriv::grad(func = function(lmd) log(besselK(x = val, nu = lmd, expon.scaled = FALSE)), x = lambda, method = "simple", method.args = list(eps = 1e-10)) + dev_cal <- function(val) numDeriv::grad(func = function(lmd) log(besselK(x = val, nu = lmd, expon.scaled = TRUE)), + x = lambda, method = "simple", method.args = list(eps = 1e-10)) E_logtau <- log(delta / kappa) + sapply(delta * kappa, dev_cal) } else - E_logtau <- log(delta / kappa) + log(besselK_ratio(delta * kappa, lmd = lambda)) + E_logtau <- log(delta / kappa) + log(besselK_ratio(delta * kappa, lambda)) } # return @@ -321,8 +322,10 @@ Estep_mvst <- function(X, nu, gamma_unscaled, mu, scatter_unscaled, alpha) { "E_invtau" = E_invtau / alpha, "E_logtau" = E_logtau + log(alpha)) - if (any(is.infinite(list_to_return$E_invtau)) || is.infinite(sum(list_to_return$E_invtau)) || - any(is.nan(list_to_return$E_invtau))) { + if (any(is.infinite(list_to_return$E_invtau)) || + any(is.nan(list_to_return$E_invtau)) || + any(is.nan(list_to_return$E_logtau)) + ) { message("Problem with the computation of E[tau], probably because of very small numbers in the evaluation of the bessel function.") browser() } @@ -332,17 +335,16 @@ Estep_mvst <- function(X, nu, gamma_unscaled, mu, scatter_unscaled, alpha) { - # https://www.researchgate.net/journal/Journal-of-Inequalities-and-Applications-1029-242X[…]mating-the-modified-Bessel-function-of-the-second-kind.pdf -besselK_ratio <- function(x, lmd) { - if (lmd < 100) - return(besselK(x = x, nu = lmd + 1, expon.scaled = TRUE) / besselK(x = x, nu = lmd, expon.scaled = TRUE)) +besselK_ratio <- function(x, nu) { + if (nu < 51) + return(besselK(x = x, nu = nu + 1, expon.scaled = TRUE) / besselK(x = x, nu = nu, expon.scaled = TRUE)) else { - lmd_i <- lmd - floor(lmd) + 10 - R_i <- besselK(x = x, nu = lmd_i + 1, expon.scaled = TRUE) / besselK(x = x, nu = lmd_i, expon.scaled = TRUE) - while (lmd_i != lmd) { - R_i <- 1/R_i + (2 * lmd_i + 2)/x - lmd_i <- lmd_i + 1 + nu_i <- nu - floor(nu) + 50 + R_i <- besselK(x = x, nu = nu_i + 1, expon.scaled = TRUE) / besselK(x = x, nu = nu_i, expon.scaled = TRUE) + while (nu_i != nu) { + R_i <- 1/R_i + (2 * nu_i + 2)/x + nu_i <- nu_i + 1 } return(R_i) } @@ -350,22 +352,36 @@ besselK_ratio <- function(x, lmd) { -log_besselK <- function(x, lmd) { - if (lmd >= 0) - stop("lmd should be negative in this function.") - lmd_i <- lmd - floor(lmd) - 1 - K_lmd_i <- besselK(x, nu = lmd_i) # K_{lmd_i}(x) - log_values <- log(K_lmd_i) - R_lmd_i <- K_lmd_i/besselK(x, nu = lmd_i - 1) # R_{lmd_i -1}(x) - log_values <- log_values - log(R_lmd_i) - lmd_i <- lmd_i - 1 - while(lmd_i != lmd) { - R_lmd_i_inv <- R_lmd_i - 2*lmd_i/x - R_lmd_i <- 1/R_lmd_i_inv # R_{i -1}(x) - log_values <- log_values - log(R_lmd_i) - lmd_i <- lmd_i - 1 +log_besselK <- function(x, nu) { + if (nu <= 10 && nu >= -1) { + return(log(besselK(x, nu))) + } else if (nu >= 0){ + nu_i <- nu - floor(nu) + 9 + K_nu_i <- besselK(x, nu = nu_i) # K_{nu_i}(x) + log_values <- log(K_nu_i) + R_nu_i <- besselK(x, nu = nu_i + 1)/K_nu_i # R_{nu_i}(x) + log_values <- log_values + log(R_nu_i) + nu_i <- nu_i + 1 + while (nu_i != nu) { + R_nu_i <- 1/R_nu_i + 2 * nu_i /x + log_values <- log_values + log(R_nu_i) + nu_i <- nu_i + 1 + } + return(log_values) + } else { + nu_i <- nu - floor(nu) - 9 + K_nu_i <- besselK(x, nu = nu_i) # K_{nu_i}(x) + log_values <- log(K_nu_i) + R_nu_i <- K_nu_i/besselK(x, nu = nu_i - 1) # R_{nu_i -1}(x) + log_values <- log_values - log(R_nu_i) + nu_i <- nu_i - 1 + while (nu_i != nu) { + R_nu_i_inv <- R_nu_i - 2*nu_i/x + R_nu_i <- 1/R_nu_i_inv # R_{i -1}(x) + log_values <- log_values - log(R_nu_i) + nu_i <- nu_i - 1 + } + return(log_values) } - return(log_values) } - diff --git a/tests/testthat/fitted_mvst_check.RData b/tests/testthat/fitted_mvst_check.RData index 89f443d781b71132904254424e9e985bc09307ca..301a406c1f710cbcbb68e06c0e533e404429a14f 100644 GIT binary patch delta 633 zcmV-<0*3v_1;hoA907-s9X@|c5HTY84yx5fmj)>6G{yIXYjNWx3m$Xa`)G@q%=k6o z)ZWWMg*4-2gv&Z99@7%WgBgS0a~fS<#Y(vYC48437d#*3J;oCfu70Csv^i?`Y@WDJ zqQ3y15E~Z%krzG7EmwLKykAY8L! zp~H$PBm(Pn9>FaA&trf3M*E$8#kt;bC6EtREV+ETqFku17^%iM3KsHQW{hZ$K~3|-|=`Q`NL>KKb~k8+Xfq2+b6gqO4cd*sEXCZR-&M6 zKB1_;WifqS!x(Km1;V7iSZ|#tbGjuhmYWkES%m0X8e-|oK*x%Udg#Ga66-%XSmM0& z;V1DBwhTo9Td+8X+q~n(wUspUAMTz~!K@~iw=yn`w4RO7XUqTq;^3Jm00Ex`&J}?c70ssI200CKA`fM}~ delta 637 zcmV-@0)qX-1;_=E907@u9X@|b?iJvt$7GAZ*fZr~-Z2!FqKeRQYakX^odD?cx~UKP zN*WAlX2IAxZ(j+a$?pLZYoG)px37q+^vsNm;Hkp(><-vc&1g3*CzP=1bfb_17c2L7 zoiW*%Xycb(VX1|aKb8~?zq@We>10U+bNI4{~jp84S7ko z@5LR{S28v91^w>SM@fI(D_72H;nLQX5IY{}fJ@OdMu+&-tf%l+uV~i^DNoT@|9yMg ztaL%v2F8;=DtKi!gz{qj?Fg;@Cc^>;&vI^CoFSMecdoLrng+kj-}PAM4jORvz1ViO zn_jiky^^N`+Rrc>VK~`!B$1wQ|IX$%=vm&4@3QXZH#MK%%V2-MUtSZHb^xi-@(e%V zp{P)_%P7)Q$C;C6atr0IgU^^El%bCVkhM*)Vcnd8d+F|?Xl|~4o#*1I+ff8$J+Q)=XytSClURmySu)n$v-OExrr_npWwrtaeo=g=b4Vkee8IMv{k zbW1yO=Hpw4#9V68lo~Z_oKJ^AGi@lK#?-1DIPFo^wH&m0aJ24500000Cm?%Y00E%| X&