tests/abc.R

In SimInf: A Framework for Data-Driven Stochastic Disease Spread Simulations

## This file is part of SimInf, a framework for stochastic
## disease spread simulations.
##
## Copyright (C) 2015 -- 2022 Stefan Widgren
##
## SimInf is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## SimInf is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <https://www.gnu.org/licenses/>.

library(SimInf)
library(tools)
source("util/check.R")

## Define a tolerance
tol <- 1e-8

##
## Generate proposals without particles from the previous generation.
##

proposals_exp <- matrix(c(0.788305135443807,    0.4089769218117,
                          0.883017404004931,    0.940467284293845,
                          0.0455564993899316,   0.528105488047004,
                          0.892419044394046,    0.551435014465824,
                          0.456614735303447,    0.956833345349878,
                          0.453334156190977,    0.677570635452867,
                          0.572633401956409,    0.102924682665616,
                          0.899824970401824,    0.24608773435466,
                          0.0420595335308462,   0.327920719282702,
                          0.954503649147227,    0.889539316063747,
                          0.6928034061566,      0.640506813768297,
                          0.994269776623696,    0.655705799115822,
                          0.708530468167737,    0.544066024711356,
                          0.59414202044718,     0.28915973729454,
                          0.147113647311926,    0.963024232536554,
                          0.902299045119435,    0.690705278422683,
                          0.795467417687178,    0.0246136845089495,
                          0.477795971091837,    0.758459537522867,
                          0.216407935833558,    0.318181007634848,
                          0.231625785352662,    0.142800022382289,
                          0.414546335814521,    0.413724326295778,
                          0.368845450924709,    0.152444747742265,
                          0.13880606344901,     0.233034099452198,
                          0.465962450252846,    0.265972640365362,
                          0.857827715342864,    0.0458311666734517,
                          0.442200074205175,    0.798924845643342,
                          0.12189925997518,     0.560947983758524,
                          0.20653138961643,     0.127531650243327,
                          0.753307864302769,    0.895045359153301,
                          0.374462775886059,    0.665115194628015,
                          0.0948406609240919,   0.383969637798145,
                          0.27438364457339,     0.814640038879588,
                          0.448516341391951,    0.810064353048801,
                          0.812389509519562,    0.794342321110889,
                          0.439831687603146,    0.754475158639252,
                          0.629221131559461,    0.710182401351631,
                          0.000624773325398564, 0.475316574098542,
                          0.220118885161355,    0.379816537722945,
                          0.612771003274247,    0.351797909243032,
                          0.111135424347594,    0.243619472719729,
                          0.66805558744818,     0.417646779678762,
                          0.788195834029466,    0.102864644257352,
                          0.434892741497606,    0.984956979984418,
                          0.893051114398986,    0.886469060787931,
                          0.175052650272846,    0.130695691565052,
                          0.653101925039664,    0.343516472261399,
                          0.656758127966896,    0.320373242488131,
                          0.187691119266674,    0.782294301316142,
                          0.0935949867125601,   0.46677904156968,
                          0.511505459900945,    0.59998895926401,
                          0.332823540316895,    0.488613033667207,
                          0.954473827499896,    0.4829023971688,
                          0.890350222121924,    0.914438186911866,
                          0.608734982321039,    0.410689776530489,
                          0.147094690939412,    0.935299803270027,
                          0.301228899974376,    0.0607205715496093,
                          0.947726940037683,    0.720596273429692,
                          0.142294295597821,    0.549284656066447,
                          0.95409123855643,     0.585483353119344,
                          0.404510281747207,    0.64789347932674,
                          0.319820617092773,    0.307720010867342,
                          0.219767631264403,    0.369488865835592,
                          0.984219203470275,    0.154202300822362,
                          0.0910439998842776,   0.141906907781959,
                          0.690007101511583,    0.619256483390927,
                          0.891394117148593,    0.672999092610553,
                          0.737077737925574,    0.521135725779459,
                          0.65983844967559,     0.821805460145697,
                          0.786281551700085,    0.979821917368099,
                          0.439431536244228,    0.31170220207423,
                          0.409474952612072,    0.0104671118315309,
                          0.18384952400811,     0.842729318886995,
                          0.231161782052368,    0.239099955651909,
                          0.0766911653336138,   0.245723678031936,
                          0.73213520552963,     0.847453165100887,
                          0.497527267085388,    0.387909029843286,
                          0.246448994148523,    0.111096461303532,
                          0.389994435245171,    0.571935313986614,
                          0.216892762808129,    0.44476800202392,
                          0.21799066872336,     0.502299563260749,
                          0.353904571849853,    0.64998515881598,
                          0.374713956611231,    0.355445380788296,
                          0.533687945455313,    0.740334360394627,
                          0.221102937823161,    0.41274611861445,
                          0.265686686849222,    0.629973053466529,
                          0.183828490786254,    0.86364411143586,
                          0.746568004135042,    0.668284649727866,
                          0.618017873261124,    0.372238060226664,
                          0.529835685854778,    0.874682342866436,
                          0.581750099780038,    0.839767764788121,
                          0.312448164913803,    0.70829032221809,
                          0.265017806086689,    0.594343194039539,
                          0.481289800489321,    0.26503273146227,
                          0.564590434776619,    0.913188223028556,
                          0.901874389499426,    0.274166621500626,
                          0.321482756407931,    0.985640884377062,
                          0.619993310188875,    0.937314089154825,
                          0.466532702324912,    0.406832593260333,
                          0.659230324206874,    0.152346616843715,
                          0.572867058217525,    0.238726026844233),
                        nrow = 100, ncol = 2, byrow = TRUE,
                        dimnames = list(NULL, c("beta", "gamma")))

