# R/CPLogLik.grads.R In MAINT.Data: Model and Analyse Interval Data

```library(sn)

#Note:  All these functions assume that vectors with non-diferent elements of symmetrical
#       matrices use a lower-triangular format

ltdind <- function(p) {
if (p==1) return(1)
i <- (p-1):1
c(1,(p*(p+1)-i*(i+1))/2+1)
}
utdind <- function(p) {i <- 1:p ;  i*(i+1)/2}

fmdind <- function(p) (1:p)^2

dvecind <- function(p)
{
ind <- 1:p
if (p>1) for (r in 2:p) ind <- c(ind,(r-1)*p+(r:p))
ind
}

dvec0ind <- function(p)
{
ind <- 2:p
if (p>2) for (r in 2:(p-1)) ind <- c(ind,(r-1)*p+((r+1):p))
ind
}

lftD <- function(vM,p)
{
tmpM <- matrix(nrow=p,ncol=p)
tmpM[row(tmpM)>=col(tmpM)] <- vM
tmpM[row(tmpM)<col(tmpM)] <- t(tmpM)[row(tmpM)<col(tmpM)]
matrix(tmpM,p^2,1)
}

lftDplus <- function(vecM,p,indices=dvecind(p),nrows=p*(p+1)/2)
matrix(vecM[indices],nrows,1)

rghtD <- function(vecM,p,ncols=p*(p+1)/2)
{
tmpM <- matrix(vecM,p,p)
tmpM[row(tmpM)>col(tmpM)] <-
tmpM[row(tmpM)>col(tmpM)] + t(tmpM)[row(tmpM)>col(tmpM)]
matrix(tmpM[row(tmpM)>=col(tmpM)],1,ncols)
}

{
Omega <- matrix(nrow=p,ncol=p)
Omega[lower.tri(Omega,diag=TRUE)] <- x
Omega[upper.tri(Omega)] <- t(Omega)[upper.tri(Omega)]
OmegaI <- Safepdsolve(Omega,maxlnk2=limlnk2,scale=TRUE)
if(is.null(OmegaI)) {
warning("Singular Omega matix found in the course of the mle estimation\n")
return(matrix(0.,nrow=nvcovpar,ncol=nvcovpar))
}

nvcovpar <- p*(p+1)/2
Jacob <- matrix(nrow=nvcovpar,ncol=nvcovpar)
ind <- 0
for (r in 1:p) for (c in r:p) {
ind <- ind+1
tmpM <- -outer(OmegaI[,r],OmegaI[,c])
if (r!=c) tmpM <- tmpM -outer(OmegaI[,c],OmegaI[,r])
Jacob[,ind] <- tmpM[lower.tri(tmpM,diag=TRUE)]
}
Jacob
}

{
logDet <- attr(OmegaInv,"log.det")
if (is.matrix(y)) {
y0 <- scale(y,center=ksi,scale=FALSE)
nOmgminusS0 <- n*Omega - t(y0)%*%y0
}
else {
y0 <- y-ksi
nOmgminusS0 <- Omega - outer(y0,y0)
}
diag(nOmgminusS0) <- diag(nOmgminusS0)/2
Omegapar <- Omega[lower.tri(Omega,diag=TRUE)]
if (is.matrix(y0)) {
sum(sapply(y0%*%eta,zeta1))*eta )
}
else  {
ksigrad <- drop( OmegaInv%*%y0 - zeta(1,y0%*%eta)*eta )
}
}

p=ifelse(is.matrix(y),ncol(y),length(y)),limlnk2)
{
nvcovpar <- p*(p+1)/2
ksi <- param[1:p]
Omega <- matrix(nrow=p,ncol=p)
Omega[lower.tri(Omega,diag=TRUE)] <- param[(p+1):(p+nvcovpar)]
Omega[upper.tri(Omega)] <- t(Omega)[upper.tri(Omega)]
omega <- sqrt(diag(Omega))
alpha <- param[p+nvcovpar+1:p]
eta <- alpha/omega
OmegaI <- Safepdsolve(Omega,maxlnk2=limlnk2,scale=TRUE)
if(is.null(OmegaI)) {
warning("Singular Omega matix found in the course of the mle estimation\n")
return(rep(0.,2*p+nvcovpar))
} else {
attr(OmegaI,"log.det") <- determinant(Omega,logarithm=TRUE)\$modulus
}
diagOmegaind <- ltdind(p)

}

b <- sqrt(2./pi)
b0 <- 2/(4-pi)
zeta1 <- function(x) zeta(1,x)

p=ifelse(is.matrix(y),ncol(y),length(y)),inoptm=TRUE,
PenF=1E12,ldRtol=log(1E-5)+1-p,c2tol=1E-6,limlnk2)
{
if ( (is.matrix(y) && p!=ncol(y)) || (!is.matrix(y) && p!=length(y)) )
{
stop("Dimension of y is not compatible with argument p\n")
}
if (!all(is.finite(param))) { return(rep(0.,length(param))) }
nvcovpar <- p*(p+1)/2

mu <- param[1:p]
Sigma <- matrix(nrow=p,ncol=p)
Sigma[lower.tri(Sigma,diag=TRUE)] <- param[(p+1):(p+nvcovpar)]
Sigma[upper.tri(Sigma)]  <- t(Sigma)[upper.tri(Sigma)]
gamma1 <- param[(p+nvcovpar+1):(2*p+nvcovpar)]
intres\$D23,intres\$D32,intres\$D33,intres\$Dtld32,intres\$Dtld33,
y,n,p,PenF,ldRtol,c2tol,limlnk2=limlnk2)
}

ldRtol=-500,c2tol=1e-6,beta0tol=1e-6,limlnk2)
{
nvcovpar <- p*(p+1)/2
sigma <- sqrt(diag(Sigma))
mu0 <- sign(gamma1)*sigma*(b0*abs(gamma1))^(1/3)
SigmaI <- Safepdsolve(Sigma,maxlnk2=limlnk2,scale=TRUE)
if (is.null(SigmaI))  {
if (inoptm)
{
warning("Non-positive definite covariance matrix found during gradient computations (which returned 0).\n")
return(rep(0.,2*p+nvcovpar))
}  else {
return(NULL)
}
} else {
attr(SigmaI,"log.det") <- determinant(Sigma,logarithm=TRUE)\$modulus
}
lRdet <- attr(SigmaI,"log.det") - sum(log(diag(Sigma)))
if (lRdet < ldRtol)
{
tmpM <- 2*SigmaI
diag(tmpM) <- diag(tmpM)/2 - 1./diag(Sigma)
} else {
}

SigmaImu0 <- drop(SigmaI%*%mu0)
beta02 <- drop(mu0%*%SigmaImu0)
beta0 <- sqrt(beta02)
p2 <- p^2

Ip <- diag(p)
D23 <- apply(kronecker(Ip,mu0)+kronecker(mu0,Ip),2,
lftDplus,p=p,nrows=nvcovpar)

if (beta0 < beta0tol)  {
Dtld32 <- matrix(0.,p,p*(p+1)/2)
Dtld33 <- matrix(0.,p,p)
} else {
mu0bar <- mu0/(sigma*beta0)
tmpsum <- matrix(0.,p,p*(p+1)/2)
for (i in 1:p)  {
tmp <- rep(0.,p2)
tmp[(i-1)*p+i] <- mu0bar[i]
tmpsum[i,] <- tmpsum[i,] + rghtD(tmp,p=p)/sigma[i]
}
Dtld32 <- beta0*tmpsum/2
Dtld33 <- (b0/(3*beta02))*sigma*diag(1./mu0bar^2)
}

ksi <- mu - mu0
Omega <- Sigma + outer(mu0,mu0)
OmegaInv <- Safepdsolve(Omega,maxlnk2=limlnk2,scale=TRUE)
if (is.null(OmegaInv)) {
if (inoptm)
{
warning("Non-positive definite matrix Omega found during gradient computations (which returned 0).\n")
return(rep(0.,2*p+nvcovpar))
}  else {
return(NULL)
}
} else {
attr(OmegaInv,"log.det") <- determinant(Omega,logarithm=TRUE)\$modulus
}
b2 <- b^2
omega <- sqrt(diag(Omega))
delta <- mu0/(b*omega)
sclMat <- 1./outer(omega,omega)
OmgbI <- OmegaInv / sclMat
OmgbIdelta <- drop(OmgbI%*%delta)
c2 <- 1. - delta %*% OmgbIdelta

if (c2 < c2tol) {

tmpM <- 2*outer(OmgbIdelta,OmgbIdelta)
diag(tmpM) <- diag(tmpM)/2
vtmpM <- tmpM[lower.tri(tmpM,diag=TRUE)]
vdOmgb <- vtmpM %*% dOmgbdOmg

dAomgm1 <- (1/3)*(sigma/(b*omega)) * (2/(4-pi))^(1/3) / (gamma1^2)^(1/3)
tmp <- -sigma/(b*omega^2) * sign(gamma1)*(2*abs(gamma1)/(4-pi))^(1/3)
dOmgdgam1 <- D23%*%Dtld33
for (j in 1:p)

if ( c2 < 0. || isTRUE(all.equal(c2,0.)) ) return(penaltygrad)
}

c1 <- sqrt((b2-(1.-b2)*beta02)/(1.+beta02))
q1 <- 1./(c1*(1.+beta02))
q2 <- q1*(2*c1-q1)/2

tmpM <- q1*q2*SigmaImu0%*%t(SigmaImu0)
D32 <- -matrix( apply(kronecker(t(SigmaImu0),q1*SigmaI-tmpM),1,
rghtD,p=p,ncols=nvcovpar),
p,nvcovpar,byrow=TRUE )
D33 <- q1*SigmaI-2*tmpM
eta <- q1*SigmaImu0
list(ksi=ksi,Omega=Omega,eta=eta,OmegaInv=OmegaInv,
D23=D23,D32=D32,D33=D33,Dtld32=Dtld32,Dtld33=Dtld33)
}

D23,D32,D33,Dtld32,Dtld33,
y,n=ifelse(is.matrix(y),nrow(y),1),
p=ifelse(is.matrix(y),ncol(y),length(y)),inoptm=FALSE,
PenF=0.,ldRtol=log(1E-5)+1-p,c2tol=1E-6,limlnk2)
{

}

{
# Note: Only three non-negative elemenst per row !!!
#      Later try to find compact representations that multiply
#      these matrices more efficiently

omega <- sqrt(diag(Omg))
if (is.null(sclMat)) sclMat <- 1./outer(omega,omega)
if (is.null(Omgb)) Omgb <- sclMat * Omg

Jacobdim <- Omgdim*(Omgdim+1)/2
Jacob <- matrix(0.,Jacobdim,Jacobdim)
diagind <- ltdind(Omgdim)
Jind <- 0
for (c in 1:(Omgdim-1)) {
Jind <- Jind + 1
for (r in (c+1):Omgdim) {
Jind <- Jind + 1
Jacob[Jind,Jind] <- sclMat[r,c]
Jacob[Jind,diagind[r]] <- -Omgb[r,c]/(2*Omg[r,r])
Jacob[Jind,diagind[c]] <- -Omgb[r,c]/(2*Omg[c,c])
}
}
Jacob
}
```

## 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.