R/IdtMaxLikNC2.R

Defines functions C2GetCovStderr PartoMatC2 CovtoParC2 VCovGC2LogLikC.grad SigmagradCnt SigmasrgradCnt Sigmagrad C2GetCov initparconf2 GC2LogLikC.grad GC2mLogLik.grad GC2mLogLik Cnf2MaxLik

Documented in C2GetCov C2GetCovStderr Cnf2MaxLik CovtoParC2 GC2LogLikC.grad GC2mLogLik GC2mLogLik.grad initparconf2 PartoMatC2 Sigmagrad SigmagradCnt SigmasrgradCnt VCovGC2LogLikC.grad

Cnf2MaxLik <- function(Data,initpar=NULL,EPS=1E-6,OptCntrl=list())
{
   #  Note  -  The Data argument should be a matrix containing the mid-points in the first columns and the log-ranges in the following columns 

	sd0_default <- 0.01

	n <- nrow(Data)			    # number of observations
	p <- ncol(Data)			    # total number of variables (mid-points + log-ranges)
	q <- p/2			    # number of interval variables
	intcomb <- q*(q+1)/2	   	    # Number of possible combinations between pairs of interval variables
	npar <- 2*intcomb + q	   	    # Total number of parameters to optimize
	if (is.null(initpar))  initpar <- initparconf2(var(Data)*(n-1)/n,n,q)	
	if (!is.null(OptCntrl$sd0)) sd0 <- OptCntrl$sd0 
	else sd0 <- sd0_default 
	OptCntrl$sd0 <- NULL 
#	parsd <- rep(sd0,npar)      	    # standard deviation hyper-parameters - used to generate random starting points
        parsd <- sd0*c(rep(1./sqrt(n),p),rep(0.1,npar-p))   # standard deviation hyper-parameters
	parlb <- rep(-Inf,npar)    	    # vector of lower bounds
	for (i in 1:q)  { 
		parlb[i*(i-1)/2+i] <- EPS 	    # diagonal elements of Choleski decomposition must be positive
		parlb[intcomb+i*(i-1)/2+i] <- EPS
	}
	res <- RepLOptim(initpar,parsd,fr=GC2mLogLik,gr=GC2mLogLik.grad,lower=parlb,X=Data,control=OptCntrl)
	list(lnLik=-res$val,SigmaSr=res$par,optres=res)
}

GC2mLogLik <- function(t0,X,ue=NULL,maxsdratio=10000)  # Minus Log-likelihood for the Gaussian model with Configuration 2
{
   # Wrapper function to be minimized in the ML estimation of a gaussian model assuming configuration 2 
   # (all correlations allowed except for mid-points and log-ranges between different interval variables) 
   #
   # Arguments:
   # 
   #  t0 - vector of parameters to be optimized. It should have as its components the non-null elements (row by row)
   #        of the lower-triangular Choleski decomposition of the covariance matrix among mid-points, followed by 
   #        the non-null elements for the log-ranges, and finally the non-null elements for the pairs mid-points
   #        log-ranges of the same variables 
   #
   #  Other arguments:
   #
   #  X  -  		Data matrix containing the mid-points in the first columns and the log-ranges in the following columns 
   #  ue -  		Matrix (with the same dimension as X) containing observation specific estimates of the X means 
   #	    		If NULL (default) the data is assumed to be centered
   #  maxsdratio -	Maxmimum allowale value for ratios of standard deviations before the covariance matrix is considered
   #                    to be numerically singular, and a (large) penalty is added to the negative log-likelihood 
   #
   #  Value: -1 times the log-likelihood 

	PenF <- 1E12					# large penalty for unfeasible parameters                             
	if ( any(!is.finite(t0)) ) return(PenF)   	# make sure that all parameters are valid real number

	p <- ncol(X)					# total number of variables (mid-points + log-ranges)
	q <- p/2					# number of interval variables

   	Sigmasr <- matrix(0.,nrow=p,ncol=p)
	intcomb <- q*(q+1)/2			# Number of possible combinations between pairs of interval variables
	for (i in 1:q) {
		for (j in 1:i) {		
			Sigmasr[i,j] <- t0[i*(i-1)/2+j]			# loadings among mid-points  
			Sigmasr[q+i,q+j] <- t0[intcomb+i*(i-1)/2+j]	# loadings among log-ranges
		}	 
		Sigmasr[q+i,i] <- t0[2*intcomb+i]	# loadings between mid-points and log-ranges of the same variables
	}
	for (i in 1:(q-1))		# non-free values of the loading matrix, required to ensure null-correlations
		for (j in (i+1):q)  {	# between mid-points and log-ranges of diferent interval variables	
			tempsum <- 0.		
			for (k in i:(j-1))  tempsum <- tempsum + Sigmasr[q+i,k]*Sigmasr[j,k]
			Sigmasr[q+i,j] <- -tempsum/Sigmasr[j,j]
		}
	Sigsd <- diag(Sigmasr)
	sdratio <- max(Sigsd)/min(Sigsd)
	if (sdratio > maxsdratio) return(PenF*(sdratio-maxsdratio))
	if (!is.null(ue)) X <- scale(X,center=ue,scale=FALSE) 

	-sum( apply(X, 1, ILogLikNC1, SigmaSrInv=forwardsolve(Sigmasr,diag(p)), const=-0.5*p*log(2*pi)-sum(log(diag(Sigmasr)))) )
}

