# R/CompoundAuxFit.R In GB2: Generalized Beta Distribution of the Second Kind: Properties, Likelihood, Estimation

```pkl.cavgb2 <- function(z,lambda){
L <- dim(lambda)[2]+1
pkl <- exp(z%*%lambda)
ck <- apply(pkl,1,sum)+1
pkl <- pkl/ck
pkl <- cbind(pkl,1-rowSums(pkl))
colnames(pkl) <- paste("comp", 1:L, sep="")
return(pkl)
# pkl: matrix I x L
}

lambda0.cavgb2 <- function(pl0, z, w=rep(1, dim(z)[1])){
# vl0 - initial values of vl (do not depend on k!!!)
vl0 <- vofp.cgb2(pl0)
sw <- sum(w)
I <- dim(z)[2]
wz <- (w*z)
swz <- apply(wz, 2, sum)
# means - the vector of mean values of the aux vars
# sum(w_k*z_{ik})/sum(w_k)
means <- swz/sw
lambda0 <- (1/(I*means))%*%t(vl0)
return(lambda0)
}

logl.cavgb2 <- function(fac, z, lambda, w=rep(1, dim(fac)[1])){
sw <- sum(w)
pkL <- pkl.cavgb2(z,lambda)
mixt <- rowSums(pkL*fac)
logcompaux <- log(mixt)
logL <- sum(w*logcompaux)/sw
return(logL)
}

scores.cavgb2 <- function(fac, z, lambda, w=rep(1, dim(fac)[1])){
# fac    : (N x L) - matrix of Gamma factors
# z      : (N x I) - matrix of auxiliary variables
# w      : (N x 1) - vector of weights
# lambda : (I x L-1)- matrix of parameters
sw <- sum(w)
evl <- exp(z%*%lambda)
ck <- apply(evl,1,sum)+1
pkl <- evl/ck #dimension N x L-1
pkL <- cbind(pkl,1-rowSums(pkl)) #dimension N x L
L <-dim(pkL)[2]
denom <- rowSums(pkL*fac)
num <- fac[,-L]
midt <- num/as.vector(denom) - 1
dlogL <- t(pkl*midt)%*%(z*w)/sw
return(dlogL)
}

ml.cavgb2 <- function (fac, z, lambda0, w = rep(1, dim(fac)[1]), maxiter = 100, fnscale=length(w))
{
dl <- dim(lambda0)
lambda0vec <- as.vector(lambda0)
fn <- function(lambdavec, fac, z, w) {
lambda <- matrix(lambdavec, nrow = dl[1], ncol = dl[2])
return(-logl.cavgb2(fac, z, lambda, w))
}
gr <- function(lambdavec, fac, z, w) {
lambda <- matrix(lambdavec, nrow = dl[1], ncol = dl[2])
sc <- -scores.cavgb2(fac, z, lambda, w)
sctr <- t(sc)
asvecsc <- as.vector(sctr)
return(asvecsc)
}
opt <- optim(lambda0vec, fn, gr, fac, z, w, method = "BFGS",
control = list(maxit = maxiter, fnscale = fnscale),
hessian = FALSE)
lambdafit <- opt\$par
lambdafitm <- matrix(lambdafit, nrow = dl[1], ncol = dl[2])
colnames(lambdafitm) <- paste("comp", 1:dl[2], sep = "")
rownames(lambdafitm) <- colnames(z)
return(list(lambdafitm, opt))
}
```

## Try the GB2 package in your browser

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

GB2 documentation built on May 2, 2019, 5:53 a.m.