R/utils_SMN.R

In PartCensReg: Estimation and Diagnostics for Partially Linear Censored Regression Models Based on Heavy-Tailed Distributions

#######################################################################################################################
## Square Matrix ##
#######################################################################################################################

matrix.sqrt <- function(A){
  sva <- svd(A)
  if (min(sva$d)>=0)
    Asqrt <- t(sva$v %*% (t(sva$u) * sqrt(sva$d)))
  else
    stop("Matrix square root is not defined")
  return(Asqrt)
}

#######################################################################################################################
## AIC for the Scale Mixtures of Normal distributions (SMN) ##
#######################################################################################################################


AICsmn<-function(lambda,logver,ff,K,N,sigma2,p,type){
  q<-ncol(K)
  B<-solve(matrix.sqrt(t(N)%*%N))
  L<-sigma2*B%*%K%*%B
  df<-sum(diag(solve(diag(q)+lambda*L)))
  lp<- logver - 0.5*lambda*(t(ff)%*%K%*%ff)
  if(type=="Normal"){AICm<- -2*lp+2*(p+1+df)}
  if(type=="T"){AICm<- -2*lp+2*(p+1+1+df)}
  if(type=="Slash"){AICm<- -2*lp+2*(p+1+1+df)}
  if(type=="NormalC"){AICm<- -2*lp+2*(p+1+2+df)}
  return(-AICm)
}

#######################################################################################################################
## BIC for the Scale Mixtures of Normal distributions (SMN) ##
#######################################################################################################################


BICsmn<-function(lambda,logver,ff,K,N,sigma2,n, p,type){
  q<-ncol(K)
  B<-solve(matrix.sqrt(t(N)%*%N))
  L<-sigma2*B%*%K%*%B
  df<-sum(diag(solve(diag(q)+lambda*L)))
  lp<- logver - 0.5*lambda*(t(ff)%*%K%*%ff)
  if(type=="Normal"){BICm<- -2*lp+log(n)*(p+1+df)}
  if(type=="T"){BICm<- -2*lp+log(n)*(p+1+1+df)}
  if(type=="Slash"){BICm<- -2*lp+log(n)*(p+1+1+df)}
  if(type=="NormalC"){BICm<- -2*lp+log(n)*(p+1+2+df)}
  return(-BICm)
}


#######################################################################################################################
## Densities of the SMN with location scale and  Censure ##
## Densities of the Normal(densN), Student's-t(densT), Slash(densSL) and Contaminated  Normal(densCN) with censure  ##
#######################################################################################################################


dNormal <- function(cc, y, mu, sigma2 = 1, cens){
  if(cens=="1"){
    densN        <- vector(mode = "numeric", length = length(y))
    densN[cc==0] <- dnorm(y[cc==0], mu[cc==0], sqrt(sigma2))
    densN[cc==1] <- pnorm((y[cc==1] - mu[cc==1])/sqrt(sigma2))
    if(length(which(densN == 0)) > 0) densN[which(densN == 0)] <- .Machine$double.xmin
  }
  if(cens=="2"){
    densN        <- vector(mode = "numeric", length = length(y))
    densN[cc==0] <- dnorm(y[cc==0], mu[cc==0], sqrt(sigma2))
    densN[cc==1] <- 1-pnorm((y[cc==1] - mu[cc==1])/sqrt(sigma2))
    if(length(which(densN == 0)) > 0) densN[which(densN == 0)] <- .Machine$double.xmin
  }
  return(densN)
}


dT <- function(cc, y, mu, sigma2 = 1, nu=4, cens){
  if(cens=="1"){
    densT        <- vector(mode = "numeric", length = length(y))
    aux          <- (y-mu)/sqrt(sigma2)
    densT[cc==0] <- dt(aux[cc==0],df=nu)/sqrt(sigma2)
    densT[cc==1] <- pt(aux[cc==1],df=nu)
    if(length(which(densT == 0)) > 0) densT[which(densT == 0)] <- .Machine$double.xmin
  }
  if(cens=="2"){
    densT        <- vector(mode = "numeric", length = length(y))
    aux          <- (y-mu)/sqrt(sigma2)
    densT[cc==0] <- dt(aux[cc==0],df=nu)/sqrt(sigma2)
    densT[cc==1] <- 1-pt(aux[cc==1],df=nu)
    if(length(which(densT == 0)) > 0) densT[which(densT == 0)] <- .Machine$double.xmin
  }
  return(densT)
}



