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

```# Gamma factors, which multiplied my the GB2 density give the component densities
fg.cgb2 <- function(x, shape1, scale, shape2, shape3, pl0, decomp="r"){
# pl0 is a vector of probabilities which sums to 1
# for ex. pl0 = c(1/3,1/3,1/3)
if (decomp=="r")
{sh <- shape3
a0 <- shape1}
if (decomp=="l")
{sh <- shape2
a0 <- -shape1}
t <- (x/scale)^a0+1
ncomp <- length(pl0)
fac <- matrix(rep(0, times=length(x)*ncomp), ncol=ncomp)
u <- qgamma(cumsum(pl0),sh)
u1 <-c(0,u[-ncomp])     # dim = ncomp
pq <- shape2+shape3
fac[,1] <-pgamma(t*u[1],pq)/pgamma(u[1],sh)
for(i in 2:ncomp){
fac[,i] <-(pgamma(t*u[i],pq) - pgamma(t*u1[i],pq))/(pgamma(u[i],sh) - pgamma(u1[i],sh))
}
dimnames(fac)[[2]] <- paste("fac",1:dim(fac)[2],sep="")
return(fac)
}

dl.cgb2 <- function(x, shape1, scale, shape2, shape3, pl0, decomp="r"){
# pl0 is a vector of probabilities which sums to 1
# for ex. pl0 = c(1/3,1/3,1/3)
L <- length(pl0)
fac <- fg.cgb2(x, shape1, scale, shape2, shape3, pl0, decomp)
dcl <- fac
for (i in 1:L){
dcl[,i] <- dgb2(x,shape1,scale,shape2,shape3)*fac[,i]
}
colnames(dcl) <- paste("comp",1:L,sep="")
return(dcl) #returns a matrix with the l-th component density in the l-th column
}

# Distribution function
library(cubature)

pl.cgb2 <- function(y, shape1, scale, shape2, shape3, pl0, decomp="r", tol=1e-05){
ncomp <- length(pl0)
Fl <- matrix(rep(1, times=length(y)*ncomp), ncol=ncomp)
dimnames(Fl)[[2]] <- paste("comp",1:dim(Fl)[2],sep="")
v <- (y<Inf)
if (length(y[v]) !=0){
for (i in 1:ncomp){
integrand <- function(x) dl.cgb2(x,shape1,scale,shape2,shape3,pl0,decomp)[i]
Fl[v,i] <- apply(as.matrix(y[v]),1,fun)
}
}
vv <- (y==Inf)
Fl[vv,] <- 1
return(Fl)
}

# Density
dcgb2 <- function(x, shape1, scale, shape2, shape3, pl0, pl, decomp="r"){
# x: value at which the density is evaluated, could be a vector or a scalar
# shape1,scale,shape2,shape3: GB2 parameters
# vector pl0 (sums to 1) of length L, the initial proportions
# vector pl gives the fitted proportions of the L components
# if pl=pl0 the function returns the GB2 density
ncomp <- length(pl0)
v <- (x<Inf)
fl <- matrix(rep(0, times=length(x)*ncomp), ncol=ncomp)
fl[v,] <- dl.cgb2(x[v],shape1,scale,shape2,shape3,pl0,decomp)
vv <- (x==Inf)          #new 7.06.12
fl[vv,] <- 0           #new 7.06.12
compd <- fl%*%pl
dimnames(compd)[[2]] <- list("compound")
return(compd)
}

# Cumulative distribution function
pcgb2 <- function(y, shape1, scale, shape2, shape3, pl0, pl, decomp="r"){
n <- length(y)
compp <- matrix(rep(1,n),ncol=1)
v <- (y<Inf)
if (length(y[v])!=0){
integrand <- function(x) dcgb2(x,shape1,scale,shape2,shape3,pl0,pl,decomp)
compp[v,1] <- apply(as.matrix(y[v]),1,fun)
}
vv <- (y==Inf)
compp[vv,] <- 1
return(compp)
}

prcgb2 <- function(y1, y2, shape1, scale, shape2, shape3, pl0, pl, decomp="r", tol=1e-08, debug=FALSE){
# Given  2 arguments y1 and y2, returns the probability that the GB2 r.v. Y is P(min(y1,y2) < Y < max(y1,y2))

integrand <- function(x) dcgb2(x,shape1,scale,shape2,shape3,pl0,pl,decomp)
M1 <- max(y1,y2)
m1 <- min(y1,y2)
if (m1 == M1) return(0)
if (M1 < Inf)  {if (!debug) {
}
else {
}
}
if (M1 == Inf) {if (!debug) {