Nothing
#############################################################################
##
## Setup for generator.
##
#############################################################################
## We store all data about the hat and squeeze in particular intervals
## in a matrix where each interval is stored in a column.
## We need the following parameter for each interval:
iv.struct <-
c(
## Data for hat and squeeze:
"x", ## left boundary of interval
"c", ## parameter for transformation
"ht.a", ## intercept for hat in transformed scale
"ht.b", ## slope for hat in transformed scale
"ht.y", ## anchor point of linear function
"sq.a", ## intercept for squeeze in transformed scale
"sq.b", ## slope for squeeze in transformed scale
"sq.y", ## anchor point of linear function
"A.ht", ## area below hat
"A.sq", ## area below squeeze
## Data required for setup only:
"type", ## type of interval
"Tfx", ## log-density at left boundary
"dTfx", ## derivative of log-density at left boundary
"d2Tfx", ## 2nd derivative of log-density at left boundary
"next" ## poor man's linked list:
## index of next interval to the right ('Inf' for last interval)
)
## Names of list elements for generator:
gen.struct <-
c(
"ivs", ## data for hat and squeeze
"lpdf", ## log-density
"A.ht.tot", ## total area below hat
"A.sq.tot", ## total area below squeeze
"env", ## environment for evaluating log-density
"iniv", ## initial intervals (for print.Tinflex)
"Acum", ## cumulated areas
"gt", ## guide table
"haved2Tf" ## whether 2nd derivative is used (for print.Tinflex)
)
## --------------------------------------------------------------------------
Tinflex.setup <- function(lpdf, dlpdf, d2lpdf=NULL, ib, cT=0, rho=1.1, max.intervals=1001) {
## ----------------------------------------------------------------------
## Setup: compute hat and squeeze for density.
## ----------------------------------------------------------------------
## lpdf ... log-density
## dlpdf ... derivative of log-density
## d2lpdf ... 2nd derivative of log-density
## ib ... interval boundaries of decomposition
## cT ... parameter for transformation (global or for each interval)
## rho ... performance parameter: requested upper bound for ratio
## between area below hat and area below squeeze
## max.intervals ... maximal numbers of intervals
## ----------------------------------------------------------------------
## return: S3 object of class "Tinflex" that contains
## ivs ... parameters for hat and squeeze for each interval
## the last interval just stores the right boundary of domain
## lpdf ... log-density
## A.ht.tot ... total area below hat
## A.sq.tot ... total area below squeeze
## env ... environment for evaluating log-density in C
## iniv ... string that describes initial intervals
## Acum ... cumulated areas
## gt ... guide table
## haved2Tf ... whether 2nd derivative is used (for print.Tinflex)
## ----------------------------------------------------------------------
## Check arguments.
if (!isTRUE(rho >= 1.0001)) {
warning ("argument 'rho' too small or invalid, using 1.0001")
rho <- 1.0001
}
if (!isTRUE(max.intervals > 51)) {
warning ("argument 'max.intervals' too small or invalid, using default")
max.intervals <- 1001
}
if (missing(lpdf) || !is.function(lpdf)) {
stop ("parameter 'lpdf' missing or invalid")
}
if (missing(dlpdf) || !is.function(dlpdf)) {
stop ("parameter 'dlpdf' missing or invalid")
}
if (! (is.function(d2lpdf) || is.null(d2lpdf) ) ) {
stop ("parameter 'd2lpdf' is invalid")
}
if (missing(ib)) {
stop ("argument 'ib' is missing, with no default")
}
if (length(ib)<2 || length(ib) > max.intervals / 2) {
stop ("argument 'ib' invalid")
}
## The boundaries must be sorted.
ib <- sort(ib)
## Check parameters for transformation.
if (! ( is.numeric(cT) && (length(cT)==1 || length(cT)==length(ib)-1 )) ) {
stop ("argument 'cT' invalid: its length must equal either 1 or number of intervals")
}
if ((! is.finite(ib[1]) && ! isTRUE(cT[1] > -1.)) ||
(! is.finite(ib[length(ib)]) && ! isTRUE(cT[length(cT)] > -1.)) ) {
stop ("(first and last) entry of argument 'cT' must be greater than -1 for unbounded domains")
}
## ........................................................................
if (is.null(d2lpdf)) {
generator <- setup.no2ndD (lpdf, dlpdf, ib, cT, rho, max.intervals)
}
else {
generator <- setup.with2ndD(lpdf, dlpdf, d2lpdf, ib, cT, rho, max.intervals)
}
## return
generator
} ## -- end of Tinflex.setup() -- ##
## --------------------------------------------------------------------------
setup.with2ndD <- function(lpdf, dlpdf, d2lpdf, ib, cT, rho, max.intervals) {
## ----------------------------------------------------------------------
## Setup: compute hat and squeeze for density.
## Version that makes use of 2nd derivative.
## ----------------------------------------------------------------------
## lpdf ... log-density
## dlpdf ... derivative of log-density
## d2lpdf ... 2nd derivative of log-density
## ib ... interval boundaries of decomposition
## cT ... parameter for transformation (global or for each interval)
## rho ... performance parameter: requested upper bound for ratio
## between area below hat and area below squeeze
## max.intervals ... maximal numbers of intervals
## ----------------------------------------------------------------------
## return: S3 object of class "Tinflex" that contains
## ivs ... parameters for hat and squeeze for each interval
## the last interval just stores the right boundary of domain
## lpdf ... log-density
## A.ht.tot ... total area below hat
## A.sq.tot ... total area below squeeze
## env ... environment for evaluating log-density in C
## iniv ... string that describes initial intervals
## Acum ... cumulated areas
## gt ... guide table
## ----------------------------------------------------------------------
## Store some of the arguments for print.Tinflex().
iniv <- paste(paste(" initial intervals ="),
paste(ib,collapse=" | "),
paste("\n cT ="),
paste(cT,collapse=" | "),
paste("\n rho =",rho,"\n"))
## Get parameter 'cT' for transformation.
if (length(cT)==1) {
cT <- rep(cT[1], length(ib))
}
else { ## Just duplicate last element.
cT <- append(cT, cT[length(cT)])
}
## Create a table for hat and squeeze.
params <- matrix(nrow=length(iv.struct), ncol=max.intervals)
rownames(params) <- iv.struct
## We need transformed densities and their derivatives on
## both boundaries of each interval.
## However, the boundary to the right of one interval is the
## boundary to the left of the consecutive interval.
## If their parameters 'c' coincide we only need to store the
## transformed density for the boundary to the left.
## If these parameters differ we insert an interval of length 0
## in order to store the transformed density for the boundary
## to the right.
n.ivs <- 0
for (i in 1:length(ib)) {
n.ivs <- n.ivs+1
params["x", n.ivs] <- ib[i]
params["c", n.ivs] <- cT[i]
params[c("Tfx","dTfx","d2Tfx"), n.ivs] <- Tfdd(lpdf, dlpdf, d2lpdf, cT[i], ib[i])
params["next", n.ivs] <- n.ivs+1
## Parameter 'cT' changes for next interval.
## Thus we insert an interval of length 0.
if (isTRUE(i < length(ib)) && !identical(cT[i], cT[i+1])) {
n.ivs <- n.ivs+1
params["x", n.ivs] <- ib[i+1]
params["c", n.ivs] <- cT[i]
params[c("Tfx","dTfx","d2Tfx"), n.ivs] <- Tfdd(lpdf, dlpdf, d2lpdf, cT[i], ib[i+1])
params["next", n.ivs] <- n.ivs+1
}
}
## Update number of intervals (after creating initial intervals).
n.ivs <- n.ivs-1
## Terminate linked list.
params["next", n.ivs+1] <- Inf
## NOTE: params[,n.ivs+1] now stores just the right boundary of the domain.
## Compute parameters for hat and squeeze for initial intervals.
for (i in 1:n.ivs) {
params["type",i] <- type.with2ndD(left=params[,i], right=params[,i+1])
params[,i] <- hat.iv(left=params[,i], right=params[,i+1])
## DEBUG:
## print(params[,i])
}
## DEBUG:
## print(params[,n.ivs+1])
## Compute total areas for initial hat and squeeze.
A.ht.tot <- sum(params["A.ht", 1:n.ivs])
A.sq.tot <- sum(params["A.sq", 1:n.ivs])
## We have to split intervals where the area
## between hat and squueze is too large.
while (! isTRUE(A.ht.tot / A.sq.tot <= rho)) {
## Compute average area.
threshold <- 0.99 * (A.ht.tot - A.sq.tot) / n.ivs
## Find all intervals where area between hat and squeeze is too large.
toolarge <- which(params["A.ht", 1:(n.ivs+1)] - params["A.sq", 1:(n.ivs+1)] >= threshold)
if (length(toolarge) <= 0) break;
if (n.ivs + length(toolarge) >= max.intervals) {
warning ("maximum number of intervals exceeded")
break;
}
## Split all intervals.
for (i in toolarge) {
## Get boundary points of interval.
left <- params[,i]
right <- params[,params["next",i]]
## Splitting point for interval (use "arc mean").
p <- arc.mean(left["x"],right["x"])
## Increment number of intervals
n.ivs <- n.ivs+1
## Create new interval on r.h.s.
params["x", n.ivs+1] <- p
params["c", n.ivs+1] <- left["c"]
params[c("Tfx","dTfx","d2Tfx"), n.ivs+1] <- Tfdd(lpdf, dlpdf, d2lpdf, left["c"], p)
params["type",n.ivs+1] <- type.with2ndD(left=params[,n.ivs+1], right=right)
params[,n.ivs+1] <- hat.iv(left=params[,n.ivs+1], right)
params["next",n.ivs+1] <- left["next"]
## Update interval on l.h.s.
params["type",i] <- type.with2ndD(left=params[,i], right=params[,n.ivs+1])
params[,i] <- hat.iv(left=params[,i], params[,n.ivs+1])
params["next",i] <- n.ivs+1
}
## Update total areas.
A.ht.tot <- sum(params["A.ht",], na.rm=TRUE)
A.sq.tot <- sum(params["A.sq",], na.rm=TRUE)
}
## Check result.
if (! is.finite(A.ht.tot)) {
stop ("Cannot create hat function. A_hat is not finite: A_hat=", A.ht.tot)
}
if (! isTRUE(A.ht.tot > 0.)) {
stop ("Cannot create hat function. A_hat is not > 0: A_hat=", A.ht.tot)
}
if (! is.finite(A.sq.tot)) {
stop ("Cannot create squeeze. A_squeeze is not finite: A_squeeze=", A.sq.tot)
}
if (! isTRUE(A.sq.tot >= 0.)) {
stop ("Cannot create squeeze. A_squeeze is not >= 0: A_squeeze=", A.sq.tot)
}
if (isTRUE(A.ht.tot < A.sq.tot))
stop ("Invalid input: A_hat < A_squeeze!")
if (!isTRUE(A.ht.tot / A.sq.tot <= rho))
warning ("ratio A_hat / A_squeeze = ", A.ht.tot / A.sq.tot,
" is larger than rho = ",rho)
## Truncate working array.
params <- params[,1:(n.ivs+1)]
## Sort array (so that the linked list becomes obsolete).
params <- params[,order(params["x",])]
## Compute guide table.
## For R version use
## guidetable <- make.guidetable(params)
## Acum <- guidetable$Acum
## gt <- guidetable$gt
## C version:
Acum <- numeric(n.ivs)
gt <- integer(n.ivs)
A.ht.tot <- .Call(C_Tinflex_RC_make_guide_table, params, Acum, gt)
## Create S3 class that contains generator.
generator <- list(
ivs=params, ## data for hat and squeeze
lpdf=lpdf, ## log-density
A.ht.tot=A.ht.tot, ## total area below hat
A.sq.tot=A.sq.tot, ## total area below hat
env=parent.frame(2), ## environment for evaluating log-density
iniv=iniv, ## initial intervals (for print.Tinflex)
Acum=Acum, ## cumulated areas
gt=gt, ## guide table
haved2Tf=TRUE ## whether 2nd derivative is used (for print.Tinflex)
)
class(generator) <- "Tinflex"
## Return generator object.
generator
} ## -- end of setup.with2ndD() -- ##
## --------------------------------------------------------------------------
setup.no2ndD <- function(lpdf, dlpdf, ib, cT=0, rho=1.1, max.intervals=1001) {
## ----------------------------------------------------------------------
## Setup: compute hat and squeeze for density.
## Version without 2nd derivative of transformed density.
## ----------------------------------------------------------------------
## lpdf ... log-density
## dlpdf ... derivative of log-density
## ib ... interval boundaries of decomposition
## cT ... parameter for transformation (global or for each interval)
## rho ... performance parameter: requested upper bound for ratio
## between area below hat and area below squeeze
## max.intervals ... maximal numbers of intervals
## ----------------------------------------------------------------------
## return: S3 object of class "Tinflex" that contains
## ivs ... parameters for hat and squeeze for each interval
## the last interval just stores the right boundary of domain
## lpdf ... log-density
## A.ht.tot ... total area below hat
## A.sq.tot ... total area below squeeze
## env ... environment for evaluating log-density in C
## iniv ... string that describes initial intervals
## Acum ... cumulated areas
## gt ... guide table
## ----------------------------------------------------------------------
## ......................................................................
## Store some of the arguments for print.Tinflex().
iniv <- paste(paste(" initial intervals ="),
paste(ib,collapse=" | "),
paste("\n cT ="),
paste(cT,collapse=" | "),
paste("\n rho =",rho,"\n"))
## Get parameter 'cT' for transformation.
if (length(cT)==1) {
cT <- rep(cT[1], length(ib))
}
else { ## Just duplicate last element.
cT <- append(cT, cT[length(cT)])
}
## Create a table for hat and squeeze.
params <- matrix(nrow=length(iv.struct), ncol=max.intervals)
rownames(params) <- iv.struct
## We need transformed densities and their derivatives on
## both boundaries of each interval.
## However, the boundary to the right of one interval is the
## boundary to the left of the consecutive interval.
## If their parameters 'c' coincide we only need to store the
## transformed density for the boundary to the left.
## If these parameters differ we insert an interval of length 0
## in order to store the transformed density for the boundary
## to the right.
n.ivs <- 0
for (i in 1:length(ib)) {
n.ivs <- n.ivs+1
params["x", n.ivs] <- ib[i]
params["c", n.ivs] <- cT[i]
params[c("Tfx","dTfx","d2Tfx"), n.ivs] <- Tfdd(lpdf, dlpdf, NULL, cT[i], ib[i])
params["next", n.ivs] <- n.ivs+1
## Parameter 'cT' changes for next interval.
## Thus we insert an interval of length 0.
if (isTRUE(i < length(ib)) && !identical(cT[i], cT[i+1])) {
n.ivs <- n.ivs+1
params["x", n.ivs] <- ib[i+1]
params["c", n.ivs] <- cT[i]
params[c("Tfx","dTfx","d2Tfx"), n.ivs] <- Tfdd(lpdf, dlpdf, NULL, cT[i], ib[i+1])
params["next", n.ivs] <- n.ivs+1
}
}
## Update number of intervals (after creating initial intervals).
n.ivs <- n.ivs-1
## Terminate linked list.
params["next", n.ivs+1] <- Inf
## NOTE: params[,n.ivs+1] now stores just the right boundary of the domain.
## Determine type of initial intervals.
for (i in 1:n.ivs) {
params["type",i] <- type <-
type.no2ndD.init(left=params[,i], right=params[,i+1], lpdf, dlpdf)
if (type == IIb_IVa__IIa_IVa || ## Case (3.3.3)
type == IIIa_IVb__IIIb_IVb ) { ## Case (4.3.3)
## we have to split the interval
## Splitting point for interval (use "arc mean").
left <- params[,i]
p <- arc.mean(left["x"],params["x",i+1])
## Increment number of intervals
n.ivs <- n.ivs+1
## Create new interval for r.h.s.
params["x", n.ivs+1] <- p
params["c", n.ivs+1] <- left["c"]
params[c("Tfx","dTfx","d2Tfx"), n.ivs+1] <- Tfdd(lpdf, dlpdf, NULL, left["c"], p)
params["next",n.ivs+1] <- left["next"]
params["type",n.ivs+1] <- if (type == IIb_IVa__IIa_IVa) { IIa_IVa } else { IIIb_IVb }
## Update interval on l.h.s.
params["next",i] <- n.ivs+1
params["type",i] <- if (type == IIb_IVa__IIa_IVa) { IIb_IVa } else { IIIa_IVb }
}
}
## Compute parameters for hat and squeeze for initial intervals.
for (i in 1:(n.ivs+1)) {
if (is.finite(params["next",i])) {
params[,i] <- hat.iv(left=params[,i], right=params[,params["next",i]])
}
}
## DEBUG:
## print(params[,1:(n.ivs+1)]) ## DEBUG
## Compute total areas for initial hat and squeeze.
A.ht.tot <- sum(params["A.ht", 1:n.ivs], na.rm=TRUE)
A.sq.tot <- sum(params["A.sq", 1:n.ivs], na.rm=TRUE)
## We have to split intervals where the area
## between hat and squueze is too large.
while (! isTRUE(A.ht.tot / A.sq.tot <= rho)) {
## Compute average area.
threshold <- 0.99 * (A.ht.tot - A.sq.tot) / n.ivs
## Find all intervals where area between hat and squeeze is too large.
toolarge <- which(params["A.ht", 1:(n.ivs+1)] - params["A.sq", 1:(n.ivs+1)] >= threshold)
if (length(toolarge) <= 0) break;
if (n.ivs + length(toolarge) >= max.intervals) {
warning ("maximum number of intervals exceeded")
break;
}
## Split all intervals.
for (i in toolarge) {
## Get boundary points of interval.
left <- params[,i]
right <- params[,params["next",i]]
## Splitting point and types of subinervals
split.iv <- type.no2ndD.split(left=left, right=right, lpdf, dlpdf)
p <- split.iv[1]
Tfp <- split.iv[2]
dTfp <- split.iv[3]
typel <- split.iv[4]
typer <- split.iv[5]
## Increment number of intervals
n.ivs <- n.ivs+1
## Create new interval on r.h.s.
params["x", n.ivs+1] <- p
params["c", n.ivs+1] <- left["c"]
params[c("Tfx","dTfx","d2Tfx"), n.ivs+1] <- c(Tfp, dTfp, NA_real_)
## or equivalently:
## params[c("Tfx","dTfx","d2Tfx"), n.ivs+1] <- Tfdd(lpdf, dlpdf, NULL, left["c"], p)
params["type",n.ivs+1] <- typer
params[ ,n.ivs+1] <- hat.iv(left=params[,n.ivs+1], right=right)
params["next",n.ivs+1] <- left["next"]
## Update interval on l.h.s.
params["type",i] <- typel
params[ ,i] <- hat.iv(left=params[,i], right=params[,n.ivs+1])
params["next",i] <- n.ivs+1
}
## Update total areas.
A.ht.tot <- sum(params["A.ht",], na.rm=TRUE)
A.sq.tot <- sum(params["A.sq",], na.rm=TRUE)
}
## DEBUG:
## print(params[,1:(n.ivs+1)])
## Check result.
if (! is.finite(A.ht.tot)) {
stop ("Cannot create hat function. A_hat is not finite: A_hat=", A.ht.tot)
}
if (! isTRUE(A.ht.tot > 0.)) {
stop ("Cannot create hat function. A_hat is not > 0: A_hat=", A.ht.tot)
}
if (! is.finite(A.sq.tot)) {
stop ("Cannot create squeeze. A_squeeze is not finite: A_squeeze=", A.sq.tot)
}
if (! isTRUE(A.sq.tot >= 0.)) {
stop ("Cannot create squeeze. A_squeeze is not >= 0: A_squeeze=", A.sq.tot)
}
if (isTRUE(A.ht.tot < A.sq.tot))
stop ("Invalid input: A_hat < A_squeeze!")
if (!isTRUE(A.ht.tot / A.sq.tot <= rho))
warning ("ratio A_hat / A_squeeze = ", A.ht.tot / A.sq.tot,
" is larger than rho = ",rho)
## Truncate working array.
params <- params[,1:(n.ivs+1)]
## Sort array (so that the linked list becomes obsolete).
params <- params[,order(params["x",])]
## Compute guide table.
## For R version use
## guidetable <- make.guidetable(params)
## Acum <- guidetable$Acum
## gt <- guidetable$gt
## C version:
Acum <- numeric(n.ivs)
gt <- integer(n.ivs)
A.ht.tot <- .Call(C_Tinflex_RC_make_guide_table, params, Acum, gt)
## Create S3 class that contains generator.
generator <- list(
ivs=params, ## data for hat and squeeze
lpdf=lpdf, ## log-density
A.ht.tot=A.ht.tot, ## total area below hat
A.sq.tot=A.sq.tot, ## total area below hat
env=parent.frame(2), ## environment for evaluating log-density
iniv=iniv, ## initial intervals (for print.Tinflex)
Acum=Acum, ## cumulated areas
gt=gt, ## guide table
haved2Tf=FALSE ## whether 2nd derivative is used (for print.Tinflex)
)
class(generator) <- "Tinflex"
## Return generator object.
return (generator)
} ## -- end of Tinflex.setup.no2ndD() -- ##
## --------------------------------------------------------------------------
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.