tests/testthat/test-rmvnorm.R

In spam: SPArse Matrix

# HEADER ####################################################
# This is file spam/tests/testthat/test-helper.R.           #
# It is part of the R package spam,                         #
#  --> https://CRAN.R-project.org/package=spam              #
#  --> https://CRAN.R-project.org/package=spam64            #
#  --> https://git.math.uzh.ch/reinhard.furrer/spam         #
# by Reinhard Furrer [aut, cre], Florian Gerber [aut],      #
#    Roman Flury [aut], Daniel Gerber [ctb],                #
#    Kaspar Moesinger [ctb]                                 #
# HEADER END ################################################

rm(list = ls())
source("helper.R")

## library("testthat")
## library("spam64", lib.loc = LIB.LOC)
## library("spam", lib.loc = "../../../lib/")


context("test-rmvnorm.R")



########
verbose <- FALSE

set.seed(14)


n <- 5
Sigma <- .25^abs(outer(1:n,1:n,'-'))
Q <- as.spam(solve(Sigma))
b <- 1:n

struct <- chol(Q, verbose=verbose)

test_that("rmvnorm", {

    set.seed(14)
    tmp1 <- rmvnorm(10,  Sigma=Sigma, verbose=verbose)
    set.seed(14)
    tmp2 <- rmvnorm.spam(10, Sigma=Sigma)
    expect_equal( tmp1, tmp2)

    set.seed(14)
    tmp1 <- rmvnorm(10,  Sigma=Q)
    set.seed(14)
    tmp2 <- rmvnorm.spam(10, Sigma=Q)
    expect_equal( tmp1, tmp2)



    set.seed(14)
    tmp1 <- rmvnorm.canonical(10, b, Q)
    set.seed(14)
    spamtest_eq( rmvnorm.canonical(10, b, Q, Lstruct=struct), tmp1 )
    set.seed(14)
    spamtest_eq( rmvnorm.prec(10, solve(Q,b), Q), tmp1 )
    set.seed(14)
    spamtest_eq( rmvnorm.prec(10, solve(Q,b), Q, Lstruct=struct), tmp1 )
    set.seed(14)
    ## cat("For rmvnorm.canonical:\n- comparing sample mean with truth:\n")
    ## for (i in 10^(1:4))
    ##     cat('    sample size n=',i,' yields  Frobenius-norm:',
    ##         norm( apply(rmvnorm.canonical(i, b, Q, Lstruct=struct), 2,mean)- solve(Q,b),'f'),'\n')
    ## cat("- comparing sample variance with truth:\n")
    ## for (i in 10^(1:4)){
    ##     cat('    sample size n=',i,' yields Frobenius-norm:',
    ##         norm( var( rmvnorm.canonical(i, b, Q=Q, Lstruct=struct))- Sigma,'f'),'\n')
    ## set.seed(14)
    ## cat("For rmvnorm.prec:\n- comparing sample mean with truth:\n")
    ## for (i in 10^(1:4))
    ##     cat('    sample size n=',i,' yields  Frobenius-norm:',
    ##         norm( apply(rmvnorm.prec(i, b, Q, Lstruct=struct), 2,mean)- b,'f'),'\n')
    ## cat("- comparing sample variance with truth:\n")
    ## for (i in 10^(1:4)){
    ##     cat('    sample size n=',i,' yields Frobenius-norm:',
    ##         norm( var( rmvnorm.prec(i, Q=Q, Lstruct=struct))- Sigma,'f'),'\n')
    ## }
    ## set.seed(14)
    ## cat("For rmvnorm.spam:\n- comparing sample mean with truth:\n")
    ## for (i in 10^(1:4))
    ##   cat('    sample size n=',i,' yields  Frobenius-norm:',
    ##       norm( apply(rmvnorm.spam(i, b, as.spam(Sigma), Lstruct=struct), 2,mean)- b,'f'),'\n')
    ## cat("- comparing sample variance with truth:\n")
    ## for (i in 10^(1:4)){
    ##   cat('    sample size n=',i,' yields Frobenius-norm:',
    ##       norm( var( rmvnorm.spam(i, b, as.spam(Sigma), Lstruct=struct))- Sigma,'f'),'\n')
    ## }
})


test_that("rmvt", {

  set.seed(14)
  tmp1 <- rmvt.spam(10,  Sigma=Sigma, df = 0)
  set.seed(14)
  tmp2 <- rmvnorm.spam(10, Sigma = Sigma)
  expect_equal( tmp1, tmp2)

  set.seed(14)
  tmp1 <- rmvt.spam(10, Sigma = Sigma, df = 4)
  set.seed(14)
  tmp2 <- rmvnorm.spam(10, Sigma = Sigma)/sqrt(rchisq(10, 4)/4)
  expect_equal( tmp1, tmp2)

  set.seed(14)
  tmp1 <- rmvt.spam(10, Sigma = Sigma, df = 4, type = "shifted")
  set.seed(14)
  tmp2 <- sweep(rmvnorm.spam(10, Sigma = Sigma)/sqrt(rchisq(10, df = 4)/4), 2, rep(0, nrow(Sigma)), "+")
  expect_equal( tmp1, tmp2)

  set.seed(14)
  tmp1 <- rmvt(10, Sigma = Sigma, df = 4, type = "shifted")
  set.seed(14)
  tmp2 <- rmvt.spam(10, Sigma = Sigma, df = 4, type = "shifted")
  expect_equal( tmp1, tmp2)
})



