inst/tests/test-mcmc-normal.R

# Part of the "structmcmc" package, https://github.com/rjbgoudie/structmcmc
#
# This software is distributed under the GPL-3 license.  It is free,
# open source, and has the attribution requirements (GPL Section 7) in
#   https://github.com/rjbgoudie/structmcmc
#
# Note that it is required that attributions are retained with each function.
#
# Copyright 2008 Robert J. B. Goudie, University of Warwick

context("MCMC Normal Fast")

test_that("Two node, very simple", {
  x1 <- c(-10, -2)
  x2 <- c(-11, -4)
  data <- cbind(x1, x2)

  initial <- bn(integer(0), integer(0))

  sampler <- BNSampler(data = data,
                       initial = initial,
                       prior = priorUniform(initial),
                       logScoreFUN = list(offline = logScoreNormalOffline,
                                          online  = logScoreNormalIncremental,
                                          prepare = logScoreNormalPrepare))

  samples <- draw(sampler, 2500, burnin = 1000, verbose = F)
  bnp <- bnpostmcmc(sampler, samples)

  set.seed(2121)
  fam <- enumerateBNSpace(2)
  scores <- logScoreNormal(fam, data)

  exact <- bnpost(bnspace     = fam,
                  logScore    = scores,
                  data        = data,
                  logScoreFUN = logScoreNormal)

  expect_that(max(ep(bnp) - ep(exact)) < 0.01, is_true())
})