attr(proposals_exp, "ancestor") <- rep(NA_integer_, 100)

set.seed(123)
proposals_obs <-
    .Call(SimInf:::SimInf_abc_proposals, ## function
          c("beta", "gamma"),            ## parameter
          c("uniform", "uniform"),       ## distribution
          c(0, 0),                       ## p1
          c(1, 1),                       ## p2
          100L,                          ## n
          NULL,                          ## x
          NULL,                          ## w
          NULL)                          ## sigma

stopifnot(identical(dim(proposals_obs), dim(proposals_exp)))
stopifnot(identical(dimnames(proposals_obs), dimnames(proposals_exp)))
stopifnot(all(abs(proposals_obs - proposals_exp) < tol))
stopifnot(identical(attr(proposals_obs, "ancestor"),
                    attr(proposals_exp, "ancestor")))

##
## Determine weights for the proposals
##

w_exp <- rep(0.01, 100)

w_obs <- .Call(SimInf:::SimInf_abc_weights, ## function
               c("uniform", "uniform"),     ## distribution
               c(0, 0),                     ## p1
               c(1, 1),                     ## p2
               NULL,                        ## x
               proposals_obs,               ## xx
               NULL,                        ## w
               NULL)                        ## sigma

stopifnot(identical(length(w_obs), length(w_exp)))
stopifnot(all(abs(w_obs - w_exp) < tol))

##
## Generate proposals with particles from a previous generation.
##