dSL <- function(cc, y, mu, sigma2 = 1, nu=2, cens){
  
  if(sum(cc)==0){
  
  if(cens=="1"){
    densSL        <- vector(mode = "numeric", length = length(y))
    aux           <- (y-mu)/sqrt(sigma2)
    densSL[cc==0] <- dSlash(aux[cc==0],0,1,nu)/sqrt(sigma2)
    if(length(which(densSL == 0)) > 0) densSL[which(densSL == 0)] <- .Machine$double.xmin
  }
  if(cens=="2"){
    densSL        <- vector(mode = "numeric", length = length(y))
    aux           <- (y-mu)/sqrt(sigma2)
    densSL[cc==0] <- dSlash(aux[cc==0],0,1,nu)/sqrt(sigma2)
    if(length(which(densSL == 0)) > 0) densSL[which(densSL == 0)] <- .Machine$double.xmin
  }
  
  }
  else {
    if(cens=="1"){
    densSL        <- vector(mode = "numeric", length = length(y))
    aux           <- (y-mu)/sqrt(sigma2)
    densSL[cc==0] <- dSlash(aux[cc==0],0,1,nu)/sqrt(sigma2)
    densSL[cc==1] <- AcumSlash(aux[cc==1],0,1,nu)
    if(length(which(densSL == 0)) > 0) densSL[which(densSL == 0)] <- .Machine$double.xmin
      }
  if(cens=="2"){
      densSL        <- vector(mode = "numeric", length = length(y))
      aux           <- (y-mu)/sqrt(sigma2)
      densSL[cc==0] <- dSlash(aux[cc==0],0,1,nu)/sqrt(sigma2)
      densSL[cc==1] <- 1-AcumSlash(aux[cc==1],0,1,nu)
      if(length(which(densSL == 0)) > 0) densSL[which(densSL == 0)] <- .Machine$double.xmin
     }
  } 
   
  return(densSL)
}



dCN <- function(cc, y, mu, sigma2 = 1, nu=c(0.1,0.1), cens){
  if(cens=="1"){
    densCN        <- vector(mode = "numeric", length = length(y))
    aux           <-(y-mu)/sqrt(sigma2)
    densCN[cc==0] <- dNormalC(aux[cc==0],0,1,nu)/sqrt(sigma2)
    densCN[cc==1] <- AcumNormalC(aux[cc==1],0,1,nu)
    if(length(which(densCN == 0)) > 0) densCN[which(densCN == 0)] <- .Machine$double.xmin
  }    
  if(cens=="2"){
    densCN        <- vector(mode = "numeric", length = length(y))
    aux           <-(y-mu)/sqrt(sigma2)
    densCN[cc==0] <- dNormalC(aux[cc==0],0,1,nu)/sqrt(sigma2)
    densCN[cc==1] <- 1-AcumNormalC(aux[cc==1],0,1,nu)
    if(length(which(densCN == 0)) > 0) densCN[which(densCN == 0)] <- .Machine$double.xmin
  }  
  return(densCN)
}


#######################################################################################################################
## Probability density function of the SMN ##
## Cumulative distribution function of the SMN ##
#######################################################################################################################

AcumPearsonVII <- function(y,mu,sigma2,nu,delta){
  Acum         <- vector(mode = "numeric", length = length(y))
  z            <- vector(mode = "numeric", length = length(y))
  sigma2a      <- sigma2*(delta/nu)
  for (i in 1:length(y)){
    z[i]       <- (y[i]-mu)/sqrt(sigma2a)
    Acum[i]    <- pt(z[i],df=nu)
  }
  return(Acum)
}


AcumSlash   <- function(y,mu,sigma2,nu){
  Acum      <- z <- vector(mode = "numeric", length = length(y))
  for (i in 1:length(y)){
    z[i]    <- (y[i]-mu)/sqrt(sigma2)
    f1      <- function(u) nu*u^(nu-1)*pnorm(z[i]*sqrt(u))
    Acum[i] <- integrate(f1,0,1)$value
  }
  return(Acum)
}


