Nothing
#' Sample covariance function
#'
#' This function provides an empirical covariance function for data
#' in a \code{\link{CMBDataFrame}} or data.frame. It assumes that data are from a stationary spherical
#' random field and the covariance depends only on a geodesic distance between locations.
#' Output is a binned covariance.
#'
#' @param cmbdf is a \code{\link{CMBDataFrame}} or \code{\link{data.frame}}
#' @param num.bins specifies the number of bins
#' @param sample.size optionally specify the size of a simple random
#' sample to take before calculating covariance. This may be useful if
#' the full covariance computation is too slow.
#' @param max.dist an optional number between 0 and pi specifying the
#' maximum geodesic distance to use for calculating covariance. Only
#' used if \code{breaks} are unspecified.
#' @param breaks optionally specify the breaks manually using a
#' vector giving the break points between cells. This vector
#' has length \code{num.bins} since the last break point is taken
#' as \code{max.dist}. If \code{equiareal = TRUE} then
#' these breaks should be \eqn{cos(r_i)} where \eqn{r_i} are radii.
#' If \code{equiareal = FALSE} then these breaks should be \eqn{r_i}.
#' @param equiareal if TRUE then the bins have equal spherical
#' area. If false then the bins have equal annular widths.
#' Default is TRUE.
#' @param calc.max.dist if TRUE then the \code{max.dist} will be
#' calculated from the locations in cmbdf. Otherwise
#' either \code{max.dist} must be specified or max.dist will
#' default to pi.
#'
#' @return
#'
#' An object of the class CMBCovariance that is a modification of \code{\link[geoR]{variog}}
#' from the package \strong{geoR} with variogram values replaced by covariances.
#'
#' The attribute "breaks" contains the break points used to create bins.
#' The result has \code{num.bins + 1} values since the first value, the sample
#' variance, is not counted as a bin.
#'
#'
#' \describe{
#' \item{u}{a vector with distances.}
#' \item{v}{a vector with estimated covariance values at distances given in u.}
#' \item{n}{number of pairs in each bin}
#' \item{sd}{standard deviation of the values in each bin}
#' \item{bins.lim}{ limits defining the interval spanned by each bin}
#' \item{ind.bin}{a logical vector indicating whether the number of pairs
#' in each bin is greater or equal to the value in the argument pairs.min}
#' \item{var.mark}{variance of the data}
#' \item{beta.ols}{parameters of the mean part of the model fitted by ordinary
#' least squares}
#' \item{output.type}{echoes the option argument}
#' \item{max.dist}{maximum distance between pairs allowed in the covariance calculations}
#' \item{n.data}{number of data}
#' \item{direction}{direction for which the covariance was computed}
#' \item{call}{the function call}
#' }
#'
#'
#'@references \strong{geoR} package, \code{\link[geoR]{variog}}, \code{\link{variogramCMB}}, \code{\link{corrCMB}}
#'
#'@examples
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # cmbdf <- sampleCMB(df, sample.size = 100000)
#' # Cov <- covCMB(cmbdf, max.dist = 0.03, num.bins = 10)
#' # Cov
#'
#' @export
covCMB <- function(cmbdf,
num.bins = 10,
sample.size,
max.dist = pi,
breaks,
equiareal = TRUE,
calc.max.dist = FALSE)
{
if (!is.CMBDataFrame(cmbdf)) {
stop("cmbdf must be a CMBDataFrame")
}
if (!missing(sample.size)) {
cmbdf <- rcosmo::sampleCMB(cmbdf, sample.size = sample.size)
}
if (is.null(coords(cmbdf)) || coords(cmbdf) != "cartesian") {
coords(cmbdf) <- "cartesian"
}
if (!all(c("x", "y", "z", "I") %in% names(cmbdf)))
{
stop("cmbdf must have columns named 'x', 'y', 'z', 'I' in that order")
}
if (missing(breaks))
{
if (calc.max.dist)
{
max.dist <- maxDist_internal2(cmbdf)
}
if (equiareal)
{
breaks <-
seq(cos(max.dist), 1, length.out = num.bins + 1)[-(num.bins + 1)]
}
else
{
breaks <- cos(rev(seq(0, max.dist, length.out = num.bins + 1)[-1]))
}
}
covs <- covCMB_internal_var(cmbdf[, c("x", "y", "z", "I")], breaks)
# Reverse order since cosine is decreasing
covs[2:nrow(covs), 1] <- rev(covs[2:nrow(covs), 1])
covs[2:nrow(covs), 2] <- rev(covs[2:nrow(covs), 2])
covs[2:nrow(covs), 3] <- rev(covs[2:nrow(covs), 3])
v <- c(0, rev(acos(breaks)))
# Drop the throw-away bin (distances greater than max.dist)
covs <- covs[-nrow(covs), ]
# Find the bin centers (break_{i+1} - break_i)/2 with break_0 = 0.
centers <-
c(0, v[1:(length(v) - 1)] + (v[2:length(v)] - v[1:(length(v) - 1)]) / 2)
pairs.min <- 2
n <- as.integer(c(covs[, 2][1], covs[, 2][-1] / 2))
indp <- (n >= pairs.min)
sd1 <- sqrt(covs[,3])
result <- list(
u = centers,
v = covs[, 1],
n = n,
sd1 = sd1,
bins.lim = v,
ind.bin = indp
)
call.fc <- match.call()
option <- "bin"
estimator.type <- "classical"
trend <- "cte"
data.var <- stats::var(cmbdf$I)
n.data <- length(cmbdf$I)
##
beta.ols <- mean(cmbdf$I)
umax <- max(centers[centers < max.dist])
result <- c(
result,
list(
var.mark = data.var,
beta.ols = beta.ols,
output.type = option,
max.dist = max.dist,
estimator.type = estimator.type,
n.data = n.data,
lambda = 1,
trend = trend,
pairs.min = pairs.min
)
)
result$nugget.tolerance <- 1e-12
result$direction <- "omnidirectional"
result$tolerance <- "none"
result$uvec <- centers
result$call <- call.fc
oldClass(result) <- "CMBCovariance"
return(result)
}
#' Sample correlation function
#'
#' This function provides an empirical correlation function for data
#' in a \code{\link{CMBDataFrame}} or data.frame. It assumes that data are from a stationary spherical
#' random field and the correlation depends only on a geodesic distance between locations.
#' Output is a binned correlation.
#'
#' @param cmbdf is a \code{\link{CMBDataFrame}} or \code{\link{data.frame}}
#' @param num.bins specifies the number of bins
#' @param sample.size optionally specify the size of a simple random
#' sample to take before calculating correlation. This may be useful if
#' the full correlation computation is too slow.
#' @param max.dist an optional number between 0 and pi specifying the
#' maximum geodesic distance to use for calculating correlation. Only
#' used if \code{breaks} are unspecified.
#' @param breaks optionally specify the breaks manually using a
#' vector giving the break points between cells. This vector
#' has length \code{num.bins} since the last break point is taken
#' as \code{max.dist}. If \code{equiareal = TRUE} then
#' these breaks should be \eqn{cos(r_i)} where \eqn{r_i} are radii.
#' If \code{equiareal = FALSE} then these breaks should be \eqn{r_i}.
#' @param equiareal if TRUE then the bins have equal spherical
#' area. If false then the bins have equal annular widths.
#' Default is TRUE.
#' @param calc.max.dist if TRUE then the \code{max.dist} will be
#' calculated from the locations in cmbdf. Otherwise
#' either \code{max.dist} must be specified or max.dist will
#' default to pi.
#'
#' @return
#'#' An object of the class CMBCorrelation that is a modification of \code{\link[geoR]{variog}}
#' from the package \strong{geoR} with variogram values replaced by correlation.
#'
#' The attribute "breaks" contains the break points used to create bins.
#' The result has \code{num.bins + 1} values since the first value at distance 0 is not
#' counted as a bin.
#'
#' \describe{
#' \item{u}{a vector with distances.}
#' \item{v}{a vector with estimated correlation values at distances given in u.}
#' \item{n}{number of pairs in each bin}
#' \item{sd}{standard deviation of the values in each bin}
#' \item{bins.lim}{ limits defining the interval spanned by each bin}
#' \item{ind.bin}{a logical vector indicating whether the number of pairs
#' in each bin is greater or equal to the value in the argument pairs.min}
#' \item{var.mark}{variance of the data}
#' \item{beta.ols}{parameters of the mean part of the model fitted by ordinary
#' least squares}
#' \item{output.type}{echoes the option argument}
#' \item{max.dist}{maximum distance between pairs allowed in the correlation calculations}
#' \item{n.data}{number of data}
#' \item{direction}{direction for which the correlation was computed}
#' \item{call}{the function call}
#' }
#'
#'
#'@references \strong{geoR} package,\code{\link[geoR]{variog}}, \code{\link{variogramCMB}}, \code{\link{covCMB}}
#'
#'@examples
#'
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # cmbdf <- sampleCMB(df, sample.size = 100000)
#' # corcmb <- corrCMB(cmbdf, max.dist = 0.03, num.bins = 10, sample.size=1000)
#' # corcmb
#'
#' @export
corrCMB <- function(cmbdf,
num.bins = 10,
sample.size,
max.dist = pi,
breaks,
equiareal = TRUE,
calc.max.dist = FALSE) {
corrCMB<- covCMB(cmbdf, num.bins,
sample.size,
max.dist,
breaks,
equiareal ,
calc.max.dist)
corrCMB$v <- corrCMB$v/corrCMB$v[1]
oldClass(corrCMB) <- "CMBCorrelation"
return(corrCMB)
}
#' Sample variogram
#'
#' This function provides an empirical variogram for data in a
#' \code{\link{CMBDataFrame}} or data.frame. It assumes that data are from a
#' stationary spherical random field and the covariance depends only on a
#' geodesic distance between locations. Output is a binned variogram.
#'
#'
#' @param cmbdf is a \code{\link{CMBDataFrame}} or \code{\link{data.frame}}
#' @param num.bins specifies the number of bins
#' @param sample.size optionally specify the size of a simple random
#' sample to take before calculating variogram. This may be useful if
#' the full covariance computation is too slow.
#' @param max.dist an optional number between 0 and pi specifying the
#' maximum geodesic distance to use for calculating covariance. Only
#' used if \code{breaks} are unspecified.
#' @param breaks optionally specify the breaks manually using a
#' vector giving the break points between cells. This vector
#' has length \code{num.bins} since the last break point is taken
#' as \code{max.dist}. If \code{equiareal = TRUE} then
#' these breaks should be \eqn{cos(r_i)} where \eqn{r_i} are radii.
#' If \code{equiareal = FALSE} then these breaks should be \eqn{r_i}.
#' @param equiareal if TRUE then the bins have equal spherical
#' area. If false then the bins have equal annular widths.
#' Default is TRUE.
#' @param calc.max.dist if TRUE then the \code{max.dist} will be
#' calculated from the locations in cmbdf. Otherwise
#' either \code{max.dist} must be specified or max.dist will
#' default to pi.
#'
#' @return
#' An object of class \code{\link[geoR]{variog}} specified in the package \strong{geoR}.
#'
#' The attribute "breaks" contains the break points used to create bins.
#' The result has \code{num.bins + 1} values since the first value at distance 0 is not
#' counted as a bin.
#'
#'\describe{
#' \item{u}{a vector with distances.}
#' \item{v}{a vector with estimated variogram values at distances given in u.}
#' \item{n}{number of pairs in each bin}
#' \item{sd}{standard deviation of the values in each bin}
#' \item{bins.lim}{ limits defining the interval spanned by each bin}
#' \item{ind.bin}{a logical vector indicating whether the number of pairs
#' in each bin is greater or equal to the value in the argument pairs.min}
#' \item{var.mark}{variance of the data}
#' \item{beta.ols}{parameters of the mean part of the model fitted by ordinary
#' least squares}
#' \item{output.type}{echoes the option argument}
#' \item{max.dist}{maximum distance between pairs allowed in the variogram calculations}
#' \item{n.data}{number of data}
#' \item{direction}{direction for which the variogram was computed}
#' \item{call}{the function call}
#' }
#'
#'
#'@references \strong{geoR} package, \code{\link[geoR]{variog}}, \code{\link{covCMB}}, \code{\link{corrCMB}}
#'
#'@examples
#'
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # cmbdf <- sampleCMB(df, sample.size = 100000)
#' # varcmb <- variogramCMB(cmbdf, max.dist = 0.1, num.bins = 30, sample.size=100)
#' # varcmb
#'
#' @export
variogramCMB <- function(cmbdf,
num.bins = 10,
sample.size,
max.dist = pi,
breaks,
equiareal = TRUE,
calc.max.dist = FALSE) {
varCMB<- covCMB(cmbdf, num.bins ,
sample.size,
max.dist ,
breaks,
equiareal ,
calc.max.dist )
varCMB$v <- varCMB$v[1]-varCMB$v
oldClass(varCMB) <- "variogram"
return(varCMB)
}
#'Plot sample variogram
#'
#'Plots sample (empirical) variogram. Uses \code{\link[geoR]{plot.variogram}} from
#'\strong{geoR} package.
#'
#'@param x An object of class variogram.
#'@param ... Extra arguments as in \code{\link[geoR]{plot.variogram}} passed to plot.default.
#'
#'
#'@return Produces a plot with the sample variogram.
#'
#'@references \strong{geoR} package, \code{\link{variogramCMB}}, \code{\link[geoR]{variog}},
#'\code{\link[geoR]{plot.variogram}}
#'
#'@examples
#'
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # cmbdf <- sampleCMB(df, sample.size = 100000)
#' # varcmb <- variogramCMB(cmbdf, max.dist = 0.1, num.bins = 30, sample.size=1000)
#' # plot(varcmb)
#'
#'@import geoR
#'
#'@name plot.variogram
#'
NULL
#'Plot sample CMBCovariance
#'
#'Plots sample (empirical) covariance function. Uses \code{\link[geoR]{plot.variogram}} from
#'\strong{geoR} package.
#'
#'@param x An object of class CMBCovariance.
#'@param ... Extra arguments as in \code{\link[geoR]{plot.variogram}} passed to plot.default.
#'
#'@return Produces a plot with the sample covariance function.
#'
#'#'@references \strong{geoR} package, \code{\link{covCMB}}, \code{\link[geoR]{variog}},
#'\code{\link[geoR]{plot.variogram}}
#'
#'@examples
#'
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # cmbdf <- sampleCMB(df, sample.size = 100000)
#' # Cov <- covCMB(cmbdf, max.dist = 0.03, num.bins = 10)
#' # plot(Cov)
#'
#'@import geoR
#'
#' @export
plot.CMBCovariance <- function (x, ...) {
x0 <- x
attributes(x0)$class <- "variogram"
graphics::plot(x0, ylab = "sample covariance", ...)
}
#'Plot sample CMBCorrelation
#'
#'Plots sample (empirical) correlation function. Uses \code{\link[geoR]{plot.variogram}} from
#'\strong{geoR} package.
#'
#'@param x An object of class CMBCorrelation.
#'@param ... Extra arguments as in \code{\link{plot.variogram}} passed to plot.default.
#'
#'@return Produces a plot with the sample correlation function.
#'
#'@references \strong{geoR} package, \code{\link{corrCMB}}, \code{\link[geoR]{variog}},
#'\code{\link{plot.variogram}}
#'
#'@import geoR
#'@examples
#'
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # cmbdf <- sampleCMB(df, sample.size = 100000)
#' # corcmb <- corrCMB(cmbdf, max.dist = 0.03, num.bins = 10, sample.size=1000)
#' # plot(corcmb)
#'
#' @export
plot.CMBCorrelation <- function (x, ...) {
x0 <- x
attributes(x0)$class <- "variogram"
graphics::plot(x0, ylab= "sample correlation", ...)
}
#' Covariance estimate via power spectra
#'
#'This function provides a covariance estimate
#'using the values of the estimated
#'power spectra.
#'
#'
#'@param PowerSpectra A data frame which first
#'column lists values of multipole
#'moments and the second column gives
#'the corresponding values of CMB power
#'spectra.
#'@param Ns A number of points in which
#'the covariance estimate is computed on
#'the interval [-1,1]
#'
#'
#'@return
#' A data frame which first column is 1-d grid starting at -1+1/Ns and
#' finishing at 1 with the step 2/Ns. The second column is the values of
#' estimated covariances on this grid.
#'
#'@references Formula (2.1) in Baran A., Terdik G. Power spectrum estimation
#' of spherical random fields based on covariances. Annales Mathematicae et
#' Informaticae 44 (2015) pp. 15–22.
#'
#' Power Spectra data are from Planck Legacy Archive
#' \url{http://pla.esac.esa.int/pla/#cosmology}
#'
#'
#'@examples
#'
#' ## Download the power spectrum first
#' # N <- 20000
#' # COM_PowerSpectra <- downloadCMBPS(link=1)
#' #
#' # Cov_est <- covPwSp(COM_PowerSpectra[,1:2], N)
#' # plot(Cov_est, type="l")
#'
#' ## Plot the covariance estimate as a function of angular distances
#' # plot(acos(Cov_est[,1]), Cov_est[,2], type ="l",
#' # xlab ="angular distance", ylab ="Estimated Covariance")
#'
#'@export
covPwSp <- function(PowerSpectra, Ns)
{
if (requireNamespace("gsl", quietly = TRUE)) {
Nl <- length(PowerSpectra[, 1])
Pls <- gsl::legendre_Pl_array(lmax = PowerSpectra[Nl, 1],
x = seq(-1 + 1 / Ns, 1, 2 / Ns))
Cov_est <- Cov_func(Pls, PowerSpectra[, 2], PowerSpectra[, 1])
df <- data.frame(t = seq(-1 + 1 / Ns, 1, 2 / Ns), Estimated_Cov = Cov_est)
return(df)
} else {
stop("Package \"gsl\" needed for this function. Please install it.")
}
}
# Helper function for covPwSp
Cov_func <- function(mat, Dfl , l) {
apply(mat[l + 1, ], function(col) {
sum(Dfl * (2 * l + 1) / (4 * pi * l * (l + 1)) * col)
}, MARGIN = 2)
}
#' Power spectra estimate via correlation
#'
#'This function provides an angular power spectra estimate
#'using the values of the sample correlations. The approach is based
#'on Lawson-Hanson algorithm for non-negative least squares.
#'
#'
#'@param corcmb An object of the class CMBCorrelation.
#'@param lmax A number of angular power spectra components to estimate
#'
#'
#'@return
#' A data frame which first column is 1-d grid of l values
#' from 0 to \code{lmax}. The second column is
#' estimated angular power spectra components on this grid.
#'
#'@references Formula (2.1) in Baran A., Terdik G. Power spectrum
#'estimation of spherical random fields based on covariances.
#'Annales Mathematicae et Informaticae 44 (2015) pp. 15–22.
#'
#'
#'@examples
#'
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # Corrf <- corrCMB(df, max.dist = 0.1, num.bins = 30,
#' # sample.size=10000)
#' # pw <- pwSpCorr(Corrf)
#'
#'@export
pwSpCorr <- function(corcmb, lmax=20*length(corcmb$u)){
if (requireNamespace("gsl", quietly = TRUE)) {
Pls <- gsl::legendre_Pl_array(lmax, x = cos(corcmb$u))
A <- as.matrix(t((2*(0:lmax)+1)/(4*pi)*Pls))
fit <- nnls::nnls(A,corcmb$v)
return(data.frame(L=(0:lmax),f_l=fit$x))}
else {
stop("Package \"gsl\" needed for this function. Please install it.")
}
}
#'Plot angular scatterplots and means
#'
#'For specified measurements from \code{\link{CMBDataFrame}} this function
#'produces scatterplots and binned means versus theta and phi angles.
#'
#'@param cmbdf A \code{\link{CMBDataFrame}} object.
#'
#'@param intensities The name of a column of \code{cmbdf},
#'containing measured values.
#'
#'@return
#' 2x2 plot. The first row shows scatterplots. The second row gives
#' barplots of the corresponding means computed over bins. The first column
#' corresponds to the values of theta and the second one is for psi.
#'
#'
#'@examples
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # df.sample <- CMBDataFrame(df, sample.size = 80000)
#' # win <- CMBWindow(theta = c(pi/4,pi/2,pi/2,pi/4), phi = c(0,0,pi/2,pi/2))
#' # cmbdf.win <- window(df.sample, new.window = win)
#' #
#' # intensities <- "I"
#' # plotAngDis(cmbdf.win, intensities)
#'
#'@export
plotAngDis <- function(cmbdf, intensities = "I")
{
coords(cmbdf) <- "spherical"
thetabreaks <-
cut(
cmbdf$theta,
breaks = graphics::hist(cmbdf$theta, plot = FALSE)$breaks,
right = FALSE
)
theta.mean <-
t(tapply(as.data.frame(cmbdf)[, intensities], thetabreaks, mean,
na.rm = TRUE))
phibreaks <-
cut(cmbdf$phi,
breaks = graphics::hist(cmbdf$phi, plot = FALSE)$breaks,
right = FALSE)
phi.mean <-
t(tapply(as.data.frame(cmbdf)[, intensities], phibreaks, mean,
na.rm = TRUE))
graphics::par(mfrow = c(2, 2), mar = c(1, 4, 1, 1) + 0.5)
graphics::plot(as.data.frame(cmbdf[, c("theta", intensities)]), type = "p", col = "red")
graphics::plot(
as.data.frame(cmbdf[, c("phi", intensities)]),
ylab = "",
type = "p",
col = "blue"
)
graphics::barplot(
theta.mean,
legend = rownames(theta.mean),
col = "red",
ylim = 1.1 * range(theta.mean),
xlab = " ",
ylab = "Mean Values",
xaxt = 'n'
)
graphics::title(xlab = "theta",
line = 0,
cex.lab = 1)
graphics::barplot(
phi.mean,
legend = rownames(phi.mean),
col = "blue",
ylim = 1.1 * range(phi.mean),
xlab = "",
ylab = "",
xaxt = 'n'
)
graphics::title(xlab = "phi",
line = 0,
cex.lab = 1)
}
#' Take a simple random sample from a \code{\link{CMBDataFrame}}
#'
#' This function returns a \code{\link{CMBDataFrame}} which size equals to sample.size,
#' whose rows comprise a simple random sample of the rows
#' from the input CMBDataFrame.
#'
#'@param cmbdf a \code{\link{CMBDataFrame}}.
#'@param sample.size the desired sample size.
#'
#'@return
#' A \code{\link{CMBDataFrame}} which size equals to sample.size,
#' whose rows comprise a simple random sample of the rows
#' from the input CMBDataFrame.
#'
#'@examples
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # plot(sampleCMB(df, sample.size = 800000))
#'
#' df <- CMBDataFrame(nside = 16, I = rnorm(12 * 16 ^ 2), ordering = "nested")
#' df.sample <- sampleCMB(df, sample.size = 100)
#' df.sample
#'
#'@export
sampleCMB <- function(cmbdf, sample.size)
{
if ( !is.CMBDataFrame(cmbdf) )
{
stop("Argument must be a CMBDataFrame")
}
srows <- sample(1:nrow(cmbdf), sample.size)
cmbdf[srows, ]
}
#' First Minkowski functional
#'
#' This function returns an area of the spherical region
#' where measured values
#' are above of the specified threshold level \eqn{alpha}.
#'
#'@param cmbdf A \code{\link{CMBDataFrame}}.
#'@param alpha A numeric threshold level.
#'@param intensities A \code{\link{CMBDataFrame}} column with measured values.
#'@return
#' The area of the exceedance region
#'
#' @references Leonenko N., Olenko A. (2014) Sojourn measures of Student
#' and Fisher-Snedecor random fields. Bernoulli, 20:1454-1483.
#'
#'@examples
#'
#' n <- 64
#' cmbdf <- CMBDataFrame(nside=n, I = rnorm(12*n^2),
#' coords = "cartesian",
#' ordering = "nested")
#' fmf(cmbdf, 0, 4)
#' fmf(cmbdf, 2, 4)
#'
#' win <- CMBWindow(theta = c(0,pi/2,pi/2), phi = c(0,0,pi/2))
#' cmbdf.win <- window(cmbdf, new.window = win)
#' fmf(cmbdf.win, 0, 4)
#'
#'@export
fmf <- function(cmbdf, alpha, intensities = "I")
{
if ( !is.CMBDataFrame(cmbdf) )
{
stop("Argument must be a CMBDataFrame")
}
pixelArea(cmbdf)*sum(cmbdf[,intensities, drop = TRUE] > alpha)
}
#' Threshold exceedance probability
#'
#' This function returns an estimated exceedance probability for the specified
#' \code{\link{CMBDataFrame}} column \code{intensities}, threshold level
#' \eqn{alpha} and \code{\link{CMBWindow}} region.
#'
#'@param cmbdf A \code{\link{CMBDataFrame}}.
#'@param win A \code{\link{CMBWindow}}
#'@param alpha A numeric threshold level.
#'@param intensities The name of the column in \code{cmbdf}
#'that contains the measured values.
#'
#'@return
#'
#'Estimated threshold exceedance probability
#'
#'@examples
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # cmbdf <- sampleCMB(df, sample.size = 1000)
#'
#' # intensities <- "I"
#' # alpha <- mean(cmbdf[,intensities, drop = TRUE])
#' # alpha
#'
#' # win1 <- CMBWindow(theta = c(0,pi/2,pi/2), phi = c(0,0,pi/2))
#' # exprob(cmbdf, win1, alpha)
#'
#'@export
exprob <- function(cmbdf, win, alpha, intensities = "I")
{
if ( !is.CMBDataFrame(cmbdf) )
{
stop("Argument must be a CMBDataFrame")
}
if ( !is.CMBWindow(win) )
{
stop("Argument must be a CMBWindow")
}
cmbdf.win <- window(cmbdf, new.window = win)
sum(cmbdf.win[,intensities, drop = TRUE] > alpha)/sum(1-is.na(cmbdf.win[,intensities, drop = TRUE]))
}
#' Quantile-Quantile plots for \code{\link{CMBWindow}}s
#'
#' This function is a modification of standard \link{qqplot} functions to work
#' with \code{\link{CMBWindow}} regions.
#'
#' \code{\link{qqplotWin}} produces a QQ plot of quantiles of observations in
#' two \code{\link{CMBWindow}}s against each other for
#' the specified \code{\link{CMBDataFrame}} column
#' \code{intensities}. The function automatically adds a diagonal line.
#'
#'@param cmbdf A \code{\link{CMBDataFrame}}.
#'@param win1 A \code{\link{CMBWindow}}
#'@param win2 A \code{\link{CMBWindow}}
#'@param intensities A \code{\link{CMBDataFrame}} column with measured values.
#'
#'@return
#'
#' A list with quantile components x and y and a QQ plot with a diagonal line
#'
#'@references \link{qqnormWin}, \link{qqnorm}, \link{qqplot}
#'
#'@examples
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # cmbdf <- sampleCMB(df, sample.size = 10000)
#'
#' # win1 <- CMBWindow(theta = c(0,pi/2,pi/2), phi = c(0,0,pi/2))
#' # win2 <- CMBWindow(theta = c(2*pi/3,3*pi/4,3*pi/4, 2*pi/3),
#' # phi = c(pi/4,pi/4,pi/3,pi/3))
#'
#' # qqplotWin(cmbdf, win1, win2)
#'
#'@export
qqplotWin <- function(cmbdf, win1, win2, intensities = "I")
{
if ( !is.CMBDataFrame(cmbdf) )
{
stop("Argument must be a CMBDataFrame")
}
if ( !is.CMBWindow(win1) | !is.CMBWindow(win2) )
{
stop("Argument must be a CMBWindow")
}
cmbdf.win1 <- window(cmbdf, new.window = win1)
cmbdf.win2 <- window(cmbdf, new.window = win2)
stats::qqplot(cmbdf.win1[,intensities, drop = TRUE], cmbdf.win2[,intensities, drop = TRUE])
graphics::abline(c(0,1))
stats::qqplot(cmbdf.win1[,intensities, drop = TRUE], cmbdf.win2[,intensities, drop = TRUE],plot.it = FALSE)
}
#' Normal QQ plot for \code{\link{CMBWindow}}
#'
#' This function is a modification of standard \link{qqnorm} functions to work
#' with \code{\link{CMBWindow}} regions.
#'
#'\code{\link{qqnormWin}} returns a normal QQ plot of for the specified
#'\code{\link{CMBDataFrame}} column \code{intensities} and \code{\link{CMBWindow}}
#'region. The function automatically adds a line of a “theoretical” normal
#'quantile-quantile plot.
#'
#'
#'@param cmbdf A \code{\link{CMBDataFrame}}.
#'@param win A \code{\link{CMBWindow}}
#'@param intensities A \code{\link{CMBDataFrame}} column with measured values.
#'
#'@return
#'
#' A list with quantile components x and y and a normal QQ plot with QQ line
#'
#'@references \link{qqnorm}, \link{qqplot}, \link{qqplotWin}
#'
#'@examples
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # cmbdf <- sampleCMB(df, sample.size = 1000)
#'
#' # win1 <- CMBWindow(theta = c(0,pi/2,pi/2), phi = c(0,0,pi/2))
#' # qqnormWin(cmbdf, win1)
#'
#'@export
qqnormWin <- function(cmbdf, win, intensities = "I")
{
if ( !is.CMBDataFrame(cmbdf) )
{
stop("Argument must be a CMBDataFrame")
}
if ( !is.CMBWindow(win) )
{
stop("Argument must be a CMBWindow")
}
cmbdf.win <- window(cmbdf, new.window = win)
stats::qqnorm(cmbdf.win[,intensities, drop = TRUE])
stats::qqline(cmbdf.win[,intensities, drop = TRUE])
stats::qqnorm(cmbdf.win[,intensities, drop = TRUE],plot.it = FALSE)
}
#' CMB Entropy
#'
#'This function returns an estimated entropy for the specified
#'\code{\link{CMBDataFrame}} column \code{intensities} and \code{\link{CMBWindow}}
#'region. The functions employs the function \link{entropy} and uses histogram
#'counts of \code{intensities} for computations. All arguments of the standard
#'\link{entropy} can be used.
#'
#'@param cmbdf A \code{\link{CMBDataFrame}}.
#'@param win A \code{\link{CMBWindow}}
#'@param intensities A \code{\link{CMBDataFrame}} column with measured values.
#'@param method A method to estimate entropy, see \link{entropy}
#'
#'@return
#'
#'Estimated Shannon entropy for observations in \code{\link{CMBWindow}}
#'
#'@references \link{entropy}
#'
#'@examples
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # cmbdf <- sampleCMB(df, sample.size = 10000)
#' # win1 <- CMBWindow(theta = c(0,pi/2,pi/2), phi = c(0,0,pi/2))
#' # entropyCMB(cmbdf, win1)
#'
#'@export
entropyCMB <- function(cmbdf, win, intensities = "I", method)
{
if ( !is.CMBDataFrame(cmbdf) )
{
stop("Argument must be a CMBDataFrame")
}
if ( !is.CMBWindow(win) )
{
stop("Argument must be a CMBWindow")
}
cmbdf.win <- window(cmbdf, new.window = win)
y <- graphics::hist(cmbdf.win[,intensities, drop = TRUE], plot = FALSE)$counts
if (missing(method)) return(entropy::entropy(y))
return(entropy::entropy(y, method = method))
}
#'Chi-squared statistic for two \code{\link{CMBWindow}}s
#'
#'This function returns the empirical chi-squared statistic for \code{intensities}
#'observations from two \code{\link{CMBWindow}}s of the specified
#'\code{\link{CMBDataFrame}}. The functions employs the function \link{chi2.empirical} and uses histogram
#'counts of \code{intensities} for computations.
#'
#'@param cmbdf A \code{\link{CMBDataFrame}}.
#'@param win1 A \code{\link{CMBWindow}}
#'@param win2 A \code{\link{CMBWindow}}
#'@param intensities A \code{\link{CMBDataFrame}} column with measured values.
#'
#'@return
#'
#'Estimated Chi-squared statistic for observations in two
#'\code{\link{CMBWindow}}s. For small sample sizes and many zero counts
#'Chi-squared statistic is inefficient.
#'
#'@references \link{chi2.empirical}
#'
#'@examples
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # cmbdf <- sampleCMB(df, sample.size = 1000)
#' #
#' # win1 <- CMBWindow(theta = c(0,pi/2,pi/2), phi = c(0,0,pi/2))
#' # win2 <- CMBWindow(theta = c(pi/2,pi,pi/2), phi = c(0,0,pi/2))
#' # plot(win1)
#' # plot(win2,col="green")
#' #
#' # chi2CMB(cmbdf, win1, win2)
#'
#'@export
chi2CMB <- function(cmbdf, win1, win2, intensities = "I")
{
if ( !is.CMBDataFrame(cmbdf) )
{
stop("Argument must be a CMBDataFrame")
}
if ( !is.CMBWindow(win1) | !is.CMBWindow(win2) )
{
stop("Argument must be a CMBWindow")
}
cmbdf.win1 <- window(cmbdf, new.window = win1)
cmbdf.win2 <- window(cmbdf, new.window = win2)
y10 <- graphics::hist(cmbdf.win1[,intensities, drop = TRUE], plot = FALSE)$counts
y20 <- graphics::hist(cmbdf.win2[,intensities, drop = TRUE], plot = FALSE)$counts
nb <- min(length(y10),length(y20))
combI <- c(cmbdf.win1[,intensities, drop = TRUE],cmbdf.win2[,intensities, drop = TRUE])
yb <- seq(min(combI), max(combI), length.out=min(nb))
y1 <- graphics::hist(cmbdf.win1[,intensities, drop = TRUE], breaks=yb, plot = FALSE)$counts
y2 <- graphics::hist(cmbdf.win2[,intensities, drop = TRUE], breaks=yb, plot = FALSE)$counts
entropy::chi2.empirical(y1, y2)
}
#' Sample Renyi function
#'
#' This function computes values of the sample Renyi function. Returns
#' the estimated values of \eqn{T(q)} for \eqn{q} taking values on a grid.
#' For large data sets could be rather time consuming.
#'
#'
#'@param cmbdf A \code{\link{CMBDataFrame}}.
#'@param q.min Left endpoint of the interval to compute the Renyi function. The default
#'value is 1.01,
#'@param q.max Right endpoint of the interval to compute the Renyi function. The default
#'value is 10
#'@param N Number of points to compute the Renyi function. The default value is 20.
#'@param k.box A dyadic decomposition level in computing the Renyi function,
#'see the references in Details. The default value is \eqn{log2(nside(cmbdf)) - 3}
#'@param intensities A CMBDataFrame column with measured values
#'
#'
#'@return Data frame which first column is the sampling grid
#'\eqn{seq(q.min, q.max, length.out = N)} of \eqn{q} values. Another column
#'consists of values of the sample Renyi function \eqn{T(q)} computed
#'on the grid using the \eqn{k.box}th level dyadic decomposition
#' of the unit ball.
#'
#'
#'@references
#'(1) Leonenko, N., and Shieh, N. 2013. Rényi function for multifractal
#'random fields. Fractals 21, Article No. 1350009.
#'
#'(2) http://mathworld.wolfram.com/RenyiEntropy.html
#'
#' @examples
#'
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' #
#' # cmbdf <- CMBDataFrame("CMB_map_smica1024.fits")
#' # win <- CMBWindow(theta = c(pi/4,pi/2,pi/2), phi = c(0,0,pi/2))
#' # cmbdf<- window(cmbdf, new.window = win)
#' # Tq <- fRen(cmbdf)
#' #
#' # plot(Tq[,1], Tq[,2], ylab =expression(D[q]), xlab = "q",
#' # main = "Sample Renyi function", pch = 20, col = "blue")
#'
#'
#' @export
fRen <- function(cmbdf,
q.min = 1.01,
q.max = 10,
N = 20,
k.box = log2(nside(cmbdf)) - 3,
intensities = "I") {
if (!is.CMBDataFrame(cmbdf))
{
stop("Argument must be a CMBDataFrame")
}
ns1 <- nside(cmbdf)
pixind <- pix.CMBDataFrame(cmbdf)
nagrpix <- setdiff(1:(12 * ns1 ^ 2), pixind)
field.comp <- rep(0, 12 * ns1 ^ 2)
field.in <- cmbdf[, intensities, drop = T]
minint <- min(field.in)
field.in <- field.in - minint
field.final <- replace(field.comp, pixind, field.in)
res.max <- log2(ns1)
npix <- 12 * 4 ^ k.box
delta <- sqrt(4 * pi / npix)
lev.diff <- 4 ^ (res.max - k.box)
if (res.max - k.box > 0) {
nagrpix <- unique(ancestor(nagrpix, res.max - k.box))
}
agrpix <- setdiff((1:npix), nagrpix)
mu <- vector(mode = "numeric", length = length(agrpix))
field.total <- 0
i <- 1
for (j in agrpix) {
pixd <- (lev.diff * (j - 1) + 1):(lev.diff * j)
field.total <- field.total + sum(field.final[pixd])
mu[i] <- sum(field.final[pixd])
i <- i + 1
}
mu <- mu / field.total
Q <- seq(q.min, q.max, length.out = N)
Tq <- vector(mode = "numeric", length = N)
ri <- 1
for (q in Q) {
Tq[ri] <- 1 / (q - 1) * log2(sum(mu ^ q)) / log2(delta)
ri <- ri + 1
}
Tqf <- data.frame(q = Q, tq = Tq)
return(Tqf)
}
#' Extreme values
#'
#' This function returns \code{n} largest extreme values for the specified
#' \code{\link{CMBDataFrame}} column \code{intensities} and
#' \code{\link{CMBWindow}} region.
#'
#'@param cmbdf A \code{\link{CMBDataFrame}}.
#'@param win A \code{\link{CMBWindow}}
#'@param n An integer value.
#'@param intensities A \code{\link{CMBDataFrame}} column with measured values.
#'
#'@return
#'
#'A \code{\link{CMBDataFrame}} with \code{n} largest extreme values
#'
#'@examples
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # cmbdf <- sampleCMB(df, sample.size = 1000)
#' #
#' # win1 <- CMBWindow(theta = c(pi/2,pi,pi/2), phi = c(0,0,pi/2))
#' # extrCMB(cmbdf, win1,5)
#' #
#' ## Ploting the window and 5 top extreme values
#' # plot(win1)
#' # plot(extrCMB(cmbdf, win1,5), col ="blue", size = 4,add = TRUE)
#'
#'@export
extrCMB <- function(cmbdf, win, n, intensities = "I")
{
if ( !is.CMBDataFrame(cmbdf) )
{
stop("Argument must be a CMBDataFrame")
}
if ( !is.CMBWindow(win) )
{
stop("Argument must be a CMBWindow")
}
cmbdf.win <- window(cmbdf, new.window = win)
cmbdf.win[order(cmbdf.win[,intensities, drop = TRUE]),][1:n,]
}
#' q-statistic
#'
#'
#'This function returns an estimated q-statistic for the specified column
#'\code{intensities} in a \code{\link{CMBDataFrame}} and the list of
#'\code{\link{CMBWindow}}s.
#'
#'
#'@param cmbdf A \code{\link{CMBDataFrame}}.
#'@param listwin A list of \code{\link{CMBWindow}}s
#'@param intensities A \code{\link{CMBDataFrame}} column with measured values.
#'
#'@details The q-statistics is used to measure spatial stratified heterogeneity
#'and takes values in [0, 1]. It gives the percent of the variance of
#'\code{intensities} explained by the stratification. 0 corresponds to no spatial
#'stratified heterogeneity, 1 to perfect spatial stratified heterogeneity.
#'
#'@return
#'
#'Estimated q-statistics for observations in a list of \code{\link{CMBWindow}}s
#'
#'@references
#'Wang, J.F, Zhang, T.L, Fu, B.J. (2016). A measure of spatial
#'stratified heterogeneity. Ecological Indicators. 67: 250–256.
#'
#'@examples
#'
#' ## Download the map first
#' # downloadCMBMap(foreground = "smica", nside = 1024)
#' # df <- CMBDataFrame("CMB_map_smica1024.fits")
#' # cmbdf <- sampleCMB(df, sample.size = 1000)
#' # win1 <- CMBWindow(theta = c(0,pi/2,pi/2), phi = c(0,0,pi/2))
#' # win2 <- CMBWindow(theta = c(pi/2,pi,pi/2), phi = c(0,0,pi/2))
#' #
#' # lw <- list(win1, win2)
#' # qstat(cmbdf, lw)
#'
#'@export
qstat <- function(cmbdf, listwin, intensities = "I")
{
if ( !is.CMBDataFrame(cmbdf) )
{
stop("Argument must be a CMBDataFrame")
}
# cmbint <- sapply(listwin,function(v1,v2){
# window(v1, new.window = v2)}, v1=cmbdf[,intensities, drop = TRUE])
cmb <- cmbdf[,intensities]
cmbint <- sapply(listwin,function(v){
window(cmb, new.window = v)})
ni <- sapply(cmbint, length)
sigmai <- sapply(cmbint, stats::var)
sigma <- stats::var(unlist(cmbint, recursive = FALSE))
1- sum((ni-1)*sigmai)/(sigma*(sum(ni)-1))
}
#' Computes values of covariance functions
#'
#'
#'This function computes the covariances given the separation distance of locations.
#'Options for different covariance functions on spheres are available. The function
#'extends the function \code{\link[geoR]{cov.spatial}} for additional
#'covariance models on spheres.
#'
#'
#'@param obj Vector of distances between pairs of spatial locations.
#'@param cov.model A type of the correlation function. Available choices are: "matern",
#'"exponential","spherical", "powered.exponential", "cauchy", "gencauchy", "pure.nugget",
#'"askey", "c2wendland", "c4wendland", "sinepower", "multiquadric". Default is "matern"
#'@param cov.pars A vector with two covariance parameters. The first parameter
#'corresponds to the variance sigma^2. The second parameter corresponds
#'to the range phi of the correlation function.
#'@param kappa A smoothness parameter of the correlation function.
#'
#'@details
#'The function returns the value of the covariance \code{C(h)} at the distance \code{h}.
#'The covariance function has the form
#'
#'\deqn{C(h) = sigma^2 * rho(h/phi).}
#'
#'The parameters of the covariance are positive numbers \code{sigma^2}, \code{phi}
#' and \code{kappa}.
#'
#'Expressions for the correlation functions which are not included in the package
#'\strong{geoR}:
#'
#'\describe{
#' \item{\strong{askey}}{
#' \deqn{rho(h/phi) = (1 - h/phi)^{kappa}, if h < phi;}
#' \deqn{0, otherwise.}}
#' \item{\strong{c2wendland}}{
#' \deqn{rho(h/phi) = (1 + kappa * h/phi) * (1 - h/phi)^{kappa}, if h < phi;}
#' \deqn{0, otherwise.}}
#' \item{\strong{c4wendland}}{
#' \deqn{rho(h/phi) = (1 + kappa * h/phi + (kappa^2 - 1) * (h/phi)^2 / 3) * (1 - h/phi)^{kappa}, if h < phi;}
#' \deqn{0, otherwise.}}
#' \item{\strong{sinepower}}{
#' \deqn{rho(h/phi) = 1 - (sin(h/(2 phi))) ^{kappa}}}
#' \item{\strong{multiquadric}}{
#' \deqn{C(h) = (1 - phi) ^{(2 * kappa)} / (1 + phi^2 - 2 * phi * cos(h))^{kappa},
#' 0<phi<1}}
#' }
#'
#'Additional information can be found in the section Details in
#'\code{\link[geoR]{cov.spatial}}.
#'
#'@return
#'
#'Values of a covariance function for the given distances.
#'
#'@references
#'\strong{geoR} package, \code{\link[geoR]{cov.spatial}}
#'
#'T. Gneiting. Strictly and non-strictly positive definite functions on spheres.
#'Bernoulli 19 (2013), no. 4, 1327-1349.
#'
#'@examples
#'
#'## Compute Askey variogram at x = pi/4
#'
#' 1 - covmodelCMB(pi/4, cov.pars = c(1, pi), kappa = 3, cov.model = "askey" )
#'
#'## Plot of the Askey covariance function
#'
#' v1.f <- function(x, ...) {covmodelCMB(x, ...)}
#' curve( v1.f(x, cov.pars = c(1, 0.03), kappa = 3, cov.model = "askey"),
#' from = 0, to = 0.1, xlab = "distance", ylab = expression(cov(h)), lty = 2,
#' main = "covariance")
#'
#'@export
covmodelCMB <- function (obj,
cov.model = "matern",
cov.pars = stop("no cov.pars argument provided"),
kappa = 0.5) {
fn.env <- sys.frame(sys.nframe())
.checkCMB.cov.model(
cov.model = cov.model,
cov.pars = cov.pars,
kappa = kappa,
env = fn.env,
output = FALSE
)
phi <- get("phi", envir = fn.env)
sigmasq <- get("sigmasq", envir = fn.env)
covs <- obj
covs[] <- 0
for (i in 1:get("ns", envir = fn.env)) {
if (phi[i] < 1e-16)
cov.model[i] <- "pure.nugget"
obj.sc <- obj / phi[i]
cov.values <- switch(
cov.model[i],
askey = ifelse(obj < phi[i], (1 - (obj.sc)) ^ (kappa[i]), 0),
c2wendland = ifelse(obj < phi[i], (1 + kappa[i] * (obj.sc)) * (1 - (obj.sc)) ^
(kappa[i]), 0),
c4wendland = ifelse(obj < phi[i], (1 + kappa[i] * (obj.sc) + (((kappa[i]) ^
2 - 1) * (obj.sc) ^ 2
) / 3) * (1 - (obj.sc)) ^ (kappa[i]), 0),
sinepower = (1 - (sin(obj.sc/ 2)) ^ (kappa[i])),
multiquadric = (1 - phi[i]) ^ (2 * kappa[i]) / (1 + (phi[i]) ^ 2 - 2 *
(phi[i]) * cos(obj)) ^ (kappa[i]),
pure.nugget = rep(0, length(obj)),
wave = (1 / obj) * (phi[i] * sin(obj.sc)),
exponential = exp(-(obj.sc)),
matern = {
if (kappa[i] == 0.5)
exp(-(obj.sc))
else
geoR::matern(u = obj,
phi = phi[i],
kappa = kappa[i])
},
gaussian = exp(-((obj.sc) ^ 2)),
spherical = ifelse(obj < phi[i], (1 - 1.5 * (obj.sc) +
0.5 * (obj.sc) ^ 3), 0),
circular = {
obj.sc[obj.sc > 1] <- 1
ifelse(obj < phi[i], (1 - (2 * ((obj.sc) *
sqrt(1 - ((obj.sc) ^ 2)) +
asin(obj.sc)
)) / pi), 0)
},
cubic = {
ifelse(obj < phi[i], (1 - (
7 * (obj.sc ^ 2) -
8.75 * (obj.sc ^ 3) +
3.5 * (obj.sc ^ 5) -
0.75 * (obj.sc ^ 7)
)), 0)
},
power = (obj) ^ phi,
powered.exponential = exp(-((obj.sc) ^ kappa[i])),
cauchy = (1 + (obj.sc) ^ 2) ^ (-kappa[i]),
gneiting = {
obj.sc <- 0.301187465825 * obj.sc
t2 <- (1 - obj.sc)
t2 <- ifelse(t2 > 0, (t2 ^ 8), 0)
(1 + 8 * obj.sc + 25 * (obj.sc ^ 2) + 32 * (obj.sc ^ 3)) * t2
},
gencauchy = (1 + (obj.sc) ^ kappa[2]) ^ (-kappa[1] / kappa[2]),
gneiting.matern = {
obj.sc <- 0.301187465825 * obj.sc / kappa[2]
t2 <- (1 - obj.sc)
t2 <- ifelse(t2 > 0, (t2 ^ 8), 0)
cov.values <-
(1 + 8 * obj.sc + 25 * (obj.sc ^ 2) + 32 * (obj.sc ^ 3)) * t2
cov.values * geoR::matern(u = obj,
phi = phi[i],
kappa = kappa[1])
},
stop("wrong or no specification of cov.model")
)
if (cov.model[i] == "power") {
A <- (max(cov.values) / sqrt(pi)) * gamma((1 + phi[i]) / 2) *
gamma(1 - (phi[i] / 2))
A <- A * max(gamma(phi[i] + (1 + (1:2)) / 2) / (gamma(1 +
phi[i]) * gamma((1 + (
1:2
)) / 2)))
cov.values <- A - cov.values
cov.values <- cov.values / max(cov.values)
}
cov.values <- ifelse(obj < 1e-16, sigmasq[i], sigmasq[i] *
cov.values)
covs <- covs + cov.values
}
if (sum(sigmasq) < 1e-16)
covs[obj < 1e-16] <- 1
if (any(!is.finite(covs)))
warning("Infinity value in covmodelCMB ")
if (any(is.na(covs)))
warning("NA value in covmodelCMB ")
if (any(is.nan(covs)))
warning("NaN value in covmodelCMB ")
return(covs)
}
#' Estimates parameters of variograms
#'
#'
#'This function estimates variogram parameters by fitting a parametric model
#'from \code{\link{covmodelCMB}} to a sample variogram. The function extends
#'\code{\link[geoR]{variofit}} from the package \strong{geoR}
#'to additional covariance models on spheres.
#'
#'@param vario An object of the class \code{variogram} obtained as an output of
#'the function \code{\link{variogramCMB}}.
#'@param ini.cov.pars A vector with initial values for the variogram parameters.
#'The first parameter corresponds to the variance sigma^2. The second parameter
#'corresponds to the range phi of the correlation function.
#'@param cov.model A type of the variogram function. Available choices are: "matern",
#'"exponential","spherical", "powered.exponential", "cauchy", "gencauchy", "pure.nugget",
#'"askey", "c2wendland", "c4wendland", "sinepower", "multiquadric". The default is "matern"
#'@param fix.nugget logical. Indicates whether the nugget variance should be regarded
#'as fixed or be estimated. The default is FALSE.
#'@param nugget A value for the nugget parameter. Regarded as a fixed values if
#'\code{fix.nugget = TRUE} or as a initial value for the minimization algorithm if
#'\code{fix.nugget = FALSE}. The default is zero.
#'@param fix.kappa logical. Indicates whether the parameter kappa should be regarded
#'as fixed or be estimated. The default is TRUE.
#'@param kappa A value for the smoothness parameter. Regarded as a fixed values if
#'\code{fix.kappa = TRUE} or as a initial value for the minimization algorithm if
#'\code{fix.kappa = FALSE}. Required not in all covariance models, see
#'\code{\link{covmodelCMB}}. The default is 0.5.
#'@param simul.number number of simulation. Used if \code{vario} has empirical variograms
#'for more than one data-set (simulations). The default is NULL
#'@param max.dist A maximum distance to fit a variogram model. The default is
#'\code{x$max.dist}.
#'@param weights Weights used in the loss function in the minimization algorithm.
#'@param limits Lower and upper limits for the model parameters used
#'in the numerical minimisation by \code{minimisation.function = "optim"}.
#'@param minimisation.function Minimization function ("optim", "nlm", "nls") to estimate
#'the parameters.
#'@param messages logical. Indicates whether or not status messages are printed on
#'the screen.
#'@param ... other minimisation parameters
#'
#'
#'@details
#'The parameter values of a variogram function from \code{\link{covmodelCMB}} are
#'found by numerical optimization using one of the functions: \code{\link{optim}},
#'\code{\link{nlm}} and \code{\link{nls}}.
#'
#'The function extends \code{\link[geoR]{variofit}} from the package \strong{geoR}
#'to additional variogram models on spheres. Available models are: "matern",
#'"exponential", "spherical", "powered.exponential", "cauchy", "gencauchy",
#'"pure.nugget", "askey", "c2wendland", "c4wendland", "sinepower", "multiquadric".
#'
#'Additionally it rescales an empirical variogram to the range \code{[0,1]} before
#'numerical optimisation and then transforms all obtained results to the original
#'scale. If \code{ini.cov.pars} are not provided then the 5x5 grid
#'\code{(seq(0,max(vario$v),l=5), seq(0,vario$max.dist,l=5))}
#'of initial values of sigma^2 and phi is used.
#'
#'
#'@return
#'An object of the class \code{variomodel} and \code{variofit}, see
#'\code{\link[geoR]{variofit}}
#'
#'@references
#'\strong{geoR} package, \code{\link[geoR]{variofit}}, \code{\link{covmodelCMB}}
#'
#'@examples
#' #
#' # df <- CMBDataFrame("../CMB_map_smica1024.fits")
#' # cmbdf <- sampleCMB(df, sample.size = 10000)
#' # varcmb <- variogramCMB(cmbdf, max.dist = 0.1, num.bins = 30)
#' # varcmb
#' #
#' # ols <- variofitCMB(varcmb, fix.nug=FALSE, wei="equal", cov.model= "matern")
#' # plot(varcmb)
#' # lines(ols, lty=2)
#' # str(ols)
#' #
#' # ols <- variofitCMB(varcmb, fix.nug = TRUE, kappa = 3, wei = "equal",
#' # cov.model = "askey")
#' # plot(varcmb, main = ols$cov.model)
#' # linesCMB(ols, lty = 2)
#' # str(ols)
#'
#'@export
variofitCMB <- function (vario, ini.cov.pars, cov.model, fix.nugget = FALSE,
nugget = 0, fix.kappa = TRUE, kappa = 0.5, simul.number = NULL,
max.dist = vario$max.dist, weights, minimisation.function,
limits = geoR::pars.limits(), messages, ...) {
cov.model <- match.arg(cov.model, choices = CMB.cov.models)
vario1 <- vario
vario1$v <- vario1$v/max(vario1$v)
if (missing(ini.cov.pars)){
ini.cov.pars <- expand.grid(seq(0,max(vario1$v),l=5), seq(0,vario1$max.dist,l=5))
}
if (cov.model %in% c("matern",
"exponential",
"spherical",
"powered.exponential",
"cauchy",
"gencauchy",
"pure.nugget")){
variofitCMB <- geoR::variofit(vario1, ini.cov.pars=ini.cov.pars, cov.model=cov.model, fix.nugget=fix.nugget,
nugget=nugget, fix.kappa=fix.kappa, kappa=kappa, simul.number=simul.number,
max.dist=max.dist, weights=weights, minimisation.function=minimisation.function,
limits = limits, messages=messages, ...)
}
if (cov.model %in% c( "askey",
"c2wendland",
"c4wendland",
"sinepower",
"multiquadric")){
variofitCMB <- variofit1(vario1, ini.cov.pars=ini.cov.pars, cov.model=cov.model, fix.nugget=fix.nugget,
nugget=nugget, fix.kappa=fix.kappa, kappa=kappa, simul.number=simul.number,
max.dist=max.dist, weights=weights, minimisation.function=minimisation.function,
limits = limits, messages=messages, ...)
}
variofitCMB$cov.pars[1] <- variofitCMB$cov.pars[1]*max(vario$v)
variofitCMB$nugget <- variofitCMB$nugget*max(vario$v)
return(variofitCMB)
}
# Helper function for variofitCMB
variofit1 <- function (vario,
ini.cov.pars,
cov.model,
fix.nugget,
nugget,
fix.kappa,
kappa,
simul.number,
max.dist,
weights,
minimisation.function,
limits,
messages,
...) {
call.fc <- match.call()
if (missing(messages))
messages.screen <-
as.logical(ifelse(is.null(getOption("geoR.messages")),
TRUE, getOption("geoR.messages")))
else
messages.screen <- messages
if (length(class(vario)) == 0 ||
all(class(vario) != "variogram"))
warning("object vario should preferably be of the geoR's class \"variogram\"")
if (!missing(ini.cov.pars)) {
if (any(class(ini.cov.pars) == "eyefit"))
cov.model <- ini.cov.pars[[1]]$cov.model
if (any(class(ini.cov.pars) == "variomodel"))
cov.model <- ini.cov.pars$cov.model
}
if (missing(cov.model))
cov.model <- "matern"
cov.model <- match.arg(cov.model, choices = CMB.cov.models)
if (cov.model == "stable")
cov.model <- "powered.exponential"
if (cov.model == "powered.exponential")
if (limits$kappa["upper"] > 2)
limits$kappa["upper"] <- 2
if (missing(weights)) {
if (vario$output.type == "cloud")
weights <- "equal"
else
weights <- "npairs"
}
else
weights <- match.arg(weights, choices = c("npairs",
"equal", "cressie"))
if (messages.screen) {
cat(paste("variofit: covariance model used is", cov.model,
"\n"))
cat(paste("variofit: weights used:", weights, "\n"))
}
if (missing(minimisation.function))
minimisation.function <- "optim"
if (any(cov.model == c("linear", "power")) &
minimisation.function ==
"nls") {
cat(
"warning: minimisation function nls can not be used with given cov.model.\n changing for \"optim\".\n"
)
minimisation.function <- "optim"
}
if (minimisation.function == "nls" & weights != "equal") {
warning(
"variofit: minimisation function nls can only be used with weights=\"equal\".\n changing for \"optim\".\n"
)
minimisation.function <- "optim"
}
if (is.matrix(vario$v) & is.null(simul.number))
stop(
"object in vario$v is a matrix. This function works for only 1 empirical variogram at once\n"
)
if (!is.null(simul.number))
vario$v <- vario$v[, simul.number]
if (mode(max.dist) != "numeric" || length(max.dist) > 1)
stop("a single numerical value must be provided in the argument max.dist")
if (max.dist == vario$max.dist)
XY <- list(u = vario$u, v = vario$v, n = vario$n)
else
XY <- list(u = vario$u[vario$u <= max.dist],
v = vario$v[vario$u <=
max.dist],
n = vario$n[vario$u <= max.dist])
if (cov.model == "pure.nugget") {
minimisation.function <- "not used"
message <-
"correlation function does not require numerical minimisation"
if (weights == "equal")
lm.wei <- rep(1, length(XY$u))
else
lm.wei <- XY$n
if (cov.model == "pure.nugget") {
if (fix.nugget) {
temp <- stats::lm((XY$v - nugget) ~ 1, weights = lm.wei)
cov.pars <- c(temp$coef, 0)
}
else {
temp <- stats::lm(XY$v ~ 1, weights = lm.wei)
nugget <- temp$coef
cov.pars <- c(0, 0)
}
}
value <- sum((temp$residuals) ^ 2)
}
else {
if (messages.screen)
cat(paste(
"variofit: minimisation function used:",
minimisation.function,
"\n"
))
umax <- max(vario$u)
vmax <- max(vario$v)
if (missing(ini.cov.pars)) {
ini.cov.pars <- as.matrix(expand.grid(c(vmax / 2, 3 *
vmax / 4, vmax), seq(0, 0.8 * umax, len = 6)))
if (!fix.nugget)
nugget <- unique(c(nugget, vmax / 10, vmax / 4, vmax / 2))
if (!fix.kappa)
kappa <- unique(c(kappa, 0.25, 0.5, 1, 1.5, 2))
if (messages.screen)
warning("initial values not provided - running the default search")
}
else {
if (any(class(ini.cov.pars) == "eyefit")) {
init <- nugget <- kappa <- NULL
for (i in 1:length(ini.cov.pars)) {
init <- drop(rbind(init, ini.cov.pars[[i]]$cov.pars))
nugget <- c(nugget, ini.cov.pars[[i]]$nugget)
if (cov.model == "gneiting.matern")
kappa <- drop(rbind(kappa, ini.cov.pars[[i]]$kappa))
else
kappa <- c(kappa, ini.cov.pars[[i]]$kappa)
}
ini.cov.pars <- init
}
if (any(class(ini.cov.pars) == "variomodel")) {
nugget <- ini.cov.pars$nugget
kappa <- ini.cov.pars$kappa
ini.cov.pars <- ini.cov.pars$cov.pars
}
}
if (is.matrix(ini.cov.pars) | is.data.frame(ini.cov.pars)) {
ini.cov.pars <- as.matrix(ini.cov.pars)
if (nrow(ini.cov.pars) == 1)
ini.cov.pars <- as.vector(ini.cov.pars)
else {
if (ncol(ini.cov.pars) != 2)
stop(
"\nini.cov.pars must be a matrix or data.frame with 2 components:
\ninitial values for sigmasq (partial sill) and phi (range parameter)\n"
)
}
}
else if (length(ini.cov.pars) > 2)
stop("\nini.cov.pars must provide initial values for sigmasq and phi\n")
if (is.matrix(ini.cov.pars) | (length(nugget) > 1) |
(length(kappa) > 1)) {
if (messages.screen)
cat("variofit: searching for best initial value ...")
ini.temp <- matrix(ini.cov.pars, ncol = 2)
grid.ini <- as.matrix(expand.grid(
sigmasq = unique(ini.temp[,
1]),
phi = unique(ini.temp[, 2]),
tausq = unique(nugget),
kappa = unique(kappa)
))
v.loss <- function(parms, u, v, n, cov.model, weights) {
sigmasq <- parms[1]
phi <- parms[2]
if (cov.model == "power")
phi <- 2 * exp(phi) / (1 + exp(phi))
tausq <- parms[3]
kappa <- parms[4]
if (cov.model == "power")
v.mod <- tausq + covmodelCMB(u,
cov.pars = c(sigmasq,
phi),
cov.model = "power",
kappa = kappa
)
else
v.mod <- (sigmasq + tausq) - covmodelCMB(u,
cov.pars = c(sigmasq, phi),
cov.model = cov.model,
kappa = kappa
)
if (weights == "equal")
loss <- sum((v - v.mod) ^ 2)
if (weights == "npairs")
loss <- sum(n * (v - v.mod) ^ 2)
if (weights == "cressie")
loss <- sum((n / (v.mod ^ 2)) * (v - v.mod) ^ 2)
return(loss)
}
grid.loss <- apply(
grid.ini,
1,
v.loss,
u = XY$u,
v = XY$v,
n = XY$n,
cov.model = cov.model,
weights = weights
)
ini.temp <- grid.ini[which(grid.loss == min(grid.loss))[1],
, drop = FALSE]
if (is.R())
rownames(ini.temp) <- "initial.value"
if (messages.screen) {
cat(" selected values:\n")
print(rbind(round(ini.temp, digits = 2), status = ifelse(
c(FALSE,
FALSE, fix.nugget, fix.kappa), "fix", "est"
)))
cat(paste("loss value:", min(grid.loss), "\n"))
}
names(ini.temp) <- NULL
ini.cov.pars <- ini.temp[1:2]
nugget <- ini.temp[3]
kappa <- ini.temp[4]
grid.ini <- NULL
}
if (ini.cov.pars[1] > 2 * vmax)
warning("unreasonable initial value for sigmasq (too high)")
if (ini.cov.pars[1] + nugget > 3 * vmax)
warning("unreasonable initial value for sigmasq + nugget (too high)")
if (vario$output.type != "cloud") {
if (ini.cov.pars[1] + nugget < 0.3 * vmax)
warning("unreasonable initial value for sigmasq + nugget (too low)")
}
if (nugget > 2 * vmax)
warning("unreasonable initial value for nugget (too high)")
if (ini.cov.pars[2] > 1.5 * umax)
warning("unreasonable initial value for phi (too high)")
if (!fix.kappa) {
if (cov.model == "powered.exponential")
Tkappa.ini <- log(kappa / (2 - kappa))
else
Tkappa.ini <- log(kappa)
}
if (minimisation.function == "nls") {
if (ini.cov.pars[2] == 0)
ini.cov.pars <- max(XY$u) / 10
if (kappa == 0)
kappa <- 0.5
if (cov.model == "power")
Tphi.ini <- log(ini.cov.pars[2] / (2 - ini.cov.pars[2]))
else
Tphi.ini <- log(ini.cov.pars[2])
XY$cov.model <- cov.model
if (fix.nugget) {
XY$nugget <- as.vector(nugget)
if (fix.kappa) {
XY$kappa <- as.vector(kappa)
res <- stats::nls((v - nugget) ~ matrix((
1 - covmodelCMB (u,
cov.pars = c(1, exp(Tphi)),
cov.model = cov.model,
kappa = kappa
)
), ncol = 1),
start = list(Tphi = Tphi.ini),
data = XY,
algorithm = "plinear",
...
)
}
else {
if (cov.model == "powered.exponential")
res <- stats::nls((v - nugget) ~ matrix((
1 - covmodelCMB(u,
cov.pars = c(1, exp(Tphi)),
cov.model = cov.model,
kappa = (2 * exp(Tkappa) /
(1 + exp(Tkappa)))
)
),
ncol = 1),
start = list(Tphi = Tphi.ini,
Tkappa = Tkappa.ini),
data = XY,
algorithm = "plinear",
...
)
else
res <- stats::nls((v - nugget) ~ matrix((
1 -
covmodelCMB(u,
cov.pars = c(1, exp(Tphi)),
cov.model = cov.model,
kappa = exp(Tkappa)
)
),
ncol = 1),
start = list(Tphi = Tphi.ini,
Tkappa = Tkappa.ini),
data = XY,
algorithm = "plinear",
...
)
kappa <- exp(stats::coef(res)["Tkappa"])
names(kappa) <- NULL
}
cov.pars <- stats::coef(res)[c(".lin", "Tphi")]
names(cov.pars) <- NULL
}
else {
if (fix.kappa) {
XY$kappa <- kappa
res <- stats::nls(
v ~ cbind(1, (
1 - covmodelCMB (
u,
cov.pars = c(1, exp(Tphi)),
cov.model = cov.model,
kappa = kappa
)
)),
start = list(Tphi = Tphi.ini),
algorithm = "plinear",
data = XY,
...
)
}
else {
if (cov.model == "powered.exponential")
res <- stats::nls(
v ~ cbind(1, (
1 - covmodelCMB (
u,
cov.pars = c(1, exp(Tphi)),
cov.model = cov.model,
kappa = (2 * exp(Tkappa) /
(1 + exp(Tkappa)))
)
)),
start = list(Tphi = Tphi.ini, Tkappa = Tkappa.ini),
algorithm = "plinear",
data = XY,
...
)
else
res <- stats::nls(
v ~ cbind(1, (
1 - covmodelCMB (
u,
cov.pars = c(1, exp(Tphi)),
cov.model = cov.model,
kappa = exp(Tkappa)
)
)),
start = list(Tphi = Tphi.ini,
Tkappa = Tkappa.ini),
algorithm = "plinear",
data = XY,
...
)
kappa <- exp(stats::coef(res)["Tkappa"])
names(kappa) <- NULL
}
nugget <- stats::coef(res)[".lin1"]
names(nugget) <- NULL
cov.pars <- stats::coef(res)[c(".lin2", "Tphi")]
names(cov.pars) <- NULL
}
if (cov.model == "power")
cov.pars[2] <-
2 * exp(cov.pars[2]) / (1 + exp(cov.pars[2]))
else
cov.pars[2] <- exp(cov.pars[2])
if (nugget < 0 | cov.pars[1] < 0) {
warning(
"\nvariofit: negative variance parameter found using the default option \"nls\".\n Try another minimisation function and/or fix some of the parameters.\n"
)
temp <- c(
sigmasq = cov.pars[1],
phi = cov.pars[2],
tausq = nugget,
kappa = kappa
)
print(rbind(round(temp, digits = 4), status = ifelse(
c(FALSE,
FALSE, fix.nugget, fix.kappa), "fix", "est"
)))
return(invisible())
}
value <- sum(stats::resid(res) ^ 2)
message <- "nls does not provides convergence message"
}
if (minimisation.function == "nlm" | minimisation.function ==
"optim") {
.global.list <- list(
u = XY$u,
v = XY$v,
n = XY$n,
fix.nugget = fix.nugget,
nugget = nugget,
fix.kappa = fix.kappa,
kappa = kappa,
cov.model = cov.model,
m.f = minimisation.function,
weights = weights
)
ini <- ini.cov.pars
if (cov.model == "power")
ini[2] <- log(ini[2] / (2 - ini[2]))
if (cov.model == "linear")
ini <- ini[1]
if (fix.nugget == FALSE)
ini <- c(ini, nugget)
if (!fix.kappa)
ini <- c(ini, Tkappa.ini)
names(ini) <- NULL
if (minimisation.function == "nlm") {
result <- stats::nlm(.loss1.vario, ini, g.l = .global.list,
...)
result$par <- result$estimate
result$value <- result$minimum
result$convergence <- result$code
if (!is.null(get(".temp.theta")))
result$par <- get(".temp.theta")
}
else {
lower.l <- sapply(limits, function(x)
x[1])
upper.l <- sapply(limits, function(x)
x[2])
if (fix.kappa == FALSE) {
if (fix.nugget) {
lower <- lower.l[c("sigmasq.lower", "phi.lower",
"kappa.lower")]
upper <- upper.l[c("sigmasq.upper", "phi.upper",
"kappa.upper")]
}
else {
lower <- lower.l[c("sigmasq.lower",
"phi.lower",
"tausq.rel.lower",
"kappa.lower")]
upper <- upper.l[c("sigmasq.upper",
"phi.upper",
"tausq.rel.upper",
"kappa.upper")]
}
}
else {
if (cov.model == "power") {
if (fix.nugget) {
lower <- lower.l[c("sigmasq.lower", "phi.lower")]
upper <- upper.l[c("sigmasq.upper", "phi.upper")]
}
else {
lower <- lower.l[c("sigmasq.lower",
"phi.lower",
"tausq.rel.lower")]
upper <- upper.l[c("sigmasq.upper",
"phi.upper",
"tausq.rel.upper")]
}
}
else {
lower <- lower.l["phi.lower"]
upper <- upper.l["phi.upper"]
}
}
result <- stats::optim(
ini,
.loss1.vario,
method = "L-BFGS-B",
hessian = TRUE,
lower = lower,
upper = upper,
g.l = .global.list,
...
)
}
value <- result$value
message <- paste(minimisation.function,
"convergence code:",
result$convergence)
if (cov.model == "linear")
result$par <- c(result$par[1], 1, result$par[-1])
cov.pars <- as.vector(result$par[1:2])
if (cov.model == "power")
cov.pars[2] <-
2 * exp(cov.pars[2]) / (1 + exp(cov.pars[2]))
if (!fix.kappa) {
if (fix.nugget)
kappa <- result$par[3]
else {
nugget <- result$par[3]
kappa <- result$par[4]
}
if (.global.list$cov.model == "powered.exponential")
kappa <- 2 * (exp(kappa)) / (1 + exp(kappa))
else
kappa <- exp(kappa)
}
else if (!fix.nugget)
nugget <- result$par[3]
}
}
if (cov.model == "matern")
# kappa <- min(0.5,kappa)
if (cov.model == "powered.exponential")
kappa <- min(1,kappa)
if (cov.model == "askey")
kappa <- max(2,kappa)
if (cov.model == "c2wendland")
kappa <- max(4,kappa)
if (cov.model == "c4wendland")
kappa <- max(6,kappa)
if (cov.model == "sinepower")
kappa <- min(2,kappa)
estimation <- list(
nugget = nugget,
cov.pars = cov.pars,
cov.model = cov.model,
kappa = kappa,
value = value,
trend = vario$trend,
beta.ols = vario$beta.ols,
practicalRange = practicalRangeCMB(
cov.model = cov.model,
phi = cov.pars[2],
kappa = kappa
),
max.dist = max.dist,
minimisation.function = minimisation.function
)
estimation$weights <- weights
if (weights == "equal")
estimation$method <- "OLS"
else
estimation$method <- "WLS"
estimation$fix.nugget <- fix.nugget
estimation$fix.kappa <- fix.kappa
estimation$lambda <- vario$lambda
estimation$message <- message
estimation$call <- call.fc
oldClass(estimation) <- c("variomodel", "variofit")
return(estimation)
}
#'Plot theoretical CMBCovariance
#'
#'Plots theoretical covariance functions from the list defined in \code{\link{covmodelCMB}}
#'
#'@param cov.model A type of the correlation function. Available choices are: "matern",
#'"exponential","spherical", "powered.exponential", "cauchy", "gencauchy", "pure.nugget",
#'"askey", "c2wendland", "c4wendland", "sinepower", "multiquadric". The default is "matern"
#'@param sigmasq The variance parameter as documented in \code{\link{covmodelCMB}}.
#'The default is 1.
#'@param phi The range parameter as documented in \code{\link{covmodelCMB}}. The default is
#'\code{pi}.
#'@param kappa A smoothness parameter of the correlation function. The default is 0.5.
#'@param from A lower range of the plotting region. The default is \code{lb =0}
#'@param to An upper range of the plotting region. The default is \code{ub= pi}.
#'@param ... optional plotting parameters.
#'
#'
#'@return Produces a plot with the theoretical covariance function.
#'
#'@references \code{\link{covmodelCMB}}
#'
#'@examples
#'
#' plotcovmodelCMB("matern", sigmasq = 5)
#' plotcovmodelCMB("askey", phi = pi/4, to = pi/2, kappa = 4)
#'
#'
#'@export
plotcovmodelCMB <- function (cov.model = "matern", sigmasq=1,
phi = pi,
kappa = 0.5,
from =0 , to = pi, ...){
cov.model <- match.arg(cov.model, choices = CMB.cov.models)
.checkCMB.cov.model(cov.model = cov.model,
cov.pars = c(sigmasq, phi),
kappa = kappa,
output = FALSE)
x <- NULL; rm(x) # Trick R CMD Check
graphics::curve(covmodelCMB(x, cov.pars = c(sigmasq, phi),
kappa = kappa, cov.model = cov.model),
from = from,
to = to,
xlab = "distance",
ylab = expression(C(h)),
lty = 2,
main = bquote(paste(.(cov.model)~"covariance function"))
)
}
#'Plot theoretical variogram
#'
#'Plots theoretical variogram functions from the list defined in \code{\link{covmodelCMB}}
#'
#'@param cov.model A type of the variogram function. Available choices are: "matern",
#'"exponential","spherical", "powered.exponential", "cauchy", "gencauchy", "pure.nugget",
#'"askey", "c2wendland", "c4wendland", "sinepower", "multiquadric". The default is "matern"
#'@param sigmasq The variance parameter as documented in \code{\link{covmodelCMB}}.
#'The default is 1.
#'@param phi The range parameter as documented in \code{\link{covmodelCMB}}. The default is
#'\code{pi}.
#'@param kappa A smoothness parameter of the variogram function. The default is 0.5.
#'@param from A lower range of the plotting region. The default is \code{lb =0}
#'@param to An upper range of the plotting region. The default is \code{ub= pi}.
#'@param ... optional plotting parameters.
#'
#'
#'@return Produces a plot with the theoretical variogram.
#'
#'@references \code{\link{covmodelCMB}}
#'
#'@examples
#'
#' plotvariogram("matern", sigmasq = 5)
#' plotvariogram("askey", phi = pi/4, to = pi/2, kappa = 4)
#'
#'
#'@export
plotvariogram <- function (cov.model = "matern", sigmasq=1,
phi = pi,
kappa = 0.5,
from =0 , to = pi, ...){
cov.model <- match.arg(cov.model, choices = CMB.cov.models)
.checkCMB.cov.model(cov.model = cov.model,
cov.pars = c(sigmasq, phi),
kappa = kappa,
output = FALSE)
x <- NULL; rm(x) # Trick R CMD Check
graphics::curve(covmodelCMB(0, cov.pars = c(sigmasq, phi),
kappa = kappa,
cov.model = cov.model)
- covmodelCMB(x, cov.pars = c(sigmasq, phi),
kappa = kappa,
cov.model = cov.model),
from = from,
to = to,
xlab = "distance",
ylab = expression(gamma(h)),
lty = 2,
main = bquote(paste(.(cov.model)~"variogram"))
)
}
#' Adds lines of fitted variograms to variogram plots
#'
#'
#'This function adds a line with the variogram model fitted by the function
#'\code{\link{variofitCMB}} to a current variogram plot. The function modifies
#'\code{\link[geoR]{lines.variomodel.variofit}} from the package \strong{geoR}
#'for additional covariance models on spheres.
#'
#'
#'@param x An object of the class \code{variofit} containing information about
#'the fitted model obtained as an output of the function \code{\link{variofitCMB}}.
#'@param max.dist A maximum distance to draw the variogram line. The default is
#'\code{x$max.dist}.
#'@param scaled logical. If TRUE the sill in the plot is 1.
#'@param ... other plotting parameters passed to \code{\link[graphics]{curve}}
#'
#'@details
#'The function adds a line with fitted variogram model to a plot. It is used
#'to compare empirical variograms against fitted models returned by
#'\code{\link{variofitCMB}}.
#'
#'Available models are: "matern", "exponential",
#'"spherical", "powered.exponential", "cauchy", "gencauchy", "pure.nugget",
#'"askey", "c2wendland", "c4wendland", "sinepower", "multiquadric".
#'
#'@return
#'A line with a fitted variogram model is added to a plot.
#'
#'@references
#'\strong{geoR} package, \code{\link[geoR]{lines.variomodel.variofit}},
#'\code{\link{covmodelCMB}}, \code{\link{variofitCMB}}
#'
#'@examples
#' ## Plot the fitted Matern variogram versus its empirical variogram
#' #
#' # df <- CMBDataFrame("../CMB_map_smica1024.fits")
#' # cmbdf <- sampleCMB(df, sample.size = 10000)
#' # varcmb <- variogramCMB(cmbdf, max.dist = 0.1, num.bins = 30)
#' # varcmb
#' # ols <- variofitCMB(varcmb, fix.nug=FALSE, wei="equal", cov.model= "matern")
#' # plot(varcmb)
#' # lines(ols, lty=2)
#' #
#' ## Plot the fitted Askey variogram versus its empirical variogram
#' #
#' # ols <- variofitCMB(vario1, ini.cov.pars = c(1, 0.03), fix.nug = TRUE,
#' # kappa = 3, wei = "equal", cov.model = "askey")
#' # plot(varcmb, main = ols$cov.model)
#' # linesCMB(ols, lty = 2)
#'
#'@export
linesCMB <- function (x, max.dist, scaled = FALSE, ...) {
my.l <- list()
if (missing(max.dist)) {
my.l$max.dist <- x$max.dist
if (is.null(my.l$max.dist))
stop("argument max.dist needed for this object")
}
else
my.l$max.dist <- max.dist
if (any(
x$cov.model == c(
"matern",
"powered.exponential",
"cauchy",
"gencauchy",
"gneiting.matern",
"askey",
"c2wendland",
"c4wendland",
"sinepower",
"multiquadric"
)
))
my.l$kappa <- x$kappa
else
kappa <- NULL
if (is.vector(x$cov.pars))
my.l$sill.total <- x$nugget + x$cov.pars[1]
else
my.l$sill.total <- x$nugget + sum(x$cov.pars[, 1])
my.l$nugget <- x$nugget
my.l$cov.pars <- x$cov.pars
my.l$cov.model <- x$cov.model
if (scaled) {
if (is.vector(x$cov.model))
my.l$cov.pars[1] <- my.l$cov.pars[1] / my.l$sill.total
else
my.l$cov.pars[, 1] <-
my.l$cov.cov.pars[, 1] / my.l$sill.total
my.l$sill.total <- 1
}
gamma.f <- function(x, my.l) {
if (any(my.l$cov.model == c("linear", "power")))
return(my.l$nugget + my.l$cov.pars[1] * (x ^ my.l$cov.pars[2]))
else
return(
my.l$sill.total -
covmodelCMB(x,
cov.model = my.l$cov.model,
kappa = my.l$kappa,
cov.pars = my.l$cov.pars
)
)
}
x <- NULL; rm(x) #Trick R CMD Check
graphics::curve(gamma.f(x, my.l = my.l),
from = 0,
to = my.l$max.dist,
add = TRUE,
...
)
return(invisible())
}
#' Practical range for covariance function
#'
#' This function computes the practical range for covariance functions on spheres.
#' The function extends \code{\link[geoR]{practicalRange}} from the
#' package \strong{geoR} to additional covariance models on spheres.
#'
#'
#'@param cov.model A type of the correlation function. Available choices are: "matern",
#'"exponential","spherical", "powered.exponential", "cauchy", "gencauchy", "pure.nugget",
#'"askey", "c2wendland", "c4wendland", "sinepower", "multiquadric".
#'@param phi The range parameter as documented in \code{\link{covmodelCMB}}
#'@param kappa A smoothness parameter of the correlation function.
#'@param correlation A correlation threshold (default is 0.05)
#'@param ... other optimisation parameters
#'
#'@details
#' The practical(effective) range for a covariance function is the distance at which
#' a covariance function first time reaches the specified value \code{correlation}. For
#' covariance functions on spheres the practical range does not exceed \eqn{pi}, the
#' distance beyond which a covariance function is not defined. For the covariance
#' functions "spherical", "askey", "c2wendland", "c4wendland" their practical ranges
#' are equal to lengths of their support.
#'
#'@return
#'Value of the practical range for the covariance function specified in \code{\link{covmodelCMB}}
#'
#'@references
#'\strong{geoR} package, \code{\link[geoR]{practicalRange}}, \code{\link{covmodelCMB}}
#'
#'@examples
#'
#'practicalRangeCMB(cov.model = "sinepower", phi = 0.1, kappa = 0.5)
#'practicalRangeCMB(cov.model = "askey", phi = 0.1, kappa = 0.5)
#'
#'@export
practicalRangeCMB <- function (cov.model,
phi,
kappa = 0.5,
correlation = 0.05,...){
cov.model <- match.arg(cov.model, choices = CMB.cov.models)
.checkCMB.cov.model(
cov.model = cov.model,
cov.pars = c(1, phi),
kappa = kappa,
output = FALSE
)
if (cov.model %in% c("circular",
"cubic",
"spherical",
"askey",
"c2wendland",
"c4wendland"))
return(min(phi, pi))
if (any(cov.model %in% c("pure.nugget")))
return(0)
if (any(cov.model %in% c("linear", "sinepower", "multiquadric")))
return(pi)
if (any(cov.model %in% c("power")))
return(Inf)
findRange <- function(range, cm, p, k, cor){
covmodelCMB (
range,
cov.model = cm,
kappa = k,
cov.pars = c(1, p)
) - cor
}
pr <- stats::uniroot(findRange,
interval = c(0, 50 * phi + 1),
cm = cov.model,
p = phi,
k = kappa,
cor = correlation,...)$root
return(min(pr, pi))
}
CMB.cov.models <- c(
"matern",
"exponential",
"spherical",
"powered.exponential",
"cauchy",
"gencauchy",
"pure.nugget",
"askey",
"c2wendland",
"c4wendland",
"sinepower",
"multiquadric"
)
# Helper function for variofitCMB
.loss1.vario <- function (theta, g.l)
{
if (g.l$cov.model == "linear")
theta <- c(theta[1], 1, theta[-1])
##
## Imposing constraints for nlm
##
if (g.l$m.f == "nlm") {
assign(".temp.theta", NULL)
if (!g.l$fix.kappa) {
if (g.l$fix.nugget) {
if (g.l$cov.model == "power")
theta.minimiser <- theta[1]
else
theta.minimiser <- theta[1:2]
Tkappa <- theta[3]
}
else{
if (g.l$cov.model == "power")
theta.minimiser <- theta[c(1:3)]
else
theta.minimiser <- theta[1:3]
Tkappa <- theta[4]
}
}
else
theta.minimiser <- theta
penalty <- 10000 * sum(0 - pmin(theta.minimiser, 0))
theta <- pmax(theta.minimiser, 0)
if (!g.l$fix.kappa)
theta <- c(theta.minimiser, Tkappa)
if (any(theta.minimiser < 0))
assign(".temp.theta", theta)
else
penalty <- 0
}
else
penalty <- 0
##
## reading parameters
##
if (!g.l$fix.kappa) {
if (g.l$fix.nugget) {
tausq <- g.l$nugget
Tkappa <- theta[3]
}
else{
tausq <- theta[3]
Tkappa <- theta[4]
}
## kappa now > 0 for both nlm() and optim()
##if(g.l$m.f == "nlm"){
if (g.l$cov.model == "powered.exponential")
kappa <- 2 * (exp(Tkappa)) / (1 + exp(Tkappa))
else
kappa <- exp(Tkappa)
##}
##else kappa <- Tkappa
}
else{
kappa <- g.l$kappa
if (g.l$fix.nugget)
tausq <- g.l$nugget
else
tausq <- theta[3]
}
##
sigmasq <- theta[1]
phi <- theta[2]
if (g.l$cov.model == "power")
phi <- 2 * exp(phi) / (1 + exp(phi))
sill.total <- sigmasq + tausq
##
## Computing values for the theoretical variogram
##
if (any(g.l$cov.model == c("linear", "power")))
gammaU <- tausq + sigmasq * (g.l$u ^ phi)
else
gammaU <-
sill.total - covmodelCMB(g.l$u,
cov.model = g.l$cov.model,
kappa = kappa,
cov.pars = c(sigmasq, phi)
)
##
## Computing loss function
##
if (g.l$weight == "equal")
loss <- sum((g.l$v - gammaU) ^ 2)
if (g.l$weights == "npairs")
loss <- sum(g.l$n * (g.l$v - gammaU) ^ 2)
if (g.l$weights == "cressie")
loss <- sum((g.l$n / (gammaU ^ 2)) * (g.l$v - gammaU) ^ 2)
if (loss > (.Machine$double.xmax ^ 0.5) |
loss == Inf | loss == -Inf | is.nan(loss))
loss <- .Machine$double.xmax ^ 0.5
return(loss + penalty)
}
# Helper function to check validity of cov.model
.checkCMB.cov.model <-
function(cov.model,
cov.pars,
kappa,
env = NULL,
output = TRUE)
{
## extracting covariance parameters
if (is.vector(cov.pars))
sigmasq <- cov.pars[1]
else
sigmasq <- cov.pars[, 1]
if (is.vector(cov.pars))
phi <- cov.pars[2]
else
phi <- cov.pars[, 2]
if (missing(kappa) || is.null(kappa))
kappa <- NA
## checking for nested models
cov.pars <- drop(cov.pars)
if (is.vector(cov.pars))
ns <- 1
else{
ns <- nrow(cov.pars)
if (length(cov.model) == 1)
cov.model <- rep(cov.model, ns)
if (length(kappa) == 1)
kappa <- rep(kappa, ns)
}
if (length(cov.model) != ns)
stop("wrong length for cov.model")
##
cov.model <- sapply(cov.model, match.arg, CMB.cov.models)
cov.model[cov.model == "stable"] <- "powered.exponential"
if (any(cov.model == c("gneiting.matern", "gencauchy"))) {
if (length(kappa) != 2 * ns)
stop(
paste(
"wrong length for kappa, ",
cov.model,
"model requires two values for the argument kappa"
)
)
}
else{
if (length(kappa) != ns)
stop('wrong length for kappa')
}
## settings for power model (do not reverse order of the next two lines!)
phi[cov.model == "linear"] <- 1
cov.model[cov.model == "linear"] <- "power"
## checking input for cov. models with extra parameter(s)
if (any(cov.model == 'gneiting.matern') && ns > 1)
stop('nested models including the gneiting.matern are not implemented')
for (i in 1:ns) {
if (any(
cov.model[i] == c(
"matern",
"powered.exponential",
"cauchy",
"gneiting.matern",
"gencauchy",
"askey",
"c2wendland",
"c4wendland",
"sinepower",
"multiquadric"
)
)) {
if (any(cov.model[i] == c("gneiting.matern", "gencauchy"))) {
if (any(is.na(kappa)) || length(kappa) != 2 * ns)
stop(
paste(
cov.model[i],
"correlation function model requires a vector with 2 parameters in the argument kappa"
)
)
}
else{
if (is.na(kappa[i]) | is.null(kappa[i]))
stop(
"for matern, powered.exponential, gencauchy, askey, c2wendland,
c4wendland, sinepower, multiquadric and cauchy covariance
functions the parameter kappa must be provided"
)
}
if ((
cov.model[i] == "matern" | cov.model[i] == "powered.exponential" |
cov.model[i] == "askey" | cov.model[i] == "c2wendland" |
cov.model[i] == "c4wendland" |
cov.model[i] == "sinepower" |
cov.model[i] == "multiquadric"
) & length(kappa) != 1 * ns)
stop("kappa must have 1 parameter for this correlation function")
if (cov.model[i] == "matern" &
kappa[i] == 0.5)
cov.model[i] == "exponential"
}
if (cov.model[i] == "power")
if (any(phi[i] >= 2) | any(phi[i] <= 0))
stop("for power model the phi parameters must be in the interval ]0,2[")
}
if (!is.null(env)) {
assign("sigmasq", sigmasq, envir = env)
assign("phi", phi, envir = env)
assign("kappa", kappa, envir = env)
assign("ns", ns, envir = env)
assign("cov.model", cov.model, envir = env)
}
if (output)
return(list(
cov.model = cov.model,
sigmasq = sigmasq,
phi = phi,
kappa = kappa,
ns = ns
))
else
return(invisible())
}
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.