proposals_exp <- matrix(c(0.40449021397453,    0.209258724620122,
                          0.36482482505885,    0.177208382897285,
                          0.110673292350425,   0.0794986914792483,
                          0.187167868494588,   0.138234511613991,
                          0.0129193391123231,  0.0535672393807684,
                          0.371875050064337,   0.134845566874378,
                          0.247248031392758,   0.126673609970299,
                          0.183642149340511,   0.0731331396573393,
                          0.20911740588746,    0.0948903274398692,
                          0.262342094567642,   0.125179104685734,
                          0.177528318171729,   0.0741054597654854,
                          0.156151696536369,   0.090429058755398,
                          0.2878964747592,     0.144758670272637,
                          0.115570346272429,   0.096959225400175,
                          0.35527286855174,    0.161557776257793,
                          0.0650153072101651,  0.00255140678233257,
                          0.372258330220992,   0.133720658379678,
                          0.402004557396214,   0.190127049086695,
                          0.288281903108544,   0.150967395091388,
                          0.177909054629067,   0.0998662267148997,
                          0.336508383404385,   0.193255960409785,
                          0.493830751357006,   0.235285023981504,
                          0.150045602591045,   0.0978408919339508,
                          0.0991368666223921,  0.114314801264312,
                          0.218399403384471,   0.0989904109886274,
                          0.185099330012481,   0.0892108392300413,
                          0.161944605578993,   0.0604990966997036,
                          0.222789306144925,   0.151835826025647,
                          0.065428159161678,   0.0629737938634062,
                          0.23615906645682,    0.102537558499362,
                          0.0121276214604985,  0.0628933335766671,
                          0.201775823307792,   0.126217268774852,
                          0.137263474185558,   0.112391434285224,
                          0.264624142608849,   0.169454608095463,
                          0.312728091948025,   0.0963945630588326,
                          0.175526458577421,   0.0686229881214495,
                          0.0593637021715585,  0.0522592803469739,
                          0.215230839090112,   0.155696377876435,
                          0.0529094844818739,  0.0320773020404419,
                          0.33619295761988,    0.179355314279367,
                          0.181477854077582,   0.10047336472656,
                          0.21069018904093,    0.14241445146947,
                          0.139211751314129,   0.0689293727497798,
                          0.352694196052565,   0.16279166759017,
                          0.226482666581105,   0.126181747015655,
                          0.269213756712608,   0.134008000696826,
                          0.0819893183892228,  0.0804205926286409,
                          0.316344402602474,   0.141277619454973,
                          0.0415392652110742,  0.064295265585175,
                          0.0578856463374477,  0.087630895538036,
                          0.200832741112366,   0.096648847157614,
                          0.148811803643498,   0.0750742362153645,
                          0.191379030672455,   0.115737154467271,
                          0.161899095545659,   0.0750730031877879,
                          0.309156788420073,   0.168264661165334,
                          0.0859624973882138,  0.030587660584528,
                          0.0603779474270183,  0.0454745455018321,
                          0.350172184371153,   0.176682327720304,
                          0.122998229064858,   0.0882698562573068,
                          0.241334749076453,   0.147564901163679,
                          0.0766590044181838,  0.0434986125841068,
                          0.348278201059334,   0.191013869028047,
                          0.257268077430999,   0.179321515928957,
                          0.223064948567959,   0.0592622522014079,
                          0.427815559979751,   0.175047967508291,
                          0.235672458510888,   0.121165312322979,
                          0.167714899755654,   0.0903647430326278,
                          0.0915963546018915,  0.0651409817705193,
                          0.218331597808469,   0.104343950311149,
                          0.280915332879878,   0.126599331161323,
                          0.106748273493958,   0.0870869332381226,
                          0.267450910188283,   0.165256775670049,
                          0.0903771196756504,  0.0498756856503441,
                          0.258634833807708,   0.134139700799639,
                          0.0295659758705841,  0.0775270276812072,
                          0.227391392219538,   0.0816484593169099,
                          0.130314347249228,   0.029628912819613,
                          0.159126993527458,   0.0856407681638652,
                          0.290396981413321,   0.143042180482158,
                          0.0373833335256739,  0.0680495387444054,
                          0.302136009179183,   0.161372082543716,
                          0.442464566486759,   0.197717305508349,
                          0.220956904210761,   0.132325899715497,
                          0.383435091126591,   0.189195961844464,
                          0.00381862106974051, 0.0320324001946579,
                          0.193629157284585,   0.131457415291685,
                          0.124570282588458,   0.114510381022667,
                          0.247574724494243,   0.146305479285253,
                          0.135945749855924,   0.0830472764106766,
                          0.221986437891198,   0.161642051632531,
                          0.101839500171702,   0.0763136426325731,
                          0.145580556982575,   0.0920176018942565,
                          0.0931858029996674,  0.0538420972542555,
                          0.224730087525381,   0.105317389453258,
                          0.225838257380999,   0.133549398107083,
                          0.306352939853714,   0.188866355047975,
                          0.247124496788588,   0.125364413806328,
                          0.21457415007365,    0.0886153625185365,
                          0.234312573547032,   0.101556928688184,
                          0.352387419464283,   0.16910808806427),
                        nrow = 100, ncol = 2, byrow = TRUE,
                        dimnames = list(NULL, c("beta", "gamma")))

attr(proposals_exp, "ancestor") <-
    c(6L, 3L, 9L, 9L, 2L, 6L, 9L, 6L, 3L, 9L, 3L, 6L, 6L, 2L, 9L, 3L,
      3L, 9L, 9L, 6L, 6L, 9L, 6L, 6L, 6L, 6L, 3L, 2L, 9L, 9L, 6L, 6L,
      2L, 3L, 3L, 3L, 3L, 2L, 3L, 3L, 3L, 2L, 3L, 3L, 3L, 9L, 2L, 3L,
      6L, 2L, 6L, 3L, 2L, 6L, 3L, 6L, 2L, 3L, 3L, 9L, 3L, 9L, 9L, 6L,
      3L, 6L, 6L, 6L, 2L, 3L, 3L, 3L, 6L, 3L, 2L, 3L, 6L, 3L, 6L, 2L,
      3L, 9L, 9L, 9L, 2L, 2L, 6L, 3L, 6L, 3L, 2L, 6L, 2L, 3L, 6L, 6L,
      3L, 3L, 9L, 3L)

sigma <- matrix(c(0.0111588476485889,  0.00441287501075721,
                  0.00441287501075721, 0.00237388735640902),
                nrow = 2, ncol = 2, byrow = TRUE,
                dimnames = list(c("beta", "gamma"), c("beta", "gamma")))

w <- c(0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1)

x <- matrix(c(0.100721582304686,  0.0719040972180665,
              0.120237128576264,  0.0883633010089397,
              0.241502250079066,  0.118690322386101,
              0.196719169151038,  0.111476795747876,
              0.0473383399657905, 0.0265292769763619,
              0.142388037405908,  0.0744408417958766,
              0.118218665011227,  0.103439392056316,
              0.246003416366875,  0.0970179655123502,
              0.277967476984486,  0.150379287777469,
              0.126839369302616, 0.0596327660605311),
            nrow = 10, ncol = 2, byrow = TRUE,
            dimnames = list(NULL, c("beta", "gamma")))

