library("mclust") library("adehabitat") ############################################ ## Classification matrix given the groups ## ############################################ clas <- function(groups,k){ n <- length(groups) z <- array(0,c(n,k),dimnames=list(1:n,paste("comp",1:k,sep=""))) for(i in 1:n) z[i,groups[i]] <- 1 return(z) } ########################################### ## Generate random pxp covariance matrix ## ########################################### covariance.matrix <- function(p){ if(p>1) { sigma <- matrix(rnorm(p*p),ncol=p) sigma <- crossprod(sigma)+ diag(rep(0.1, p)) } else sigma <- matrix(rchisq(1,1),ncol=p) return(sigma) } ############################### ## polinomial shape modifier ## ############################### qu <- function(r,beta,d){ qua <- 1 + beta/(8*d*(d+2))*(r^2 - 2*(d+2)*r + d*(d+2)) return(qua) } ############################# ## Weighted Sample Moments ## ############################# WSM <- function(X,weights=NULL,eps=0.001){ X <- as.matrix(X) d <- ncol(X) n <- nrow(X) if(is.null(weights)) weights = rep(1,n) betamax <- min(4*d,4/5*d*(d+2)) mu <- apply(X,2,weighted.mean,w=weights) Sigma <- cov.wt(x=X, wt = weights, cor = FALSE, center = TRUE, method = "ML")$cov beta <- max(eps,min(sum(weights/sum(weights)*mahalanobis(X, center = mu, cov = Sigma, inverted = FALSE)^2)-d*(d+2),(betamax-eps))) return( list( d = d, betamax = betamax, mu = mu, Sigma = Sigma, beta = beta ) ) } ################################################### ## Multivariate Lepto-Kurtic Normal Distribution ## ################################################### MD <- function(x,mu,Sigma,beta){ d <- length(mu) A <- 1+beta/(8*d*(d+2))*(mahalanobis(x=x, center=mu, cov=Sigma, inverted=FALSE)^2 - 2*(d+2)*mahalanobis(x=x, center=mu, cov=Sigma, inverted=FALSE) + d*(d+2)) B <- (2*pi)^(-d/2)*det(Sigma)^(-1/2)*exp(-1/2*mahalanobis(x=x, center=mu, cov=Sigma, inverted=FALSE)) A*B } ##################################################################### ## M-step of the EM algorithm to fit Mixtures of LKN distributions ## ##################################################################### Mstep <- function(X,weights=NULL,method="BFGS",startmu=NULL,startSigma=NULL,startbeta=NULL,maxit=1000,reltol=1e-15,trace=0,eps=0.0001){ X <- as.matrix(X) d <- ncol(X) n <- nrow(X) betamax <- min(4*d,4/5*d*(d+2)) if(is.null(weights)) weights <- rep(1,n) nj <- sum(weights) if(is.null(startmu)) startmu <- apply(X,2,weighted.mean,w=weights) if(is.null(startSigma)) startSigma <- cov.wt(x=X, wt = weights, cor = FALSE, center = TRUE, method = "ML")$cov #(n-1)/n*cov(X) if(is.null(startbeta)) startbeta <- max(eps,min(sum(weights/sum(weights)*mahalanobis(X, center = startmu, cov = startSigma, inverted = FALSE)^2)-d*(d+2),(betamax-eps))) # initialization initial.values <- numeric(d+d*(d+1)/2+1) initial.values[1:d] <- startmu startU <- base::chol(startSigma, pivot=FALSE) cont <- 0 for(i in 1:d) for(j in i:d){ cont <- cont + 1 initial.values[d+cont] <- startU[i,j] } initial.values[d+d*(d+1)/2+1] <- log(startbeta/(betamax-startbeta)) # objective function f <- function(par,X,weights,d,betamax,eps){ # -- # # mu # # -- # mu <- array(0,c(d),dimnames=list(paste("X.",1:d,sep=""))) for(i in 1:d) mu[i] <- par[i] # ----- # # Sigma # # ----- # U <- array(0,c(d,d),dimnames=list(paste("X.",1:d,sep=""),paste("X.",1:d,sep=""))) cont <- 0 for(i in 1:d) for(j in i:d){ cont <- cont + 1 U[i,j] <- par[d+cont] } Sigma <- crossprod(x=U,y=U) # ---- # # beta # # ---- # beta <- (betamax-eps)*exp(par[d+d*(d+1)/2+1])/(1+exp(par[d+d*(d+1)/2+1])) # -------------- # # -loglikelihood # # -------------- # sum(weights*log((1+beta/(8*d*(d+2))*(mahalanobis(x=X,center=mu,cov=Sigma,inverted=FALSE)^2-2*(d+2)*mahalanobis(x=X,center=mu,cov=Sigma,inverted=FALSE)+d*(d+2)))*((2*pi)^(-d/2)*det(Sigma)^(-1/2)*exp(-1/2*mahalanobis(x=X,center=mu,cov=Sigma,inverted=FALSE))))) } grr <- function(par,X,weights,d,betamax,eps){ n <- length(weights) nj <- sum(weights) # -- # # mu # # -- # mu <- array(0,c(d),dimnames=list(paste("X.",1:d,sep=""))) for(i in 1:d) mu[i] <- par[i] # ----- # # Sigma # # ----- # U <- array(0,c(d,d),dimnames=list(paste("X.",1:d,sep=""),paste("X.",1:d,sep=""))) cont <- 0 for(i in 1:d) for(j in i:d){ cont <- cont + 1 U[i,j] <- par[d+cont] } Sigma <- crossprod(x=U,y=U) # ---- # # beta # # ---- # beta <- (betamax-eps)*exp(par[d+d*(d+1)/2+1])/(1+exp(par[d+d*(d+1)/2+1])) # ----------- # # Derivatives # # ----------- # invSigma <- solve(Sigma) invU <- solve(t(U)) v <- beta/(2*d*(d+2))*(mahalanobis(x = X, center = mu, cov = Sigma) - (d+2))/qu(r=mahalanobis(x = X, center = mu, cov = Sigma), beta=beta, d=d) dmu <- colSums(weights*(1-v)*(X-matrix(rep(mu,n),nrow=n,byrow=T)) %*% invSigma) dLambda <- -nj*invU + invU %*% crossprod((weights*(1-v))^(1/2)*(X-matrix(rep(mu,n),nrow=n,byrow=T))) %*% invSigma dLa <- numeric(d*(d+1)/2) cont <- 0 for(i in 1:d) for(j in i:d){ cont <- cont + 1 dLambda[i,j] -> dLa[cont] } dga <- 1/(1+beta/(betamax-beta))*weights %*% (1-1/(qu(r=mahalanobis(x = X, center = mu, cov = Sigma), beta=beta, d=d))) return(c(dmu,dLa,dga)) } res <- optim(par=initial.values, fn=f, gr=grr, X=X, weights=weights, d=d, betamax=betamax, eps=eps, method=method, control=list(maxit=maxit,reltol=reltol,trace=trace,fnscale=-1)) loglik <- res$value est <- res$par # -- # # mu # # -- # estmu <- array(0,c(d),dimnames=list(paste("X.",1:d,sep=""))) for(i in 1:d) estmu[i] <- est[i] # ----- # # Sigma # # ----- # estU <- array(0,c(d,d),dimnames=list(paste("X.",1:d,sep=""),paste("X.",1:d,sep=""))) cont <- 0 for(i in 1:d) for(j in i:d){ cont <- cont + 1 estU[i,j] <- est[d+cont] } estSigma <- crossprod(x=estU,y=estU) # ---- # # beta # # ---- # estbeta <- (betamax-eps)*exp(est[d+d*(d+1)/2+1])/(1+exp(est[d+d*(d+1)/2+1])) return( list( res = res, loglik = loglik, mu = estmu, U = estU, Sigma = estSigma, beta = estbeta ) ) } Mstepwithout <- function(X,weights=NULL,method="BFGS",startmu=NULL,startSigma=NULL,startbeta=NULL,maxit=1000,reltol=1e-15,trace=0,eps=0.0001){ X <- as.matrix(X) d <- ncol(X) n <- nrow(X) betamax <- min(4*d,4/5*d*(d+2)) if(is.null(weights)) weights <- rep(1,n) nj <- sum(weights) if(is.null(startmu)) startmu <- apply(X,2,weighted.mean,w=weights) if(is.null(startSigma)) startSigma <- cov.wt(x=X, wt = weights, cor = FALSE, center = TRUE, method = "ML")$cov #(n-1)/n*cov(X) if(is.null(startbeta)) startbeta <- max(eps,min(sum(weights/sum(weights)*mahalanobis(X, center = startmu, cov = startSigma, inverted = FALSE)^2)-d*(d+2),(betamax-eps))) # initialization initial.values <- numeric(d+d*(d+1)/2+1) #for(i in 1:d) # initial.values[i] <- startmu[i] initial.values[1:d] <- startmu startU <- base::chol(startSigma, pivot=FALSE) cont <- 0 for(i in 1:d) for(j in i:d){ cont <- cont + 1 initial.values[d+cont] <- startU[i,j] } initial.values[d+d*(d+1)/2+1] <- log(startbeta/(betamax-startbeta)) # objective function f <- function(par,X,weights,d,betamax,eps){ # -- # # mu # # -- # mu <- array(0,c(d),dimnames=list(paste("X.",1:d,sep=""))) for(i in 1:d) mu[i] <- par[i] # ----- # # Sigma # # ----- # U <- array(0,c(d,d),dimnames=list(paste("X.",1:d,sep=""),paste("X.",1:d,sep=""))) cont <- 0 for(i in 1:d) for(j in i:d){ cont <- cont + 1 U[i,j] <- par[d+cont] } Sigma <- crossprod(x=U,y=U) # ---- # # beta # # ---- # beta <- (betamax-eps)*exp(par[d+d*(d+1)/2+1])/(1+exp(par[d+d*(d+1)/2+1])) # -------------- # # -loglikelihood # # -------------- # -sum(weights*log((1+beta/(8*d*(d+2))*(mahalanobis(x=X,center=mu,cov=Sigma,inverted=FALSE)^2-2*(d+2)*mahalanobis(x=X,center=mu,cov=Sigma,inverted=FALSE)+d*(d+2)))*((2*pi)^(-d/2)*det(Sigma)^(-1/2)*exp(-1/2*mahalanobis(x=X,center=mu,cov=Sigma,inverted=FALSE))))) } res <- optim(par=initial.values, fn=f, X=X, weights=weights, d=d, betamax=betamax, eps=eps, method=method, control=list(maxit=maxit,reltol=reltol,trace=trace)) loglik <- -res$value est <- res$par # -- # # mu # # -- # estmu <- array(0,c(d),dimnames=list(paste("X.",1:d,sep=""))) for(i in 1:d) estmu[i] <- est[i] # ----- # # Sigma # # ----- # estU <- array(0,c(d,d),dimnames=list(paste("X.",1:d,sep=""),paste("X.",1:d,sep=""))) cont <- 0 for(i in 1:d) for(j in i:d){ cont <- cont + 1 estU[i,j] <- est[d+cont] } estSigma <- crossprod(x=estU,y=estU) # ---- # # beta # # ---- # estbeta <- (betamax-eps)*exp(est[d+d*(d+1)/2+1])/(1+exp(est[d+d*(d+1)/2+1])) return( list( loglik = loglik, mu = estmu, U = estU, Sigma = estSigma, beta = estbeta ) ) } ############################################################################## ## EM algorithm (second version of the paper) ## ## for Mixtures of LKN distributions (with partial derivatives for optim()) ## ############################################################################## EMLKN <- function( X, # matrix of data k, # number of groups initialization="kmeans", # initialization procedure: "random.soft", "random.hard", "manual", or "mclust" eps=0.0001, # the maximum kurtosis is betamax - eps method="BFGS", # used for optim in each M-step gr=TRUE, start.z=NULL, # (n x k)-matrix of soft or hard classification: it is used only if initialization="manual" iter.max=200, # maximum number of iterations in the EM-algorithm reltol=1e-15, trace=0, loglikplot=FALSE#, # if TRUE, the log-likelihood values against the iterations are plotted ){ #call = match.call() if (is.vector(X)) X <- matrix(X,ncol=1) if (is.data.frame(X)) X <- as.matrix(X) n <- nrow(X) # sample size d <- ncol(X) # number of variables if (is.null(X)) stop('Hey, we need some data, please! X is null') #if (!is.matrix(X)) stop('X needs to be in matrix form') if (!is.numeric(X)) stop('X is required to be numeric') if (n == 1) stop('nrow(X) is equal to 1') if (any(is.na(X))) stop('No NAs allowed.') if (is.null(k)) stop('k is NULL') k <- as.integer(ceiling(k)) if (!is.integer(k)) stop('k is not a integer') if (any(k < 1)) stop('k is not a positive integer') # parameters prior <- numeric(k) # proportion of each group mu <- array(0,c(d,k),dimnames=list(paste("X.",1:d,sep=""),paste("group ",1:k,sep=""))) U <- array(0,c(d,d,k),dimnames=list(paste("X.",1:d,sep=""),paste("X.",1:d,sep=""),paste("group ",1:k,sep=""))) Sigma <- array(0,c(d,d,k),dimnames=list(paste("X.",1:d,sep=""),paste("X.",1:d,sep=""),paste("group ",1:k,sep=""))) beta <- array(0,c(k),dimnames=list(paste("group ",1:k,sep=""))) kurtosi <- d*(d+2) betamax <- min(4*d,4/5*d*(d+2)) PX <- array(0,c(n,k),dimnames=list(1:n,paste("group ",1:k,sep=""))) # ------------------------- # # posteriors initialization # # ------------------------- # if(initialization=="random.soft"){ z <- array(runif(n*k),c(n,k)) # soft posterior probabilities (no-normalized) (n x k) z <- z/rowSums(z) # soft posterior probabilities (n x k) } if(initialization=="random.hard"){ z <- t(rmultinom(n, size = 1, prob=rep(1/k,k))) # hard posterior probabilities (n x k) } if(initialization=="manual"){ # z.start can be both soft and hard initialization z <- start.z # soft or hard classification } if(initialization=="mclust"){ if(d==1) mclustfit <- Mclust(data=X, G=k, modelNames="V") if(d>1) mclustfit <- Mclust(data=X, G=k, modelNames="VVV") z <- mclustfit$z } if(initialization=="kmeans"){ if(k>1){ km <- kmeans(X,k) groups <- km$cluster z <- clas(groups=groups,k=k) } if(k==1){ groups <- rep(1,n) z <- matrix(groups,nrow=n,ncol=k) } } startmu <- startSigma <- startbeta <- NULL # ------------ # # EM algorithm # # ------------ # # Preliminary definition of convergence criterions check <- 0 iteration <- 1 loglik <- NULL aloglik <- NULL aloglik <- c(0,0) a <- NULL a <- c(0,0) if(gr==TRUE) cat("EM algorithm (M-step with exact partial derivatives): ") if(gr==FALSE) cat("EM algorithm (M-step with numeric partial derivatives): ") while(check<1){ # ++++++ # # M-step # # ++++++ # # ------- # # Weights # # ------- # prior <- colMeans(z) # ---------------- # # other parameters # # ---------------- # for(j in 1:k){ if(iteration > 1){ startmu <- mu[,j] startSigma <- Sigma[,,j] startbeta <- beta[j] } if(gr) res <- Mstep(X=X,weights=z[,j],method=method,startmu=startmu,startSigma=startSigma,startbeta=startbeta,maxit=iter.max,reltol=reltol,trace=trace) if(!gr) res <- Mstepwithout(X=X,weights=z[,j],method=method,startmu=startmu,startSigma=startSigma,startbeta=startbeta,maxit=iter.max,reltol=reltol,trace=trace) mu[,j] <- res$mu U[,,j] <- res$U Sigma[,,j] <- res$Sigma beta[j] <- res$beta } # ------- # # density # # ------- # for(j in 1:k) PX[,j] <- MD(x=X,mu=mu[,j],Sigma=Sigma[,,j],beta=beta[j]) # ------------------------------------- # # Global - Observed-data log-likelihood # # ------------------------------------- # # model-based clustering llvalues <- sum(log(rowSums(matrix(rep(prior,n),n,k,byrow=TRUE)*PX))) loglik[iteration] <- llvalues # --------------------------- # # Aitken's Stopping Criterion # # --------------------------- # if(iteration>2 & k > 1){ if(abs(loglik[iteration-1]-loglik[iteration-2])>0){ a[iteration-1] <- (loglik[iteration]-loglik[iteration-1])/(loglik[iteration-1]-loglik[iteration-2]) aloglik[iteration] <- loglik[iteration-1]+(1/(1-a[iteration-1])*(loglik[iteration]-loglik[iteration-1])) if(abs(aloglik[iteration]-loglik[iteration])10^(-322))*z # vincolo per evitare i NaN nel calcolo di tau*log(tau) hard.z <- (matrix(rep(apply(z,1,max),k),n,k,byrow=F)==z)*1 ECM <- sum(hard.z*log(z.const)) ICL <- BIC+ECM result <- list( X = X, k = k, d = d, n = n, prior = prior, mu = mu, U = U, Sigma = Sigma, beta = beta, betarel = beta/betamax, iter.stop = iteration, z = z, group = group, loglik = finalloglik, AIC = AIC, AIC3 = AIC3, AICc = AICc, AICu = AICu, AWE = AWE, CAIC = CAIC, BIC = BIC, ICL = ICL, # alla McNicholas call = match.call() ) class(result) <- "MM" return(result) alarm() } ################################################ ## Select k for Mixtures of MLN distributions ## ################################################ selectk <- function( X, maxk = 1, # maximum number of components to be tried maxR = 1, # number of initial random points to be tried initialization = "mclust", eps = 0.0001, # minimum value for the excess of kurtosis in WSM method = "BFGS", # used for optim in each M-step start.z = NULL, # (n x k)-matrix of soft or hard classification: it is used only if initialization="manual" iter.max = 200, # maximum number of iterations in the EM-algorithm reltol = 1e-15, trace = 0, loglikplot = FALSE, # if TRUE, the log-likelihood values against the iterations are plotted verbose = TRUE ){ gridk <- 1:maxk gridR <- 1:maxR numk <- length(gridk) numR <- length(gridR) IC <- array(NA,c(numk,numR,3),dimnames=list( paste("# of groups=",gridk,sep=""), paste("# of initial random starts=",gridR,sep=""), c("loglik","BIC","ICL")) ) par <- list() par1 <- list() par2 <- list() for(i in gridk){ par[[i]] <- list() par1[[i]] <- list() par2[[i]] <- list() for(h in gridR){ if(verbose){ cat("\n") cat(paste("Model with k=",gridk[i]," at iteration ",gridR[h],sep="")) cat("\n") } # ---------------------- # # with exact derivatives # # ---------------------- # tryCatch({ par1[[i]][[h]] <- EMLKN( X=X, # matrix of data k=gridk[i], # number of groups initialization=initialization, # initialization procedure: "random.soft", "random.hard", "kmeans", or "mclust" gr=TRUE, eps=eps, # minimum value for the excess of kurtosis in WSM method=method, # used for optim in each M-step start.z=start.z, # (n x k)-matrix of soft or hard classification: it is used only if initialization="manual" iter.max=iter.max, # maximum number of iterations in the EM-algorithm reltol=reltol, trace=trace, loglikplot=loglikplot # if TRUE, the log-likelihood values against the iterations are plotted ) } ,error = function(e) cat(paste("\n Model not estimated.\n",e)), par1[[i]][[h]] <- 0 ) # ------------------------ # # with numeric derivatives # # ------------------------ # tryCatch({ par2[[i]][[h]] <- EMLKN( X=X, # matrix of data k=gridk[i], # number of groups initialization=initialization, # initialization procedure: "random.soft", "random.hard", "kmeans", or "mclust" gr=FALSE, eps=eps, # minimum value for the excess of kurtosis in WSM method=method, # used for optim in each M-step start.z=start.z, # (n x k)-matrix of soft or hard classification: it is used only if initialization="manual" iter.max=iter.max, # maximum number of iterations in the EM-algorithm reltol=reltol, trace=trace, loglikplot=loglikplot # if TRUE, the log-likelihood values against the iterations are plotted ) } ,error = function(e) cat(paste("\n Model not estimated.\n",e)), par2[[i]][[h]] <- 0 ) # ------------- # # best solution # # ------------- # if(is.list(par1[[i]][[h]]) & is.list(par2[[i]][[h]])){ if(par1[[i]][[h]]$loglik >= par2[[i]][[h]]$loglik) par[[i]][[h]] <- par1[[i]][[h]] if(par1[[i]][[h]]$loglik < par2[[i]][[h]]$loglik) par[[i]][[h]] <- par2[[i]][[h]] } if(!is.list(par1[[i]][[h]])) par[[i]][[h]] <- par2[[i]][[h]] if(!is.list(par2[[i]][[h]])) par[[i]][[h]] <- par1[[i]][[h]] IC[i,h,1] <- par[[i]][[h]]$loglik IC[i,h,2] <- par[[i]][[h]]$BIC IC[i,h,3] <- par[[i]][[h]]$ICL } } # --- # # BIC # # --- # BestIndBIC <- which.max(IC[,,2]) BestIndBIC <- arrayInd(BestIndBIC, .dim=c(numk,numR)) BestBIC <- IC[BestIndBIC[1],BestIndBIC[2],2] bestBIC <- par[[BestIndBIC[1]]][[BestIndBIC[2]]] # --- # # ICL # # --- # BestIndICL <- which.max(IC[,,3]) BestIndICL <- arrayInd(BestIndICL, .dim=c(numk,numR)) BestICL <- IC[BestIndICL[1],BestIndICL[2],3] bestICL <- par[[BestIndICL[1]]][[BestIndICL[2]]] bestk <- array(0,c(2),dimnames=list(c("BIC","ICL"))) bestk <- c(gridk[BestIndBIC[1]],gridk[BestIndICL[1]]) if(verbose){ cat("\n\n") cat("# ----------------------- #","\n") cat("# Model Selection Results #","\n") cat("# ----------------------- #","\n\n") cat("Best BIC value of",BestBIC,"obtained for k =",gridk[BestIndBIC[1]],"group(s)","\n\n") cat("Best ICL value of",BestICL,"obtained for k =",gridk[BestIndICL[1]],"group(s)","\n\n") } return( structure( list( call = match.call(), X = X, fittedmodels = par, bestBICpar = par[[gridk[BestIndBIC[1]]]][[gridR[BestIndBIC[2]]]], bestICLpar = par[[gridk[BestIndICL[1]]]][[gridR[BestIndICL[2]]]], best = bestk, BIC = IC[,,2], ICL = IC[,,3], bestBIC = bestBIC, bestICL = bestICL ), class = "kbest" ) ) } ################################### ## Univariate LeptoKurtic-Normal ## ################################### univlknorm <- function(x,mu=0,sd=1,beta=0){ z <- (x-mu)/sd A <- 1+beta/24*(z^4-6*z^2+3) B <- dnorm(x,mean=mu,sd=sd) return(A*B) } ############################################################## ## Distribution function of a Univariate LeptoKurtic-Normal ## ############################################################## punivlknorm <- function(q,mu=0,sd=1,beta=0){ ploum <- function(y) univlknorm(y,mu,sd,beta) res <- sapply(q, function(i) integrate(ploum, -Inf, i)$value) return(res) } ########################################################## ## Quantile function of a Univariate LeptoKurtic-Normal ## ########################################################## qunivlknorm <- function(p,mu=0,sd=1,beta=0){ n <- length(p) res <- numeric(n) for(i in 1:n){ f <- function(par,p) (punivlknorm(q=par,mu=mu,sd=sd,beta=beta)-p)^2 res[i] <- optim(par=qnorm(p=p[i],mean=mu,sd=sd),fn=f,p=p[i],method="BFGS")$par } res } ############################################################ ## Random Generation from a Univariate LeptoKurtic-Normal ## ############################################################ runivlknorm <- function(n,mu=0,sd=1,beta=0){ values <- runif(n) sapply(1:n, function(i) qunivlknorm(p=values[i],mu=mu,sd=sd,beta=beta)) } ################################## ## Bivariate LeptoKurtic-Normal ## ################################## bivlknorm <- function(x,y,mu,Sigma,beta){ rho <- Sigma[1,2] muX <- mu[1] muY <- mu[2] sigmaX <- sqrt(Sigma[1,1]) sigmaY <- sqrt(Sigma[2,2]) mah <- 1/(1-rho^2)*( (x-muX)^2/sigmaX^2 - 2*rho*(x-muX)*(y-muY)/(sigmaX*sigmaY) + (y-muY)^2/sigmaY^2 ) A <- 1+beta/64*(mah^2-8*mah+8) B <- 1/(2*pi*sigmaX*sigmaY*sqrt(1-rho^2))*exp(-1/2*mah) A*B } library(mvtnorm) library(pracma) ############################################################################################ ## RANDOM GENERATOR ## ############################################################################################ ##################################### ## Density function of the R tilde ## ##################################### dRtilde <- function(x,d,beta=0) { sapply(1:length(x), function(i) 2*x[i]*dchisq(x[i]^2,df=d)*(1+beta/(8*d*(d+2))*(x[i]^4-2*(d+2)*x[i]^2+d*(d+2)))) } ########################################## ## Distribution function of the R tilde ## ########################################## pRtilde <- function(q,d,beta=0) { ploum <- function(y) dRtilde(y,d,beta) res <- sapply(q, function(i) integrate(ploum, 0, i)$value) return(res) } ####################################### ## Quantile function of a Modified R ## ####################################### qRtilde <- function (p,d,beta=0){ m <- length(p) res <- numeric(m) f <- function(par,p) (pRtilde(q=par,d=d,beta=beta)-p)^2 res <- sapply(1:m, function(i) optimize(f = f, p = p[i], interval = c(0,10))$minimum ) return(res) } ######################################### ## Random Generation from a Modified R ## ######################################### rRtilde <- function(m,d,beta=0){ values <- runif(m,min=0,max=1) res <- sapply(1:m, function(i) qRtilde(p=values[i], d=d, beta=beta)) } ################################################################### ## Random Generator from a Multivariate Leptokurtic Normal (MLN) ## ################################################################### rMLN <- function(m, d, mu=rep(0,d), Sigma=diag(d), beta=0){ X <- matrix(,m,d) for(i in 1:m){ A <- t(chol(Sigma)) R <- rRtilde(1,d=d,beta=beta) G <- t(rmvnorm(1,mean=rep(0,d))) U <- G/Norm(G) X[i,] <- mu + A %*% U * R } return(X) }