require('fields')

context("test-rgrf.R: locations")


tau <- 0
nx <- 49
ny <- nx
field <- rgrf(1, nx=nx, ny=ny, tau=tau, Covariance="cov.nug", theta=1)
quilt.plot(cbind(attr(field,"locs"),z=field), nx=nx, ny=ny)
points(attr(field,"locs"))

expect_error(rgrf(1, nx=3, tau=.5, Covariance="cov.nug", theta=1))

tau <- 0.499
nx <- 4
ny <- nx
field <- rgrf(1, nx=nx, ny=ny, tau=tau, Covariance="cov.nug", theta=c(.2, 1, 1.5))
plot(attr(field,"locs"))
quilt.plot(cbind(attr(field,"locs"),z=field), nx=nx, ny=ny)
points(attr(field,"locs"))

# unit grid:
nx <- 12
ny <- 7
field <- rgrf(1, nx=nx, ny=ny, xlim=c(0,nx), ylim=c(0,ny),  Covariance="cov.nug", theta=1)
quilt.plot(cbind(attr(field,"locs"),z=field), nx=nx, ny=ny)
points(attr(field,"locs"))

# effective tau:
tau <- 0.25
nx <- 12
ny <- 7
field <- rgrf(1, nx=nx, ny=ny, xlim=c(0,nx), ylim=c(0,ny), tau=tau,  Covariance="cov.nug", theta=1)
quilt.plot(cbind(attr(field,"locs"),z=field), nx=nx, ny=ny)
points(attr(field,"locs"))
attr(field,"eff")

tau <- 0.25
nx <- 12
ny <- 7
field <- rgrf(1, nx=nx, ny=ny, xlim=c(0,nx), ylim=c(0,ny), tau=tau,  Covariance="cov.nug", theta=1)
quilt.plot(cbind(round(attr(field,"locs")+.499)-.49,z=field), nx=nx, ny=ny)
points(attr(field,"locs"))
attr(field,"eff")

# min distance checking:  (within 0,1)
tau <- 0.499
nx <- 40
ny <- nx
field <- rgrf(1, nx=nx, ny=ny, tau=tau, Covariance="cov.nug", theta=c(.2, 1, 1.5))

tt <- rdist(attr(field,"locs"))
etau <- attr(field,"effective.tau")[1]

expect_lt( (1/nx-2*etau), min( tt[tt>0]))
expect_lt( max( tt[tt>0]), sqrt(2)*((nx-1)/nx+2*etau))



context("test-rgrf.R: mean")
tau <- 0.499
nx <- 20
ny <- nx
field <- rgrf(1, nx=nx, ny=ny, tau=tau, Covariance="cov.nug", theta=c(.2, 1, 1.5))

locs <- attr(field,"locs")
field <- rgrf(1, nx=nx, ny=ny, tau=tau, X=cbind(1,locs), beta=c(4,1,1), Covariance="cov.nug", theta=c(.2, 1, 1.5))

quilt.plot(cbind(attr(field,"locs"),z=field))



context("test-rgrf.R: samples")


nx <- 12
ny <- nx
set.seed(12)
field <- rgrf(1, nx=nx, ny=ny, Covariance="cov.sph", theta=c(1,.2,.2))

distmat <- rdist( attr(field,"locs"))

set.seed(12)
out <- rmvnorm(1, Sigma=cov.sph(distmat,theta=c(1,.2,.2)))

expect_equal( norm(field-out), 0)


expect_error( rgrf(1, nx=nx, ny=ny, Covariance="cov.sphd", theta=c(1,.2,.2)))
expect_error( rgrf(1, nx=nx, ny=ny, Covariance="cov.sphd", theta=c(1,.2,.2)))
expect_error( rgrf(1, nx=nx, ny=ny, Covariance="cov.sphd", theta=c(1,.2,.2)))
expect_error( rgrf(1, nx=nx, ny=ny, Covariance="cov.sphd", theta=c(1,.2,.2)))




if (F) {
    n <- 2000
    locs <- cbind(runif(n), runif(n))

    # rdist just beats the shit out of it...
#    microbenchmark::microbenchmark(
#                        a=rdist(locs),
#                        b=as.matrix(dist(locs)),
#                        c=nearest.dist(locs, upper=NULL, delta=.5),
#                        d=nearest.dist(locs, upper=NULL, delta=sqrt(2)),
#                        times=10)

    }