set.seed(123)
proposals_obs <-
    .Call(SimInf:::SimInf_abc_proposals, ## function
          c("beta", "gamma"),            ## parameter
          c("uniform", "uniform"),       ## distribution
          c(0, 0),                       ## p1
          c(1, 1),                       ## p2
          100L,                          ## n
          x,                             ## x
          w,                             ## w
          sigma)                         ## sigma

stopifnot(identical(dim(proposals_obs), dim(proposals_exp)))
stopifnot(identical(dimnames(proposals_obs), dimnames(proposals_exp)))
stopifnot(all(abs(proposals_obs - proposals_exp) < tol))
stopifnot(identical(attr(proposals_obs, "ancestor"),
                    attr(proposals_exp, "ancestor")))

##
## Determine particle weights with particles from a previous
## generation.
##

w <- c(0.0887036508782177, 0.166359113528729, 0.0837072024175504,
       0.0744001632010883, 0.135889593166249, 0.0855856292642815,
       0.0982428712389774, 0.108701885332761, 0.0751630561589793,
       0.0832468348131658)

sigma <- matrix(c(0.00492214698840636, 0.00214848981388495,
                  0.00214848981388495, 0.00148131946857759),
                nrow = 2, ncol = 2, byrow = TRUE,
                dimnames = list(c("beta", "gamma"), c("beta", "gamma")))

x <- matrix(c(0.181656674010934, 0.088313908707635,
              0.13237129296898,  0.0653391855550409,
              0.13237129296898,  0.0653391855550409,
              0.13237129296898,  0.0653391855550409,
              0.164534995786719, 0.102338963629853,
              0.253254773828567, 0.0945404905257672,
              0.164534995786719, 0.102338963629853,
              0.181656674010934, 0.088313908707635,
              0.164534995786719, 0.102338963629853,
              0.253254773828567, 0.0945404905257672),
            nrow = 10, ncol = 2, byrow = TRUE,
            dimnames = list(NULL, c("beta", "gamma")))

xx <- matrix(c(0.207941308990465,  0.110336421309656,
               0.123214899415275,  0.0707943094184964,
               0.106876498007091,  0.0521419601364841,
               0.0817667562688388, 0.0384318357782007,
               0.13190369503059,   0.0784629457261605,
               0.231164237650949,  0.0977418120716123,
               0.187665128906516,  0.0899812114257998,
               0.17464412045458,   0.0826845340300653,
               0.178604586453451,  0.120632356669452,
               0.208363762292826, 0.105931947092043),
             nrow = 10, ncol = 2, byrow = TRUE,
             dimnames = list(NULL, c("beta", "gamma")))

w_exp <- c(0.0923963866851073, 0.0836785669165191, 0.100309975935471,
           0.139249040595048, 0.08154757535186, 0.122518476357057,
           0.0807798908697413, 0.0785544941015971, 0.13056926565214,
           0.0903963275354594)

w_obs <- .Call(SimInf:::SimInf_abc_weights, ## function
               c("uniform", "uniform"),     ## distribution
               c(0, 0),                     ## p1
               c(1, 1),                     ## p2
               x,                           ## x
               xx,                          ## xx
               w,                           ## w
               sigma)                       ## sigma

stopifnot(identical(length(w_obs), length(w_exp)))
stopifnot(all(abs(w_obs - w_exp) < tol))

##
## Check that an invalid 'distribution' is detected.
##
res <- assertError(
    .Call(SimInf:::SimInf_abc_weights, ## function
          "a",                         ## distribution (invalid)
          0,                           ## p1
          0,                           ## p2
          x,                           ## x
          xx,                          ## xx
          w,                           ## w
          sigma))                      ## sigma
check_error(res, "Unknown distribution.")