test_that("Steven Hill's example", {
  dat <- matrix(c(
    0.814723686393179000000000000000,0.905791937075619000000000000000,
    1.741378976923020000000000000000,
    0.913375856139019000000000000000,0.632359246225410000000000000000,
    1.548217823880720000000000000000,
    0.278498218867048000000000000000,0.546881519204984000000000000000,
    0.870515543567796000000000000000,
    0.964888535199277000000000000000,0.157613081677548000000000000000,
    1.169740976362910000000000000000,
    0.957166948242946000000000000000,0.485375648722841000000000000000,
    1.467085801589190000000000000000,
    0.141886338627215000000000000000,0.421761282626275000000000000000,
    0.588110253173491000000000000000,
    0.792207329559554000000000000000,0.959492426392903000000000000000,
    1.768585726443530000000000000000,
    0.035711678574189600000000000000,0.849129305868777000000000000000,
    0.929843676763850000000000000000,
    0.678735154857773000000000000000,0.757740130578333000000000000000,
    1.454937624492120000000000000000,
    0.392227019534168000000000000000,0.655477890177557000000000000000,
    1.053265047476410000000000000000,
    0.706046088019609000000000000000,0.031832846377420700000000000000,
    0.776891537813086000000000000000,
    0.046171390631153900000000000000,0.097131781235847500000000000000,
    0.162790113715064000000000000000,
    0.694828622975817000000000000000,0.317099480060861000000000000000,
    1.024012667332370000000000000000,
    0.034446080502908800000000000000,0.438744359656398000000000000000,
    0.493386047438713000000000000000,
    0.765516788149002000000000000000,0.795199901137063000000000000000,
    1.565539415544490000000000000000,
    0.489764395788231000000000000000,0.445586200710899000000000000000,
    0.941949261129447000000000000000,
    0.709364830858073000000000000000,0.754686681982361000000000000000,
    1.511154042379210000000000000000,
    0.679702676853675000000000000000,0.655098003973841000000000000000,
    1.382607407839010000000000000000,
    0.118997681558377000000000000000,0.498364051982143000000000000000,
    0.646122163294443000000000000000,
    0.340385726666133000000000000000,0.585267750979777000000000000000,
    0.928642454793268000000000000000,
    0.751267059305653000000000000000,0.255095115459269000000000000000,
    1.018101170433540000000000000000,
    0.699076722656686000000000000000,0.890903252535798000000000000000,
    1.607637903753590000000000000000,
    0.547215529963803000000000000000,0.138624442828679000000000000000,
    0.726899674802380000000000000000,
    0.257508254123736000000000000000,0.840717255983663000000000000000,
    1.098995681989980000000000000000,
    0.814284826068816000000000000000,0.243524968724989000000000000000,
    1.059960984876700000000000000000,
    0.349983765984809000000000000000,0.196595250431208000000000000000,
    0.555028517889152000000000000000,
    0.616044676146639000000000000000,0.473288848902729000000000000000,
    1.121789298797190000000000000000,
    0.830828627896291000000000000000,0.585264091152724000000000000000,
    1.452678838331950000000000000000,
    0.917193663829810000000000000000,0.285839018820374000000000000000,
    1.235419980807000000000000000000,
    0.753729094278495000000000000000,0.380445846975357000000000000000,
    1.156721126575400000000000000000,
    0.075854289563063600000000000000,0.053950118666607100000000000000,
    0.157154852843988000000000000000,
    0.779167230102011000000000000000,0.934010684229183000000000000000,
    1.727993954611580000000000000000,
    0.568823660872193000000000000000,0.469390641058206000000000000000,
    1.075448942284110000000000000000,
    0.337122644398882000000000000000,0.162182308193243000000000000000,
    0.508752703343751000000000000000,
    0.311215042044805000000000000000,0.528533135506213000000000000000,
    0.874086949219283000000000000000,
    0.601981941401637000000000000000,0.262971284540144000000000000000,
    0.874128783728644000000000000000,
    0.689214503140008000000000000000,0.748151592823709000000000000000,
    1.455790325788230000000000000000,
    0.083821377996932600000000000000,0.228976968716819000000000000000,
    0.344079274750236000000000000000,
    0.152378018969223000000000000000,0.825816977489547000000000000000,
    1.017206368216340000000000000000,
    0.996134716626885000000000000000,0.078175528753183700000000000000,
    1.078366533823360000000000000000,
    0.106652770180584000000000000000,0.961898080855054000000000000000,
    1.115020149584070000000000000000,
    0.774910464711502000000000000000,0.817303220653433000000000000000,
    1.630999319295360000000000000000,
    0.084435845510910300000000000000,0.399782649098896000000000000000,
    0.508558076229965000000000000000,
    0.800068480224308000000000000000,0.431413827463545000000000000000,
    1.253275237116900000000000000000,
    0.181847028302852000000000000000,0.263802916521990000000000000000,
    0.467989132296333000000000000000,
    0.136068558708664000000000000000,0.869292207640089000000000000000,
    1.020678239949580000000000000000,
    0.549860201836332000000000000000,0.144954798223727000000000000000,
    0.720240432829115000000000000000,
    0.622055131485066000000000000000,0.350952380892271000000000000000,
    0.998546090585942000000000000000,
    0.401808033751942000000000000000,0.075966691690841900000000000000,
    0.518656110858897000000000000000,
    0.123318934835166000000000000000,0.183907788282417000000000000000,
    0.346968293961755000000000000000),
    ncol = 3, byrow = T)

  fam <- enumerateBNSpace(3)
  scores <- logScoreNormal(fam, dat)

  priors <- 1
  scores <- scores - max(scores)
  expected <- exp(scores)*priors/sum(exp(scores)*priors)

  numberOfBurnIn <- 100
  numberOfSamples <- 5000

  expectedTable <- data.frame(expected = round(expected * numberOfSamples))
  row.names(expectedTable) <- lapply(fam, function(network){
    paste(network, sep = "", collapse = ",")})

  initial <- bn(integer(0), integer(0), integer(0))

  set.seed(2221)
  sampler <- BNSampler(data = dat,
                       initial = initial,
                       prior = priorUniform(initial),
                       logScoreFUN = list(offline = logScoreNormalOffline,
                                          online  = logScoreNormalIncremental,
                                          prepare = logScoreNormalPrepare))
  samples <- lapply(seq_len(numberOfBurnIn), sampler)
  samples <- lapply(seq_len(numberOfSamples), sampler)

  outTable <- table(factor(unlist(lapply(samples,function(l){
    paste(l,sep = "",collapse = ",")}))))

  expect_that(as.vector(outTable["integer(0),c(1,3),1"]),
            is_within(894, 125))
  expect_that(as.vector(outTable["3,c(1,3),integer(0)"]),
              is_within(3428, 150))
  expect_that(as.vector(outTable["2:3,integer(0),2"]),
              is_within(162, 125))
})
rjbgoudie/structmcmc documentation built on Nov. 3, 2020, 3:41 a.m.