AcumNormalC   <- function(y,mu,sigma2,nu){
  Acum        <- vector(mode = "numeric", length = length(y))
  for (i in 1:length(y)){
    eta       <- nu[1]
    gama      <- nu[2]
    Acum[i]   <- eta*pnorm(y[i],mu,sqrt(sigma2/gama)) + (1-eta)*pnorm(y[i],mu,sqrt(sigma2))
  }
  return(Acum)
}



dPearsonVII <- function(y,mu,sigma2,nu,delta){
  f <- z <- vector(mode = "numeric", length = length(y))
  sigma2a <- sigma2*(delta/nu)
  for (i in 1:length(y)){
    z[i] <- (y[i]-mu)/sqrt(sigma2a)
    f[i] <- dt(z[i],df=nu)/sqrt(sigma2a)
  }
  return(f)
}


dSlash <- function(y,mu,sigma2,nu){
  resp <- z <- vector(mode = "numeric", length = length(y))
  for (i in 1:length(y)){
    z[i]    <- (y[i]-mu)/sqrt(sigma2)
    f2      <- function(u) nu*u^(nu-0.5)*dnorm(z[i]*sqrt(u))/sqrt(sigma2)
    resp[i] <- integrate(f2,0,1)$value
  }
  return(resp)
}


dNormalC    <- function(y,mu,sigma2,nu){
  Acum      <- vector(mode = "numeric", length = length(y))
  for (i in 1:length(y)){
    eta     <- nu[1]
    gama    <- nu[2]
    Acum[i] <- eta*dnorm(y[i],mu,sqrt(sigma2/gama)) + (1-eta)*dnorm(y[i],mu,sqrt(sigma2))
  }
  return(Acum)
}

#######################################################################################################################
## Expectations ##
#######################################################################################################################

NCensurEsperUY <-  function(y,mu,sigma2,nu,delta,type){
  EUY0   <- EUY1 <- EUY2 <- c()
  d      <- (y-mu)^2/sigma2
  n      <- length(y)
  if(type=="T"){
    EUY0 <- (nu+1)/(nu+d)
    EUY1 <- y*(nu+1)/(nu+d)
    EUY2 <- (y^2)*(nu+1)/(nu+d)
  }
  if(type=="Normal"){
    EUY0 <- rep(1,n)
    EUY1 <- y
    EUY2 <- (y^2)
  }
  if(type=="PearsonVII"){
    d    <- (y-mu)^2/sigma2
    EUY0 <- (nu+1)/(delta+d)
    EUY1 <- y*(nu+1)/(delta+d)
    EUY2 <- (y^2)*(nu+1)/(delta+d)
  }
  if(type=="Slash"){
    Num  <- GamaInc(nu+1.5,0.5*d)*(0.5*d)^(-nu-1.5)
    Den  <- GamaInc(nu+0.5,0.5*d)*(0.5*d)^(-nu-0.5)
    EUY0 <- Num/Den
    EUY1 <- y*Num/Den
    EUY2 <- (y^2)*Num/Den
  }
  
  if(type=="NormalC"){
    Num  <- 1-nu[1]+nu[1]*(nu[2])^(1.5)*exp(0.5*d*(1-nu[2]))
    Den  <- 1-nu[1]+nu[1]*(nu[2])^(0.5)*exp(0.5*d*(1-nu[2]))
    EUY0 <- Num/Den
    EUY1 <- y*Num/Den
    EUY2 <- (y^2)*Num/Den
  }
  return(list(EUY0=EUY0,EUY1=EUY1,EUY2=EUY2))
}

TN <- function(mu,sigma2,a,b){
  sigma2 <- as.vector(sigma2)
  lim1 <- (a-mu)/sqrt(sigma2)
  lim2 <- (b-mu)/sqrt(sigma2)
  EX <-   (dnorm(lim2)- dnorm(lim1))/(pnorm(lim2)- pnorm(lim1))
  EX2 <- 1-(lim2*dnorm(lim2)- lim1*dnorm(lim1))/(pnorm(lim2)- pnorm(lim1))
  EX3 <- 2*EX - (lim2^2*dnorm(lim2)- lim1^2*dnorm(lim1))/(pnorm(lim2)- pnorm(lim1))
  EX4 <- 3*EX2 - (lim2^3*dnorm(lim2)- lim1^3*dnorm(lim1))/(pnorm(lim2)- pnorm(lim1))
  EY <- mu + sqrt(sigma2)*EX
  EY2 <- mu^2 + 2*mu*sqrt(sigma2)*EX + sigma2*(EX2)
  EY3 <- mu^3 + 3*mu^2*sqrt(sigma2)*EX + 3*mu*sigma2*(EX2) + EX3
  EY4 <- mu^4 + 4*mu^3*sqrt(sigma2)*EX + 6*mu^2*sigma2*(EX2) + 4*mu*sigma2^1.5*EX3 + EX4*sigma2^2
  return(list(EY=EY, EY2=EY2, EY3=EY3, EY4=EY4))
}

