gcjc <- function(formula,
data,
category,
par,
config = 1,
zlimit = Inf,
fixed = list(),
equal.noise = TRUE,
opt = c("nlminb", "optim"),
lower=-Inf,
upper=Inf,
control=list())
{
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset", "na.action"), names(mf), 0L)
mf <- mf[c(1L, m)]
mf$drop.unused.levels <- TRUE
mf[[1L]] <- as.name("model.frame")
mf <- eval.parent(mf)
Terms <- attr(mf, "terms")
response <- model.response(mf)
x <- model.matrix(Terms, mf)
xint <- match("(Intercept)", colnames(x), nomatch=0L)
if(xint > 0) x <- x[, -xint, drop=FALSE]
dimen <- ncol(x)
if(missing(category)) category <- response
if(nrow(data) != length(category))
stop("nrow(data) and length(category) are different")
#biases <- as.vector(rowMeans(x))
if(missing(par)) par <- list(noise=NULL,bias=NULL)
if(!("gcjcStruct" %in% class(par)))
{
if(is.null(par$noise)) par$noise <- 10
if(is.null(par$bias)){
means <- sapply(split(as.data.frame(x), category), colMeans)
par$bias <- as.vector(rowMeans(means))
}
par <- gcjcStruct(noise=par$noise,bias=par$bias, config)
names(par) <- attr(Terms,"term.labels")
}
if(missing(fixed)) fixed <- list(noise=NULL, bias=NULL)
if(is.null(fixed$noise)) fixed$noise <- c(FALSE, TRUE)
if(is.null(fixed$bias)) fixed$bias <- c(FALSE, FALSE)
if(length(fixed$noise) == 1) fixed$noise <- rep(fixed$noise, 2)
if(equal.noise) fixed$noise[2] <- TRUE
if(length(fixed$bias) == 1) fixed$bias <- rep(fixed$bias, 2)
fixed <- list(
list(noise=fixed$noise[1], coeffs=c(TRUE, TRUE), bias=fixed$bias[1]),
list(noise=fixed$noise[2], coeffs=c(TRUE, TRUE), bias=fixed$bias[2])
)
names(fixed) <- names(par)
fixidx <- unlist(fixed)
new_par <- par
initpar <- unlist(new_par)[!fixidx]
lims <- apply(x,2,function(x){extendrange(range(x))})
if(missing(lower)) lower <- unlist(gcjcStruct(noise=.001,bias=lims[1,],config))[!fixidx]
if(missing(upper)) upper <- unlist(gcjcStruct(noise=500,bias=lims[2,],config))[!fixidx]
#make sure that the lower and upper bouds are indeed lower and upper
lims <- apply(rbind(lower,upper),2,range)
lower <- lims[1,]
upper <- lims[2,]
.optfun <- function(par, response, x, zlimit, skeleton, fixed, equal.noise){
ulpar <- unlist(skeleton)
ulpar[!fixed] <- par
params <- relist(ulpar, skeleton=skeleton)
if(equal.noise)
for(i in 1:length(params)){ params[[i]]$noise <- params[[1L]]$noise }
class(params) <- c("gcjcStruct","list")
-logLik.gcjcStruct(params,response,x,zlimit)
}
# Search for the parameters
opt_method <- match.arg(opt)
if(opt_method == "nlminb"){
optRes <- nlminb(start=initpar, objective = .optfun,
gradient = NULL, hessian = NULL,
response=response, x=x, zlimit=zlimit, skeleton=new_par, fixed=fixidx, equal.noise=equal.noise,
lower=lower, upper=upper, control=control)
optRes$value <- optRes$objective
} else {
optRes <- optim(par=initpar, fn = .optfun,
response=response, x=x, zlimit=zlimit, skeleton=new_par, fixed=fixidx, equal.noise=equal.noise,
lower=lower, upper=upper, method="L-BFGS-B", control=control)
}
if(optRes$convergence){
msg <- paste(opt_method, " problem, convergence error code = ",
optRes$convergence, "\n message = ", optRes$message,
sep='')
warning(msg)
}
# Format the output
fit <- NULL
fit$terms <- Terms
cl <- match.call()
cl[[1L]] <- as.name("goo")
fit$call <- cl
fit$contrasts <- attr(x, "contrasts")
fit$model <- mf
fit$category <- category
fit$convergence <- optRes$convergence
fit$initpar <- par
#format the fitted parameters
params <- unlist(par)
params[!fixidx] <- optRes$par
params <- relist(params, skeleton=par)
if(equal.noise)
for(i in 1:length(params)){ params[[i]]$noise <- params[[1L]]$noise }
fit$par <- params
fit$logLik <- -optRes$value
#fit$RMSE <- RMSE.gcjcStruct(object=params, response=response, x=x, zlimit=zlimit)
fit$convergence <- optRes$convergence
attr(fit$logLik, "df") <- length(optRes$par)
class(fit$logLik) <- "logLik"
class(fit) <- "gcjc"
fit
}
gcjcStruct <- function(noise, bias, config=c(1,2,3,4))
{
if(length(noise) == 1)
noise <- rep(noise, length.out=2)
if(length(bias) != 2)
stop("length(bias) must be 2")
signs <- list(c(-1,1,-1,1), c(1,-1,-1,1), c(1,-1,1,-1), c(-1,1,1,-1))[[config]]
par <- list(
list(noise=noise[1], coeffs= signs[1]*c(1,0), bias=signs[2]*bias[1]),
list(noise=noise[2], coeffs= signs[3]*c(0,1), bias=signs[4]*bias[2])
)
class(par) <- c("gcjcStruct", "list")
par
}
print.gcjc <- function(x, digits = getOption("digits"), ...)
{
cat('\nNoise:\n');
print(sapply(x$par, "[[", "noise"))
cat('\nBoundary Parameters :\n');
print(unlist(coef(x)))
cat('\nNegative Loglikelihood:\n');
cat(format(-x$logLik, digits=digits),
" (df=",format(attr(x$logLik,"df")),")\n",sep="")
cat('\nAIC score:\n');
print(AIC(x), digits=digits)
}
extractAIC.gcjc <- function(fit, scale, k = 2, ...)
{
loglik <- fit$logLik
edf <- attr(loglik, "df")
c(edf, -2 * loglik + k * edf)
}
coef.gcjc <- function(object, ...)
{
.coef <- function(obj)
{
val <- sum(obj$bias * obj$coeffs * -1)
names(val) <- "(Intercept)"
val
}
lapply(object$par, .coef)
}
logLik.gcjcStruct <-
function(object,
response,
x,
zlimit = Inf, ...)
{
dimen <- ncol(as.matrix(x))
x <- as.matrix(cbind(x,1))
if(nrow(x) != length(response))
stop("nrow(x) and length(response) are different")
g <- as.factor(response)
lev <- levels(g)
#define a fuction to calculate probability
.prob <- function(obj, x, zlimit)
{
a <- obj$coeffs
if(is.null(obj$angle)) a <- obj$coeffs
else a <- angle2cart(obj$angle)
z_coefs <- c(a, obj$bias) / obj$noise
z <- as.matrix(x) %*% as.matrix(z_coefs)
# Truncate the large z-scores
z[z < -zlimit] <- -zlimit
z[z > zlimit] <- zlimit
pnorm(z, lower.tail=FALSE, log.p=FALSE)
}
#do the rest
pr <- sapply(object, .prob, x, zlimit)
prA <- apply(pr, 1, prod)
pr <- c( prA[g==lev[1]], 1 - prA[g!=lev[1]] )
res <- sum(log(pr))
attr(res, "df") <- dimen + 1
class(res) <- "logLik"
res
}
logLik.gcjc <- function(object, ...)
{
val <- object$logLik
class(val) <- "logLik"
val
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.