Nothing
"initialize.x" <-
function(model, x1,x2,xrest=NULL) {
## Extract model components
n <- model$n
logp.position <- model$logp.position
mask.x <- model$mask.x
## x1, x2 and logp at max
x1.max <- as.double(rep(NA,n))
x2.max <- as.double(rep(NA,n))
logp.max <- as.double(rep(NA,n))
print(paste("i in", length(x1)))
## Loop over grid
for(i in 1:length(x1)) {
print(i)
for(j in 1:length(x2)) {
x <- cbind(rep(x1[i],n),rep(x2[j],n),xrest)
logp <- logp.position(x)
## Keep point with max logp that satisfies the mask
keep <- mask.x(x[ , 1:2]) & (is.na(logp.max) | (logp > logp.max))
logp.max[keep] <- logp[keep]
x1.max[keep] <- x1[i]
x2.max[keep] <- x2[j]
}
}
## The full x at the maximum
x <- cbind(x1.max,x2.max,xrest)
## Replace fixed points
##x[model$fixed.x,] <- model$start.x[model$fixed.x,]
x
}
## "position.logp" <-
## function (model, x1, x2, xrest = NULL, subset = 1:model$n, initialize.x = TRUE,
## start = NULL, end = NULL, prob = 0.8)
## {
## n <- model$n
## logp.position <- model$logp.position
## mask.x <- model$mask.x
## logp <- array(0, c(length(x1), length(x2), length(subset)))
## mask <- array(0, c(length(x1), length(x2), length(subset)))
## print(paste("i in ", length(x1)))
## for (i in 1:length(x1)) {
## print(i)
## for (j in 1:length(x2)) {
## x <- cbind(rep(x1[i], n), rep(x2[j], n), xrest)
## logp[i, j, ] <- logp.position(x)[subset]
## mask[i, j, ] <- mask.x(x)[subset]
## }
## }
## res <- list(x = x1, y = x2, logp = logp, mask = mask)
## if (initialize.x) {
## #logp.init <- function(d.model, logp, ) {
## #nq = 6
## xy <- expand.grid(x = x1, y = x2)
## logp.quantile <- logp > 0
## for (i in 2:(n - 1)) logp.quantile[, , i] <-
## logp[, , i] > quantile(logp[, , i], prob, na.rm = TRUE)
## logp <- logp.quantile
## xx <- matrix(0, nrow = n, ncol = 3)
## for (i in 2:(n - 1)) {
## logp.quantile[,,i] <- ((logp[, , i - 1] * logp[,
## , i]) & (logp[, , i] * logp[, , i + 1]) == 1) * mask[,,i]
## #this <- this * mask[,,i]
## xx[i, 1:2] <- apply(xy[as.logical(logp.quantile[,,i] ), ], 2, mean, na.rm = TRUE)
## }
## #browser()
## res$logp.quantile <- logp.quantile
## X <- xx
## if (!is.null(start)) X[1,1:2] <- start
## if (!is.null(end)) X[nrow(X), 1:2] <- end
## require(zoo)
## res$naX <- X
## X <- as.matrix(na.approx(X, na.rm = FALSE))
## res[["X"]] <- X
## }
## class(res) <- c("diag", class(res))
## res
## }
"position.logp" <-
function (model, x1, x2, xrest = NULL, subset = 1:model$n, initialize.x = TRUE,
start = NULL, end = NULL, prob = 0.8, winoffset = 5)
{
n <- model$n
if (n < (2 * winoffset - 1)) {
stop("too few twilights (", n, ") for winoffset value of ", winoffset, "\n try reducing winoffset value")
}
logp.position <- model$logp.position
mask.x <- model$mask.x
logp <- array(0, c(length(x1), length(x2), length(subset)))
mask <- array(0, c(length(x1), length(x2), length(subset)))
print(paste("i in ", length(x1)))
for (i in 1:length(x1)) {
cat(i, "\n")
for (j in 1:length(x2)) {
x <- cbind(rep(x1[i], n), rep(x2[j], n), xrest)
logp[i, j, ] <- logp.position(x)[subset]
mask[i, j, ] <- mask.x(x)[subset]
}
}
res <- list(x = x1, y = x2, logp = logp, mask = mask)
if (initialize.x) {
xy <- expand.grid(x = x1, y = x2)
logp.quantile <- logp > 0
for (i in 2:(n - 1)) {
logp.quantile[, , i] <- logp[, ,i] > quantile(logp[, , i], prob, na.rm = TRUE)
}
logp <- logp.quantile
xx <- cbind(matrix(as.numeric(NA), nrow = n, ncol = 2), xrest)
for (i in winoffset:(n - winoffset)) {
mm <- mask[,,i] * 0
for (ik in (i-winoffset+1):(i+winoffset -1)) {
mm <- mm + logp[,,ik]
}
#logp.quantile[,,i] <- (!mm == 0) & mask[,,i]
xy1 <- xy[which.max(as.vector(mm) * mask[,,i]), ]
if (nrow(xy1) > 0) {
xx[i, 1:2] <- as.matrix(xy1)[ceiling(nrow(xy1)/2), ,drop = TRUE]
} else {
xx[i,1:2] <- rep(NA, 2)
}
}
res$logp.quantile <- logp.quantile
X <- xx
if (!is.null(start)) {
X[1:(winoffset-1), 1:2] <- matrix(start, nrow = winoffset - 1, ncol = 2, byrow = TRUE)
}
if (!is.null(end)) {
X[(nrow(X) - winoffset+2):nrow(X), 1:2] <- matrix(end, nrow = winoffset - 1, ncol = 2, byrow = TRUE)
}
## require(zoo) removed by namespace addition MDS 2011-10-06
res$naX <- X
X <- as.matrix(na.approx(X, na.rm = FALSE))
res[["X"]] <- X
}
class(res) <- c("diag", class(res))
res
}
"light.quantile" <-
function(model,chain,day,seg,probl=c(0.025,0.5,0.975)) {
n <- length(day)
## Parameters for this segment
xs <- t(chain$x[seg,,])
## Matrix of quantiles
qs <- matrix(0,n,length(probl))
for(i in 1:n) {
## Predictions for this time over all xs
lgt <- model$light.predict(xs,day[i])
qs[i,] <- quantile(lgt,prob=probl)
}
qs
}
"show.segment" <- function (model, chain, segment, day, light, k, n = 50, ...)
{
segment <- factor(segment)
day.seg <- day[segment == k]
plot(day.seg, light[segment == k], ...)
text(mean(day.seg), mean(light[segment == k])-50, k, font = 2)
ds <- seq(min(day.seg), max(day.seg), length = n)
qs <- light.quantile(model, chain, ds, k)
matlines(ds, qs, lty = 1, col = "blue", lwd = 2)
}
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.