R Under development (unstable) (2024-11-03 r87286 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > #' > #' Header for all (concatenated) test files > #' > #' Require spatstat.univar > #' Obtain environment variable controlling tests. > #' > #' $Revision: 1.5 $ $Date: 2020/04/30 05:31:37 $ > > require(spatstat.univar) Loading required package: spatstat.univar spatstat.univar 3.1-1 > FULLTEST <- (nchar(Sys.getenv("SPATSTAT_TEST", unset="")) > 0) > ALWAYS <- TRUE > cat(paste("--------- Executing", + if(FULLTEST) "** ALL **" else "**RESTRICTED** subset of", + "test code -----------\n")) --------- Executing **RESTRICTED** subset of test code ----------- > # tests/weightedstats.R > # $Revision: 1.2 $ $Date: 2023/11/05 01:40:53 $ > > local({ + if(ALWAYS) { # depends on hardware + ## whist() + ## check agreement between C and interpreted code for whist() + set.seed(98123) + x <- runif(1000) + w <- sample(1:5, 1000, replace=TRUE) + b <- seq(0,1,length=101) + aC <- whist(x,b,w, method="C") + aR <- whist(x,b,w, method="interpreted") + if(!all(aC == aR)) + stop("Algorithms for whist disagree") + } + if(FULLTEST) { + ## cases of 'unnormdensity()' + x <- rnorm(20) + d0 <- unnormdensity(x, weights=rep(0, 20)) + dneg <- unnormdensity(x, weights=c(-runif(19), 0)) + } + + }) > > #' > #' tests/parzen.R > #' > #' Tests of the Parzen-Rosenblatt estimator > #' (fixed bandwidth, no boundary correction) > #' > #' $Revision: 1.1 $ $Date: 2023/10/22 02:39:49 $ > > local({ + if(FULLTEST) { + #' code in kernels.R + kernames <- c("gaussian", "rectangular", "triangular", + "epanechnikov", "biweight", "cosine", "optcosine") + X <- rnorm(20) + U <- runif(20) + for(ker in kernames) { + dX <- dkernel(X, ker) + fX <- pkernel(X, ker) + qU <- qkernel(U, ker) + m0 <- kernel.moment(0, 0, ker) + m1 <- kernel.moment(1, 0, ker) + m2 <- kernel.moment(2, 0, ker) + m3 <- kernel.moment(3, 0, ker) + fa <- kernel.factor(ker) + sq <- kernel.squint(ker) + } + } + }) > > local({ + if(ALWAYS) { + ## unnormdensity + x <- rnorm(20) + d0 <- unnormdensity(x, weights=rep(0, 20)) + dneg <- unnormdensity(x, weights=c(-runif(19), 0)) + } + }) > # > # tests/NAinCov.R > # > # Testing the response to the presence of NA's in covariates > # > # $Revision: 1.10 $ $Date: 2024/09/30 23:13:54 $ > > if(FULLTEST) { + local({ + #' quantile.ewcdf + f <- ewcdf(runif(100), runif(100)) + qf <- quantile(f, probs=c(0.1, NA, 0.8)) + #' quantile.density + f <- density(runif(100)) + qf <- quantile(f, probs=c(0.1, NA, 0.8)) + }) + } > > #' tests/direct.R > #' > #' Check output of densityBC() by comparing with > #' kernel estimates computed directly in R code. > #' > #' $Revision: 1.2 $ $Date: 2023/10/22 02:39:35 $ > > INTERACTIVE <- FALSE > > #' ...... direct implementation ........................ > #' > #' Biweight boundary kernel (for estimation at point r) > > bdry.bwt <- function(x,r,h=1){ + u <- x/h + k <- (15/(16*h))*(1-u^2)^2*(u^2<1) + p <- r/h + p[p>1] <- 1 + a0 <- (3*p^5 - 10*p^3 + 15*p + 8)/16 + a1 <- (5*p^6 - 15*p^4 + 15*p^2 -5)/32 + a2 <- (15*p^7 - 42*p^5 + 35*p^3 + 8)/112 + bk <- (a2-a1*u)*k*(u > #' Standard biweight kernel > > bwt <- function(x,r,h=1){ + # r is ignored + u <- x/h + k <- (15/(16*h))*(1-u^2)^2*(u^2<1) + k + } > > #' Kernel estimate > > kernR <- function(x, h, kname="bwt", from, to, nr=200) { + ker <- get(kname, mode="function") + fhat <- numeric(0) + rvalues <- seq(from, to, length=nr) + for (r in rvalues) + fhat <- c(fhat,mean(ker(r-x, r, h=h))) + return(list(x=rvalues, y=fhat)) + } > > #' ...................... RUN EXAMPLE ....................... > > sim.dat <- rexp(500) > > from <- 0.01 > to <- 5 > nr <- 200 > > if(INTERACTIVE) opa <- par(ask=TRUE) > > cat("--- Fixed bandwidth ---\n") --- Fixed bandwidth --- > > fhatR <- kernR(sim.dat, 0.4, "bwt", from=from, to=to, nr=nr) > fhatC <- densityBC(sim.dat, "biweight", h=0.4, from=from, to=to, n=nr) > > stopifnot(length(fhatR$x) == length(fhatC$x)) > stopifnot(all(fhatR$x == fhatC$x)) > rvalues <- fhatR$x > ftrue <- dexp(rvalues) > > if(INTERACTIVE) { + plot(c(from, to), c(0, 1.1), type="n", xlab="r", ylab="density", + main="Fixed bandwidth") + lines(rvalues, fhatR$y, col=1) + lines(rvalues, fhatC$x, col=2) + lines(rvalues, ftrue, col=3) + } > > cat("Range of absolute discrepancies between estimates (fixed bandwidth):\n") Range of absolute discrepancies between estimates (fixed bandwidth): > print(range(fhatC$y - fhatR$y)) [1] -6.661338e-16 6.661338e-16 > > cat("\n--- Boundary kernel ---\n") --- Boundary kernel --- > > fhat.R <- kernR(sim.dat, 0.4, "bdry.bwt", from=from, to=to, nr=nr) > fhat.C <- densityBC(sim.dat, "biweight", h=0.4, from=from, to=to, n=nr, + zerocor="bdrykern") > > stopifnot(length(fhat.R$x) == length(fhat.C$x)) > stopifnot(all(fhat.R$x == fhat.C$x)) > > rvalues <- fhat.R$x > ftrue <- dexp(rvalues) > > if(INTERACTIVE) { + plot(c(from, to), c(0, 1.1), type="n", xlab="r", ylab="density", + main="Boundary kernel") + lines(rvalues, fhat.R$y, col=1) + lines(rvalues, fhat.C$y, col=2) + lines(rvalues, ftrue, col=3) + } > > cat("Range of absolute discrepancies between estimates (boundary kernel):\n") Range of absolute discrepancies between estimates (boundary kernel): > print(range(fhat.C$y - fhat.R$y)) [1] -6.661338e-16 3.330669e-16 > > rel.discrep <- (fhat.C$y - fhat.R$y)/ftrue > rel.discrep <- rel.discrep[ftrue > 0.1] > cat("Range of relative discrepancies between estimates (boundary kernel):\n") Range of relative discrepancies between estimates (boundary kernel): > print(range(rel.discrep)) [1] -8.865342e-16 6.911665e-16 > if(max(abs(rel.discrep)) > 0.01) + stop("Relative discrepancies between C and R code exceed 1 percent") > > if(INTERACTIVE) par(opa) > > #' tests/kermom.R > #' > #' Check R function kernel.moment() against C function kermom () > #' > #' $Revision: 1.2 $ $Date: 2024/10/31 10:36:27 $ > #' > > moo <- 1 > sdee <- 0.5 > xx <- seq(moo - 4 * sdee, moo + 4 * sdee, length=512) > > kernames <- c("gaussian", "rectangular", "triangular", + "epanechnikov", "biweight", "cosine", "optcosine") > > eps <- sqrt(.Machine$double.eps) > > for(m in 0:2) { + cat("Incomplete moment of order", m, fill=TRUE) + for(ker in kernames) { + Rvalues <- kernel.moment(m, xx, ker, mean=moo, sd=sdee) + Cvalues <- kermom(m, xx, ker, mean=moo, sd=sdee) + discrep <- max(abs(Rvalues-Cvalues)) + if(discrep > eps) + stop(paste("kernel.moment and kermom disagree", + "for m =", m, "for kernel", sQuote(ker), + "\n\tDiscrepancy", discrep)) + cat("Kernel", sQuote(ker), "\tdiscrepancy", discrep, fill=TRUE) + } + } Incomplete moment of order 0 Kernel 'gaussian' discrepancy 0 Kernel 'rectangular' discrepancy 0 Kernel 'triangular' discrepancy 2.220446e-16 Kernel 'epanechnikov' discrepancy 1.110223e-16 Kernel 'biweight' discrepancy 1.110223e-16 Kernel 'cosine' discrepancy 0 Kernel 'optcosine' discrepancy 0 Incomplete moment of order 1 Kernel 'gaussian' discrepancy 0 Kernel 'rectangular' discrepancy 0 Kernel 'triangular' discrepancy 3.330669e-16 Kernel 'epanechnikov' discrepancy 2.220446e-16 Kernel 'biweight' discrepancy 2.220446e-16 Kernel 'cosine' discrepancy 2.333721e-17 Kernel 'optcosine' discrepancy 0 Incomplete moment of order 2 Kernel 'gaussian' discrepancy 0 Kernel 'rectangular' discrepancy 2.220446e-16 Kernel 'triangular' discrepancy 4.440892e-16 Kernel 'epanechnikov' discrepancy 4.440892e-16 Kernel 'biweight' discrepancy 4.440892e-16 Kernel 'cosine' discrepancy 2.220446e-16 Kernel 'optcosine' discrepancy 5.938278e-17 > > proc.time() user system elapsed 0.35 0.07 0.42