GC2mLogLik.grad <- function(t0,X,ue=NULL)
{
    # gradient of GC2mLogLik 

	PenF <- 1E12					# large penalty for unfeasible parameters                             
	if ( any(!is.finite(t0)) ) return(PenF)   	# make sure that all parameters are valid real number

	p <- ncol(X)					# total number of variables (mid-points + log-ranges)
	n <- nrow(X)					# total number of observations
	q <- p/2					# number of interval variables

   	Sigmasr <- matrix(0.,nrow=p,ncol=p)
	intcomb <- q*(q+1)/2			# Number of possible combinations between pairs of interval variables
	difintcomb <- q*(q-1)/2			# Number of possible combinations between pairs of different interval variables
	for (i in 1:q) {
		for (j in 1:i) {		
			Sigmasr[i,j] <- t0[i*(i-1)/2+j]		# loadings among mid-points  
			Sigmasr[q+i,q+j] <- t0[intcomb+i*(i-1)/2+j]	# loadings among log-ranges
		}	 
		Sigmasr[q+i,i] <- t0[2*intcomb+i]	# loadings between mid-points and log-ranges of the same variables
	}
	for (i in 1:(q-1))		# non-free values of the loading matrix, required to ensure null-correlations
		for (j in (i+1):q)  {	# between mid-points and log-ranges of diferent interval variables	
			tempsum <- 0.		
			for (k in i:(j-1))  tempsum <- tempsum + Sigmasr[q+i,k]*Sigmasr[j,k]
			Sigmasr[q+i,j] <- -tempsum/Sigmasr[j,j]
		}
        if (is.null(ue)) {
		if (!is.matrix(X)) X <- as.matrix(X)
		V = t(X) %*% X / n
	}
	else  {
		Xdev = as.matrix(X-ue)
        	V = t(Xdev) %*% Xdev / n
	}
	SigmaInv <- chol2inv(t(Sigmasr))
	fixpgrad = matrix(0.,nrow=2*q,ncol=2*q)
	for (a in 1:(q-1))  for (b in q:(a+1)) {
		fixpgrad[q+a,b] = fixpgrad[q+a,b] + fL.grad(dfx=Sigmagrad,L=Sigmasr,l1=q+a,l2=b,totald=TRUE,L2=SigmaInv,n=n,V=V) 
		if (b > a+1) for (k in (a+1):(b-1)) fixpgrad[q+a,k] = fixpgrad[q+a,k] - Sigmasr[b,k]/Sigmasr[b,b] * fixpgrad[q+a,b]
	}

 	grad = array(dim=2*intcomb+q)
	for (i in 1:q)  { 
		for (j in 1:i) {		
			grad[i*(i-1)/2+j] = fL.grad(dfx=Sigmagrad,L=Sigmasr,l1=i,l2=j,totald=TRUE,L2=SigmaInv,n=n,V=V) 
			if (j<i) for (a in 1:(i-1)) grad[i*(i-1)/2+j] = grad[i*(i-1)/2+j] - Sigmasr[q+a,j]/Sigmasr[i,i] * fixpgrad[q+a,i]
			else if (j>1) for (a in 1:(j-1)) for (k in a:(j-1))
				grad[j*(j-1)/2+j] = grad[j*(j-1)/2+j] + Sigmasr[q+a,k]*Sigmasr[j,k]/Sigmasr[j,j]^2 * fixpgrad[q+a,j]
			grad[intcomb+i*(i-1)/2+j] = fL.grad(dfx=Sigmagrad,L=Sigmasr,l1=q+i,l2=q+j,totald=TRUE,L2=SigmaInv,n=n,V=V) 
		}
		grad[2*intcomb+i] = fL.grad(dfx=Sigmagrad,L=Sigmasr,l1=q+i,l2=i,totald=TRUE,L2=SigmaInv,n=n,V=V)
		if (i<q) for (b in (i+1):q) grad[2*intcomb+i] = grad[2*intcomb+i] - Sigmasr[b,i]/Sigmasr[b,b] * fixpgrad[q+i,b]
	}
	-matrix(grad,1,length(grad))
}