##
## Check that it works to use a gamma distribution for the proposals.
##
set.seed(123)
proposals_exp <- matrix(
    c(121.276958026798, 135.226111526571, 50.0608928839561,
      99.026703921349, 155.215472878808, 109.737500944335,
      60.0091438721518, 49.2206033058395, 136.474675481633,
      106.413871945226, 107.754150099189, 98.4420970793161,
      78.6402240451908, 139.098313398995, 122.148636952499,
      93.29139417767, 80.9864194166459, 88.4630763817778,
      66.0081318785976, 92.3853209979804, 50.1249564939848,
      62.4717320740012, 63.158652859414, 85.3912889151501,
      86.1229540420162, 91.0531624988248, 122.010318305364,
      117.410881625831, 112.839950617046, 93.1013355716318,
      104.037295209279, 61.1054179733021, 88.6996330858316,
      90.0144916130854, 135.879819760401, 63.5369008383957,
      107.848052867766, 69.6426105934016, 92.447768799244,
      82.1477056461401, 93.6832384540696, 103.698087374237,
      65.4192775797346, 53.2608207538562, 64.8808971531191,
      101.772343053721, 107.0615900602, 80.1481734870256,
      172.820347168567, 58.3237080596067, 104.585705127255,
      109.316983192402, 96.6407236074985, 125.552635673176,
      168.694971762368, 80.4681820354317, 69.2642980395984,
      74.3983787452227, 115.046380711507, 94.8090527724932,
      61.1003325370593, 72.4062993620932, 95.1777472059826,
      107.246241967593, 83.9189626963106, 72.7904320016698,
      105.501405100422, 131.814488410219, 108.886651704555,
      85.2196923468244, 101.467307888719, 80.5279999858835,
      102.500688357744, 76.6323000623859, 139.006679836589,
      81.2494406295684, 148.110470354515, 87.8741136774812,
      100.7456878299, 157.734303060981, 96.1650133533107,
      134.712273479757, 78.5646757768523, 93.6172210533339,
      109.443334449738, 111.932724563851, 88.0353929988547,
      78.0941759541126, 110.556760626146, 86.9859186312889,
      111.683664814931, 104.508903382769, 98.2850776999855,
      76.2783750257764, 126.957683883533, 123.141669446248,
      108.346657967421, 105.767089522976, 114.27236205513,
      161.436490159227),
    nrow = 100, ncol = 1, dimnames = list(NULL, "beta"))

attr(proposals_exp, "ancestor") <- rep(NA_integer_, 100)

proposals_obs <- .Call(SimInf:::SimInf_abc_proposals, ## function
                       "beta",                        ## parameter
                       "gamma",                       ## distribution
                       10,                            ## p1
                       0.1,                           ## p2
                       100L,                          ## n
                       NULL,                          ## x
                       NULL,                          ## w
                       NULL)                          ## sigma

stopifnot(all(abs(proposals_obs - proposals_exp) < tol))

proposals_exp <- matrix(
    c(195.005336846255, 140.901856658208, 154.907184516442,
      141.727518910812, 38.8627963471009, 97.4310586052267,
      77.9025517251395, 161.553738349774, 14.6661883074001,
      138.538188048932, 135.791097877793, 26.8124231064739,
      58.2166821923673, 88.9042366085132, 125.392952230297,
      36.4609319468952, 87.3801617715897, 44.450678744187,
      121.338336949295, 81.7396336435478, 85.8797448305749,
      65.9594817591945, 119.222178557239, 88.7666717784036,
      161.215764012589, 144.907147873647, 93.1307406434064,
      128.658939113042, 106.441515784153, 97.5774119114729,
      47.0868669053793, 38.6724597152766, 60.0168467240897,
      127.171904393379, 25.6809402954593, 107.766033910909,
      47.0098212565299, 66.8680237581495, 73.3515474282523,
      106.457640635664, 110.921501906587, 205.691203703484,
      156.995403049369, 156.271004316669, 62.4625040304671,
      82.9020590751595, 121.916803758159, 96.378331880152,
      185.349514517733, 80.596742844214, 92.3981388259838,
      94.0004606629752, 66.0515508010587, 37.3968637643038,
      38.1447180927324, 118.130218920069, 145.276845109855,
      67.7079233919813, 91.3783382695343, 22.210128875797,
      44.2076870110467, 195.260355686186, 74.6184146825141,
      177.573584431022, 140.835447969089, 126.947780191398,
      61.5761331369285, 144.305281196236, 39.3401179774552,
      112.535377741002, 146.891864943471, 3.69622962315792,
      39.4318184166903, 33.317670867482, 29.30667571148,
      144.360016400905, 140.576403248937, 108.256040447012,
      4.96692849268332, 75.4542720813571, 134.799250355623,
      141.590525114084, 31.6742831157865, 122.63568014853,
      143.086340668081, 138.252218077257, 48.2821523236046,
      57.4456206901338, 36.0037611028634, 121.776246615544,
      112.778325036825, 61.3014719001834, 60.6066480472759,
      121.09400043126, 103.265776758772, 60.5899765586154,
      90.830668022743, 105.691068362283, 88.0208565279985,
      69.5282553651992),
    nrow = 100, ncol = 1, dimnames = list(NULL, "beta"))

attr(proposals_exp, "ancestor") <-
    c(79L, 41L, 89L, 95L, 4L, 51L, 89L, 54L, 44L, 95L, 44L, 66L, 57L,
      10L, 89L, 23L, 4L, 32L, 95L, 89L, 69L, 63L, 99L, 66L, 69L, 54L,
      60L, 29L, 13L, 95L, 89L, 69L, 1L, 48L, 76L, 23L, 32L, 23L, 13L,
      41L, 41L, 35L, 16L, 13L, 23L, 48L, 26L, 86L, 4L, 44L, 79L, 13L,
      57L, 19L, 13L, 76L, 89L, 66L, 10L, 38L, 26L, 82L, 44L, 82L, 82L,
      79L, 44L, 76L, 63L, 73L, 1L, 48L, 23L, 38L, 60L, 35L, 10L, 23L,
      66L, 41L, 79L, 10L, 44L, 99L, 89L, 89L, 16L, 13L, 66L, 35L, 66L,
      32L, 19L, 79L, 10L, 48L, 51L, 60L, 32L, 48L)

