"pu"<- function(quant, p, m, n) { subintervals <- 20000 parameter <- c(m, p, m + n - p) if(p > m) { p <- parameter[1] m <- parameter[2] n <- parameter[3] } if(p == 1) prob <- pbeta(quant, n/2, m/2) if(p == 2) prob <- pbeta(sqrt(quant), n - 1, m) if(p >= 3) { r <- p %/% 2 tol <- 1e-08 if(p %% 2 == 0) M <- max( - log(qbeta(tol^(1/r), n + 1 - 2 * (1:r), m)) ) * r else M <- max( - log(qbeta(tol^(1/(r + 1)), n + 1 - 2 * (1:r), m)), -0.5 * log(qbeta(tol^(1/(r + 1)), (n - p + 1)/2, m/2))) * (r + 1) step <- M/subintervals x <- seq(0, M, by = step) f <- fft(dbeta(exp( - x), n - 1, m) * exp( - x)) if(r >= 2) for(i in 2:r) f <- fft(dbeta(exp( - x), n + 1 - 2 * i, m) * exp( - x)) * f if(p %% 2 == 1) f <- fft(dbeta(exp(-2 * x), (n - p + 1)/2, m/2) * 2 * exp(-2 * x)) * f f <- Re(fft(f, inverse = T)) cte <- sum(f) * step f <- f/cte prob <- sum(f[x >= - log(sqrt(quant))]) * step } prob } "qu"<- function(prob, p, m, n) { subintervals <- 20000 parameter <- c(m, p, m + n - p) if(p > m) { p <- parameter[1] m <- parameter[2] n <- parameter[3] } if(p == 1) crit.point.u <- qbeta(prob, n/2, m/2) if(p == 2) crit.point.u <- qbeta(prob, n - 1, m)^2 if(p >= 3) { r <- p %/% 2 tol <- 1e-08 if(p %% 2 == 0) M <- max( - log(qbeta(tol^(1/r), n + 1 - 2 * (1:r), m)) ) * r else M <- max( - log(qbeta(tol^(1/(r + 1)), n + 1 - 2 * (1:r), m)), -0.5 * log(qbeta(tol^(1/(r + 1)), (n - p + 1)/2, m/2))) * (r + 1) step <- M/subintervals x <- seq(0, M, by = step) f <- fft(dbeta(exp( - x), n - 1, m) * exp( - x)) if(r >= 2) for(i in 2:r) f <- fft(dbeta(exp( - x), n + 1 - 2 * i, m) * exp( - x)) * f if(p %% 2 == 1) f <- fft(dbeta(exp(-2 * x), (n - p + 1)/2, m/2) * 2 * exp(-2 * x)) * f f <- Re(fft(f, inverse = T)) cte <- sum(f) * step f <- f/cte * step cumprob <- cumsum(f) index <- max((1:length(x))[cumprob <= (1 - prob)]) crit.point.x <- x[index + 1] - ((x[index + 1] - x[index]) * ( cumprob[index + 1] - (1 - prob)))/(cumprob[index + 1] - cumprob[index]) crit.point.u <- exp(-2 * crit.point.x) } Cfactor <- ( - (n - (p - m + 1)/2) * log(crit.point.u))/qchisq(1 - prob, p * m) list(critical.point = crit.point.u, Cfactor = Cfactor) } "constante"<- function(p, r, cte0 = 0.01,cte1=100) { gammap <- gamma(p/2) fonction <- function(x, cte, gammap, p) { y <- 2^(-p/2+1)*x^(p - 1) *exp( - x^2/2)/gammap ind <- (x <= cte) x[ind] <- 3*x[ind]^2/cte^2 -3*x[ind]^4/(cte^4) + x[ind]^6/ (cte^6) x[!ind] <- 1 x * y } epsilon <- 0.001 precision <- 1 fonction0_integrate(fonction,0,Inf,cte=cte0,gammap=gammap,p=p)$integral-r fonction1_integrate(fonction,0,Inf,cte=cte1,gammap=gammap,p=p)$integral-r while(precision > epsilon) { cte2 <- (cte0 + cte1)/2 fonction2_integrate(fonction,0,Inf,cte=cte2,gammap=gammap,p=p)$integral-r if(fonction2<0) cte1_cte2 else cte0_cte2 precision <- 100*abs(cte1 - cte0)/cte0 } cte0 } "ale2"<-function(a1, b1, a2, b2, nr) { va2 <- as.vector(a2) vb2 <- as.vector(b2) nr1 <- nr + 1 div <- seq(0, 1, length = nr1) lambda <- runif(nr, div[ - nr1], div[-1]) out1 <- a1 + outer(b1 - a1, lambda) out2 <- va2 + outer(vb2 - va2, lambda) list(mu = out1, sigma = out2) } "M.scale"<-function(x, b1, cte, s0) { epsilon <- 0.001 precision <- 1 while(precision > epsilon) { s1 <- s0 + (mean(rho(x/s0, cte)) - b1)/mean((psi(x/s0, cte) * x)/ s0^2) o1 <- min(abs(x)/s1) while((s1 < 0) | (o1 > cte)) { s1 <- (s0 + s1)/2 o1 <- min(abs(x)/s1) } while(abs(mean(rho(x/s1, cte)) - b1) - abs(mean(rho(x/s0, cte)) - b1) > 0) { s1 <- (s0 + s1)/2 } precision <- abs(s1 - s0)/s0 s0 <- s1 } s0 } "delta12"<-function(x, mu, C, cte, s, k) { p <- ncol(x) MDstand <- sqrt(mahalanobis(x, mu, C))/s um <- u(MDstand, cte) mu <- apply(um * x, 2, sum)/sum(um) sigma <- (p * (t(x) - mu) %*% (um * t(t(x) - mu)))/sum(um * MDstand^2 - rho(MDstand, cte) + k) list(mu = mu, sigma = sigma) } "det"<-function(A) { prod(eigen(A)$values) } "rho"<-function(d, cte) { ind <- (abs(d) <= cte) d[ind] <- d[ind]^2/2 - d[ind]^4/(2 * cte^2) + d[ind]^6/(6 * cte^4) d[!ind] <- cte^2/6 d } "psi"<-function(d, cte) { ind <- (abs(d) <= cte) d[ind] <- d[ind] * ((d[ind]/cte)^2 - 1)^2 d[!ind] <- 0 d } "u"<-function(d, cte) { ind <- (d == 0) d[ind] <- 1 d[!ind] <- psi(d[!ind], cte)/d[!ind] d } "s.estimate"<- function(x, r, nr, Nsamp) { nr1 <- nr + 1 s.tilde <- 10^9 p <- ncol(x) n <- nrow(x) cte <- constante( p, r, cte0 = .01,cte1=100) k <- (r * cte^2)/6 for(l in (1:Nsamp)) { rang <- 0 while(rang != p) { J <- sample(n, p + 1) muJ0 <- apply(x[J, ], 2, mean) sigmaJ0 <- var(x[J, ]) rang <- qr(sigmaJ0)$rank } if(l > 1) { musigmaJ <- ale2(as.vector(muJ0), as.vector(mu.tilde), sigmaJ0, sigma.tilde, nr) muJ <- cbind(as.vector(muJ0), musigmaJ$mu) sigmaJ <- cbind(as.vector(sigmaJ0), musigmaJ$sigma) for(i in (1:nr1)) { muJi <- muJ[, i] sigmaJi <- matrix(sigmaJ[, i], ncol = p, nrow = p) CJi <- (det(sigmaJi))^(-1/p) * sigmaJi MD.2 <- sqrt(mahalanobis(x, muJi, CJi)) if(mean(rho(MD.2/s.tilde, cte)) - k < 0) { mu.tilde <- muJi C.tilde <- CJi MD <- MD.2 s.tilde <- M.scale(MD, k, cte, s.tilde) sigma.tilde <- s.tilde^2 * C.tilde print("mu=") print(mu.tilde) print("V=") print(sigma.tilde) print("determinant=") print(s.tilde^(2 * p)) delta <- delta12(x, mu.tilde, C.tilde, cte, s.tilde, k) condition <- 0 val.m <- 0 while((condition) == 0 & (val.m <= 10)) { mu.2tilde <- mu.tilde * (1 - 2^( - val.m)) + delta$mu * 2^( - val.m) sigma.2tilde <- sigma.tilde * (1 - 2^( - val.m)) + delta$sigma * 2^( - val.m) C.2tilde <- (det(sigma.2tilde))^(-1/p) * sigma.2tilde MD.2 <- sqrt(mahalanobis(x, mu.2tilde, C.2tilde)) if(mean(rho(MD.2/s.tilde, cte)) - k < 0) { condition <- 1 mu.tilde <- mu.2tilde C.tilde <- C.2tilde MD <- MD.2 s.tilde <- M.scale(MD, k, cte, s.tilde) sigma.tilde <- s.tilde^2 * C.tilde print("mu=") print(mu.tilde) print("V=") print(sigma.tilde) print("determinant=") print(s.tilde^(2 * p)) } else val.m <- val.m + 1 } } } } else { mu.tilde <- muJ0 C.tilde <- det(sigmaJ0)^(-1/p) * sigmaJ0 MD <- sqrt(mahalanobis(x, mu.tilde, C.tilde)) s.tilde <- M.scale(MD, k, cte, s.tilde) sigma.tilde <- s.tilde^2 * C.tilde print("mu=") print(mu.tilde) print("V=") print(sigma.tilde) print("determinant=") print(s.tilde^(2 * p)) delta <- delta12(x, mu.tilde, C.tilde, cte, s.tilde, k) condition <- 0 val.m <- 0 while((condition == 0) & (val.m <= 10)) { mu.2tilde <- mu.tilde * (1 - 2^( - val.m)) + delta$mu * 2^( - val.m) sigma.2tilde <- sigma.tilde * (1 - 2^( - val.m) ) + delta$sigma * 2^( - val.m) C.2tilde <- (det(sigma.2tilde))^(-1/p) * sigma.2tilde MD.2 <- sqrt(mahalanobis(x, mu.2tilde, C.2tilde )) if(mean(rho(MD.2/s.tilde, cte)) - k < 0) { condition <- 1 mu.tilde <- mu.2tilde C.tilde <- C.2tilde MD <- MD.2 s.tilde <- M.scale(MD, k, cte, s.tilde) sigma.tilde <- s.tilde^2 * C.tilde print("mu=") print(mu.tilde) print("V=") print(sigma.tilde) print("determinant=") print(s.tilde^(2 * p)) } else val.m <- val.m + 1 } } } for(i in 1:3) { delta <- delta12(x, mu.tilde, C.tilde, cte, s.tilde, k) condition <- 0 val.m <- 0 while((condition == 0) & (val.m <= 10)) { mu.2tilde <- mu.tilde * (1 - 2^( - val.m)) + delta$mu * 2^( - val.m) sigma.2tilde <- sigma.tilde * (1 - 2^( - val.m)) + delta$sigma * 2^( - val.m) C.2tilde <- (det(sigma.2tilde))^(-1/p) * sigma.2tilde MD.2 <- sqrt(mahalanobis(x, mu.2tilde, C.2tilde)) if(mean(rho(MD.2/s.tilde, cte)) - k < 0) { condition <- 1 mu.tilde <- mu.2tilde C.tilde <- C.2tilde MD <- MD.2 s.tilde <- M.scale(MD, k, cte, s.tilde) sigma.tilde <- s.tilde^2 * C.tilde print("mu=") print(mu.tilde) print("V=") print(sigma.tilde) print("determinant=") print(s.tilde^(2 * p)) } else val.m <- val.m + 1 } } list(distance.mahalanobis = MD/s.tilde, determinant = s.tilde^(2 * p),mu = mu.tilde, V = sigma.tilde) } qqbeta_function(x) { x_as.matrix(x) p_ncol(x) n_nrow(x) a_p/2 b_(n-p-1)/2 alpha_0.5*(a-1)/a beta_0.5*(b-1)/b x_sort(mahalanobis(x,apply(x,2,"mean"),var(x))) y_qbeta(((1:n)-alpha)/(n-alpha-beta+1),a,b)*(n-1)^2/n plot(x,y,xlab="Squared distances",ylab="Beta quantiles") } asymp_function(p, r) { cte <- constante(p, r) b0 <- (r * cte^2)/6 pi <- 3.14159 f1 <- function(x, cte, gammap, p) { y <- (2 * pi^(p/2))/gammap * x^(p - 1) * (2 * pi)^( - p/2) * exp( - x^2/2) ind <- ((x) <= cte) x[ind] <- x[ind] - (2 * x[ind]^3)/(cte^2) + x[ind]^5/(cte^4) x[!ind] <- 0 (x^2 * y)/p } f2 <- function(x, cte, gammap, p) { y <- (2 * pi^(p/2))/gammap * x^(p - 1) * (2 * pi)^( - p/2) * exp( - x^2/2) ind <- ((x) <= cte) x1 <- x x2 <- x x1[ind] <- x[ind] - (2 * x[ind]^3)/(cte^2) + x[ind]^5/(cte^4) x2[ind] <- 1 - (6 * x[ind]^2)/(cte^2) + (5 * x[ind]^4)/(cte^4) x1[!ind] <- 0 x2[!ind] <- 0 (x2/p + ((1 - 1/p) * x1)/x) * y } f3 <- function(x, cte, gammap, p) { y <- (2 * pi^(p/2))/gammap * x^(p - 1) * (2 * pi)^( - p/2) * exp( - x^2/2) ind <- ((x) <= cte) x1 <- x x1[ind] <- x[ind] - (2 * x[ind]^3)/(cte^2) + x[ind]^5/(cte^4) x1[!ind] <- 0 (x1^2 * x^2) * y } f4 <- function(x, cte, gammap, p) { y <- (2 * pi^(p/2))/gammap * x^(p - 1) * (2 * pi)^( - p/2) * exp( - x^2/2) ind <- ((x) <= cte) x1 <- x x2 <- x x1[ind] <- x[ind] - (2 * x[ind]^3)/(cte^2) + x[ind]^5/(cte^4) x2[ind] <- 1 - (6 * x[ind]^2)/(cte^2) + (5 * x[ind]^4)/(cte^4) x1[!ind] <- 0 x2[!ind] <- 0 (x2 * x^2 + (p + 1) * x1 * x) * y } f5 <- function(x, b0, cte, gammap, p) { y <- (2 * pi^(p/2))/gammap * x^(p - 1) * (2 * pi)^( - p/2) * exp( - x^2/2) ind <- ((x) <= cte) x[ind] <- x[ind]^2/2 - x[ind]^4/(2 * cte^2) + x[ind]^6/(6 * cte^ 4) x[!ind] <- cte^2/6 (x - b0)^2 * y } f6 <- function(x, cte, gammap, p) { y <- (2 * pi^(p/2))/gammap * x^(p - 1) * (2 * pi)^( - p/2) * exp( - x^2/2) ind <- ((x) <= cte) x1 <- x x1[ind] <- x[ind] - (2 * x[ind]^3)/(cte^2) + x[ind]^5/(cte^4) x1[!ind] <- 0 x1 * x * y } gammap <- gamma(p/2) e1 <- integrate(f1, 0, Inf, cte = cte, gammap = gammap, p = p)$integral e2 <- integrate(f2, 0, Inf, cte = cte, gammap = gammap, p = p)$integral e3 <- integrate(f3, 0, Inf, cte = cte, gammap = gammap, p = p)$integral e4 <- integrate(f4, 0, Inf, cte = cte, gammap = gammap, p = p)$integral e5 <- integrate(f5, 0, Inf, b0 = b0, cte = cte, gammap = gammap, p = p)$ integral e6 <- integrate(f6, 0, Inf, cte = cte, gammap = gammap, p = p)$integral lambda <- e1/e2^2 sigma1 <- (p * (p + 2) * e3)/e4^2 sigma2 <- - (2/p) * sigma1 + ((4 * e5)/e6^2) list(lambda = round(lambda, digits = 3), sigma1 = round(sigma1, digits = 3), sigma2 = round(sigma2, digits = 3)) }