GC2LogLikC.grad <- function(t0,X,ue=NULL)   
# Gradient of observation contributions for Log-likelihood of Gaussian model with Configuration 2
{
	PenF <- 1E12					# large penalty for unfeasible parameters                             
	if ( any(!is.finite(t0)) ) return(PenF)   	# make sure that all parameters are valid real number

	p <- ncol(X)					# total number of variables (mid-points + log-ranges)
	n <- nrow(X)					# total number of observations
	q <- p/2					# number of interval variables

   	Sigmasr <- matrix(0.,nrow=p,ncol=p)
	intcomb <- q*(q+1)/2			# Number of possible combinations between pairs of interval variables
	difintcomb <- q*(q-1)/2			# Number of possible combinations between pairs of different interval variables
	for (i in 1:q) {
		for (j in 1:i) {		
			Sigmasr[i,j] <- t0[i*(i-1)/2+j]		# loadings among mid-points  
			Sigmasr[q+i,q+j] <- t0[intcomb+i*(i-1)/2+j]	# loadings among log-ranges
		}	 
		Sigmasr[q+i,i] <- t0[2*intcomb+i]	# loadings between mid-points and log-ranges of the same variables
	}
	for (i in 1:(q-1))		# non-free values of the loading matrix, required to ensure null-correlations
		for (j in (i+1):q)  {	# between mid-points and log-ranges of diferent interval variables	
			tempsum <- 0.		
			for (k in i:(j-1))  tempsum <- tempsum + Sigmasr[q+i,k]*Sigmasr[j,k]
			Sigmasr[q+i,j] <- -tempsum/Sigmasr[j,j]
		}
        if (is.null(ue)) { if (!is.matrix(X)) X <- as.matrix(X) }
	else  X <- as.matrix(X-ue)
        Vcnt <- apply(X,1,function(v) outer(v,v)/n)
	dim(Vcnt) <- c(p,p,n)
	SigmaInv <- chol2inv(t(Sigmasr))
	fixpgradc = array(0.,dim=c(p,p,n))
	for (a in 1:(q-1))  for (b in q:(a+1)) {
		fixpgradc[q+a,b,] = fixpgradc[q+a,b,] + apply(Vcnt,3,SigmasrgradCnt,L=Sigmasr,l1=q+a,l2=b,L2=SigmaInv,n=n) 
		if (b > a+1) for (k in (a+1):(b-1)) fixpgradc[q+a,k,] = fixpgradc[q+a,k,] - Sigmasr[b,k]/Sigmasr[b,b] * fixpgradc[q+a,b,]
	}
 	gradc = matrix(nrow=2*intcomb+q,ncol=n)
	for (i in 1:q)  { 
		for (j in 1:i) {		
			gradc[i*(i-1)/2+j,] =  apply(Vcnt,3,SigmasrgradCnt,L=Sigmasr,l1=i,l2=j,L2=SigmaInv,n=n) 
			if (j<i) for (a in 1:(i-1)) gradc[i*(i-1)/2+j,] = gradc[i*(i-1)/2+j,] - Sigmasr[q+a,j]/Sigmasr[i,i] * fixpgradc[q+a,i,]
			else if (j>1) for (a in 1:(j-1)) for (k in a:(j-1))
				gradc[j*(j-1)/2+j,] = gradc[j*(j-1)/2+j,] + Sigmasr[q+a,k]*Sigmasr[j,k]/Sigmasr[j,j]^2 * fixpgradc[q+a,j,]
			gradc[intcomb+i*(i-1)/2+j,] =  apply(Vcnt,3,SigmasrgradCnt,L=Sigmasr,l1=q+i,l2=q+j,L2=SigmaInv,n=n) 
		}
		gradc[2*intcomb+i,] = apply(Vcnt,3,SigmasrgradCnt,L=Sigmasr,l1=q+i,l2=i,L2=SigmaInv,n=n) 
		if (i<q) for (b in (i+1):q) gradc[2*intcomb+i,] = gradc[2*intcomb+i,] - Sigmasr[b,i]/Sigmasr[b,b] * fixpgradc[q+i,b,]
	}
	t(gradc)
}

