Nothing
gof <-
function(observed, simulated, weights=NULL) {
if (class(observed) == "numeric") {
observed <- matrix(observed,nrow=1)
}
if (class(observed) != "matrix" | class(simulated) != "gof.preprocess") {
stop("simulated must be a gof.preprocess object.")
}
if (ncol(observed) != simulated$Variates) {
stop("Dimensionality of function parameters do not match.")
}
observations = nrow(observed)
if (is.null(weights)) {
if (length(observed) > 1) {
weights = matrix(0, simulated$Variates+1, simulated$Variates+1)
diag(weights) <- 1 / (simulated$Variates+1)
} else {
weights = 1
}
}
foldedObservations <- sapply(1:observations, function(i) abs(observed[i,]-simulated$Expectation))
p <- .Call("doTest", as.matrix(foldedObservations), as.matrix(simulated$FoldedSimulations),
as.matrix(simulated$CenteredComparisons), as.matrix(simulated$ComparisonCrease),
as.matrix(weights))
# Comparisons in R instead of C
# tailComparisons.observed = sapply(1:observations, function (i) gof.compare(abs(observed[i,]-simulated$Expectation), simulated$FoldedSimulations) )
# tailComparisons.observed.centered = abs(tailComparisons.observed - simulated$ComparisonCrease)
# testStatistic.observed = as.numeric(sapply(1:observations, function(i) t(tailComparisons.observed.centered[,i]) %*% weights %*% tailComparisons.observed.centered[,i]))
# testStatistic.simulated = as.numeric(sapply(1:simulated$Iterations, function(i) t(simulated$CenteredComparisons[,i]) %*% weights %*% simulated$CenteredComparisons[,i]))
# edf.observed = sapply(1:observations, function (i) sum(testStatistic.simulated <= testStatistic.observed[i]) / simulated$Iterations )
# p <- 1 - abs(1 - 2 * edf.observed)
mhd <- sapply(1:observations, function(i) t(observed[i,] - simulated$Expectation) %*% simulated$InverseCovariance %*% (observed[i,] - simulated$Expectation) )
p.mhd <- sapply(1:observations, function (i) 1 - abs(1 - 2 * sum(mhd[i] >= simulated$MHDistances)/simulated$Iterations) )
return = list(p=p, #v.obs=testStatistic.observed, v.sim=testStatistic.simulated,
Observation=observed, Preprocess=simulated, MHDistance=mhd, p.MHD = p.mhd)
class(return) <- "gof.result"
return
}
print.gof.result <- function (x, ...) {
with(x, cat("snopgof results\nVariates:", Preprocess$Variates ,"\nSimulations:", Preprocess$Iterations ,"\nP-Value:", p,"\nP-Value(MHD):", p.MHD))
}
plot.gof.result <- function (x, standardize=3, violin=TRUE,
main="SNOPGOF Goodness of Fit Test", ylim=NULL, xlim=NULL,
xlab=NULL, ylab=NULL, key=NULL, perc=.05, ...) {
itns <- x$Preprocess$Iterations
vars <- x$Preprocess$Variates
sims <- x$Preprocess$Original
obs <- x$Observation
n.obs <- nrow(obs)
if (standardize==3) {
sims.median <- apply(sims, 2, median)
sims.min <- apply(sims, 2, min)
sims.max <- apply(sims, 2, max)
sims <- sapply(1:ncol(sims), function(i) (sims[,i] - sims.median[i])/(sims.max[i] - sims.min[i] ) )
obs <- matrix(sapply(1:ncol(sims), function(i) (obs[,i] - sims.median[i])/(sims.max[i] - sims.min[i] ) ), nrow=n.obs )
} else if (standardize==2) {
sims.min <- apply(sims, 2, min)
sims.max <- apply(sims, 2, max)
sims <- sapply(1:ncol(sims), function(i) (sims[,i] - sims.min[i])/(sims.max[i] - sims.min[i] ) )
obs <- matrix(sapply(1:ncol(sims), function(i) (obs[,i] - sims.min[i])/(sims.max[i] - sims.min[i] ) ), nrow=n.obs )
} else if (standardize==1) {
sims.mean <- apply(sims, 2, mean)
sims <- sapply(1:ncol(sims), function(i) (sims[,i] - sims.mean[i]) )
obs <- matrix(sapply(1:ncol(sims), function(i) (obs[,i] - sims.mean[i]) ), nrow=n.obs )
}
if (is.null(ylim)) {
ylim = c(min(obs, sims), max(obs, sims))
}
if (is.null(xlim)) {
xlim = c(0, ncol(obs)+1)
}
if (is.null(xlab)) {
xlab= paste( paste("p:", round(x$p, 3), collapse = " "), paste("p(MHD):", round(x$p.MHD, 3), collapse = " "), collapse = "\n")
}
if (is.null(ylab)) {
ylab = "Statistic Values"
}
xAxis <- (1:vars)
plot(obs[1,]~xAxis, col="white", type="p", ylim=ylim, xlim=xlim, main=main, xlab=xlab, ylab=ylab, axes=FALSE, ...)
if (!is.null(key)) {
if (length(key) != ncol(obs)) {
stop("Key length does not match the number of variates.")
}
axis(1, at=xAxis, lab=key)
} else {
axis(1, at=xAxis, lab=paste("v", xAxis, sep=""))
}
ind.lower = round(itns * perc/2)
ind.upper = round(itns * (1-perc/2))
yperc.lower = sapply(1:ncol(sims), function(i) sort(sims[,i])[ind.lower] )
yperc.upper = sapply(1:ncol(sims), function(i) sort(sims[,i])[ind.upper] )
lines(yperc.lower~xAxis, lty=3, col = "gray", lwd=3)
lines(yperc.upper~xAxis, lty=3, col = "gray", lwd=3)
if (violin) {
require(vioplot)
for (i in 1:ncol(sims)) {
vioplot(sims[,i], at=xAxis[i], add=TRUE, col="gray", wex=.75, ...)
}
} else {
boxplot(as.numeric(sim)~rep(1:vars, each=itns), add=TRUE, ...)
}
for(i in 1:nrow(obs)) {
lines(obs[i,]~xAxis, col="red", type="l", lwd=1, ...)
lines(obs[i,]~xAxis, col="red", type="p", lwd=3, pch=19, ...)
text(xAxis, obs[i,], labels=round(x$Observation[i,],3), pos=4)
}
}
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.