gqc <- function(formula,
data,
category,
par = list(),
zlimit = Inf,
fixed = list(),
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(dimen < 2) stop("too few number of dimentions")
if(missing(category)) category <- response
if(nrow(data) != length(category))
stop("nrow(data) and length(category) are different")
#Set up the initial parameters, find one if missing
if(missing(par)) par <- NULL
if(is.null(par$pnoise)) par$pnoise <- 10
if(is.null(par$cnoise)) par$cnoise <- 100
if(is.null(par$coeff) | is.null(par$bias))
{
mc <- mcovs.formula(formula, data=data, pooled=FALSE)
par <- qdb(means=mc$means, covs=mc$covs)
}
ncoeffs <- length(par$coeffs)
if(missing(fixed)) fixed <- list(pnoise=FALSE, cnoise=FALSE, coeffs=FALSE, bias=FALSE)
if(is.null(fixed$pnoise)) fixed$pnoise <- FALSE
if(is.null(fixed$cnoise)) fixed$cnoise <- FALSE
if(is.null(fixed$coeffs)) fixed$coeffs <- rep(FALSE, ncoeffs)
if(length(fixed$coeffs) != ncoeffs)
fixed$coeffs <- rep(fixed$coeffs, ncoeffs)[1:ncoeffs]
if(is.null(fixed$bias)) fixed$bias <- FALSE
fixidx <- unlist(fixed[c("pnoise","cnoise","coeffs","bias")])
if(all(fixidx)) stop("no free parameters to fit.")
initpar <- unlist(par)[!fixidx]
if(missing(lower)){
lower = c(.1, .1, rep(-Inf, length(unlist(par))-2))
lower = lower[!fixidx]
}
if(missing(upper)){
upper = c(5000, 5000, rep( Inf, length(unlist(par))-2))
upper = upper[!fixidx]
}
.optfun <- function(par, response, x, zlimit, skeleton, fixed){
params <- unlist(skeleton)
params[!fixed] <- par
params <- relist(params, skeleton=skeleton)
-logLik(params, response, x, zlimit)
}
#Fit the model
opt <- match.arg(opt)
if(opt == "nlminb"){
optRes <- nlminb(start=initpar, objective=.optfun,
gradient = NULL, hessian = NULL,
response = response, x = x, zlimit=zlimit, skeleton=par, fixed=fixidx,
lower = lower, upper = upper, control = control)
optRes$value <- optRes$objective
} else {
optRes <- optim(par=initpar, fn=.optfun,
response = response, x = x, zlimit=zlimit, skeleton=par, fixed = fixidx,
method = "L-BFGS-B", lower = lower, upper = upper, control = control)
}
if(optRes$convergence){
msg <- paste(opt, " problem, convergence error code = ",
optRes$convergence, "\n message = ", optRes$message,
sep='')
warning(msg)
}
fit <- NULL
fit$terms <- Terms
cl <- match.call()
cl[[1L]] <- as.name("gqc")
fit$call <- cl
fit$contrasts <- attr(x, "contrasts")
fit$model <- mf
fit$category <- category
fit$initpar <- par
#format the fitted parameters
params <- unlist(par)
params[!fixidx] <- optRes$par
fit$par <- relist(params, skeleton=par)
fit$logLik <- -optRes$value
attr(fit$logLik, "df") <- length(optRes$par) - 1
class(fit$logLik) <- "logLik"
class(fit) <- "gqc"
fit
}
gqcStruct <- function(pnoise, cnoise, coeffs, bias)
{
par <- NULL
par$pnoise <- pnoise
par$cnoise <- cnoise
par$coeffs <- coeffs
par$bias <- bias
class(par) <- c("gqcStruct", "list")
par
}
print.gqc <- function(x, digits = max(5, getOption("digits") - 3), ...)
{
cat('\nPerceptual Noise:\n');
print(x$par$pnoise, digits=digits);
cat('\nCriterial Noise:\n');
print(x$par$cnoise, digits=digits);
cat('\nBoundary Parameters:\n');
print(c(x$par$coeff, x$par$bias), digits=digits);
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)
}
logLik.gqcStruct <- function(object, response, x, zlimit = Inf, ...)
{
x <- as.matrix(x)
dimen <- ncol(x)
nr <- nrow(x)
if(nrow(x) != length(response)) stop("nrow(x) and length(response) are different")
lev <- levels(g <- as.factor(response))
ncoeff <- length(coeffs <- object$coeffs)
idx <- 1:dimen
if(dimen < 2 | dimen > 3) stop("Unsupported number of dimenensions.")
if((sum(idx) + dimen) != ncoeff) stop("Inappropriate number of coeffs")
A <- diag(coeffs[idx], nrow=dimen, ncol=dimen)
A[lower.tri(A)] <- A[upper.tri(A)] <- (.5 * coeffs[(dimen+1):(ncoeff-dimen)])
b <- coeffs[(ncoeff-dimen+1):(ncoeff)]
c0 <- as.vector(object$bias)
pnoise <- object$pnoise
cnoise <- object$cnoise
# Compute mean of h(x) for each x (Eqn.11, Ashby p.462).
## obtain index for product terms
comb <- cbind(rbind(idx,idx),combn(idx,2))
tmp1 <- pnoise*sum(diag(A))
meanhxs <- tmp1 + cbind(x[,comb[1,]] * x[,comb[2,]], x, 1) %*% c(coeffs, c0)
# Compute variance of h(x) for each x (Eqn.12, Ashby p.462)
# Take advantage of fact that perceptual noise matrix is diagonal with
# pnoise the same in all three dimenensions
bvec <- matrix(rep(b,nr), ncol=nr, nrow=dimen, byrow=FALSE)
varhxs <- 2*tmp1*tmp1 + pnoise*colSums((bvec + 2*A %*% t(x))^2)
#Compute z-scores for each data point
z <- drop( meanhxs / sqrt(varhxs+cnoise) )
#Truncate the large z-scores
z[z < -zlimit] <- -zlimit
z[z > zlimit] <- zlimit
log_A_probs <- pnorm(z[g == lev[1]], lower.tail=FALSE, log.p=TRUE)
log_B_probs <- pnorm(z[g != lev[1]], lower.tail=TRUE, log.p=TRUE)
res <- sum(log_A_probs, log_B_probs)
attr(res, "df") <- length(unlist(object)) - 1
class(res) <- "logLik"
res
}
logLik.gqc <- function(object, ...)
{
val <- object$logLik
class(val) <- "logLik"
val
}
extractAIC.gqc <- function(fit, scale, k = 2, ...)
{
loglik <- fit$logLik
edf <- attr(loglik, "df")
c(edf, -2 * loglik + k * edf)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.