GamaInc <- function(a,x){
  res <- vector(mode = "numeric", length = length(x))
  f <-function(t) exp(-t)*t^(a-1)
  for  (i in 1:length(x)){
    res[i] <- integrate(f,0,x[i])$value
  }
  return(res)
}


E_phi <- function(r,a,nu,delta,type=type,cens=cens){
  n <- length(a)
  b <- rep(Inf,n)
  b1<- rep(-Inf,n)
  if(setequal(a,b)== TRUE | setequal(a,b1)== TRUE){
    resp <- rep(0,n)
  }
  else{
    if(type=="Normal"){
      resp <- dnorm(a)
    }
    if(type=="T"){
      Aux0 <- gamma(0.5*(nu+2*r))
      Aux1 <- gamma(nu/2)*sqrt(2*pi)
      Aux2 <- Aux0/Aux1
      Aux3 <- (0.5*nu)^(nu/2)
      Aux4 <- (0.5*(a^2+nu))^(-0.5*(nu+2*r))
      resp <- Aux2*Aux3*Aux4
    }
    if(type=="PearsonVII"){
      Aux0 <- gamma(0.5*(nu+2*r))
      Aux1 <- gamma(nu/2)*sqrt(2*pi)
      Aux2 <- Aux0/Aux1
      Aux3 <- (0.5*delta)^(nu/2)
      Aux4 <- (0.5*(a^2+delta))^(-0.5*(nu+2*r))
      resp <- Aux2*Aux3*Aux4
    }
    if(type=="Slash"){
      Aux0 <- nu/sqrt(2*pi)
      Aux1 <- (0.5*a^2)^(-(nu+r))
      Aux2 <- GamaInc(nu+r,0.5*a^2)
      resp <- Aux0*Aux1*Aux2
    }
    if(type=="NormalC"){
      Aux0 <- nu[1]*nu[2]^(r)*dnorm(a*sqrt(nu[2]))
      Aux1 <- (1-nu[1])*dnorm(a)
      resp <- Aux0 + Aux1
    }
  }
  return(resp)
}

E_Phi <- function(r,a,nu,delta,type=type){
  n <- length(a)
  if(type=="Normal"){
    resp <- pnorm(a)
  }
  if(type=="T"){
    Aux0 <- gamma(0.5*(nu+(2*r)))
    Aux1 <- gamma(nu/2)
    Aux2 <- Aux0/Aux1
    Aux3 <- (0.5*nu)^(-r)
    Aux4 <- AcumPearsonVII(a,0,1,nu+(2*r),nu)
    resp <- Aux2*Aux3*Aux4
  }
  if(type=="PearsonVII"){
    Aux0 <- gamma(0.5*(nu+(2*r)))
    Aux1 <- gamma(nu/2)
    Aux2 <- Aux0/Aux1
    Aux3 <- (0.5*delta)^(-r)
    Aux4 <- AcumPearsonVII(a,0,1,nu+(2*r),delta)
    resp <- Aux2*Aux3*Aux4
  }
  if(type=="Slash"){
    Aux0 <- nu/(nu+r)
    Aux1 <- AcumSlash(a,0,1,nu+r)
    resp <- Aux0*Aux1
  }
  if(type=="NormalC"){
    Aux0 <- nu[2]^(r)*AcumNormalC(a,0,1,nu)
    Aux1 <- (1-nu[2]^(r))*(1-nu[1])*pnorm(a)
    resp <- Aux0 + Aux1
  }
  return(resp)
}