initparconf2 <- function(S,n,q)	
{					 
# 	Find initial parameter estimates in configuration 2 (all correlations allowed except for mid-points and log-ranges 
#	between different interval variables) by replacing values in the Choleski decomposition of the within covariance 
#	matrix in order to inforce the required null-correlations

	EPS <- 1E-9
	intcomb <- q*(q+1)/2			# Number of possible combinations between pairs of interval variables
	initpar <- array(dim=2*intcomb+q)
	Segv <- eigen(S,symmetric=TRUE,only.values=TRUE)$values	
	if (Segv[2*q]<EPS) Sigmasr <- diag(1,2*q)			
	else Sigmasr <- t(chol(S))	
	cnt <- 1	
	for (i in 1:q) {
		for (j in 1:i)  {
			initpar[cnt] <- Sigmasr[i,j] 
			initpar[intcomb+cnt] <- Sigmasr[q+i,q+j]
			cnt <- cnt + 1
		}
		initpar[2*intcomb+i] <- Sigmasr[q+i,i]
	}
	initpar
} 

C2GetCov <- function(par,OStdv,q,tol=1E-8)

#  Gets the covariance matrix from the optimal Configuration 2 optimization applied to scaled data

#  Arguments

#  par    - The optimal parameter values from the (scaled data) Configuration 2 optimization
#  OStdv  - Matrix with the outer product of the standard deviations
#  q      - Number of integer variables
  
{
	L <- matrix(0.,2*q,2*q)
	intcomb <- q*(q+1)/2		# Number of possible combinations between pairs of interval variables
	cnt <- 1	
	for (i in 1:q) {
		for (j in 1:i)  {
			L[i,j] <- par[cnt]  
			L[q+i,q+j] <- par[intcomb+cnt]  
			cnt <- cnt + 1
		}
		L[q+i,i] <- par[2*intcomb+i]  
	}
	for (i in 1:q)  {					# non-free values of the loading matrix, required to 
		if (i>1) for (j in 1:(i-1)) L[q+i,j] <- 0. 	# ensure null-correlations between mid-points and 	
		if (i<q) for (j in (i+1):q)  {			# log-ranges of diferent interval variables	
			tempsum <- 0.		
			for (k in i:(j-1))  tempsum <- tempsum + L[q+i,k]*L[j,k]
			L[q+i,j] <- -tempsum/L[j,j]
		}
	}
	StdSigma <- L %*% t(L)
	StdSigma[abs(StdSigma)<tol]  <- 0.
	OStdv * StdSigma
}

Sigmagrad <- function(j1,j2,n,V,L2)
{
   grad <- n * ( sum(outer(L2[,j1],L2[,j2])*V) - L2[j1,j2] )
   if (j1==j2)  grad <- grad / 2
   grad #  return(grad) 
}

SigmasrgradCnt <- function(VCnt,L2,L,l1,l2,n)  fL.grad(dfx=SigmagradCnt,L=L,l1=l1,l2=l2,totald=TRUE,L2=L2,n=n,VCnt=VCnt)

SigmagradCnt <- function(VCnt,InvMat,j1,j2,n)
{
    grad <- n*sum(outer(InvMat[,j1],InvMat[,j2])*VCnt) - InvMat[j1,j2] 
    if (j1==j2)  grad <- grad / 2
    grad #  return(grad) 
}