x <- proposals_obs
sigma <- 2 * cov(x)
set.seed(123)
proposals_obs <- .Call(SimInf:::SimInf_abc_proposals, ## function
                       "beta",                        ## parameter
                       "gamma",                       ## distribution
                       10,                            ## p1
                       0.1,                           ## p2
                       100L,                          ## n
                       x,                             ## x
                       rep(0.01, 100),                ## w
                       sigma)                         ## sigma

stopifnot(all(abs(proposals_obs - proposals_exp) < tol))

w_exp <-
    c(0.00651534674717087, 0.00984796740883172, 0.00837633901987569,
      0.00974951167168201, 0.00278220798043043, 0.0155575235280236,
      0.0152955767286177, 0.00782971162134828, 1.08450996490866e-05,
      0.0101374234187459, 0.0104874043706198, 0.000471362147479918,
      0.010039155098466, 0.015946517384311, 0.0119218523425067,
      0.00212726267952287, 0.0159435097489079, 0.00463101738461317,
      0.0125123478192166, 0.0156975911661324, 0.0159153946288713,
      0.0127043058985561, 0.0128227314478057, 0.0159472657242274,
      0.00785514716156604, 0.00938364135901829, 0.0158323412720233,
      0.0114557332641912, 0.0146082203934316, 0.0155457556084599,
      0.00562342343268247, 0.00272692808985884, 0.0107118965247735,
      0.0116665764386191, 0.000371884012831849, 0.0144388327309886,
      0.00559364440093777, 0.012971812046599, 0.0145495783097439,
      0.0146061930791147, 0.0140152650272004, 0.00660209791824104,
      0.00819402858322751, 0.00825617118399543, 0.0115789747602285,
      0.0157804696048089, 0.0124276183858865, 0.0156377738844014,
      0.00665373930369284, 0.0155987993282588, 0.0158641257942813,
      0.015788585386551, 0.0127319040306528, 0.00237129878282316,
      0.00257663423008185, 0.0129827571322452, 0.00934249530924648,
      0.0132093694191666, 0.0159003122141, 0.000160671246114197,
      0.00454270566789663, 0.00651439471946911, 0.0147873168565003,
      0.00690693996069261, 0.00985594685087353, 0.0116985669592717,
      0.011271619423175, 0.00945125415014257, 0.00292329993088291,
      0.0137900519504907, 0.00916621317183316, 2.13578824730994e-10,
      0.0029508041752752, 0.00141656611045132, 0.000754289856796076,
      0.00944507326311541, 0.00988715841702764, 0.0143747856285377,
      2.5354782960133e-09, 0.0149314073605162, 0.0106171613156254,
      0.00976575058511817, 0.00111263535838425, 0.0123224815498804,
      0.00959055004797866, 0.0101731901822058, 0.00609012404850771,
      0.00974356928499434, 0.00201338985696298, 0.0124481970836683,
      0.0137557517117914, 0.0111747332836431, 0.0109263230245733,
      0.0125481621414489, 0.0149882002226354, 0.0109203059677858,
      0.0159157334572017, 0.0147015558158957, 0.0159478292853504,
      0.0136910878586874)

x <- x[attr(proposals_obs, "ancestor"), , drop = FALSE]
xx <- proposals_obs
w_obs <- .Call(SimInf:::SimInf_abc_weights, ## function
               "gamma",                    ## distribution
               10,                          ## p1
               0.1,                         ## p2
               x,                           ## x
               xx,                          ## xx
               rep(0.01, 100),              ## w
               sigma)                       ## sigma

stopifnot(all(abs(w_obs - w_exp) < tol))

