Nothing
###########################################
#### AUTHOR: Arnost Komarek ####
#### (2003) ####
#### ####
#### FILE: convertCDA.R ####
#### ####
#### FUNCTIONS: c2a ####
#### a2c ####
#### derivative.expAD ####
#### find.c ####
#### give.c ####
#### derivative.cc3 ####
###########################################
### ========================================
### c2a: Compute a coefficients from c's
### ========================================
c2a <- function(ccoef, which.zero = which.max(ccoef), toler = 1e-6){
ccoef[ccoef < toler] <- toler
c.zero <- ccoef[which.zero]
acoef <- log(ccoef/c.zero)
return(acoef)
}
### ========================================
### a2c: Function to compute c's from a's
### =========================================
a2c <- function(acoef){
ccoef <- exp(acoef)
sum.exp.a <- sum(ccoef)
if (is.nan(sum.exp.a)) return(NULL)
ccoef <- ccoef/sum.exp.a
return(ccoef)
}
### =======================================================================================
### derivative.expAD: Function to compute derivatives of non-zero exp(a)'s w.r.t. exp(d)'s
### =======================================================================================
##
## * there are g - 1 non-zero a's
## * there are g - 3 d's
##
## * g - 3 a's are equal to d's
## * 2 a's are a function of the rest
## * 1 a is equal to zero
##
## INPUT: knots ....... vector of knots
## sdspline .... standard deviation of a basis G-spline
## last.three ... vector with indeces of the three a's which are to be computed from the rest
## a[last.three[1]] = 0
## a[last.three[2]] = first function(d's)
## a[last.three[3]] = second function(d's)
## all ......... do I want the full matrix or only two columns w.r.t.
## to the two a's
## OUTPUT: Omega.....
## if all == TRUE, matrix (g - 2) x g (there is one zero column)
## all == FALSE, matrix (g - 2) x 2
## (the first row is always an intercept)
derivative.expAD <- function(knots, sdspline, last.three, all = TRUE){
g <- length(knots)
if (g < 4){
stop("Too short input 'knots' vector ")
}
if (length(last.three) != 3){
stop("Incorrect 'last.three' parameter ")
}
if (sum(last.three %in% (1:g)) != 3){
stop("Incorrect 'last.three' parameter ")
}
if (last.three[1] == last.three[2] || last.three[1] == last.three[3] || last.three[2] == last.three[3]){
stop("Incorrect 'last.three' parameter ")
}
if (sdspline >= 1 || sdspline <= 0){
stop("Incorrect 'sdspline' parameter ")
}
which.zero <- last.three[1]
l1 <- last.three[2]
l2 <- last.three[3]
s02 <- sdspline * sdspline
if (abs(knots[l1]) < 1e-4){
stop("Zero reference knot in derivative.expAD ")
}
kn2.kn1 <- knots[l2] - knots[l1]
jsmm <- 1 - s02 + knots[l1]*knots[l2]
## Knots with removed the two ones corresponding to the two special a's
knotsMin2 <- knots[-c(l1, l2)]
## Index of the zero a in the shorter (by 2) a's sequence
which.zero2 <- ifelse(which.zero < min(l1, l2),
which.zero,
ifelse(which.zero < max(l1, l2),
which.zero - 1,
which.zero - 2))
## Compute the two columns of the resulting matrix corresponding to the two a's
vec2b <- -(1/kn2.kn1) * (knotsMin2 - knots[l1])
vec2c <- (1/jsmm) * (1 - s02 + knots[l1] * knotsMin2)
vec2 <- (vec2b & vec2c)
vec1a <- knots[l2]/knots[l1]
vec1d <- (1/knots[l1]) * knotsMin2
vec1 <- -vec1a * vec2 - vec1d
int1 <- vec1[which.zero2]
int2 <- vec2[which.zero2]
slope1 <- vec1[-which.zero2]
slope2 <- vec2[-which.zero2]
slope <- cbind(slope1, slope2)
intercept <- c(int1, int2)
Omega <- rbind(intercept, slope)
if (all){
leftMat <- rbind(matrix(0, nrow = 1, ncol = g - 3), diag(g - 3))
zeroCol <- matrix(c(1, rep(0, g - 3)), nrow = g - 2, ncol = 1)
k <- 1
kk <- 1
for (j in 1:g){
if (j == which.zero){
Omega <- cbind(Omega, zeroCol)
}
else{
if (j == l1 || j == l2){
Omega <- cbind(Omega, Omega[,k])
k <- k + 1
}
else{
Omega <- cbind(Omega, leftMat[,kk])
kk <- kk + 1
}
}
}
## the first two columns are now obscure -> remove them
Omega <- Omega[,-c(1,2)]
}
return(Omega)
}
### ===================================================================
### find.c: Find mixture proportions that approximate
### given distribution (dist) by a G-spline mixture with knots
### and standard deviation sdspline
### ===================================================================
## RETURNS: vector with c coeff.
## or NULL if problems to find them
find.c <- function(knots, sdspline, dist = "dnorm"){
nsplines <- length(knots)
in.knots <- list(x=as.numeric(knots))
right.side <- do.call(dist, in.knots)
right.side[abs(right.side) < 1e-10] <- 0
mus <- matrix(rep(knots, rep(nsplines, nsplines)), nrow = nsplines)
knotsmat <- matrix(rep(knots, nsplines), nrow = nsplines)
Cmat <- dnorm(knotsmat, mean=mus, sd=sdspline)
Cmat <- qr(Cmat, tol = 1e-07)
if (Cmat$rank == ncol(Cmat$qr)){
ccoef <- solve(Cmat, right.side)
tempmin <- min(ccoef[ccoef > 0])
ccoef[ccoef <= 0] <- tempmin ## this is not theoretically possible but numerically it is
}
else
ccoef <- NULL
return(ccoef)
}
### =================================================================
### give.c: Give a vector of all c's satisfying the three constrains
### if remaining (g-3) c's are given.
### =================================================================
## INPUT: knots ......... knots
## sdspline ...... standard deviation of a G-spline
## c.rest ..... remaining g-3 c coefficients
## last.three... indeces of the three c coefficients which are a function of the rest ones
## RETURN: a vector of all c's
give.c <- function(knots, sdspline, last.three, c.rest)
{
g <- length(knots)
if (length(c.rest) != g - 3)
stop("Incorrect dimension of the input vector ")
if (length(last.three) != 3) stop("Incorrect 'last.three' parameter ")
if (sum(last.three %in% 1:g) != 3) stop("Incorrect 'last.three' parameter ")
if (length(unique(last.three)) != 3) stop("Incorrect 'last.three' parameter ")
## Matrix to compute last c's from the first g - 3 ones
Omega <- derivative.cc3(knots, sdspline, last.three, all = TRUE)
## compute all c's
Omega0 <- Omega[1,]
Omega1 <- matrix(Omega[2:(g - 2),], nrow = g - 3)
c.all <- (t(Omega1) %*% c.rest) + Omega0
return(as.numeric(c.all))
}
### ==========================================================
### derivative.cc3: Derivatives of all c's w.r.t. (g - 3) c's
### ==========================================================
## INPUT: knots ....... knots
## sdspline .... standard deviation of a G-spline
## last.three... indeces of the three c coefficients which are a function of the rest ones
## all ....... if TRUE, matrix to compute all c's from the first three ones is returned
## if FALSE, matrix to compute last three c's from the first three ones is returned
## RETURN: Matrix where the first row is an intercept
## and remaining (g-3) rows is dc/da
derivative.cc3 <- function(knots, sdspline, last.three, all = TRUE){
g <- length(knots)
if (g < 4) stop("Too short input 'knots' vector ")
if (length(last.three) != 3) stop("Incorrect 'last.three' parameter ")
if (sum(last.three %in% 1:g) != 3) stop("Incorrect 'last.three' parameter ")
if (length(unique(last.three)) != 3) stop("Incorrect 'last.three' parameter ")
last.three <- last.three[order(last.three)]
s02 = sdspline * sdspline
## Indices of the three c coefficients
l0 <- last.three[1]
l1 <- last.three[2]
l2 <- last.three[3]
## Compute first the matrix used to compute last three c's from the first g - 3 ones
Omega = matrix(0, nrow = g - 2, ncol = 3)
## 1st row (intercept)
Omega[1, 1] = (1 - s02 + knots[l2]*knots[l1])/((knots[l2] - knots[l0])*(knots[l1] -knots[l0]))
Omega[1, 2] = -(1 - s02 + knots[l2]*knots[l0])/((knots[l2] - knots[l1])*(knots[l1] -knots[l0]))
Omega[1, 3] = 1 - Omega[1, 1] - Omega[1, 2]
## the rest (loop over rows)
i <- 2
for (j in 1:g){
if (j == l0 || j == l1 || j == l2) next;
Omega[i, 1] = -((knots[l2] - knots[j])*(knots[l1] - knots[j]))/((knots[l2] - knots[l0])*(knots[l1] -knots[l0]));
Omega[i, 2] = ((knots[l2] - knots[j])*(knots[l0] - knots[j]))/((knots[l2] - knots[l1])*(knots[l1] -knots[l0]));
Omega[i, 3] = -(1 + Omega[i, 1] + Omega[i, 2])
i <- i + 1
}
## Add the components to compute all c's from the first g - 3 ones
if (all){
leftmat <- rbind(rep(0, g - 3), diag(g - 3))
OmegaAll <- NULL
k <- 1; kk <- 1
for (j in 1:g){
if (j == l0 || j == l1 || j == l2){
OmegaAll <- cbind(OmegaAll, Omega[,k])
k <- k + 1
}
else{
OmegaAll <- cbind(OmegaAll, leftmat[, kk])
kk <- kk + 1
}
}
Omega <- OmegaAll
}
return(Omega)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.