#VCovGC2LogLikC.grad <- function(t0,X,ue=NULL)   
VCovGC2LogLikC.grad <- function(t0,X,limlnk2,ue=NULL)   
# Variance-Covariance Gradient of observation contributions for Log-likelihood of Gaussian model with Configuration 2
{
	PenF <- 1E12					# large penalty for unfeasible parameters                             
	if ( any(!is.finite(t0)) ) return(PenF)   	# make sure that all parameters are valid real number

	p <- ncol(X)					# total number of variables (mid-points + log-ranges)
	n <- nrow(X)					# total number of observations
	q <- p/2					# number of interval variables

   	Sigma <- matrix(0.,nrow=p,ncol=p)
	intcomb <- q*(q+1)/2			# Number of possible combinations between pairs of interval variables
	difintcomb <- q*(q-1)/2			# Number of possible combinations between pairs of different interval variables
	for (i in 1:q) {
		for (j in 1:i) {		
			Sigma[i,j] <- Sigma[j,i] <- t0[i*(i-1)/2+j]			# covariances among mid-points  
			Sigma[q+i,q+j] <- Sigma[q+j,q+i] <- t0[intcomb+i*(i-1)/2+j]	# covariances among log-ranges
		}	 
		Sigma[q+i,i] <- Sigma[i,q+i] <- t0[2*intcomb+i]		# covariances between mid-points and log-ranges of the same variables
	}
        if (is.null(ue)) { if (!is.matrix(X)) X <- as.matrix(X) }
	else  X <- as.matrix(X-ue)
        Vcnt <- apply(X,1,function(v) outer(v,v)/n)
	dim(Vcnt) <- c(p,p,n)
#	SigmaInv <- pdwt.solve(Sigma)
	SigmaInv <- Safepdsolve(Sigma,maxlnk2=limlnk2,scale=TRUE)
 	gradc = matrix(nrow=2*intcomb+q,ncol=n)
	for (i in 1:q)  { 
		for (j in 1:i) {		
			gradc[i*(i-1)/2+j,] =  apply(Vcnt,3,SigmagradCnt,j1=i,j2=j,InvMat=SigmaInv,n=n) 
			gradc[intcomb+i*(i-1)/2+j,] =  apply(Vcnt,3,SigmagradCnt,j1=q+i,j2=q+j,InvMat=SigmaInv,n=n) 
		}
		gradc[2*intcomb+i,] = apply(Vcnt,3,SigmagradCnt,j1=q+i,j2=i,InvMat=SigmaInv,n=n) 
	}
	t(gradc)
}

CovtoParC2 <- function(Sigma,q,Sqrt)	
{					 
# 	Gets the vectorized form of covariance matrix according to Configuration 2

	intcomb <- q*(q+1)/2	   	    # Number of possible combinations between pairs of interval variables
	par <- array(dim=2*intcomb+q)
	if (Sqrt==TRUE) mat <- t(chol(Sigma))
	else mat <- Sigma
	cnt <- 1	
	for (i in 1:q) {
		for (j in 1:i)  {
			par[cnt] <- mat[i,j] 
			par[intcomb+cnt] <- mat[q+i,q+j]
			cnt <- cnt + 1
		}
		par[2*intcomb+i] <- mat[q+i,i]
	}
	par
} 

PartoMatC2 <- function(par,q,initval)

#  Converts a vectorized form of parameters (or standard errors) according to Configuration 2, to ist matrix form
{
	mat <- matrix(initval,2*q,2*q)
	intcomb <- q*(q+1)/2		# Number of possible combinations between pairs of interval variables
	cnt <- 1	
	for (i in 1:q) {
		for (j in 1:i)  {
			mat[i,j] <- mat[j,i] <- par[cnt]  
			mat[q+i,q+j] <- mat[q+j,q+i] <- par[intcomb+cnt]  
			cnt <- cnt + 1
		}
		mat[q+i,i] <- mat[i,q+i] <- par[2*intcomb+i]  
	}
	mat
}


#C2GetCovStderr <- function(mleSigmaC2,Data,q,ue=NULL)
C2GetCovStderr <- function(mleSigmaC2,Data,q,limlnk2,ue=NULL)

#  Arguments

#  mleSigmaC2    - The maximum likelihood estimators of the covariance matrix under configuration C2
#  Data      	 - Data frame with the MidPoints followed by the Log-Ranges of the Interval Data set.
#                  When ue is set to NULL (default) Data is assumed to be centered.
#  q             - Number of Integer Varaibles
#  ue		 - Vector with the MidPoints and Log-Ranges means (or NULL if Data is centered)
  
{
	if (!is.null(ue)) Data <- scale(Data,center=ue,scale=FALSE)
	par <- CovtoParC2(mleSigmaC2,q,Sqrt=FALSE)
#	gradc <- VCovGC2LogLikC.grad(par,Data)
	gradc <- VCovGC2LogLikC.grad(par,Data,limlnk2=limlnk2)
	VCov <- mleVCov(gradc)
	stderr <- sqrt(diag(VCov))
	PartoMatC2(stderr,q,NA)
}

Try the MAINT.Data package in your browser

Any scripts or data that you put into this service are public.

MAINT.Data documentation built on Sept. 21, 2021, 5:11 p.m.