##
## Check that it works to use a normal distribution for the proposals.
##
set.seed(123)
proposals_exp <- matrix(
    c(10.0800554282007, 10.1190206626991, 9.83104443359712,
      10.1239495885998, 9.98910340276845, 9.9882758038212,
      10.0183082613838, 10.1280554878202, 9.82727293711329,
      10.1690184353646, 10.0503812447155, 10.2528336550704,
      10.0549096735636, 10.0238212920794, 9.89511068564135,
      10.1294763254584, 10.0825539843874, 9.99443139882799,
      9.92156177824118, 9.92664967774085, 9.97841346059533,
      9.96650872422866, 9.89143008646971, 9.99145767355659,
      10.1070610538271, 9.98546064490073, 9.883445515211,
      9.91814842774869, 10.0684936077925, 9.96799435807232,
      9.86884775886032, 9.9400391672056, 9.98705893109928,
      10.0886736147004, 9.98486040376642, 10.0329791200758,
      9.67726771702325, 9.92282082287538, 10.0286548567354,
      9.87794880176339, 10.0434550377113, 10.0800176865835,
      9.9836069031357, 10.1242918774937, 9.90656149419448,
      10.0393708652216, 10.0403631455697, 9.91135632835832,
      9.86810623964634, 10.0028843908516, 9.95678702065287,
      10.1689872518646, 10.1228392782136, 10.0276023478266,
      9.89510244955269, 9.94791306561671, 10.1623202521731,
      9.89299317706557, 10.1685887244205, 9.9758310232261,
      9.95317995213595, 9.92270217718191, 10.2149919335756,
      9.86656463717927, 10.0495870479997, 10.1233976240157,
      10.0634362124806, 10.0412022274758, 10.0793585307679,
      9.98475893669352, 9.97711041848692, 9.90992082494195,
      9.92649738442157, 9.85723142164066, 10.0619283534849,
      9.99938017377304, 9.9314293154307, 9.97206664723628,
      9.92172697247167, 9.92210027602678, 9.96251999068716,
      9.96806061911706, 10.0084543768051, 9.92315263970603,
      9.93740890869283, 9.90991291447621, 10.0663728669675,
      10.030027911812, 10.0074856823862, 10.0206372695351,
      9.95110771647057, 9.93720483421767, 9.99530832735717,
      10.0162618115322, 10.1292305914973, 9.953644349848,
      10.0305463227482, 9.99160112867597, 10.0410363449162,
      10.0183678240529),
    nrow = 100, ncol = 1, dimnames = list(NULL, "beta"))

attr(proposals_exp, "ancestor") <- rep(NA_integer_, 100)

proposals_obs <- .Call(SimInf:::SimInf_abc_proposals, ## function
                       "beta",                        ## parameter
                       "normal",                      ## distribution
                       10,                            ## p1
                       0.1,                           ## p2
                       100L,                          ## n
                       NULL,                          ## x
                       NULL,                          ## w
                       NULL)                          ## sigma

stopifnot(all(abs(proposals_obs - proposals_exp) < tol))

proposals_exp <- matrix(
    c(10.249416561216, 10.207608519523, 10.1616680108471,
      10.1805772834956, 9.91479236877152, 9.93191420178626,
      9.89396473893688, 10.1527586225798, 9.99011930878491,
      10.1694897128394, 10.4112048126883, 9.96355706675316,
      10.1239136939972, 10.1081469462128, 10.059063060151,
      9.79861661562053, 10.0834608890643, 9.88213968735834,
      10.1096951751754, 9.90730418809612, 9.99937601512637,
      10.1134158987824, 10.0582441676788, 10.1789383298931,
      10.2612784240681, 10.0948874715195, 9.96999644377345,
      10.1234876364007, 10.1515596606107, 10.0270913480332,
      9.78683534485534, 9.83526193230731, 9.4443095703895,
      9.86708730357155, 10.0748322929171, 9.80619868825253,
      10.0465057153763, 9.891036435502, 9.90432555379842,
      10.036523818283, 10.0878646938157, 10.103383114257,
      10.2275566574943, 10.3509405541753, 10.3247897200789,
      9.8890099552063, 9.92093009237337, 10.0927567574681,
      9.85583874175354, 10.4240471022101, 10.2193240271061,
      9.89270707977462, 10.1083088711292, 10.1511512754972,
      9.82209594339388, 9.91412880041124, 10.1275946027116,
      10.1281885528548, 9.65709650665204, 10.1057285028862,
      10.1167480535188, 9.7579238538185, 9.82260437758685,
      10.178553405421, 10.1985406280894, 10.1170661109179,
      9.98934757012056, 10.0128174121219, 10.1531997009741,
      10.2185910734405, 10.0208748331058, 9.96138254760106,
      10.1691045430748, 9.64557445557514, 9.80894477380542,
      9.79653874264101, 9.74811484536115, 10.014341440229,
      10.2877829963427, 10.0482092018837, 9.88761215054977,
      9.98008282535699, 10.040112708211, 10.2913085470995,
      10.0492472140046, 10.0701110581701, 10.1205733590029,
      10.1037677411139, 9.97300364487065, 9.98122755264019,
      9.9955103553499, 9.93582994235593, 10.2624138028386,
      9.94072074019406, 9.90278375356397, 9.99246700023077,
      10.1580742187487, 9.84336310352742, 9.90896822969643,
      10.0136618792768),
    nrow = 100, ncol = 1, dimnames = list(NULL, "beta"))