CensEsperUY1 <- function(mu,sigma2,nu,delta,Lim1,Lim2,type=type,cens=cens){
  Lim11      <- (Lim1-mu)/sqrt(sigma2)
  Lim21      <- (Lim2-mu)/sqrt(sigma2)
  n          <- length(Lim11)
  if(type=="Normal"){
    EU     <-  1
    FNIb   <- pnorm(Lim21)
    FNIa   <- pnorm(Lim11)
  }
  if(type=="T"){
    EU     <-  1
    FNIb   <- pt(Lim21,nu)
    FNIa   <- pt(Lim11,nu)
  }
  if(type=="PearsonVII"){
    FNIb   <- AcumPearsonVII(Lim21,0,1,nu,delta)
    FNIa   <- AcumPearsonVII(Lim11,0,1,nu,delta)
  }
  if(type=="Slash"){
    FNIb   <- AcumSlash(Lim21,0,1,nu)
    FNIa   <- AcumSlash(Lim11,0,1,nu)
  }
  if(type=="NormalC")
  {
    EU     <- (nu[1]*nu[2]) + (1-nu[1])
    FNIb   <- AcumNormalC(Lim21,0,1,nu)
    FNIa   <- AcumNormalC(Lim11,0,1,nu)
  }
  if (cens=="1"){
    Aux11  <- rep(0,n)
  }else{
    Aux11  <- Lim11
  }
  if (cens=="2"){
    Aux22  <- rep(0,n)
  }else{
    Aux22  <- Lim21
  }
  Kaux     <- (FNIb-FNIa)
  if(length(which(Kaux == 0)) > 0) {Kaux[which(Kaux == 0)] <- .Machine$double.xmin}
  K     <- 1/Kaux
  EUX0  <- K*(E_Phi(1,Lim21, nu,delta,type)       - E_Phi(1,Lim11, nu,delta,type))
  EUX1  <- K*(E_phi(0.5,Lim11,nu,delta,type,cens) - E_phi(0.5,Lim21,nu,delta,type,cens))
  EUX2  <- K*(E_Phi(0,Lim21, nu,delta,type)       - E_Phi(0,Lim11, nu,delta,type) + Aux11*E_phi(0.5,Lim11,nu,delta,type,cens) - Aux22*E_phi(0.5,Lim21,nu,delta,type,cens))
  EUX20 <- K*(E_Phi(2,Lim21, nu,delta,type)       - E_Phi(2,Lim11, nu,delta,type))
  EUX21 <- K*(E_phi(1.5,Lim11,nu,delta,type,cens) - E_phi(1.5,Lim21,nu,delta,type,cens))
  EUY0  <- EUX0
  EUY1  <- mu*EUX0 + sqrt(sigma2)*EUX1
  EUY2  <- EUX0*mu^2 + 2*mu*sqrt(sigma2)*EUX1 + sigma2*EUX2
  EUY20 <- EUX20
  EUY21 <- mu*EUX20 + sqrt(sigma2)*EUX21
  return(list(EUY0=EUY0,EUY1=EUY1,EUY2=EUY2, EUY20=EUY20, EUY21=EUY21))
}

rmix <- function(n, p1, rF1, arg1, arg2){
  x1 <- vector(mode = "numeric", length = n)
  for (i in 1:n){
    u <- runif(1)
    if (u < p1) x1[i] <- do.call("rF1", c(list(1), arg1))
    if (u > p1) x1[i] <- do.call("rF1", c(list(1), arg2))
  }
  return(x1)
}

#######################################################################################################################
## Generate of the SMN ##
#######################################################################################################################

geraPearsonVII <- function(n,mu,sigma2,nu,delta){
  y <- vector(mode = "numeric", length =n)
  for (i in 1:n){
    y[i] <- mu[i] + (rgamma(1, 0.5*nu, 0.5*delta))^(-1/2)*rnorm(1,0,sqrt(sigma2))
  }
  return(y)
}

geraS <- function(n,mu, sigma2,nu){
  y <- vector(mode = "numeric", length =n)
  u2 <- rbeta(n,nu,1) 
  for (i in 1:n){
    y[i] <- mu[i] + (u2[i])^(-1/2)*rnorm(1,0, sqrt(sigma2))
  }
  return(y)
}

geraNC <- function(n,mu,sigma2,nu){
  y <- vector(mode = "numeric", length =n)
  p1<- nu[1]
  gama<-nu[2]
  u<-runif(1)
  for(i in 1:n){
    if(u<p1) y[i]<- rnorm(1,mu[i],sqrt(sigma2/gama))
    if(u>p1) y[i]<- rnorm(1,mu[i],sqrt(sigma2))
  }
  return(y)
}