Nothing
awstestprop <- function(dy,
hmax,
theta = 1,
family = "Gaussian",
lkern = "Triangle",
aws = TRUE,
memory = FALSE,
shape = 2,
homogeneous = TRUE,
varadapt = FALSE,
ladjust = 1,
spmin = 0.25,
seed = 1,
minlevel = 1e-6,
maxz = 25,
diffz = .5,
maxni = FALSE,
verbose = FALSE) {
if (length(dy) > 3) {
cat(
"maximum array dimension is 3\n contents of first argument will be interpreted as array of
parameters\n"
)
nnn <- length(dy)
} else {
nnn <- prod(dy)
}
if (minlevel < 5 / nnn) {
minlevel <- 5 / nnn
cat("minlevel reset to",
minlevel,
"due to insufficient size of test sample\n")
}
set.seed(seed)
par(
mfrow = c(1, 1),
mar = c(3, 3, 3, 3),
mgp = c(2, 1, 0)
)
if (length(dy) <= 3) {
y <- array(switch(
family,
"Gaussian" = rnorm(nnn),
"Poisson" = rpois(nnn, theta),
"Exponential" = rexp(nnn, 1),
"Bernoulli" = rbinom(nnn, 1, theta),
"Volatility" = rnorm(nnn),
"Variance" = rchisq(nnn, shape) / shape,
"NCchi" = sqrt(rchisq(nnn, shape, theta ^ 2))
), dy)
} else {
ddy <- if (!is.null(dim(dy)))
dim(dy)
else
length(dy)
y <- array(switch(
family,
"Gaussian" = rnorm(nnn, dy - mean(dy)),
"Poisson" = rpois(nnn, dy - mean(dy) + theta),
"Exponential" = rexp(nnn, dy - mean(dy) + 1),
"Bernoulli" = rbinom(nnn, 1, dy - mean(dy) + theta),
"Volatility" = rnorm(nnn, dy - mean(dy)),
"Variance" = rchisq(nnn, dy - mean(dy) + shape) /
(dy - mean(dy) + shape),
"NCchi" = sqrt(rchisq(nnn, shape, (
dy - mean(dy) + theta
) ^ 2))
), ddy)
dy <- ddy
}
cat("minlevel ",
minlevel,
"due to insufficient size of test sample\n")
if (family == "NCchi") {
varstats <-
sofmchi(shape / 2) # precompute table of mean, sd and var for
#
# NCchi for noncentral chi with shape=degrees of freedom and theta =NCP
#
}
z <- seq(0, maxz, diffz)
nz <- length(z)
elevel <- trunc(log(1e-6, 10))
levels <- as.vector(outer(c(.5, .2, .1), 10 ^ (-0:elevel)))
levels <- levels[levels >= minlevel]
wghts <- switch(length(dy), c(0, 0), c(1, 0), c(1, 1))
cpar <-
setawsdefaults(dy,
mean(y),
family,
lkern,
"Uniform",
aws,
memory,
ladjust,
hmax,
shape,
wghts)
lambda <- cpar$lambda
hmax <- cpar$hmax
shape <- cpar$shape
d <- cpar$d
n <- length(y)
if (!homogeneous & family == "Gaussian") {
sigma2 <- array(rchisq(prod(dy), shape) / shape, dy)
} else
sigma2 <- 1
zfamily <- awsfamily(family, y, sigma2, shape, 0, lambda, cpar)
cpar <- zfamily$cpar
lambda <- zfamily$lambda
sigma2 <- zfamily$sigma2
h0 <- zfamily$h0
y <- zfamily$y
lkern <- cpar$lkern
rm(zfamily)
n1 <- switch(d, n, dy[1], dy[1])
n2 <- switch(d, 1, dy[2], dy[2])
n3 <- switch(d, 1, 1, dy[3])
maxvol <- cpar$maxvol
# k <- cpar$k
k <- 1
kstar <- cpar$kstar
h <- numeric(kstar)
if (k > 1)
h[1:(k - 1)] <- 1 + (0:(k - 2)) * .001
exceedence <-
exceedencena <-
matrix(0, nz, kstar) # this is used to store exceedence probabilities for adaptive and nonadaptive estimates
zobj <- zobj0 <- list(ai = y, bi = rep(1, n))
bi <- rep(1, n)
yhat <- y / shape
lambda0 <- 1e50
total <- cumsum(1.25 ^ (1:kstar)) / sum(1.25 ^ (1:kstar))
#
# get initial conditions for a comparison
#
if (family == "Bernoulli")
y0 <- (10 * y + 1) / 12
if (family == "Poisson")
y0 <- y + .1
#
# this corresponds to the regularization used to avoid Inf distances
#
kldistnorm1 <- function(th1, y, df) {
L <- df / 2
m1 <-
sqrt(pi / 2) * gamma(L + 1 / 2) / gamma(1.5) / gamma(L) * gsl::hyperg_1F1(-0.5, L, -th1 ^
2 / 2, give = FALSE, strict = TRUE)
(m1 - y) ^ 2 / 2 / (2 * L + th1 ^ 2 - m1 ^ 2)
}
KLdist0 <- switch(
family,
"Gaussian" = y ^ 2 / 2,
"Poisson" = (theta - y0 + y0 * (log(y0) - log(theta))),
"Exponential" = (log(y) - 1 + 1 / y),
"Bernoulli" = (y0 * log(y0 / theta) +
(1 - y0) * log((1 - y0) / (1 - theta))),
"Volatility" = (log(y) - 1 + 1 / y) / 2,
"Variance" = shape / 2 * (log(y) - 1 + 1 / y),
"NCchi" = kldistnorm1(theta, y, shape)
)
exceedence0 <- .Fortran(C_exceed,
as.double(KLdist0),
as.integer(length(KLdist0)),
as.double(z),
as.integer(nz),
exprob = double(nz)
)$exprob
#
# now iterate
#
t0 <- Sys.time()
cat("using lambda=", lambda, "\n")
while (k <= kstar) {
t1 <- Sys.time()
hakt0 <-
gethani(1, 1.25 * hmax, lkern, 1.25 ^ (k - 1), wghts, 1e-4)
hakt <- gethani(1, 1.25 * hmax, lkern, 1.25 ^ k, wghts, 1e-4)
cat("step", k, "hakt", hakt, "\n")
if (lkern == 5) {
# assume hmax was given in FWHM units (Gaussian kernel will be truncated at 4)
hakt <- hakt * 0.42445 * 4
}
h[k] <- hakt
dlw <- (2 * trunc(hakt / c(1, wghts)) + 1)[1:d]
#
# get nonadaptive estimate
#
if (!homogeneous & family == "Gaussian") {
zobj0 <- .Fortran(C_chaws,
as.double(y),
as.double(sigma2),
as.integer(1:n),#position
as.integer(n1),
as.integer(n2),
as.integer(n3),
hakt = as.double(hakt),
as.double(1e40),
double(n),
bi = double(n),
bi2 = double(n),
double(n),
ai = double(n),
as.integer(cpar$mcode),
as.integer(lkern),
as.double(spmin),
double(prod(dlw)),
as.double(wghts)
)[c("bi", "bi2", "ai", "hakt")]
} else {
zobj0 <- .Fortran(C_caws,
as.double(y),
as.integer(1:n),# position for full mask
as.integer(n1),
as.integer(n2),
as.integer(n3),
hakt = as.double(hakt),
as.double(1e40),
double(n),
bi = double(n),
bi2 = double(n),
double(n),
ai = double(n),
as.integer(cpar$mcode),
as.integer(lkern),
as.double(spmin),
double(prod(dlw)),
as.double(wghts)
)[c("bi", "bi2", "ai", "hakt")]
}
if (family %in% c("Bernoulli", "Poisson"))
zobj0 <- regularize(zobj0, family)
yhat0 <- zobj0$ai / zobj0$bi
dim(yhat0) <- dy
bi <- zobj0$bi
ih <- as.integer(hakt)
ind1 <- (ih + 1):(dy[1] - ih)
if (length(dy) > 1)
ind2 <- (ih + 1):(dy[2] - ih)
if (length(dy) > 2)
ind3 <- (ih + 1):(dy[3] - ih)
yhat0 <-
switch(length(dy), yhat0[ind1], yhat0[ind1, ind2], yhat0[ind1, ind2, ind3])
ni <- max(zobj0$bi)
KLdist0 <- switch(
family,
"Gaussian" = yhat0 ^ 2 / 2,
"Poisson" = (theta - yhat0 + yhat0 * (log(yhat0) -
log(theta))),
"Exponential" = (log(yhat0) - 1 + 1 / yhat0),
"Bernoulli" = (yhat0 * log(yhat0 / theta) +
(1 - yhat0) * log((1 - yhat0) / (1 -
theta))),
"Volatility" = (log(yhat0) - 1 + 1 / yhat0) / 2,
"Variance" = shape / 2 * (log(yhat0) - 1 + 1 /
yhat0),
"NCchi" = kldistnorm1(theta, yhat0, shape)
)
exceedencena[, k] <- .Fortran(C_exceed,
as.double(KLdist0),
as.integer(length(KLdist0)),
as.double(z / ni),
as.integer(nz),
exprob = double(nz)
)$exprob
#
# get adaptive estimate
#
if (!homogeneous & family == "Gaussian") {
zobj <- .Fortran(C_chaws,
as.double(y),
as.double(sigma2),
as.integer(1:n),# position for full mask
as.integer(n1),
as.integer(n2),
as.integer(n3),
hakt = as.double(hakt),
as.double(lambda0),
as.double(yhat),
bi = as.double(bi),
bi2 = double(n),
double(n),
ai = double(n),
as.integer(cpar$mcode),
as.integer(lkern),
as.double(spmin),
double(prod(dlw)),
as.double(wghts)
)[c("bi", "bi2", "ai", "hakt")]
} else {
if (cpar$mcode != 6) {
zobj <- .Fortran(C_caws,
as.double(y),
as.integer(1:n),# position for full mask
as.integer(n1),
as.integer(n2),
as.integer(n3),
hakt = as.double(hakt),
as.double(lambda0),
as.double(yhat),
bi = as.double(bi),
bi2 = double(n),
double(n),
ai = double(n),
as.integer(cpar$mcode),
as.integer(lkern),
as.double(spmin),
double(prod(dlw)),
as.double(wghts)
)[c("bi", "bi2", "ai", "hakt")]
} else {
zobj <- .Fortran(C_caws6,
as.double(y),
as.integer(1:n),# position for full mask
as.integer(n1),
as.integer(n2),
as.integer(n3),
hakt = as.double(hakt),
as.double(lambda0),
as.double(yhat),
as.double(fncchiv(yhat, varstats) / 2),
bi = as.double(bi),
bi2 = double(n),
double(n),
ai = double(n),
as.integer(lkern),
as.double(spmin),
double(prod(dlw)),
as.double(wghts)
)[c("bi", "bi2", "ai", "hakt")]
}
}
if (family %in% c("Bernoulli", "Poisson"))
zobj <- regularize(zobj, family)
dim(zobj$ai) <- dy
yhat <- zobj$ai / zobj$bi
dim(yhat) <- dy
if (varadapt)
bi <- bi ^ 2 / zobj$bi2
if (maxni)
bi <- pmax(bi, zobj$bi)
else
bi <- zobj$bi
lambda0 <- lambda
yhat0 <-
switch(length(dy), yhat[ind1], yhat[ind1, ind2], yhat[ind1, ind2, ind3])
KLdist1 <- switch(
family,
"Gaussian" = yhat0 ^ 2 / 2,
"Poisson" = (theta - yhat0 + yhat0 * (log(yhat0) -
log(theta))),
"Exponential" = (log(yhat0) - 1 + 1 / yhat0),
"Bernoulli" = (yhat0 * log(yhat0 / theta) +
(1 - yhat0) * log((1 - yhat0) / (1 -
theta))),
"Volatility" = (log(yhat0) - 1 + 1 / yhat0) / 2,
"Variance" = shape / 2 * (log(yhat0) - 1 + 1 /
yhat0),
"NCchi" = kldistnorm1(theta, yhat0, shape)
)
exceedence[, k] <- .Fortran(C_exceed,
as.double(KLdist1),
as.integer(length(KLdist1)),
as.double(z / ni),
as.integer(nz),
exprob = double(nz)
)$exprob
contour(
z,
0:k,
cbind(exceedence0, exceedence[, 1:k]),
levels = levels,
ylab = "step",
xlab = "z",
main = paste(family, length(dy), "-dim. ladj=", ladjust, " Exceed. Prob.")
)
contour(
z,
0:k,
cbind(exceedence0, exceedencena[, 1:k]),
levels = levels,
ylab = "step",
xlab = "z",
add = TRUE,
col = 2,
lty = 3
)
yaxp <- par("yaxp")
at <-
unique(as.integer(seq(yaxp[1], yaxp[2], length = yaxp[3] + 1)))
at <- at[at > 0 & at <= k]
axis(4, at = at, labels = as.character(signif(h[at], 3)))
mtext("bandwidth", 4, 1.8)
if (max(total) > 0) {
cat("Progress:", signif(total[k], 2) * 100, "% . ", sep = "")
}
t2 <- Sys.time()
tpar <-
if (family %in% c("Bernoulli", "Poisson"))
paste("theta=", theta)
else
""
cat(
family,
"(dim:",
length(dy),
tpar,
") ni=",
ni,
" Time: Step",
format(signif(difftime(t2, t1), 3)),
"Total",
format(signif(difftime(t2, t0), 3)),
"\n"
)
k <- k + 1
gc()
}
if (family %in% c("Bernoulli", "Poisson"))
y <- y0
list(
h = h,
z = z,
prob = exceedence,
probna = exceedencena,
y = if (verbose)
y
else
NULL ,
theta = if (verbose)
yhat
else
NULL,
levels = levels,
family = family
)
}
pawstestprop <- function(dy,
hmax,
theta = 1,
family = "Gaussian",
lkern = "Triangle",
aws = TRUE,
patchsize=1,
shape = 2,
ladjust = 1,
spmin = 0.25,
seed = 1,
minlevel = 1e-6,
maxz = 25,
diffz = .5,
maxni = FALSE,
verbose = FALSE) {
varadapt <- FALSE
patchadapt <- "max"
if (length(dy) > 3) {
cat(
"maximum array dimension is 3\n contents of first argument will be interpreted as array of
parameters\n"
)
nnn <- length(dy)
} else {
nnn <- prod(dy)
}
if (minlevel < 5 / nnn) {
minlevel <- 5 / nnn
cat("minlevel reset to",
minlevel,
"due to insufficient size of test sample\n")
}
set.seed(seed)
par(
mfrow = c(1, 1),
mar = c(3, 3, 3, 3),
mgp = c(2, 1, 0)
)
if (length(dy) <= 3) {
y <- array(switch(
family,
"Gaussian" = rnorm(nnn),
"Poisson" = rpois(nnn, theta),
"Exponential" = rexp(nnn, 1),
"Bernoulli" = rbinom(nnn, 1, theta),
"Volatility" = rnorm(nnn),
"Variance" = rchisq(nnn, shape) / shape,
"NCchi" = sqrt(rchisq(nnn, shape, theta ^ 2))
), dy)
} else {
ddy <- if (!is.null(dim(dy)))
dim(dy)
else
length(dy)
y <- array(switch(
family,
"Gaussian" = rnorm(nnn, dy - mean(dy)),
"Poisson" = rpois(nnn, dy - mean(dy) + theta),
"Exponential" = rexp(nnn, dy - mean(dy) + 1),
"Bernoulli" = rbinom(nnn, 1, dy - mean(dy) + theta),
"Volatility" = rnorm(nnn, dy - mean(dy)),
"Variance" = rchisq(nnn, dy - mean(dy) + shape) /
(dy - mean(dy) + shape),
"NCchi" = sqrt(rchisq(nnn, shape, (
dy - mean(dy) + theta
) ^ 2))
), ddy)
dy <- ddy
}
cat("minlevel ",
minlevel,
"due to insufficient size of test sample\n")
if (family == "NCchi") {
varstats <-
sofmchi(shape / 2) # precompute table of mean, sd and var for
#
# NCchi for noncentral chi with shape=degrees of freedom and theta =NCP
#
}
z <- seq(0, maxz, diffz)
nz <- length(z)
elevel <- trunc(log(1e-6, 10))
levels <- as.vector(outer(c(.5, .2, .1), 10 ^ (-0:elevel)))
levels <- levels[levels >= minlevel]
wghts <- switch(length(dy), c(0, 0), c(1, 0), c(1, 1))
cpar <-
setawsdefaults(dy,
mean(y),
family,
lkern,
"Uniform",
aws,
FALSE,
ladjust,
hmax,
shape,
wghts)
lambda <- cpar$lambda
# additional adjustments for taking the maximum of s_{ij} over patches
# adjusted using simulations such that for homogeneous structures
# loss by adaptation < 1% in MAE and <.1 in PSNR
if(length(dy)==1){
ladjust <- switch(patchsize,.55,.55,.5,.5)*ladjust
} else if(length(dy)==2){
ladjust <- switch(patchsize,1.1,1.1,1.2)*ladjust
} else {
ladjust <- switch(patchsize,1.44,1.8)*ladjust
}
hmax <- cpar$hmax
shape <- cpar$shape
d <- cpar$d
n <- length(y)
sigma2 <- 1
zfamily <- awsfamily(family, y, sigma2, shape, 0, lambda, cpar)
cpar <- zfamily$cpar
lambda <- zfamily$lambda
sigma2 <- zfamily$sigma2
h0 <- zfamily$h0
y <- zfamily$y
lkern <- cpar$lkern
rm(zfamily)
n1 <- switch(d, n, dy[1], dy[1])
n2 <- switch(d, 1, dy[2], dy[2])
n3 <- switch(d, 1, 1, dy[3])
maxvol <- cpar$maxvol
# k <- cpar$k
k <- 1
kstar <- cpar$kstar
h <- numeric(kstar)
if (k > 1)
h[1:(k - 1)] <- 1 + (0:(k - 2)) * .001
exceedence <-
exceedencena <-
matrix(0, nz, kstar) # this is used to store exceedence probabilities for adaptive and nonadaptive estimates
zobj <- zobj0 <- list(ai = y, bi = rep(1, n))
bi <- rep(1, n)
yhat <- y / shape
lambda0 <- 1e50
total <- cumsum(1.25 ^ (1:kstar)) / sum(1.25 ^ (1:kstar))
#
# get initial conditions for a comparison
#
if (family == "Bernoulli")
y0 <- (10 * y + 1) / 12
if (family == "Poisson")
y0 <- y + .1
#
# this corresponds to the regularization used to avoid Inf distances
#
kldistnorm1 <- function(th1, y, df) {
L <- df / 2
m1 <-
sqrt(pi / 2) * gamma(L + 1 / 2) / gamma(1.5) / gamma(L) * gsl::hyperg_1F1(-0.5, L, -th1 ^
2 / 2, give = FALSE, strict = TRUE)
(m1 - y) ^ 2 / 2 / (2 * L + th1 ^ 2 - m1 ^ 2)
}
KLdist0 <- switch(
family,
"Gaussian" = y ^ 2 / 2,
"Poisson" = (theta - y0 + y0 * (log(y0) - log(theta))),
"Exponential" = (log(y) - 1 + 1 / y),
"Bernoulli" = (y0 * log(y0 / theta) +
(1 - y0) * log((1 - y0) / (1 - theta))),
"Volatility" = (log(y) - 1 + 1 / y) / 2,
"Variance" = shape / 2 * (log(y) - 1 + 1 / y),
"NCchi" = kldistnorm1(theta, y, shape)
)
exceedence0 <- .Fortran(C_exceed,
as.double(KLdist0),
as.integer(length(KLdist0)),
as.double(z),
as.integer(nz),
exprob = double(nz)
)$exprob
#
# define patch area
#
patchsize <- min(patchsize,5-length(dy))
np1 <- if (patchsize > 0)
2 * patchsize + 1
else
7
## patchsize == 0 includes immediate neighbors only
np2 <- if (n2 > 1)
2 * patchsize + 1
else
1
np3 <- if (n3 > 1)
2 * patchsize + 1
else
1
t0 <- Sys.time()
cat("using lambda=", lambda, "\n")
mask <- array(TRUE,dy)
#
# now iterate
#
while (k <= kstar) {
t1 <- Sys.time()
hakt0 <-
gethani(1, 1.25 * hmax, lkern, 1.25 ^ (k - 1), wghts, 1e-4)
hakt <- gethani(1, 1.25 * hmax, lkern, 1.25 ^ k, wghts, 1e-4)
cat("step", k, "hakt", hakt, "\n")
if (lkern == 5) {
# assume hmax was given in FWHM units (Gaussian kernel will be truncated at 4)
hakt <- hakt * 0.42445 * 4
}
h[k] <- hakt
dlw <- (2 * trunc(hakt / c(1, wghts)) + 1)[1:d]
#
# define interior mask
#
ih <- trunc(h)+patchsize
if(nnn==1) mask[c(1:ih,n1+1-1:ih)] <- FALSE
if(nnn==2) mask[c(1:ih,n1+1-1:ih),
c(1:ih,n2+1-1:ih)] <- FALSE
if(nnn==3) mask[c(1:ih,n1+1-1:ih),
c(1:ih,n2+1-1:ih),
c(1:ih,n3+1-1:ih)] <- FALSE
#
# get nonadaptive estimate
#
zobj0 <- .Fortran(C_caws,
as.double(y),
as.integer(1:n),# position for full mask
as.integer(n1),
as.integer(n2),
as.integer(n3),
hakt = as.double(hakt),
as.double(1e40),
double(n),
bi = double(n),
bi2 = double(n),
double(n),
ai = double(n),
as.integer(cpar$mcode),
as.integer(lkern),
as.double(spmin),
double(prod(dlw)),
as.double(wghts)
)[c("bi", "bi2", "ai", "hakt")]
if (family %in% c("Bernoulli", "Poisson"))
zobj0 <- regularize(zobj0, family)
yhat0 <- zobj0$ai / zobj0$bi
dim(yhat0) <- dy
bi <- zobj0$bi
ih <- as.integer(hakt)
ind1 <- (ih + 1):(dy[1] - ih)
if (length(dy) > 1)
ind2 <- (ih + 1):(dy[2] - ih)
if (length(dy) > 2)
ind3 <- (ih + 1):(dy[3] - ih)
yhat0 <-
switch(length(dy), yhat0[ind1], yhat0[ind1, ind2], yhat0[ind1, ind2, ind3])
ni <- max(zobj0$bi)
KLdist0 <- switch(
family,
"Gaussian" = yhat0 ^ 2 / 2,
"Poisson" = (theta - yhat0 + yhat0 * (log(yhat0) -
log(theta))),
"Exponential" = (log(yhat0) - 1 + 1 / yhat0),
"Bernoulli" = (yhat0 * log(yhat0 / theta) +
(1 - yhat0) * log((1 - yhat0) / (1 -
theta))),
"Volatility" = (log(yhat0) - 1 + 1 / yhat0) / 2,
"Variance" = shape / 2 * (log(yhat0) - 1 + 1 /
yhat0),
"NCchi" = kldistnorm1(theta, yhat0, shape)
)
exceedencena[, k] <- .Fortran(C_exceedm,
as.double(KLdist0),
as.integer(length(KLdist0)),
as.double(z / ni),
as.integer(nz),
exprob = double(nz),
as.integer(mask)
)$exprob
#
# get adaptive estimate
#
zobj <- .Fortran(C_pcaws,
as.double(y),
as.double(1:n),# full mask
as.integer(n1),
as.integer(n2),
as.integer(n3),
hakt = as.double(hakt),
as.double(lambda0),
as.double(yhat),
as.double(zobj$bi),
bi2 = double(n),
bi = double(n), #biout
theta = double(n),
as.integer(cpar$mcode),
as.integer(lkern),
as.double(spmin),
double(prod(dlw)),
as.double(wghts),
as.integer(patchsize))[c("bi", "bi2", "theta", "hakt")]
if (family %in% c("Bernoulli", "Poisson"))
zobj <- regularize(zobj, family)
yhat <- zobj$theta
dim(yhat) <- dy
if (varadapt)
bi <- bi ^ 2 / zobj$bi2
if (maxni)
bi <- pmax(bi, zobj$bi)
else
bi <- zobj$bi
lambda0 <- lambda
yhat0 <-
switch(length(dy), yhat[ind1], yhat[ind1, ind2], yhat[ind1, ind2, ind3])
KLdist1 <- switch(
family,
"Gaussian" = yhat0 ^ 2 / 2,
"Poisson" = (theta - yhat0 + yhat0 * (log(yhat0) -
log(theta))),
"Exponential" = (log(yhat0) - 1 + 1 / yhat0),
"Bernoulli" = (yhat0 * log(yhat0 / theta) +
(1 - yhat0) * log((1 - yhat0) / (1 -
theta))),
"Volatility" = (log(yhat0) - 1 + 1 / yhat0) / 2,
"Variance" = shape / 2 * (log(yhat0) - 1 + 1 /
yhat0),
"NCchi" = kldistnorm1(theta, yhat0, shape)
)
exceedence[, k] <- .Fortran(C_exceedm,
as.double(KLdist1),
as.integer(length(KLdist1)),
as.double(z / ni),
as.integer(nz),
exprob = double(nz),
as.integer(mask)
)$exprob
contour(
z,
0:k,
cbind(exceedence0, exceedence[, 1:k]),
levels = levels,
ylab = "step",
xlab = "z",
main = paste(family, length(dy), "-dim. ladj=", ladjust, " Exceed. Prob.")
)
contour(
z,
0:k,
cbind(exceedence0, exceedencena[, 1:k]),
levels = levels,
ylab = "step",
xlab = "z",
add = TRUE,
col = 2,
lty = 3
)
yaxp <- par("yaxp")
at <-
unique(as.integer(seq(yaxp[1], yaxp[2], length = yaxp[3] + 1)))
at <- at[at > 0 & at <= k]
axis(4, at = at, labels = as.character(signif(h[at], 3)))
mtext("bandwidth", 4, 1.8)
if (max(total) > 0) {
cat("Progress:", signif(total[k], 2) * 100, "% . ", sep = "")
}
t2 <- Sys.time()
tpar <-
if (family %in% c("Bernoulli", "Poisson"))
paste("theta=", theta)
else
""
cat(
family,
"(dim:",
length(dy),
tpar,
") ni=",
ni,
" mean(bi/ni)=",
mean(bi[mask])/ni,
" Time: Step",
format(signif(difftime(t2, t1), 3)),
"Total",
format(signif(difftime(t2, t0), 3)),
"\n"
)
k <- k + 1
gc()
}
if (family %in% c("Bernoulli", "Poisson"))
y <- y0
list(
h = h,
z = z,
prob = exceedence,
probna = exceedencena,
y = if (verbose)
y
else
NULL ,
theta = if (verbose)
yhat
else
NULL,
levels = levels,
family = family
)
}
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.