attr(proposals_exp, "ancestor") <-
    c(79L, 41L, 89L, 95L, 4L, 51L, 89L, 54L, 44L, 95L, 44L, 66L, 57L,
      10L, 89L, 23L, 4L, 32L, 95L, 89L, 69L, 63L, 99L, 66L, 69L, 54L,
      60L, 29L, 13L, 95L, 89L, 69L, 79L, 1L, 48L, 76L, 23L, 32L, 23L,
      13L, 41L, 41L, 35L, 16L, 13L, 23L, 48L, 26L, 86L, 4L, 44L, 79L,
      13L, 57L, 19L, 13L, 76L, 89L, 38L, 66L, 10L, 38L, 26L, 82L, 44L,
      82L, 82L, 79L, 44L, 76L, 63L, 73L, 1L, 48L, 23L, 38L, 60L, 35L,
      10L, 23L, 66L, 41L, 79L, 10L, 44L, 99L, 89L, 89L, 16L, 13L, 66L,
      35L, 66L, 32L, 19L, 79L, 10L, 48L, 51L, 60L)

x <- proposals_obs
sigma <- 2 * cov(x)
set.seed(123)
proposals_obs <- .Call(SimInf:::SimInf_abc_proposals, ## function
                       "beta",                        ## parameter
                       "normal",                      ## distribution
                       10,                            ## p1
                       0.1,                           ## p2
                       100L,                          ## n
                       x,                             ## x
                       rep(0.01, 100),                ## w
                       sigma)                         ## sigma

stopifnot(all(abs(proposals_obs - proposals_exp) < tol))

w_exp <-
    c(0.000142318649743074, 0.00579600739030807, 0.00521010488029535,
      0.0165539935379971, 0.00172895512773806, 0.0140048724453806,
      0.00888657462900812, 0.00749296302718984, 0.00512353136012425,
      0.0164045467559143, 0.0034692201875978, 0.00368085124683484,
      0.0135665386787454, 0.0115382622999379, 0.0109993286318955,
      0.0213240392591689, 0.0120365872070376, 0.0152173041076391,
      0.013686467402762, 0.0096297749493388, 0.00905993715711376,
      0.00842245374843249, 0.0123677530593766, 0.0160396979463799,
      0.00647004757375468, 0.0106877125112904, 0.0130034348006057,
      0.0125367112467912, 0.010223394143662, 0.00731372071440856,
      0.00320829921000589, 0.00123451323762445, 0.000203430658368934,
      0.00203082608980474, 0.00337865836650964, 0.00474672289135411,
      0.00373028711065374, 0.0152326925857242, 0.0159847854575654,
      0.0121396094574977, 0.0119484439807612, 0.0114489276857043,
      0.00135077302681736, 0.00772632809892925, 0.00170492307554806,
      0.0173489562342237, 0.0149667748648686, 0.0074986382601814,
      0.0181847424556568, 0.00289572448163313, 0.0153933717284334,
      0.0163302426547517, 0.0120472217505778, 0.0161782108795572,
      0.0163845462220557, 0.0057351981162584, 0.006504438466244,
      0.00715726804137168, 0.00554456349043892, 0.0135771081805756,
      0.012437515912019, 0.0126350733628201, 0.00710539406351158,
      0.0020376794150082, 0.0159574787319664, 0.00470310241470761,
      0.0123331881474831, 0.00816167932410601, 0.0157587662712092,
      0.00217622297865884, 0.00187329469995823, 0.0124872852168028,
      0.0119426597105867, 0.00666924543834568, 0.0213177343576396,
      0.0150601089567754, 0.00410673993607958, 0.0118248989047837,
      0.0214930835596806, 0.00363721209302184, 0.00108469621550827,
      0.0104620062458924, 0.00613617592913788, 0.0213476067899275,
      0.00938821961969724, 0.0122272869750957, 0.00761776670233406,
      0.00862798022767443, 0.00380709274766305, 0.00973035169542611,
      0.00552821133025288, 0.0126413903108567, 0.0130108587964331,
      0.0140347426335291, 0.0159479734612046, 0.0097851155898044,
      0.0167082434400863, 0.0181775098321298, 0.0140845927002193,
      0.0115004838257942)

x <- x[attr(proposals_obs, "ancestor"), , drop = FALSE]
xx <- proposals_obs
w_obs <- .Call(SimInf:::SimInf_abc_weights, ## function
               "normal",                    ## distribution
               10,                          ## p1
               0.1,                         ## p2
               x,                           ## x
               xx,                          ## xx
               rep(0.01, 100),              ## w
               sigma)                       ## sigma

stopifnot(all(abs(w_obs - w_exp) < tol))

res <- assertError(
    .Call(SimInf:::SimInf_abc_weights, ## function
          "normal",                    ## distribution
          10,                          ## p1
          -0.1,                        ## p2
          x,                           ## x
          xx,                          ## xx
          rep(0.01, 100),              ## w
          sigma))                      ## sigma
check_error(res, "Invalid weight detected (non-finite or < 0.0).")