diff --git a/DESCRIPTION b/DESCRIPTION index 48870c7..97c04d9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,18 +1,16 @@ Package: Bclim Type: Package Title: Bayesian Palaeoclimate Reconstruction from Pollen -Version: 2.2 -Date: 2013-01-23 +Version: 2.3.1 +Date: 2014-04-15 Author: Andrew Parnell, Thinh Doan and James Sweeney Maintainer: Andrew Parnell -Depends: mclust, hdrcde, statmod, MASS +Depends: MASS, mclust +Imports: hdrcde, statmod Suggests: Bchron, doMC, foreach -Description: This package takes pollen/time data from lake cores and - produces a Bayesian posterior distribution of palaeoclimate - from that location after fitting a non-linear non-Gaussian - state-space model. This latest version fixes a bug where - volatilities were estimated incorrectly +Description: This package takes pollen/time data from lake cores and produces a Bayesian posterior distribution of palaeoclimate from that location after fitting a non-linear non-Gaussian state-space model. License: GPL (>= 2) -Packaged: 2013-01-23 13:15:04 UTC; aparnell +Packaged: 2014-04-17 14:26:45 UTC; aparnell +NeedsCompilation: yes Repository: CRAN -Date/Publication: 2013-01-23 16:20:13 +Date/Publication: 2014-04-18 08:39:17 diff --git a/MD5 b/MD5 index 45cfbfe..c0e6b82 100644 --- a/MD5 +++ b/MD5 @@ -1,31 +1,27 @@ -3055d527f9894b4cb6ee8d6b238ac748 *DESCRIPTION -8e419122bd6a5b07a3a495d573138afc *NAMESPACE -5f2912b0a548a38cd6d4eb0cf7b7e477 *R/Bclim-internal.R +a8f48f85af68c7fb7e284f820a54aa36 *DESCRIPTION +2f113365dfbb6c37cf8748a3a09bac07 *NAMESPACE c6f874826f8e6697c8eec0b90413f529 *R/Bclim.R -f776eb835eb11d4c792b0cf2cca720e5 *R/BclimCompile.R -4cff32608adc1207131105a1982713e9 *R/BclimInterp.R -af57a23af5cd03c0f897ade4349cb3d2 *R/BclimLayer.R -041a88b7f1a4e3db237178fe3c5422f8 *R/BclimMCMC.R -c382207be84fe453dc5d4b7c35c4da95 *R/BclimMixPar.R -89588b14f7687294d938dd3f7ca1b1cc *R/BclimMixSer.R +276e9d9aea394b90e58fbfcf9f45253c *R/BclimCompile.R +fa6543be87cf722c53fe6a046e918762 *R/BclimInterp.R +0595934946efef63e52463392895e476 *R/BclimLayer.R +5d5c2210cce82ea17fcc11839bbc97fa *R/BclimMCMC.R +b5af690c5e7e0c3d19da9a480d23a767 *R/BclimMixPar.R +58f1410cb5c2c35da61f24e2b5266f83 *R/BclimMixSer.R 856245e1980b34b2eea6bb3f130c9c08 *R/BclimRun.R -f48b41d0a322f1bfa36fac12444f3bd4 *R/NIGB.R -3c8d4d5c7857fb8cf3a1eb8a841642f2 *R/NIGp3.R -9b83e343f86b627f055a360f78511deb *R/plotBclim.R -d385f12dbe1227314f75a2b58b53e2ea *R/plotBclimVol.R -eba7fbd97a896a794b26ebcee5e784fc *man/Bclim-package.Rd -3b3ae1db0dd7153c87186a432aeff7c6 *man/BclimCompile.Rd +254beb9b64f7819c7ff3846767f8a44e *R/plotBclim.R +9a1092b29526425879e14855de9144e6 *R/plotBclimVol.R +90f8e07927a8651a6deeb9eb5a87fdea *R/zzz.R +a68bbc80ffe7b1c1e5d50110ac6d6c06 *man/Bclim-package.Rd +63baf8f9f1b0106d60434f345faaf428 *man/BclimCompile.Rd 844027ba79272149b55c70e30dde004a *man/BclimInterp.Rd -cb66a992d9bdf6fdd8f0c9449f573f4f *man/BclimLayer.Rd -fd2b950219625b2cc6fed00790569963 *man/BclimMCMC.Rd -bc6b60f2b815ecf01afe6d3379879482 *man/BclimMixPar.Rd -d3e0c28f0ba477c4b69f042d18564b1f *man/BclimMixSer.Rd -d098a650d48d745dfcd4ef089baa1cb4 *man/BclimRun.Rd -2a5fdad1ccd8655ef2a5a9457b876578 *man/NIGB.Rd -78acf6ecd9f34501621bb24312632701 *man/NIGp3.Rd -f46732a1c4585ce1457b92572090dde7 *man/plotBclim.Rd -d869628c17b2764282187b609858549e *man/plotBclimVol.Rd -5b72d41aaa4cd831235669faa9f3e617 *src/BclimMCMC3D.c +4132c98098e3d0e0749d6f09b99ba77e *man/BclimLayer.Rd +80aea2ad9cdc6706ea63c9b56fae4463 *man/BclimMCMC.Rd +bfdbdd437f2f9bf65b3021933d9bb81d *man/BclimMixPar.Rd +afb588a381b69654837f58ecb4128afd *man/BclimMixSer.Rd +2ab97fc9ac914c86d95970d3fd267141 *man/BclimRun.Rd +9efb7750817a00efb8d54f706abb1ae7 *man/plotBclim.Rd +466695715b643ea20742560c271a83d1 *man/plotBclimVol.Rd +46000bef71469cd91818dbbd68818493 *src/BclimMCMC3D.c efa94aa42c370a9851dca166219f9c0f *src/PalaeoRecon3D.c -effde32cdc7ea252698c65b8e9b48b58 *src/use.c -1c8929a0e47b44cb06b0586a95adbbac *src/use.h +384ae2930ede0489a28287fb2fa5eb59 *src/use.c +6dab0cf9547153c2ecf4bc9825dd9fca *src/use.h diff --git a/NAMESPACE b/NAMESPACE index 991e55f..20367e2 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,4 +2,7 @@ exportPattern("^[[:alpha:]]+") useDynLib(Bclim) -export(Bclim, BclimLayer, BclimMixPar, BclimMixSer, BclimMCMC, BclimCompile, BclimInterp, plotBclim, plotBclimVol) \ No newline at end of file +import(MASS) +import(mclust) +importFrom(statmod,rinvgauss) +importFrom(hdrcde,hdr) \ No newline at end of file diff --git a/R/Bclim-internal.R b/R/Bclim-internal.R deleted file mode 100644 index ef3140e..0000000 --- a/R/Bclim-internal.R +++ /dev/null @@ -1,134 +0,0 @@ -.Random.seed <- -c(403L, 10L, 1579219043L, 450026753L, 1134948464L, -205808626L, -1257716857L, 1589496843L, -132782790L, -661769156L, 1864729471L, --119582507L, 1831445740L, -144064958L, -1530687075L, 1706542103L, -2074745214L, -286441672L, -588351461L, -1880491655L, 1016176584L, -732890102L, 1769100849L, -747900829L, -1826926542L, 185871748L, -1997706183L, 290674333L, 1278049812L, 947717306L, -709988091L, --313199905L, 96247654L, 263124720L, -475787757L, 1973864369L, --640675232L, -1279479810L, -276127863L, 2040510651L, -2133740406L, -1497450348L, -2147281681L, -1056538843L, -2136675876L, -2077403214L, --1953800531L, -631941529L, 477898126L, 387076712L, -1266512277L, --281942007L, 126260888L, 814895558L, 1409646881L, 1131680179L, --227971806L, -1006287084L, 458047383L, -1051828339L, -1130765660L, -1752427658L, -1701281035L, 1064510799L, -88378794L, 1308774656L, -1622762563L, 1762469857L, -401569968L, -1589610066L, -244861095L, -1258732907L, 22963418L, 1168633244L, 913412511L, -927190027L, -1253941900L, 763488098L, -2022533315L, 1173827511L, 1218083934L, -477983128L, -1224223493L, -626745191L, -245657304L, 2044372950L, --360102831L, -5574909L, 196585106L, 1714873380L, -617157913L, --1779368131L, -2009727052L, 1155097626L, -1966702747L, 30827007L, -1664363654L, 1105702928L, -1169689933L, -862391983L, 1607553280L, --1706896418L, 615907369L, -2064853029L, -1111030550L, 594469516L, --1864172081L, -1938771579L, 919177788L, 1014558482L, 1538151437L, --1251745977L, -459449042L, -1066705784L, 408291723L, -118025879L, --379462856L, -88341018L, 990208385L, 1332573459L, 1440608962L, --1436131980L, 686611831L, 1546839661L, -1114347260L, 2086829098L, --449064555L, 710051375L, 1172480822L, -359867552L, -804331101L, --1962024255L, -1588987216L, -1675020850L, 268626233L, -1719892661L, --2056073478L, -2046836740L, -1710571457L, 429881621L, -1605664212L, -369069826L, 1872518109L, -1543776553L, 1561829822L, -1186561032L, --332859429L, -184890823L, 1659521928L, -1230800074L, -2076448399L, -512574371L, 1840993010L, 754692676L, -154526457L, -1329077027L, -567759316L, 2005169402L, 1773497413L, 473136287L, 383589926L, -930257584L, 1847265107L, 1539785329L, -1326641248L, 218179390L, --238522807L, -129135877L, -1539687734L, -1769706452L, 1248228399L, --132326171L, -1698769508L, 587339762L, -1800748435L, 237684647L, --1258735154L, 2025590568L, -1205418709L, -1867536567L, 738159704L, --504342266L, 1369531617L, -628276109L, 1212158178L, -1034078124L, --305504937L, -1421841075L, -1677905692L, 951413962L, -1834172107L, -2013611919L, -1903094634L, 1124654144L, 1805913347L, -939097823L, -514116112L, 290197998L, -957921895L, 110469675L, -2067967206L, -244689372L, 1606983647L, 719022005L, -800194228L, 1894746018L, -1838722813L, 1864204279L, -1017214690L, 1484403672L, 245159227L, -290697945L, -347689368L, -1337380970L, -1761846511L, 642086083L, -605617362L, -1986647708L, 71257255L, -1894075907L, 734551028L, --38054694L, 1558033957L, -879088833L, -149370682L, 156861776L, --88789645L, -1516669039L, 185240000L, -1919589424L, 2060326242L, --281985304L, -1397105460L, -517588484L, 1700503026L, 445586560L, --204842508L, -1156963704L, -1383130310L, -2059934384L, -596810644L, -1911434644L, 228550738L, -2026578464L, -548481844L, 2060820128L, --1039519358L, -603029272L, 701000428L, -431645380L, -694572942L, -596214384L, 1096047140L, -1656174536L, 723794202L, -990845744L, -2089589052L, 1420178324L, -1710590078L, -194532416L, 2084189756L, -1712849296L, 1551345698L, 237726984L, -650278772L, 157110044L, -663870610L, 689720704L, 1171271700L, 1645111080L, -1868855750L, -65685968L, -1377509012L, 959273908L, -1782812142L, -59394080L, -1078925772L, 331012320L, -1071721438L, -1249739992L, 1200007692L, --445787332L, 1622140914L, 1938268912L, -1609455036L, -530934184L, --1049413830L, 497317872L, 1936804540L, 1784014100L, 703147202L, -382798944L, 1569283452L, -70869040L, 984320290L, 1392832168L, --1073209460L, 304383164L, 207237938L, -18393536L, 128955572L, --920171448L, -91025990L, -2013726192L, 854592748L, 1019736404L, -156605202L, -946944608L, 913765068L, 708315744L, -2048530494L, -924050216L, 712655276L, 965162300L, -547667086L, -2130262736L, -692051748L, 1063085368L, -463730982L, 455514000L, 1603560188L, -1595997076L, 641909826L, -569586432L, 30113852L, 827297104L, --923765598L, -1700103352L, -1740558580L, -261837092L, -1917843630L, -1656571520L, 1350299732L, 88556520L, -1289704710L, -808759984L, -340891436L, -30897292L, -1677556846L, 1276562016L, -1225666676L, --1491314848L, -820839454L, 720997288L, -83153716L, -1210004996L, -1053953970L, -2031804944L, -672186940L, -1870712488L, 1981859194L, -469343664L, 1352582076L, 1736796244L, 302571458L, -1942874080L, --1410797380L, 374020688L, -2081482782L, 1628790376L, -42541492L, --850564868L, -488132622L, 1533011200L, -1193534604L, 711972872L, -2127817146L, 353668560L, 1675875308L, -1798804588L, -724907054L, --1983025312L, -2035222196L, 1114579232L, -871108734L, -189116696L, --1214203668L, 435994172L, 1314661234L, -824850832L, 1962601508L, -109174584L, -629210982L, -2060286256L, 985244860L, -1092318572L, -1086610050L, -1423945536L, 1651270972L, -1958016880L, 1210498594L, -201325064L, 1728648460L, 1838047900L, -2081548782L, 282587520L, --1255906796L, -1519190488L, -2039725766L, 363322320L, 1920241260L, -382106420L, 53949074L, 374834528L, 129447244L, 1769382880L, -1732681822L, --1358717016L, -960303220L, -440206660L, -465691790L, -1589673744L, --1098904508L, 1477799256L, -1314041798L, 1533395824L, 931746876L, -816868116L, 393366978L, 264211680L, 1816164220L, -1951164720L, --1237902686L, -1450290136L, 1306837260L, 1229147068L, 1805092146L, -1350268864L, 291759668L, -319493176L, 1439900858L, 559073168L, -680241644L, -1200548780L, 1731107474L, 63957792L, 985588172L, --3187744L, 1617483458L, 653366824L, -1346818388L, 1267599036L, -666108402L, -76605136L, -1994471260L, -1612416328L, -924737702L, --420921968L, -510087812L, -738766572L, -2072236734L, 1578431232L, -1888791484L, 2167632L, -1559976670L, -1329206968L, -2106022846L, --1076126633L, -1770784863L, 1917683814L, -2135666588L, -2023724139L, -1406180543L, 1912147032L, -1004849034L, 1317546115L, -1907816251L, -288097762L, 568421456L, -1329248903L, -518604117L, 1985012252L, -1446133066L, 1203348799L, -1155980407L, -732292194L, -1959396788L, -1881316381L, 1185363367L, -2050307888L, -1368794514L, -1007909605L, --1600061827L, -1515647894L, 1077852424L, 80620977L, -1516959197L, --1235898588L, 1453004754L, -1843186649L, 1104866929L, 196695222L, -1153242196L, -1576358683L, 1256069359L, -1967354840L, 1839675526L, -230179507L, 65094677L, 314463858L, 419156672L, 1267611433L, -2098350693L, --2095934868L, 914310778L, 2146547631L, -181169479L, -563381042L, --426633700L, -1940640371L, 2018802679L, 715669248L, -922413730L, -642956843L, 90517805L, -1069412966L, 1069767512L, -1622302015L, -1009583475L, -975317356L, -1160633566L, -880366025L, 899595009L, -211386886L, 1382588484L, -1434406539L, -996436769L, -254810056L, --478171114L, 166420515L, -1826832411L, -1279866366L, 364992624L, --900895207L, -1642485L, -995497220L, 1376647978L, 1834069343L, -1348445417L, -1408650242L, 52809708L, 1885117245L, -1171237113L, --228487440L, 357969550L, 1959618747L, 73604253L, -1689137078L, -395522664L, -1669530031L, -1668627133L, -994036668L, 721499890L, -1573764807L, -496810863L, 813087254L, -1171953548L, 13773573L, -1030193487L, -279592120L, 1074574054L, 1613997075L, -1327637643L, --847103150L, -761058400L, 1724169865L, 2068433595L, 83931404L, --2131178342L, 1072484879L, 1401632345L, 813419246L, -903429764L, --379981331L, -565707113L, 763681120L, -1506845378L, -278557173L, -873560717L, -773179462L, 1214170616L, -1553489503L, -1553037485L, -1918355060L, 1016465410L, 1168142999L, 689904737L, 333723814L, --1100994396L, 1010295125L, 393708031L, 262961304L, 113370550L, --105689533L, 1824093061L, -1099392990L, -40297712L, 615757497L, --863696021L, 330293852L, 202555274L, 256634751L, -1358232887L, -391162718L, 1982353420L, -1667252259L, -1582232857L, -213900528L, -1675311662L, 804659419L, -1920606787L, -967529686L, -65084472L, --621747983L, -2056544669L, -1045337500L, -1138701934L, -902388121L, -269919793L, -238609546L, 1442115732L, 1602223781L, 1298557871L, -1243023464L, -1478943674L, 812944852L) -.onAttach <- -function(libname, pkgname) { - Bclimver <- read.dcf(file=system.file("DESCRIPTION", package=pkgname), - fields="Version") - packageStartupMessage(paste(pkgname, Bclimver)) - packageStartupMessage("Welcome to Bclim. Type help(Bclim) to get started.") - packageStartupMessage("See http://mathsci.ucd.ie/~parnell_a/Bclim.html for updates, bugs and a tutorial.") -} diff --git a/R/BclimCompile.R b/R/BclimCompile.R index db5a9a8..8b053a3 100644 --- a/R/BclimCompile.R +++ b/R/BclimCompile.R @@ -2,10 +2,7 @@ BclimCompile <- function(Layers,Mixtures,MCMC,Interpolations,core.name="Core") { clim.dims <- c("GDD5","MTCO","AET/PET x 1000") - results <- list(time.grid=Interpolations$time.grid,core.name=core.name,clim.interp=Interpolations$clim.interp, - vol.interp=Interpolations$vol.interp,MDP=Layers,ScMean=Mixtures$ScMean,ScVar=Mixtures$ScVar, - clim.dims=clim.dims,n=Mixtures$n,m=Mixtures$m,n.samp=Mixtures$n.samp,Chronsfile=MCMC$chron.loc, - nchron=MCMC$nchron,chron.store=MCMC$chron.store) + results <- list(time.grid=Interpolations$time.grid,core.name=core.name,clim.interp=Interpolations$clim.interp,v.interp=Interpolations$v.interp,MDP=Layers,ScMean=Mixtures$ScMean,ScVar=Mixtures$ScVar,clim.dims=clim.dims,n=Mixtures$n,m=Mixtures$m,n.samp=Mixtures$n.samp,Chronsfile=MCMC$chron.loc,nchron=MCMC$nchron,chron.store=MCMC$chron.store) class(results) <- "Bclim" return(results) } diff --git a/R/BclimInterp.R b/R/BclimInterp.R index 4f30959..615d8a3 100644 --- a/R/BclimInterp.R +++ b/R/BclimInterp.R @@ -1,161 +1,205 @@ -BclimInterp <- -function(Bclim.data, Bclim.res, tgrid=seq(0,14,by=0.1)) { - chron <- Bclim.res$chron.store # 1000-by-150 - alpha <- sqrt(Bclim.res$phi/Bclim.res$eta) - delta <- sqrt(Bclim.res$eta*Bclim.res$phi) - clim <- Bclim.res$c.store # 1000-by-150-by-3 - vol <- Bclim.res$v.store # 1000-by-149-by-3 - - n.interp <- dim(clim)[1] #1000 - n.sections <- dim(clim)[2] #150 - - L <- n.sections + length(tgrid) #150+141 - t.all.mat <- matrix(NA, nrow = n.interp, ncol = L) #1000-by-(150+141) - c.all.mat <- matrix(NA, n.interp, L) #1000-by-(150+141) - v.all.mat <- matrix(NA, n.interp, L-1) #1000-by-(150+140) - clim.interp <- array(NA, dim = c(n.interp, length(tgrid), Bclim.data$m)) # r=3 clim dim - vol.interp <- array(NA, dim = c(n.interp, length(tgrid)-1, Bclim.data$m)) +BclimInterp <- function(Bclim.data, Bclim.res, tgrid=seq(0,14,by=0.1)) { + # Get relevant parts from data given in arguments + chron <- Bclim.res$chron.store # n.s by n where n.s is the numer of stored MCMC iterations + alpha <- sqrt(Bclim.res$phi2.store/Bclim.res$phi1.store) # Need this parameterisation for bridging both are n.s by m + delta <- sqrt(Bclim.res$phi1.store*Bclim.res$phi2.store) + clim <- Bclim.res$c.store # n.s by n by m + v <- Bclim.res$v.store # n.s by n-1 by m + eps = 1e-5 # A small number - for (k in 1:Bclim.data$m){ #loop through clim dim - cat("Perform interpolation for climate dimension", k, "\n") - for (i in 1:n.interp) { #loop through each sample - t.all <- c() - v.all <- c() #IGB - c.all <- c() #NIGB - - if(i%%10==0) { - cat("\r") - #cat("Completed:",round(100*iter/iterations,2),"%") - cat("Completed:",format(round(100*i/n.interp,2), nsmall = 2),"%") - if(i t.start & tgrid < t.end], t.end) - temp <- NIGB(delta[k]*chron.diff.i[1], vol[i, 1, k], t.b.points, c.start, c.end) - t.all <- c(t.all, t.b.points[-length(t.b.points)]) #the last t point is the starting point for the next bridge - v.all <- c(v.all, temp$IGB) - c.all <- cbind(c.all, temp$NIGB[-length(temp$NIGB)]) #the last clim point is the starting point for the next bridge - j.tmp <- 1 - } else if(i.tmp>0) { #MIRROR RESULT OF EXTRAPOLATION (NIG PROCESS), I.E. EXTRAPOLATE BACKWARD - t.b.points <- c(tgrid[1], tgrid[tgrid > tgrid[1] & tgrid < tgrid[i.tmp+1]], tgrid[i.tmp+1]) - temp <- NIGp3(c.start, alpha[k], delta[k]*chron.diff.i[i.tmp], length(t.b.points), (t.start-tgrid[1]) ) - temp$IGB <- temp$IGB - t.all <- c(t.all, t.b.points[-length(t.b.points)] ) - tmp.vstart <- v.all <- c(v.all, rev(temp$IGB) ) - c.all <- c(c.all, rev(temp$NIGB)[-length(temp$NIGB)] ) - j.tmp <- 0 - } else - { - j.tmp <- length( (chron[i,] < tgrid[1])[(chron[i,] < tgrid[1])==TRUE] ) - t.b.points <- c(chron[i, j.tmp], tgrid[tgrid < chron[i, j.tmp+1]], chron[i, j.tmp+1]) - temp <- NIGB(delta[k]*chron.diff.i[j.tmp], vol[i, j.tmp, k], t.b.points, c.start, c.end) - if(chron[i, j.tmp]==tgrid[1]) { - t.all <- c(t.all, t.b.points[-length(t.b.points)]) - } else - { - t.all <- c(t.all, t.b.points[-c(1, length(t.b.points))]) - } - v.all <- c(v.all, temp$IGB) - c.all <- c(c.all, temp$NIGB[-length(temp$NIGB)]) - } - - j.tmp2 <- length( (chron[i, ] <= max(tgrid))[(chron[i, ] <= max(tgrid))==T] ) + # Create some arrays to hold the interpolated climates and squared volatilities v + # Note that the storage here is longer as we have to interpolate onto a finer grid which includes the chronologies + clim.interp <- array(NA, dim = c(n.s, n.g, m)) + v.interp.all <- array(NA, dim = c(n.s , n.g-1+n, m)) + v.interp <- array(NA, dim = c(n.s, n.g-1, m)) + + # Extrapolation function + NIGextrap = function(c.s,t.s,t.out,phi1,phi2,future=FALSE) { + #c.s=curr.clim[1];t.s=curr.chron[1];t.out=t.select;phi1=Bclim.res$phi1[1];phi2=Bclim.res$phi2[1] + t.diff = abs(diff(sort(c(t.out,t.s)))) + mu = phi1*t.diff # Re-parameterise on to the R version of the IG distribution + lambda = phi2*phi1*t.diff^2 + v.out = statmod::rinvgauss(length(t.diff),mu,lambda) + if(future==TRUE) { + # extrapolating into the future + return(list(NIG=rev(cumsum(rev(rnorm(length(t.diff),mean=0,sd=sqrt(v.out)))))+c.s,IG=v.out)) + } else { + return(list(NIG=cumsum(rnorm(length(t.diff),mean=0,sd=sqrt(v.out)))+c.s,IG=v.out)) + } + } + + # Bridging function - taken from Ribeiro 2008 + NIGB <- function(delta, IG.sum, tb.points, c.start, c.end){ + total.points <- length(tb.points) + NIG.bridge <- rep(NA,total.points) + NIG.bridge[1] <- c.start + NIG.bridge[total.points] <- c.end + IG.increment <- rep(NA, total.points-2) + z <- IG.sum + eps = 1e-5 + + #curr.cstart = c.start + + l <- 2 + while(l < total.points){ - for (j in (j.tmp+1):(j.tmp2-1) ) { - t.start <- chron[i, j] - t.end <- chron[i, j+1] - c.start <- clim[i, j, k] - c.end <- clim[i, j+1, k] + if((tb.points[l]-tb.points[l-1]) max(tgrid)) {t.start <- max(tgrid)-1E-6} - if(t.end > max(tgrid)) {t.end <- max(tgrid)-1E-6} - #if use 14, problem t.all.mat[i,] %in% t.grids (too many t points at 14) + # Reparameterise + mu <- (tb.points[total.points]-tb.points[l]) / (tb.points[l]-tb.points[l-1]) + #lambda <- (1/delta^2 * (tb.points[total.points]-tb.points[l])^2) / z + lambda <- (delta^2 * (tb.points[total.points]-tb.points[l])^2) / z - t.b.points <- c(t.start, tgrid[tgrid > t.start & tgrid < t.end], t.end) - temp <- NIGB(delta[k]*chron.diff.i[j], vol[i, j, k], t.b.points, c.start, c.end) - t.all <- c(t.all, t.b.points[-length(t.b.points)]) - v.all <- c(v.all, temp$IGB) - if(!is.null(tmp.vstart)) { - v.all[1:length(tmp.vstart)] <- v.all[1:length(tmp.vstart)] + v.all[length(tmp.vstart)+1] - tmp.vstart <- NULL + if(z==0) { + print(c(delta,IG.sum,tb.points,c.start,c.end)) + stop("Problem in Inverse Gaussian Bridge - IG sum is zero") } - c.all <- c(c.all, temp$NIGB[-length(temp$NIGB)]) - if(any(is.na(c.all))) stop("NAs in interp!") - } - - t.start <- chron[i, j.tmp2] - c.start <- clim[i, j.tmp2, k] - - if(t.start < max(tgrid)){ - t.b.points <- c(t.start, tgrid[tgrid > t.start & tgrid < max(tgrid)], max(tgrid)) - tmp1 <- NIGp3(c.start, alpha[k], delta[k]*chron.diff.i[j.tmp2], length(t.b.points), (max(tgrid)-t.start)) - tmp1$IGB <- tmp1$IGB + temp$IGB[length(temp$IGB)] - t.all <- c(t.all, t.b.points) - v.all <- c(v.all, tmp1$IGB) - c.all <- c(c.all, tmp1$NIGB) - } else - { - t.all <- c(t.all, t.b.points[length(t.b.points)]) - c.all <- c(c.all, temp$NIGB[length(temp$NIGB)]) + # Compute the roots of the chi-square random variate + s1 <- mu + (mu^2*q)/(2*lambda) - (mu/(2*lambda)) * sqrt(4*mu*lambda*q + mu^2*q^2) + if(lambdamax(curr.chron)) { + t.select.index = which(tgrid>max(curr.chron)) + t.all.select.index = which(tgrid.all>max(curr.chron)) + t.select = tgrid[t.select.index] + + # Perform extrapolation + clim.temp = NIGextrap(curr.clim[length(curr.clim)],curr.chron[length(curr.clim)],t.select,Bclim.res$phi1.store[i,j],Bclim.res$phi2.store[i,j]) + clim.interp[i,t.select.index,j] = clim.temp$NIG + v.interp.all[i,t.all.select.index-1,j] = clim.temp$IG + } + + # Now look at how many gaps their are + for(k in 1:(n-1)) { + # Find out if in this section there are any tgrid points inside + t.select.index = which(tgrid>=curr.chron[k] & tgrid=curr.chron[k] & tgrid.all0) { + # Select which bits of the grid are inbetween the current section of the chronology + t.select = tgrid[t.select.index] + + # Now interpolate using the NIGB code + clim.temp <- NIGB(delta[i,j], curr.v[k], c(curr.chron[k],t.select,curr.chron[k+1]), curr.clim[k], curr.clim[k+1]) + + clim.interp[i,t.select.index,j] = clim.temp$NIGB[-c(1,length(clim.temp$NIGB))] + v.interp.all[i,t.all.select.index,j] = clim.temp$IGB + } + + # If there's no grid points just store the v value from the original data + if(length(t.select.index)==0 & length(t.all.select.index)>0) { + v.interp.all[i,t.all.select.index,j] = curr.v[k] } - } - } + } + + # Now sum over the interpolated v.interp.all values so as to get the correct interpolated v values + tgrid.select.lower = match(tgrid,tgrid.all) + tgrid.select.upper = c((tgrid.select.lower-1)[2:length(tgrid.select.lower)],length(tgrid.all)) + for(l in 1:(n.g-1)) { + v.interp[i,l,j] = sum(v.interp.all[i,tgrid.select.lower[l]:tgrid.select.upper[l],j]) + } + + if(any(is.na(clim.interp[i,,j]))) stop("Some interpolated climates have been given NAs") + if(any(is.na(v.interp[i,,j]))) stop("Some interpolated v values have been given NAs") + + } # End of samples loop + } # End of climate dimension loop - for(k in 1:Bclim.data$m) { - clim.interp[,, k] <- clim.interp[,, k] * sqrt(Bclim.data$ScVar)[k] + Bclim.data$ScMean[k] + for(j in 1:Bclim.data$m) { + clim.interp[,, j] <- clim.interp[,, j] * sqrt(Bclim.data$ScVar)[j] + Bclim.data$ScMean[j] + v.interp[,, j] <- sqrt(v.interp[,, j]) * sqrt(Bclim.data$ScVar)[j] } - + cat("Completed! \n") - return(list(clim.interp = clim.interp, vol.interp = vol.interp, time.grid=tgrid)) + return(list(clim.interp = clim.interp, v.interp = v.interp, time.grid=tgrid)) } diff --git a/R/BclimLayer.R b/R/BclimLayer.R index c8c6d74..f075f5e 100644 --- a/R/BclimLayer.R +++ b/R/BclimLayer.R @@ -49,6 +49,5 @@ function(pollen.loc,required.data3D,nsamples=1000){ #naming dimensions dimnames(samples)=list(NULL,NULL,c("GDD5","MTCO","AET/PET")) - #return(list(Posteriors=post,samples=samples,nslices=nslices,nsamples=nsamples)) return(samples) } diff --git a/R/BclimMCMC.R b/R/BclimMCMC.R index 721c6a1..02bcb17 100644 --- a/R/BclimMCMC.R +++ b/R/BclimMCMC.R @@ -1,13 +1,15 @@ BclimMCMC <- -function(Bclimdata,chron.loc,nchron=10000,control.mcmc=list(iterations=100000,burnin=20000,thinby=40,report=100),control.chains=list(v.mh.sd=0.8,vstart=rinvgauss(Bclimdata$n,2,1),Kstart=sample(1:Bclimdata$G,Bclimdata$n,replace=TRUE)),control.priors=list(eta=rep(2.66,Bclimdata$m),phi=rep(15.33,Bclimdata$m))) { +function(Bclimdata,chron.loc,nchron=10000,control.mcmc=list(iterations=100000,burnin=20000,thinby=40,report=100),control.chains=list(v.mh.sd=2,phi1.mh.sd=1,phi2.mh.sd=10,vstart=statmod::rinvgauss(Bclimdata$n-1,2,1),Zstart=sample(1:Bclimdata$G,Bclimdata$n,replace=TRUE),phi1start=rep(3,Bclimdata$m),phi2start=rep(20,Bclimdata$m)),control.priors=list(phi1dlmean=rep(1.275,Bclimdata$m),phi1dlsd=rep(0.076,Bclimdata$m),phi2dlmean=rep(4.231,Bclimdata$m),phi2dlsd=rep(0.271,Bclimdata$m))) { # Create output matrices remaining <- (control.mcmc$iterations-control.mcmc$burnin)/control.mcmc$thinby if(remaining!=as.integer(remaining)) stop("Iterations minus burnin divided by thinby must be an integer") vout <- rep(0,length=Bclimdata$m*(Bclimdata$n-1)*remaining) +zout <- rep(0,length=Bclimdata$m*Bclimdata$n*remaining) chronout <- rep(0,length=Bclimdata$n*remaining) cout <- rep(0,length=Bclimdata$m*(Bclimdata$n)*remaining) +phi1out = phi2out = rep(0,length=Bclimdata$m*remaining) # Re-dim the precisions matrix Bclimprec <- Bclimdata$tau.mat[,1,] @@ -31,24 +33,39 @@ out <- .C("BclimMCMC3D", as.integer(control.mcmc$report), as.double(control.chains$v.mh.sd), as.double(control.chains$vstart), - as.integer(control.chains$Kstart), - as.double(control.priors$eta), - as.double(control.priors$phi), + as.integer(control.chains$Zstart), + as.double(control.chains$phi1start), + as.double(control.chains$phi2start), as.double(vout), + as.integer(zout), as.double(chronout), - as.double(cout)) + as.double(cout), + as.double(phi1out), + as.double(phi2out), + as.double(control.priors$phi1dlmean), + as.double(control.priors$phi1dlsd), + as.double(control.priors$phi2dlmean), + as.double(control.priors$phi2dlsd), + as.double(control.chains$phi1.mh.sd), + as.double(control.chains$phi2.mh.sd) + ) #,PACKAGE="Bclim") +#browser() vout <- array(NA,dim=c(remaining,Bclimdata$n-1,Bclimdata$m)) cout <- array(NA,dim=c(remaining,Bclimdata$n,Bclimdata$m)) for(i in 1:remaining) { for(j in 1:Bclimdata$m) { vout[i,,j] <- out[[18]][seq(1,Bclimdata$n-1)+(j-1)*(Bclimdata$n-1)+(i-1)*(Bclimdata$n-1)*Bclimdata$m] - cout[i,,j] <- out[[20]][seq(1,Bclimdata$n)+(j-1)*(Bclimdata$n)+(i-1)*(Bclimdata$n)*Bclimdata$m] + cout[i,,j] <- out[[21]][seq(1,Bclimdata$n)+(j-1)*(Bclimdata$n)+(i-1)*(Bclimdata$n)*Bclimdata$m] } } -chronout <- matrix(out[[19]],ncol=Bclimdata$n,nrow=remaining,byrow=TRUE) +chronout <- matrix(out[[20]],ncol=Bclimdata$n,nrow=remaining,byrow=TRUE) +zout <- matrix(out[[19]],ncol=Bclimdata$n,nrow=remaining,byrow=TRUE) +phi1out = matrix(out[[22]],ncol=Bclimdata$m,nrow=remaining,byrow=TRUE) +phi2out = matrix(out[[23]],ncol=Bclimdata$m,nrow=remaining,byrow=TRUE) -return(list(v.store=vout,chron.store=chronout,c.store=cout,eta=control.priors$eta,phi=control.priors$phi,chron.loc=chron.loc,nchron=nchron)) + +return(list(v.store=vout,chron.store=chronout,c.store=cout,z.store=zout,zout=zout,chron.loc=chron.loc,nchron=nchron,phi1.store=phi1out,phi2.store=phi2out)) } diff --git a/R/BclimMixPar.R b/R/BclimMixPar.R index 6de0a7f..9cff26f 100644 --- a/R/BclimMixPar.R +++ b/R/BclimMixPar.R @@ -1,7 +1,7 @@ BclimMixPar <- function(MDP,G=10,num.cores=4,mixwarnings=FALSE) { - registerDoMC(cores=num.cores) +doMC::registerDoMC(cores=num.cores) # Calculate n.samp = number of samples, n = number of layers, m = number of climate dimensions n.samp <- dim(MDP)[1] @@ -13,7 +13,7 @@ ScVar <- rep(1,m) MDP2 <- MDP for(i in 1:m) { ScMean[i] <- mean(MDP[,,i]) - ScVar[i] <- max(var(MDP[,,i])) + ScVar[i] <- median(diag(var(MDP[,,i]))) MDP2[,,i] <- (MDP[,,i]-ScMean[i])/sqrt(ScVar[i]) } @@ -24,35 +24,18 @@ mu.mat <- array(NA,dim=c(n,m,G)) tau.mat <- array(NA,dim=c(n,m,G)) p.mat <- matrix(NA,nrow=n,ncol=G) -ans.all <- foreach(i = icount(n)) %dopar% { +ans.all <- foreach::foreach(i = icount(n)) %dopar% { cat("\r") cat("Completed:",format(round(100*i/n,2), nsmall = 2),"%") - Mclust(MDP2[,i,],G=G,modelNames="EII",warn=mixwarnings) + mclust::Mclust(MDP2[,i,],G=G,modelNames="EII",warn=mixwarnings) } -# ans.all <- list() -# for(i in 1:n) { -# cat("\r",i) -# ans.all[[i]] <- Mclust(MDP2[,i,],G=G,modelNames="EII",warn=mixwarnings) -# } for(i in 1:n) { - - #len <- G-ncol(ans$parameters$mean) - #if(length(len)==0) len <- 1 - # if(len==0) { - mu.mat[i,,] <- ans.all[[i]]$parameters$mean - for(g in 1:G) { - tau.mat[i,,g] <- 1/diag(ans.all[[i]]$parameters$variance$sigma[,,g]) - } - p.mat[i,] <- ans.all[[i]]$parameters$pro - # } else { - # ans <- Mclust(MDP2[,i,],G=nummix,modelName=ifelse(r>1,"EII","E"),warn=TRUE) - # MixMeans[,i,] <- ans$parameters$mean - # MixVar[i,] <- rep(ans$parameters$variance$sigmasq,nummix) - # MixPro[i,] <- ans$parameters$pro - #} - - #rm(ans) + mu.mat[i,,] <- ans.all[[i]]$parameters$mean + for(g in 1:G) { + tau.mat[i,,g] <- 1/diag(ans.all[[i]]$parameters$variance$sigma[,,g]) + } + p.mat[i,] <- ans.all[[i]]$parameters$pro # End of loop through n layers } diff --git a/R/BclimMixSer.R b/R/BclimMixSer.R index b9224fb..1572437 100644 --- a/R/BclimMixSer.R +++ b/R/BclimMixSer.R @@ -1,8 +1,6 @@ BclimMixSer <- function(MDP,G=10,mixwarnings=FALSE) { - -# registerDoMC(cores=num.cores) - + # Calculate n.samp = number of samples, n = number of layers, m = number of climate dimensions n.samp <- dim(MDP)[1] n <- dim(MDP)[2] @@ -13,7 +11,7 @@ function(MDP,G=10,mixwarnings=FALSE) { MDP2 <- MDP for(i in 1:m) { ScMean[i] <- mean(MDP[,,i]) - ScVar[i] <- max(var(MDP[,,i])) + ScVar[i] <- median(diag(var(MDP[,,i]))) MDP2[,,i] <- (MDP[,,i]-ScMean[i])/sqrt(ScVar[i]) } @@ -24,38 +22,23 @@ function(MDP,G=10,mixwarnings=FALSE) { tau.mat <- array(NA,dim=c(n,m,G)) p.mat <- matrix(NA,nrow=n,ncol=G) -# ans.all <- foreach(i = icount(n)) %dopar% { -# cat(i,"\r") -# Mclust(MDP2[,i,],G=G,modelName="EII",warn=mixwarnings) -# } ans.all <- list() for(i in 1:n) { cat("\r") cat("Completed:",format(round(100*i/n,2), nsmall = 2),"%") - ans.all[[i]] <- Mclust(MDP2[,i,],G=G,modelNames="EII",warn=mixwarnings) + ans.all[[i]] <- mclust::Mclust(MDP2[,i,],G=G,modelNames="EII",warn=mixwarnings) } for(i in 1:n) { - - #len <- G-ncol(ans$parameters$mean) - #if(length(len)==0) len <- 1 - # if(len==0) { mu.mat[i,,] <- ans.all[[i]]$parameters$mean for(g in 1:G) { tau.mat[i,,g] <- 1/diag(ans.all[[i]]$parameters$variance$sigma[,,g]) } p.mat[i,] <- ans.all[[i]]$parameters$pro - # } else { - # ans <- Mclust(MDP2[,i,],G=nummix,modelName=ifelse(r>1,"EII","E"),warn=TRUE) - # MixMeans[,i,] <- ans$parameters$mean - # MixVar[i,] <- rep(ans$parameters$variance$sigmasq,nummix) - # MixPro[i,] <- ans$parameters$pro - #} - - #rm(ans) - # End of loop through n layers + # End of loop through n layers } + # Now output everything to a nice neat list Bclimdata <- list(MDP=MDP2,n=n,m=m,n.samp=n.samp,ScMean=ScMean,ScVar=ScVar,G=G,mu.mat=mu.mat,tau.mat=tau.mat,p.mat=p.mat) cat("\n") diff --git a/R/NIGB.R b/R/NIGB.R deleted file mode 100644 index b1b083d..0000000 --- a/R/NIGB.R +++ /dev/null @@ -1,46 +0,0 @@ -NIGB <- -function(delta, IG.sum, tb.points, c.start, c.end){ - total.points <- length(tb.points) - NIG.bridge <- rep(NaN, nrow=total.points) - NIG.bridge[1] <- c.start - NIG.bridge[total.points] <- c.end - IG.increment <- rep(NaN, total.points-2) - z <- IG.sum - - l <- 2 - while(l < total.points){ - # Generate chi-square random variate - q <- rnorm(1,0,1)^2 - - # Reparameterise - mu <- (tb.points[total.points]-tb.points[l]) / (tb.points[l]-tb.points[l-1]) - lambda <- (1/delta^2 * (tb.points[total.points]-tb.points[l])^2) / z - - if(z==0) { - print(c(delta,IG.sum,tb.points,c.start,c.end)) - stop("Problem in Inverse Gaussian Bridge") - } - - - # Compute the roots of the chi-square random variate - s1 <- mu + (mu^2*q)/(2*lambda) - (mu/(2*lambda)) * sqrt(4*mu*lambda*q + mu^2*q^2) - if(lambda<1e-5) s1 <- mu - s2 <- mu^2 / s1 - - # Acceptance/rejection of root - s <- ifelse(runif(1,0,1) < mu*(1+s1)/((1+mu)*(mu+s1)), s1, s2) - - - # Compute the IG incrrement - IG.increment[l-1] <- z / (1 + s) - - # Rescale the sum of left-over distance of the IG bridge - z <- z - IG.increment[l-1] - - #Compute the current value of the NIG bridge - NIG.bridge[l] <- (c.start*z + c.end*(IG.sum-z)) / IG.sum + - rnorm(1, 0, 1) * (IG.sum-z)*z / IG.sum - l <- l +1 - } - return(list(IGB = c(IG.increment, (IG.sum - sum(IG.increment))), NIGB = t(NIG.bridge))) -} diff --git a/R/NIGp3.R b/R/NIGp3.R deleted file mode 100644 index 3283bb4..0000000 --- a/R/NIGp3.R +++ /dev/null @@ -1,18 +0,0 @@ -NIGp3 <- -function(begin, alpha, delta, totalCount, distance) -{ - dt <- distance/(totalCount-1) - ig_dt <- rinvgauss(totalCount-1, mu = 1*dt/alpha, lambda=1*dt^2) - NIGprocess <- B_dt <- rep(0,totalCount) - NIGprocess[1] <- begin - - for (i in 2:totalCount) { - z <- rnorm(1,0,1) - B_dt[i] <- B_dt[i-1] + sqrt(ig_dt[i-1])*z #Time change of a standard Brownian motion - NIGprocess[i] <- sum(NIGprocess[i-1], 1*B_dt[i]) #delta*B_dt[i] - - #See page 9 of http://iriaf.univ-poitiers.fr/colloque2011/article/v1s1a3.pdf - #A two factor Levy model for stochastic mortality - Viou Ainou - } - return(list(IGB = cumsum(ig_dt), NIGB = NIGprocess)) -} diff --git a/R/plotBclim.R b/R/plotBclim.R index fa814f2..d38af08 100644 --- a/R/plotBclim.R +++ b/R/plotBclim.R @@ -1,17 +1,26 @@ plotBclim <- function(x,dim=1,title=NULL,presentleft=TRUE,blob=TRUE,MDPcol="blue",denscol="red",MDPtransp=0.1,denstransp=0.5,leg=TRUE,legloc="topleft",...) { - #dim=2;title=NULL;presentleft=TRUE;blob=TRUE;MDPcol="blue";denscol="red";MDPtransp=0.1;denstransp=0.5;leg=TRUE;legloc="topleft" - if(class(x)!="Bclim") stop("Needs a Bclim output object") - # Set up plot - par(mar=c(4,4,3,1)) - xrange <- range(c(0,x$time.grid)) - if(!presentleft) xrange <- rev(xrange) - yrange <- range(c(0,x$MDP[,,dim])) - mytitle <- title - if(is.null(title)) mytitle <- paste(x$core.name,": ",x$clim.dims[dim],sep="") - plot(1,1,type="n",xlim=xrange,ylim=yrange,xlab="Age (k cal years BP)",ylab=x$clim.dims[dim],las=1,bty="n",main=mytitle) - #plot(1,1,type="n",xlim=xrange,ylim=yrange,xlab="Age (thousands of years before present)",ylab=x$clim.dims[dim],las=1,bty="n",main=mytitle) + if(class(x)!="Bclim") stop("Needs a Bclim output object") + # Set up plot + par(mar=c(4,4,3,1)) + xrange <- range(c(0,x$time.grid)) + if(!presentleft) xrange <- rev(xrange) + yrange <- range(c(0,x$MDP[,,dim])) + mytitle <- title + if(is.null(title)) mytitle <- paste(x$core.name,": ",x$clim.dims[dim],sep="") + + if(dim==1) { + plot(1,1,type="n",xlim=xrange,ylim=yrange,xlab="Age (k cal years BP)",ylab=expression(paste("GDD5 (",degree,"C days)",sep="")),las=1,bty="n",main=mytitle,xaxt='n') + } + if(dim==2) { + plot(1,1,type="n",xlim=xrange,ylim=yrange,xlab="Age (k cal years BP)",ylab=expression(paste("MTCO (",degree,"C)",sep="")),las=1,bty="n",main=mytitle,xaxt='n') + } + if(dim==3) { + plot(1,1,type="n",xlim=xrange,ylim=yrange,xlab="Age (k cal years BP)",ylab=x$clim.dims[3],las=1,bty="n",main=mytitle,xaxt='n') + } + + axis(side=1,at=pretty(x$time.grid,n=10)) rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],col="lightgray",border="NA") grid(col="white") @@ -34,7 +43,7 @@ function(x,dim=1,title=NULL,presentleft=TRUE,blob=TRUE,MDPcol="blue",denscol="re chrontemp <- chron[1:num,i] if(var(chrontemp)>0) { - tmp <- kde2d(chrontemp,MDPtemp) + tmp <- MASS::kde2d(chrontemp,MDPtemp) # Standardise and find 95% limit z <- tmp$z/sum(tmp$z) @@ -61,7 +70,7 @@ function(x,dim=1,title=NULL,presentleft=TRUE,blob=TRUE,MDPcol="blue",denscol="re HDR50 <- matrix(NA,nrow=length(x$time.grid),ncol=50) for(i in 1:length(x$time.grid)) { if(sd(x$clim.interp[,i,dim])>0) { - TempHDR <- hdr(x$clim.interp[,i,dim],h=bw.nrd0(x$clim.interp[,i,dim]),prob=c(1,50,75,95))$hdr + TempHDR <- hdrcde::hdr(x$clim.interp[,i,dim],h=bw.nrd0(x$clim.interp[,i,dim]),prob=c(1,50,75,95))$hdr Temp95 <- TempHDR[1,!is.na(TempHDR[1,])] HDR95[i,1:length(Temp95)] <- Temp95 @@ -101,7 +110,7 @@ function(x,dim=1,title=NULL,presentleft=TRUE,blob=TRUE,MDPcol="blue",denscol="re lines(c(x$time.grid[i],x$time.grid[i]),c(HDR75[i,2*j-1],HDR75[i,2*j]),col="blue",lwd=2) } for(j in 1:(ncol(HDR50)/2)) { - lines(c(x$time.grid[i],x$time.grid),c(HDR50[i,2*j-1],HDR50[i,2*j]),col="blue",lwd=3) + lines(c(x$time.grid[i],x$time.grid[i]),c(HDR50[i,2*j-1],HDR50[i,2*j]),col="blue",lwd=3) } } } diff --git a/R/plotBclimVol.R b/R/plotBclimVol.R index f6c67ca..258eecd 100644 --- a/R/plotBclimVol.R +++ b/R/plotBclimVol.R @@ -1,13 +1,12 @@ plotBclimVol <- -function(x,dim=1,title=NULL,presentleft=TRUE,denscol="red",denstransp=0.5,leg=TRUE,dolog10=FALSE,med=FALSE,legloc="topleft",...) { - #dim=1;title=NULL;presentleft=TRUE;denscol="red";denstransp=0.5;leg=TRUE;legloc="topleft" +function(x,dim=1,title=NULL,presentleft=TRUE,denscol="red",denstransp=0.5,leg=TRUE,mean=TRUE,legloc="topleft",...) { if(class(x)!="Bclim") stop("Needs a Bclim output object") # Create HDRs for volatility errorbar <- matrix(NA,nrow=length(x$time.grid)-1,ncol=3) - if(dolog10) for(i in 1:(length(x$time.grid)-1)) errorbar[i,] <- quantile(log10(abs(x$vol.interp[,i,dim]*sqrt(x$ScVar[dim]))),probs=c(0.025,0.5,0.975)) - if(!dolog10) for(i in 1:(length(x$time.grid)-1)) errorbar[i,] <- quantile(abs(x$vol.interp[,i,dim]*sqrt(x$ScVar[dim])),probs=c(0.025,0.5,0.975)) + for(i in 1:(length(x$time.grid)-1)) errorbar[i,] <- quantile(x$v.interp[,i,dim],probs=c(0.025,0.5,0.975)) + for(i in 1:(length(x$time.grid)-1)) errorbar[i,2] <- mean(x$v.interp[,i,dim]) # Sort out colours tmp <- col2rgb(denscol) @@ -21,20 +20,35 @@ function(x,dim=1,title=NULL,presentleft=TRUE,denscol="red",denstransp=0.5,leg=TR yrange <- range(c(0,as.vector(errorbar))) mytitle <- title if(is.null(title)) mytitle <- paste(x$core.name,": ",x$clim.dims[dim],sep="") - plot(1,1,type="n",xlim=xrange,ylim=yrange,xlab="Age (k cal years BP)",ylab=paste(x$clim.dims[dim],ifelse(dolog10,"log10 volatility","volatility")),las=1,bty="n",main=mytitle) + + if(dim==1) { + plot(1,1,type="n",xlim=xrange,ylim=yrange,xaxt='n',xlab="Age (k cal years BP)",ylab=expression(paste("GDD5 volatility (",degree,"C days)",sep="")),las=1,bty="n",main=mytitle) + } + if(dim==2) { + plot(1,1,type="n",xlim=xrange,ylim=yrange,xaxt='n',xlab="Age (k cal years BP)",ylab=expression(paste("MTCO volatility (",degree,"C)",sep="")),las=1,bty="n",main=mytitle) + } + if(dim==3) { + plot(1,1,type="n",xlim=xrange,ylim=yrange,xaxt='n',xlab="Age (k cal years BP)",ylab=x$clim.dims[dim],las=1,bty="n",main=mytitle) + } + + axis(side=1,at=pretty(x$time.grid,n=10)) rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],col="lightgray",border="NA") grid(col="white") # Draw lines - #lines(x$time.grid[-1],errorbar[,1],col="red") - #lines(x$time.grid[-1],errorbar[,2],col="blue") - #lines(x$time.grid[-1],errorbar[,3],col="green") - polygon(c(x$time.grid[-1],rev(x$time.grid[-1])),c(errorbar[,1],rev(errorbar[,3])),col=mycol2,border=mycol2) - if(med) lines(x$time.grid[-1],apply(errorbar,1,"median"),col=denscol) + tgridshift = x$time.grid[-1]-c(diff(x$time.grid)/2) + lines(c(min(x$time.grid),max(x$time.grid)),c(min(yrange),min(yrange))) + for(i in 1:length(x$time.grid[-1])) { + lines(c(tgridshift[i],tgridshift[i]),c(errorbar[i,1],errorbar[i,3]),col=denscol) + if(mean) points(tgridshift[i],errorbar[i,2],col=denscol) + } + for(i in 1:length(x$time.grid)) lines(c(x$time.grid[i],x$time.grid[i]),c(min(yrange),0.05*max(yrange))) # Finally draw a legend if(leg==TRUE) { - legend(legloc,legend=c("95% credibility interval"),fill=c(denscol),bty="n") + if(!mean) legend(legloc,legend=c("95% credibility interval",'Time Interval'),lty=c(1,1),col=c(mycol2,'black'),pch=c(-1,-1),bty="n") + if(mean) legend(legloc,legend=c("95% credibility interval","Mean",'Time Interval'),lty=c(1,-1,1),col=c(mycol2,denscol,'black'),pch=c(-1,1,-1),bty="n") + #fill=c(denscol,NULL) } # End of function diff --git a/R/zzz.R b/R/zzz.R new file mode 100755 index 0000000..c6f9f5d --- /dev/null +++ b/R/zzz.R @@ -0,0 +1,7 @@ +.onAttach <- function(libname, pkgname) { + Bclimver <- read.dcf(file=system.file("DESCRIPTION", package=pkgname), + fields="Version") + packageStartupMessage(paste(pkgname, Bclimver)) + packageStartupMessage("Welcome to Bclim. Type help(Bclim) to get started.") + packageStartupMessage("See http://mathsci.ucd.ie/~parnell_a/Rpack/Bclim.htm for updates, bugs and a tutorial.") +} diff --git a/man/Bclim-package.Rd b/man/Bclim-package.Rd index 29bcead..45e809d 100644 --- a/man/Bclim-package.Rd +++ b/man/Bclim-package.Rd @@ -12,8 +12,8 @@ This package takes pollen core and chronological data to produce palaeoclimate e \tabular{ll}{ Package: \tab Bclim\cr Type: \tab Package\cr -Version: \tab 2.2\cr -Date: \tab 2013-01-23\cr +Version: \tab 2.3.1\cr +Date: \tab 2013-07-29\cr License: \tab GPL >=2 \cr } See the function \code{\link{BclimRun}} for a full description of the workings of the package and an example script for getting it to work. diff --git a/man/BclimCompile.Rd b/man/BclimCompile.Rd index dd1520a..260c039 100644 --- a/man/BclimCompile.Rd +++ b/man/BclimCompile.Rd @@ -81,8 +81,10 @@ load('required.data3D.RData') ## note that all of these functions have further options you can change with step1 <- BclimLayer(pollen.loc,required.data3D=required.data3D) -step2 <- BclimMixSer(step1) # See also the parallelised version BclimMixPar if you have doMC and foreach installed -step3 <- BclimMCMC(step2,chron.loc) # You should probably do some convergence checking after this step +step2 <- BclimMixSer(step1) +# See also the parallelised version BclimMixPar if you have doMC and foreach installed +step3 <- BclimMCMC(step2,chron.loc) +# You should probably do some convergence checking after this step step4 <- BclimInterp(step2,step3) results <- BclimCompile(step1,step2,step3,step4,core.name="Sluggan Moss") diff --git a/man/BclimLayer.Rd b/man/BclimLayer.Rd index f32436f..ca208dc 100644 --- a/man/BclimLayer.Rd +++ b/man/BclimLayer.Rd @@ -64,8 +64,10 @@ load('required.data3D.RData') ## note that all of these functions have further options you can change with step1 <- BclimLayer(pollen.loc,required.data3D=required.data3D) -step2 <- BclimMixSer(step1) # See also the parallelised version BclimMixPar if you have doMC and foreach installed -step3 <- BclimMCMC(step2,chron.loc) # You should probably do some convergence checking after this step +step2 <- BclimMixSer(step1) +# See also the parallelised version BclimMixPar if you have doMC and foreach installed +step3 <- BclimMCMC(step2,chron.loc) +# You should probably do some convergence checking after this step step4 <- BclimInterp(step2,step3) results <- BclimCompile(step1,step2,step3,step4,core.name="Sluggan Moss") diff --git a/man/BclimMCMC.Rd b/man/BclimMCMC.Rd index a544897..aa22649 100644 --- a/man/BclimMCMC.Rd +++ b/man/BclimMCMC.Rd @@ -8,7 +8,13 @@ Bclim Markov chain Monte Carlo run A Bclim Markov chain Monte Carlo (MCMC) run to determine the volatilities and climate from a set of marginal data posteriors approximated as mixtures. } \usage{ -BclimMCMC(Bclimdata, chron.loc, nchron = 10000, control.mcmc = list(iterations = 1e+05, burnin = 20000, thinby = 40, report = 100), control.chains=list(v.mh.sd=0.8,vstart=rinvgauss(Bclimdata$n,2,1), Kstart=sample(1:Bclimdata$G,Bclimdata$n,replace=TRUE)), control.priors = list(eta = rep(2.66, Bclimdata$m), phi = rep(15.33, Bclimdata$m))) +BclimMCMC(Bclimdata,chron.loc,nchron=10000,control.mcmc=list(iterations=100000, +burnin=20000,thinby=40,report=100),control.chains=list(v.mh.sd=2,phi1.mh.sd=1, +phi2.mh.sd=10,vstart=statmod::rinvgauss(Bclimdata$n-1,2,1), +Zstart=sample(1:Bclimdata$G,Bclimdata$n,replace=TRUE),phi1start=rep(3,Bclimdata$m), +phi2start=rep(20,Bclimdata$m)),control.priors=list(phi1dlmean=rep(1.275,Bclimdata$m), +phi1dlsd=rep(0.076,Bclimdata$m),phi2dlmean=rep(4.231,Bclimdata$m), +phi2dlsd=rep(0.271,Bclimdata$m))) } \arguments{ \item{Bclimdata}{ @@ -24,10 +30,10 @@ The number of chronologies in the chronologies file A list containing elements that control the MCMC, including the number of iterations, the size of the burn-in period, the amount to thinby, and how often for the algorithm to report its progress. } \item{control.chains}{ -A list containing elements that control the starting values of the parameters (vstart and Kstart) and the Metropolis-Hastings proposal standard deviation for v. +A list containing elements that control the starting values of the parameters (vstart, Zstart, phi1start and phi2) and the Metropolis-Hastings proposal standard deviation for v, phi1 and phi2. } \item{control.priors}{ -A list containing the prior parameters for the volatilities, givem by eta and phi, both of which should be 3-vectors. +A list containing the prior parameters for the volatilities, given by phi1 and phi2, both of which should be the log-mean and log-sd of the . } } \details{ @@ -38,8 +44,9 @@ A list object wiht the following elements \item{v.store }{Samples of the posterior estimated volatilities} \item{chron.store }{Samples of the used chronologies} \item{c.store }{Samples of the posterior estimated climates} - \item{eta }{Values used for the IG prior on v for each climate dimension} - \item{phi }{Values used for the IG prior on v for each climate dimension} + \item{z.store }{Samples of the posterior mixture indices} + \item{phi1 }{Values used for the IG prior on v for each climate dimension} + \item{phi2 }{Values used for the IG prior on v for each climate dimension} \item{chron.loc }{A character string giving the location of the chronology file} \item{nchron }{The number of chronologies in the chronology file} } @@ -59,15 +66,15 @@ The main Bclim function is \code{\link{BclimRun}}. See also the other 3 stages: # Set the working directory using setwd (not shown) # Download and load in the response surfaces: -url1 <- 'http://mathsci.ucd.ie/~parnell_a/required.data3D.RData' +url1 <- 'http://mathsci.ucd.ie/~parnell_a/media/requireddata3D.RData' download.file(url1,'required_data3D.RData') # and now the pollen -url2 <- 'http://mathsci.ucd.ie/~parnell_a/SlugganPollen.txt' +url2 <- 'http://mathsci.ucd.ie/~parnell_a/media/SlugganPollen.txt' download.file(url2,'SlugganPollen.txt') # and finally the chronologies -url3 <- 'http://mathsci.ucd.ie/~parnell_a/Sluggan_2chrons.txt' +url3 <- 'http://mathsci.ucd.ie/~parnell_a/media/Sluggan_2chrons.txt' download.file(url3,'Slugganchrons.txt') # Create variables which state the locations of the pollen and chronologies @@ -75,12 +82,14 @@ pollen.loc <- paste(getwd(),'/SlugganPollen.txt',sep='') chron.loc <- paste(getwd(),'/Slugganchrons.txt',sep='') # Load in the response surfaces -load('required.data3D.RData') +load('required_data3D.RData') ## note that all of these functions have further options you can change with step1 <- BclimLayer(pollen.loc,required.data3D=required.data3D) -step2 <- BclimMixSer(step1) # See also the parallelised version BclimMixPar if you have doMC and foreach installed -step3 <- BclimMCMC(step2,chron.loc) # You should probably do some convergence checking after this step +step2 <- BclimMixSer(step1) +# See also the parallelised version BclimMixPar if you have doMC and foreach installed +step3 <- BclimMCMC(step2,chron.loc) +# You should probably do some convergence checking after this step step4 <- BclimInterp(step2,step3) results <- BclimCompile(step1,step2,step3,step4,core.name="Sluggan Moss") diff --git a/man/BclimMixPar.Rd b/man/BclimMixPar.Rd index e6793d4..9765c4d 100644 --- a/man/BclimMixPar.Rd +++ b/man/BclimMixPar.Rd @@ -74,8 +74,10 @@ load('required.data3D.RData') ## note that all of these functions have further options you can change with step1 <- BclimLayer(pollen.loc,required.data3D=required.data3D) -step2 <- BclimMixPar(step1) # See also the serial version BclimMixSer if you cannot install doMC and foreach -step3 <- BclimMCMC(step2,chron.loc) # You should probably do some convergence checking after this step +step2 <- BclimMixPar(step1) +# See also the serial version BclimMixSer if you cannot install doMC and foreach +step3 <- BclimMCMC(step2,chron.loc) +# You should probably do some convergence checking after this step step4 <- BclimInterp(step2,step3) results <- BclimCompile(step1,step2,step3,step4,core.name="Sluggan Moss") diff --git a/man/BclimMixSer.Rd b/man/BclimMixSer.Rd index 0651a1a..dae5271 100644 --- a/man/BclimMixSer.Rd +++ b/man/BclimMixSer.Rd @@ -71,8 +71,10 @@ load('required.data3D.RData') ## note that all of these functions have further options you can change with step1 <- BclimLayer(pollen.loc,required.data3D=required.data3D) -step2 <- BclimMixSer(step1) # See also the parallelised version BclimMixPar if you have doMC and foreach installed -step3 <- BclimMCMC(step2,chron.loc) # You should probably do some convergence checking after this step +step2 <- BclimMixSer(step1) +# See also the parallelised version BclimMixPar if you have doMC and foreach installed +step3 <- BclimMCMC(step2,chron.loc) +# You should probably do some convergence checking after this step step4 <- BclimInterp(step2,step3) results <- BclimCompile(step1,step2,step3,step4,core.name="Sluggan Moss") diff --git a/man/BclimRun.Rd b/man/BclimRun.Rd index 3dacae1..410e155 100644 --- a/man/BclimRun.Rd +++ b/man/BclimRun.Rd @@ -7,7 +7,8 @@ Run all stages of Bclim together This function is intended for beggining/light users of Bclim who wish to create palaeoclimate reconstructions from their own pollen data. It collects together the functions \code{\link{BclimLayer}}, \code{\link{BclimMixSer}}, \code{\link{BclimMCMC}}, \code{\link{BclimInterp}} and \code{\link{BclimCompile}}. } \usage{ -BclimRun(pollen.loc, chron.loc, core.name = "Core", time.grid = seq(0, 14, length = 100), required.data3D = NULL, nchrons = 10000, parallel = FALSE, save.as.you.go = TRUE) +BclimRun(pollen.loc, chron.loc, core.name = "Core", time.grid = seq(0, 14, length = 100), + required.data3D = NULL, nchrons = 10000, parallel = FALSE, save.as.you.go = TRUE) } \arguments{ \item{pollen.loc}{ @@ -90,7 +91,8 @@ chron.loc <- paste(getwd(),'/Slugganchrons.txt',sep='') load('required.data3D.RData') # Now run Bclim the simple way -RunSluggan <- BclimRun(pollen.loc,chron.loc,core.name="Sluggan Moss",time.grid=seq(0,14,length=100),nchrons=10000,required.data3D=required.data3D,parallel=FALSE) +RunSluggan <- BclimRun(pollen.loc,chron.loc,core.name="Sluggan Moss",time.grid=seq(0,14,length=100), + nchrons=10000,required.data3D=required.data3D,parallel=FALSE) # Create some plots of climate plotBclim(RunSluggan,dim=1) diff --git a/man/NIGB.Rd b/man/NIGB.Rd deleted file mode 100644 index 89106b5..0000000 --- a/man/NIGB.Rd +++ /dev/null @@ -1,46 +0,0 @@ -\name{NIGB} -\alias{NIGB} -\title{ -Normal-Inverse Gaussian bridging function -} -\description{ -Normal-Inverse Gaussian bridging function -} -\usage{ -NIGB(delta, IG.sum, tb.points, c.start, c.end) -} -\arguments{ - \item{delta}{ -IG parameter 1 -} - \item{IG.sum}{ -Sum of IG values -} - \item{tb.points}{ -Time points over which to interpolate -} - \item{c.start}{ -Climate start value -} - \item{c.end}{ -Climate end value -} -} -\details{ -Not meant for use outside of \code{\link{BclimInterp}}. -} -\value{ -Interpolated NIG values returned to \code{\link{BclimInterp}}. -} -\author{ -Thinh K Doan -} -\seealso{ -\code{\link{Bclim}} -} -\examples{ -## See other Bclim functions -} -\keyword{ model } -\keyword{ multivariate } -\keyword{ smooth } diff --git a/man/NIGp3.Rd b/man/NIGp3.Rd deleted file mode 100644 index 975525b..0000000 --- a/man/NIGp3.Rd +++ /dev/null @@ -1,46 +0,0 @@ -\name{NIGp3} -\alias{NIGp3} -\title{ -Normal inverse Gaussian process -} -\description{ -Normal inverse Gaussian process. Not meant for use outside of other Bclim functions. -} -\usage{ -NIGp3(begin, alpha, delta, totalCount, distance) -} -\arguments{ - \item{begin}{ -Starting value -} - \item{alpha}{ -IG parameter -} - \item{delta}{ -IG parameter -} - \item{totalCount}{ -Total Count -} - \item{distance}{ -Distance -} -} -\details{ -Not meant for us outside of other Bclim functions. See \code{\link{Bclim}} for information. -} -\value{ -A set of values from the required NIG process -} -\author{ -Thinh K Doan -} -\seealso{ -\code{\link{Bclim}} -} -\examples{ -## See other Bclim functions -} -\keyword{ model } -\keyword{ multivariate } -\keyword{ smooth } diff --git a/man/plotBclim.Rd b/man/plotBclim.Rd index 354f1b4..00d3265 100644 --- a/man/plotBclim.Rd +++ b/man/plotBclim.Rd @@ -7,7 +7,8 @@ Plots of posterior Bclim climate histories Create plots of climate histories from a Bclim run } \usage{ -plotBclim(x, dim = 1, title = NULL, presentleft = TRUE, blob = TRUE, MDPcol = "blue", denscol = "red", MDPtransp = 0.1, denstransp = 0.5, leg = TRUE, legloc = "topleft", ...) +plotBclim(x, dim = 1, title = NULL, presentleft = TRUE, blob = TRUE, MDPcol = "blue", + denscol = "red", MDPtransp = 0.1, denstransp = 0.5, leg = TRUE, legloc = "topleft", ...) } \arguments{ \item{x}{ diff --git a/man/plotBclimVol.Rd b/man/plotBclimVol.Rd index 77938c6..2bf4b6c 100644 --- a/man/plotBclimVol.Rd +++ b/man/plotBclimVol.Rd @@ -7,7 +7,8 @@ Plotting of interpolated climate volatilities This function plots volatilities over time for the selected climate dimension after a Bclim run. } \usage{ -plotBclimVol(x, dim = 1, title = NULL, presentleft = TRUE, denscol = "red", denstransp = 0.5, leg = TRUE, dolog10 = FALSE, med = FALSE, legloc = "topleft", ...) +plotBclimVol(x, dim = 1, title = NULL, presentleft = TRUE, denscol = "red", + denstransp = 0.5, leg = TRUE, mean = TRUE, legloc = "topleft", ...) } \arguments{ \item{x}{ @@ -31,11 +32,8 @@ The transparency (between 0 and 1) of the 95\% volatility ranges \item{leg}{ Whether to include a legend or not (default TRUE) } - \item{dolog10}{ -Whether to plot the log of the volatilities instead -} - \item{med}{ -Whether to include a line indicating the median volatility + \item{mean}{ +Whether to include a line indicating the mean volatility } \item{legloc}{ Where to place the legend (e.g. 'topleft', 'bottomright', etc) diff --git a/src/BclimMCMC3D.c b/src/BclimMCMC3D.c index 050a9a6..64389ee 100644 --- a/src/BclimMCMC3D.c +++ b/src/BclimMCMC3D.c @@ -17,7 +17,7 @@ /////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////// -void BclimMCMC3D(int *G,int *n, int *m,int *nchrons, double *MixPro, double *MixMeans, double *MixPrec, char **ChronsFILE, int *iterations, int *burnin, int *thinby, int *reportevery, double *vmhsd, double *vstart, int *Kstart, double *muval, double *phival, double *vstore, double *chronstore, double *cstore) +void BclimMCMC3D(int *G,int *n, int *m,int *nchrons, double *MixPro, double *MixMeans, double *MixPrec, char **ChronsFILE, int *iterations, int *burnin, int *thinby, int *reportevery, double *vmhsd, double *vstart, int *Zstart, double *phi1start, double *phi2start, double *vstore, int *zstore, double *chronstore, double *cstore, double *phi1store, double *phi2store, double *phi1dlmean, double *phi1dlsd, double *phi2dlmean, double *phi2dlsd, double *phi1mhsd, double *phi2mhsd) { // G is number of mixture groups, n is number of layers, m is number of climate dimensions (always 3 here), nchrons is number of chronologies @@ -25,10 +25,10 @@ void BclimMCMC3D(int *G,int *n, int *m,int *nchrons, double *MixPro, double *Mix // ChronsFILE is the path to the chronologies file // iterations is number of iterations, burnin is size of burnin, thinby is thinning amount, reportevery is how often to report // vmhsd is the standard deviation of the truncated random walk for proposing new values of v -// vstart and Kstart are starting values for v and K respectively -// mu and phi are 3-vectors of the values used for mu and phi +// vstart and Zstart are starting values for v and Z respectively +// phi1 and phi2 are 3-vectors of the values used for phi1 and phi2 // vstore is variance output, chronstore are used chronologies, cstore is used for climates - + // Declare indicators int i,j,k,l,j2,j3,iter; @@ -65,13 +65,14 @@ for(i=0;i<*n;i++) { for(i=0;i<*n;i++) for(j=0;j<*G;j++) MyMixPro[i][j] = MixPro[j*(*n)+i]; // Set up matrices -double **v,*mu,*phi,*mixpro,*choldiag,*chollower,*rannormal,*cholzeros,*mvnvar; -int *K,*mixseq; +double **v,*phi1,*phi2,*mixpro,*choldiag,*chollower,*rannormal,*cholzeros,*mvnvar; +int *Z,*Zstar,*mixseq; v = (double **)calloc(*n-1, sizeof(double *)); for(i=0;i<*n-1;i++) v[i] = (double *)calloc(*m, sizeof(double)); -mu = (double *)calloc(*m, sizeof(double)); -phi = (double *)calloc(*m, sizeof(double)); -K = (int *)calloc(*n, sizeof(int)); +phi1 = (double *)calloc(*m, sizeof(double)); +phi2 = (double *)calloc(*m, sizeof(double)); +Z = (int *)calloc(*n, sizeof(int)); +Zstar = (int *)calloc(*n, sizeof(int)); mixpro = (double *)calloc(*G, sizeof(double)); for(i=0;i<*G;i++) mixpro[i]=1/((double)*G); mixseq = (int *)calloc(*G, sizeof(int)); @@ -83,9 +84,10 @@ cholzeros = (double *)calloc(*n, sizeof(double)); mvnvar = (double *)calloc(*n, sizeof(double)); if(v==NULL) error("Can't allocate memory"); -if(mu==NULL) error("Can't allocate memory"); -if(phi==NULL) error("Can't allocate memory"); -if(K==NULL) error("Can't allocate memory"); +if(phi1==NULL) error("Can't allocate memory"); +if(phi2==NULL) error("Can't allocate memory"); +if(Z==NULL) error("Can't allocate memory"); +if(Zstar==NULL) error("Can't allocate memory"); if(mixseq==NULL) error("Can't allocate memory"); if(choldiag==NULL) error("Can't allocate memory"); if(chollower==NULL) error("Can't allocate memory"); @@ -95,29 +97,40 @@ if(mvnvar==NULL) error("Can't allocate memory"); // Give starting values for(i=0;i<*n-1;i++) for(j=0;j<*m;j++) v[i][j] = vstart[i]; -for(i=0;i<*m;i++) mu[i] = muval[i]; -for(i=0;i<*m;i++) phi[i] = phival[i]; -for(i=0;i<*n;i++) K[i] = Kstart[i]-1; +for(i=0;i<*m;i++) phi1[i] = phi1start[i]; +for(i=0;i<*m;i++) phi2[i] = phi2start[i]; +for(i=0;i<*n;i++) Z[i] = Zstart[i]-1; // Set up the matrices I'll need; -double *mujk,*Djk,*VDmu,*Vz,*mujkstar,*Djkstar,*VDmustar,*Dmu,*Dmustar; -mujk = (double *)calloc(*n, sizeof(double)); -Djk = (double *)calloc(*n, sizeof(double)); -VDmu = (double *)calloc(*n, sizeof(double)); -Vz = (double *)calloc(*n, sizeof(double)); -mujkstar = (double *)calloc(*n, sizeof(double)); -Djkstar = (double *)calloc(*n, sizeof(double)); -VDmustar = (double *)calloc(*n, sizeof(double)); -Dmu = (double *)calloc(*n, sizeof(double)); -Dmustar = (double *)calloc(*n, sizeof(double)); - -if(mujk==NULL) error("Can't allocate memory"); -if(Djk==NULL) error("Can't allocate memory"); -if(VDmu==NULL) error("Can't allocate memory"); -if(Vz==NULL) error("Can't allocate memory"); -if(mujkstar==NULL) error("Can't allocate memory"); -if(Djkstar==NULL) error("Can't allocate memory"); -if(VDmustar==NULL) error("Can't allocate memory"); +double **M, **R, **DM, **Mstar, **Rstar, **DMstar; +M = (double **)calloc(*n, sizeof(double *)); +for(i=0;i<*n;i++) M[i] = (double *)calloc(*m, sizeof(double)); // n rows, m columns - hopefully! +R = (double **)calloc(*n, sizeof(double *)); +for(i=0;i<*n;i++) R[i] = (double *)calloc(*m, sizeof(double)); +DM = (double **)calloc(*n, sizeof(double *)); +for(i=0;i<*n;i++) DM[i] = (double *)calloc(*m, sizeof(double)); +Mstar = (double **)calloc(*n, sizeof(double *)); +for(i=0;i<*n;i++) Mstar[i] = (double *)calloc(*m, sizeof(double)); +Rstar = (double **)calloc(*n, sizeof(double *)); +for(i=0;i<*n;i++) Rstar[i] = (double *)calloc(*m, sizeof(double)); +DMstar = (double **)calloc(*n, sizeof(double *)); +for(i=0;i<*n;i++) DMstar[i] = (double *)calloc(*m, sizeof(double)); + +// Some final vectors I need +double *MDR,*MDM,*MDRstar,*MDMstar,*D,*Wupper,*Wdiag; +D = (double *)calloc(*n, sizeof(double)); +MDR = (double *)calloc(*m, sizeof(double)); +MDM = (double *)calloc(*m, sizeof(double)); +MDRstar = (double *)calloc(*m, sizeof(double)); +MDMstar = (double *)calloc(*m, sizeof(double)); + +if(M==NULL) error("Can't allocate memory"); +if(D==NULL) error("Can't allocate memory"); +if(R==NULL) error("Can't allocate memory"); +if(DM==NULL) error("Can't allocate memory"); +if(Mstar==NULL) error("Can't allocate memory"); +if(Rstar==NULL) error("Can't allocate memory"); +if(DMstar==NULL) error("Can't allocate memory"); // Set up somewhere to store output double **allchrons; @@ -143,27 +156,57 @@ currchron = (double *)calloc(*n, sizeof(double)); if(diffchron==NULL) error("Can't allocate memory"); if(currchron==NULL) error("Can't allocate memory"); +// Proposal vector for Z; +double *eqprob; +eqprob = (double *)calloc(*G, sizeof(double)); +for(k=0;k<*G;k++) eqprob[k] = (double)1 / *G; + // Set up MCMC variables I'll need -double *triupper,*tridiag,*trilower,vsqrtstar,vsqrtstarrat,temp,*z,tzVz,logdetratio,expnum,expden,igratio,logRv,*ones,*Vone,oneVone,exp1bit,exp2bit,*tempv,u1,u2,u3,sumv,sumxv,*tempc; -int accept,countchron=0,countv=0,countmu=0,countphi=0,currentpos,countK=0,countc=0; +double *tempv,*tempDM,*tempR,*tempDMstar,*tempRstar,*tempvstar,*triupper,*tridiag,*trilower,*logdetratio,vstar,vstarrat,logRv,logRz,*tempc,phi1star,phi2star,phi1starrat,phi2starrat,logRphi1,logRphi2,*B,*S,tBS; +tempv = (double *)calloc(*n-1, sizeof(double)); +tempvstar = (double *)calloc(*n-1, sizeof(double)); +tempDM = (double *)calloc(*n, sizeof(double)); +tempR = (double *)calloc(*n, sizeof(double)); +tempDMstar = (double *)calloc(*n, sizeof(double)); +tempRstar = (double *)calloc(*n, sizeof(double)); +logdetratio = (double *)calloc(*m, sizeof(double)); triupper = (double *)calloc(*n, sizeof(double)); tridiag = (double *)calloc(*n, sizeof(double)); trilower = (double *)calloc(*n, sizeof(double)); -tempv = (double *)calloc(*n-1, sizeof(double)); -z = (double *)calloc(*n, sizeof(double)); -ones = (double *)calloc(*n, sizeof(double)); -Vone = (double *)calloc(*n, sizeof(double)); tempc = (double *)calloc(*n, sizeof(double)); +B = (double *)calloc(*n, sizeof(double)); +S = (double *)calloc(*n, sizeof(double)); +int accept,countchron=0,countv=0,countz=0,currentpos,countc=0,countphi1=0,countphi2=0; // Start off the chronologies for(i=0;i<*n;i++) currchron[i] = allchrons[0][i]; diff(currchron,n,diffchron); +// Create the matrices I need later +for(j=0;j<*m;j++) { + // Create M,D and DM and MDM, etc + MDM[j] = 0; + MDMstar[j] = 0; + for(i=0;i<*n;i++) { + M[i][j] = MyMixMean[i][j][Z[i]]; + D[i] = MyMixPrec[i][Z[i]]; + DM[i][j] = D[i]*M[i][j]; + Mstar[i][j] = MyMixMean[i][j][Z[i]]; + DMstar[i][j] = D[i]*M[i][j]; + MDM[j] += M[i][j]*DM[i][j]; + MDMstar[j] += M[i][j]*DM[i][j]; + } +} + // Start iterations loop double progress=0; for(iter=0;iter<*iterations;iter++) { -//for(i=0;i<1;i++) { + + // Check if someone's pressed escape + R_CheckUserInterrupt(); + + // Report progress if(iter%*reportevery==0) { progress = (double) 100*iter/ *iterations; Rprintf("\r"); @@ -172,307 +215,217 @@ for(iter=0;iter<*iterations;iter++) { Rprintf("\r"); R_FlushConsole(); } - - // Check if someone's pressed escape - R_CheckUserInterrupt(); - - // Loop through climate dimensions; - for(j=0;j<*m;j++) { - // Create mujk,Djk and VDmu - for(i=0;i<*n;i++) { - mujk[i] = MyMixMean[i][j][K[i]]; - Djk[i] = MyMixPrec[i][K[i]]; - Dmu[i] = Djk[i]*mujk[i]; - //Rprintf("%lf %i \n",Djk[i],K[i]-1); - //Rprintf("%lf \n",Djk[i]); - } - - - for(l=0;l<*n-1;l++) tempv[l] = v[l][j]; - //for(l=0;l<*n;l++) Rprintf("tempv=%lf \n",tempv[l]); - //for(l=0;l<*n;l++) Rprintf("Djk=%lf \n",Djk[l]); - maketri(tempv,*n,Djk,triupper,tridiag,trilower); - trisolve (*n, trilower,tridiag, triupper, Dmu, VDmu); - - // Sample a new v - for(i=0;i<*n-1;i++) { - - vsqrtstar = truncatedwalk(sqrt(v[i][j]), *vmhsd, 0, LARGE); - vsqrtstarrat = truncatedrat(sqrt(v[i][j]), *vmhsd, 0, LARGE, vsqrtstar); - - for(l=0;l<*n;l++) z[l] = 0.0; - z[i] = 1.0; - z[i+1] = -1.0; - trisolve (*n, trilower,tridiag, triupper, z, Vz); - tzVz = Vz[i]-Vz[i+1]; - - logdetratio = -0.5*log(1+(1/pow(vsqrtstar,2)-1/v[i][j])*tzVz); - expnum = -pow(VDmu[i]-VDmu[i+1],2); - expden = 2*((1/(1/pow(vsqrtstar,2)-1/v[i][j]))+tzVz); - igratio = dlinvgauss2(pow(vsqrtstar,2),mu[j]*diffchron[i],phi[j]*diffchron[i])-dlinvgauss2(v[i][j],mu[j]*diffchron[i],phi[j]*diffchron[i]); - logRv = logdetratio + expnum/expden + igratio + 0.5*(log(v[i][j])-log(pow(vsqrtstar,2))); - - accept = (int)UpdateMCMC(logRv,1,0,vsqrtstarrat); - if(accept==1) { - v[i][j] = pow(vsqrtstar,2); - for(l=0;l<*n-1;l++) tempv[l] = v[l][j]; - maketri(tempv,*n,Djk,triupper,tridiag,trilower); - trisolve (*n, trilower,tridiag, triupper, Dmu, VDmu); - } - - - } // End of i loop for v update -/* - -if(i==2) { - Rprintf("i=%i \n",i+1); - Rprintf("v="); - for(l=0;l<*n-1;l++) Rprintf("%lf, ",tempv[l]); - Rprintf("\n"); - Rprintf("vstar=%lf \n",pow(vsqrtstar,2)); - Rprintf("vstarrat=%lf \n",vsqrtstarrat); - for(l=0;l<*n;l++) Rprintf("lower = %lf \n",trilower[l]); - for(l=0;l<*n;l++) Rprintf("diag = %lf \n",tridiag[l]); - for(l=0;l<*n;l++) Rprintf("upper = %lf \n",triupper[l]); - Rprintf("tzVz=%lf \n",tzVz); - Rprintf("logdetratio=%lf \n",logdetratio); - Rprintf("expnum=%lf \n",expnum); - Rprintf("expden=%lf \n",expden); - Rprintf("mu[j]=%lf \n",mu[j]); - Rprintf("phi[j]=%lf \n",phi[j]); - Rprintf("diffchron[i]=%lf \n",diffchron[i]); - Rprintf("ig1=%lf \n",dlinvgauss2(pow(vsqrtstar,2),mu[j]*diffchron[i],phi[j]*diffchron[i])); - Rprintf("ig2=%lf \n",dlinvgauss2(v[i][j],mu[j]*diffchron[i],phi[j]*diffchron[i])); - Rprintf("igratio=%lf \n",igratio); - Rprintf("logRv=%lf \n",logRv); - error("Stop! \n"); - } - - // Create some useful things for updating mu and phi - sumv=0; sumxv=0; - for(i=0;i<(*n-1);i++) { - sumv += v[i][j]; - sumxv += pow(diffchron[i],2)/v[i][j]; - } - - u1 = sumv; - u2 = diffchron[*n-1]-diffchron[0]-*phipar2; - u3 = sumxv; - - // Update mu - mustar = truncatedwalk(mu[j], *mumhsd, 0, LARGE); - mustarrat = truncatedrat(mu[j], *mumhsd, 0, LARGE, mustar); - logRmu = ((((double)*n-1)/2)+*mupar1-1)*(log(mustar)-log(mu[j]))-0.5*(phi[j]*u1*(1/mustar-1/mu[j]))-0.5*(phi[j]*u3+2*(*mupar2))*(mustar-mu[j]); - mu[j] = UpdateMCMC(logRmu,mustar,mu[j],mustarrat); - - // Update phi - //phi[j] = rgamma(((double)*n-1)/2+*phipar1,1/(u1/(2*mu[j])-u2+u3*mu[j]/2)); -*/ - } // End of j loop updating v mu and phi - - // Store everything - if((iter%*thinby==0) & (iter>=*burnin)) { + // Store stuff + if((iter%*thinby==0) & (iter>=*burnin)) { - currentpos = (int)((iter-*burnin)/(*thinby)); // Sample a chronology for(i=0;i<*n;i++) currchron[i] = allchrons[(int)currentpos][i]; diff(currchron,n,diffchron); - // Sample K - for(i=0;i<*n;i++) { - for(l=0;l<*G;l++) mixpro[l] = MyMixPro[i][l]; - K[i] = sample(mixseq, *G, mixpro); - } - + + // Store the used chronology for(i=0;i<*n;i++) chronstore[countchron+i] = currchron[i]; countchron += *n; - for(j=0;j<*m;j++) { - for(i=0;i<*n-1;i++) vstore[countv+i] = v[i][j]; - countv += *n-1; + // Store current z values + for(i=0;i<*n;i++) zstore[countz+i] = Z[i]; + countz += *n; + + // Store phi1 and phi2 + for(j=0;j<*m;j++) { + phi1store[countphi1+j] = phi1[j]; + phi2store[countphi2+j] = phi2[j]; + } + countphi1 += *m; + countphi2 += *m; + + for(j=0;j<*m;j++) { + + for(i=0;i<*n-1;i++) vstore[countv+i] = v[i][j]; + countv += *n-1; - // Calculate some climates - for(i=0;i<*n;i++) { - mujk[i] = MyMixMean[i][j][K[i]]; - Djk[i] = MyMixPrec[i][K[i]]; - Dmu[i] = Djk[i]*mujk[i]; - } + // Generate n independent random normals + for(l=0;l<*n;l++) rannormal[l] = rnorm(0,1); + + for(l=0;l<*n-1;l++) tempv[l] = v[l][j]; + for(l=0;l<*n;l++) tempDM[l] = DM[l][j]; + for(l=0;l<*n;l++) tempR[l] = R[l][j]; - for(l=0;l<*n-1;l++) tempv[l] = v[l][j]; - maketri(tempv,*n,Djk,triupper,tridiag,trilower); + // Create the tri-diagonal martix Q, stored in triupper, tridiag and trilower + maketri(tempv,*n,D,triupper,tridiag,trilower); // Create cholesky decomposition of the tridiag CholTriDiag(tridiag,triupper,*n,choldiag,chollower); - // Generate n independent random normals - for (l=0;l<*n;l++) rannormal[l] = rnorm(0,1); - // Solve the cholesky decomposition - trisolve (*n, chollower,choldiag, cholzeros, rannormal, mvnvar); + trisolve(*n, chollower,choldiag, cholzeros, rannormal, mvnvar); // Get the mean - trisolve (*n, trilower,tridiag, triupper, Dmu, VDmu); + trisolve(*n, trilower,tridiag, triupper, tempDM, tempR); + + // Output the climates + for(i=0;i<*n;i++) tempc[i] = tempR[i] + mvnvar[i]; - // Output the climates - for(i=0;i<*n;i++) tempc[i] = VDmu[i] + mvnvar[i]; - // Store the climates in the right place for(i=0;i<*n;i++) cstore[countc+i] = tempc[i]; countc += *n; - - } + // Need to be careful with c as now we're storing them in a weird order + + // End of j loop through climate dimensions + } + - // End of storage if statement + // End of storage loop } - -// End of iterations loop -} - -// Change the RNG state -PutRNGstate(); - -Rprintf("\r"); -R_FlushConsole(); -Rprintf("Completed: 100.00 %%"); -Rprintf("\n"); -R_FlushConsole(); -// End of function -} - -/* + // Loop through climate dimensions to update v + for(j=0;j<*m;j++) { -error("Stopped here!"); -Rprintf("muval[1]=%lf \n",muval[0]); -Rprintf("phival[1]=%lf \n",phival[0]); -Rprintf("vstart[1]=%lf \n",vstart[0]); -Rprintf("dlinvgauss2(vstart[1],mu[1],phi[1])=%lf \n",dlinvgauss2(vstart[0],muval[0],phival[0])); -error("Stop! \n"); -for(i=0;i<*n;i++) Rprintf("currchron[i]=%lf \n",currchron[i]); -for(i=0;i<*n-1;i++) Rprintf("diffchron[i]=%lf \n",diffchron[i]); - - -double *temppro,*tempphi,*temptheta,*diffthetas,sigmasqnew,alphanew,deltanew,alpharat; -int accept,countth=0,countsig=0,counta=0,countd=0,countchron=0; -diffthetas = (double *)calloc(*m-1, sizeof(double)); -temppro = (double *)calloc(*nummix, sizeof(double)); -tempphi = (double *)calloc(*m-1, sizeof(double)); -temptheta = (double *)calloc(*m, sizeof(double)); -if(diffthetas==NULL) error("Can't allocate memory"); -if(temppro==NULL) error("Can't allocate memory"); -*/ - + for(l=0;l<*n-1;l++) tempv[l] = v[l][j]; + for(l=0;l<*n;l++) tempDM[l] = DM[l][j]; + for(l=0;l<*n;l++) tempR[l] = R[l][j]; - // if(v[148][0]==100) { - // for(l=0;l<*n;l++) Rprintf("trilower=%lf tridiag=%lf triupper=%lf Dmu=%lf VDmu=%lf \n",trilower[l],tridiag[l],triupper[l],Dmu[l],VDmu[l]); - // error("BAD"); - //} - - /*if(iter=20) { - Rprintf("\n v[i][j] = \n"); - for(i=0;i<*n-1;i++) { - for(j=0;j<*m;j++) Rprintf("%lf ",v[i][j]); - Rprintf("\n"); - } - error("Got to here"); - }*/ - /*if(i==148 & j==0) { - Rprintf("\n v=%lf, vstar=%lf VDmu[i]%lf VDmu[i+1]=%lf \n",v[i][j],vstar,VDmu[i],VDmu[i+1]); - //error("Got to here"); - }*/ -/* -Rprintf("mu[j]=%lf, mustar=%lf, mumhsd=%lf, logRmu=%lf \n",mu[j],mustar,*mumhsd,logRmu); - error("Good?\n"); + // Create the tri-diagonal martix Q, stored in triupper, tridiag and trilower + maketri(tempv,*n,D,triupper,tridiag,trilower); + // Find R from solution given in appendix to paper + trisolve(*n,trilower,tridiag,triupper,tempDM,tempR); + // Finally create MDR + MDR[j] = 0; + for(l=0;l<*n;l++) MDR[j] += tempDM[l]*tempR[l]; + + // Sample a new v for(i=0;i<*n-1;i++) { - for(j=0;j<*m;j++) Rprintf("%lf ",v[i][j]); - Rprintf("\n"); - } - if(iter=500) error("V's ok?\n"); - for(i=0;i<(*n-1);i++) { - Rprintf("%i v=%lf \n",i,v[i][j]); - } - Rprintf("mu[j]=%lf, mustar=%lf, mumhsd=%lf, logRmu=%lf \n",mu[j],mustar,*mumhsd,logRmu); - Rprintf("j=%i, phi[j]=%lf \n",j,phi[j]); - Rprintf("j=%i, mu[j]=%lf, mustar=%lf, mumhsd=%lf, logRmu=%lf, bit1=%lf, bit2=%lf, bit3=%lf \n",j,mu[j],mustar,*mumhsd,logRmu,((((double)*n-1)/2)+*mupar1-1)*(log(mustar)-log(mu[j])),0.5*(phi[j]*u1*(1/mustar-1/mu[j])),0.5*(phi[j]*u3+2*(*mupar2))*(mustar-mu[j])); - Rprintf("phi=%lf, par1=%lf, par2=%lf \n",phi[j],((double)*n-1)/2+*phipar1,1/(u1/(2*mu[j])-u2+u3*mu[j]/2)); -Rprintf("Got to here \n"); - for(l=0;l<*n;l++) Rprintf("trilower=%lf tridiag=%lf triupper=%lf Dmu=%lf VDmu=%lf \n",trilower[l],tridiag[l],triupper[l],Dmu[l],VDmu[l]); - // for(l=0;l<*n;l++) Rprintf("tempv=%lf *n=%i Djk=%lf diffchron=%lf \n",tempv[l],*n,Djk[l],diffchron[l]); - for(i=0;i<*n-1;i++) { - for(l=0;l<*m;l++)Rprintf("v=%lf ",v[i][l]); - Rprintf("\n"); - } - -*/ -/* - // Update K + + // For some reason I found it easier to accept new values on the volatility scale + vstar = truncatedwalk(v[i][j], *vmhsd, 0, LARGE); + vstarrat = truncatedrat(v[i][j], *vmhsd, 0, LARGE, vstar); + for(l=0;l<*n-1;l++) tempvstar[l] = tempv[l]; + tempvstar[i] = vstar; + + // Find log ratio of determinant + for(l=0;l<*n;l++) B[l] = 0.0; + B[i] = 1.0; + B[i+1] = -1.0; + trisolve (*n, trilower,tridiag, triupper, B, S); + tBS = S[i]-S[i+1]; + logdetratio[j] = log(1+(1/vstar-1/v[i][j])*tBS); + + // Calculate Rstar + maketri(tempvstar,*n,D,triupper,tridiag,trilower); + trisolve(*n,trilower,tridiag,triupper,tempDM,tempRstar); - for(i=0;i<*n;i++) { - temp = sample(mixseq, *G, mixpro); - kstar = mixseq[(int)temp]; - for(l=0;l<*n;l++) Kstar[l] = K[l]; - Kstar[i]=kstar; + //Calculate MDRstar + MDRstar[j] = 0; + for(l=0;l<*n;l++) MDRstar[j] += tempDM[l]*tempRstar[l]; + logRv = -0.5*log(vstar/tempv[i]) -0.5*logdetratio[j] + 0.5*MDRstar[j] - 0.5*MDR[j] + dlinvgauss2(vstar,phi1[j]*diffchron[i],phi2[j]*diffchron[i]) - dlinvgauss2(tempv[i],phi1[j]*diffchron[i],phi2[j]*diffchron[i]); + + //Rprintf("logRv=%lf \n",logRv); + + //Rprintf("i=%i, j=%i, v=%lf, vstar=%lf, vstarrat=%lf logRv=%lf logdetratio[j]=%lf MDRstar[j]=%lf MDR[j]=%lf dligstar=%lf dlig=%lf\n",i,j,v[i][j],vstar,vstarrat,logRv,logdetratio[j],MDRstar[j],MDR[j],dlinvgauss2(vstar,phi1[j]*diffchron[i],phi2[j]*diffchron[i]),dlinvgauss2(tempv[i],phi1[j]*diffchron[i],phi2[j]*diffchron[i])); + + accept = (int)UpdateMCMC(logRv,1,0,vstarrat); + if(accept==1) { + v[i][j] = vstar; + tempv[i] = vstar; + maketri(tempv,*n,D,triupper,tridiag,trilower); + for(l=0;l<*n;l++) tempR[l] = tempRstar[l]; + for(l=0;l<*n;l++) R[l][j] = tempRstar[l]; + MDR[j] = MDRstar[j]; + } + + } // End of i loop for v update - logRk = 0; + // End of loop through climate dimensions + } + + + // Sample a new z + for(i=0;i<*n;i++) { + for(l=0;l<*n;l++) Zstar[l] = Z[l]; + Zstar[i] = sample(mixseq, *G, eqprob); + for(j=0;j<*m;j++) { - for(l=0;l<*n;l++) { - mujk[l] = MyMixMean[l][j][(int)K[l]]; - - Djk[l] = MyMixPrec[l][(int)K[l]]; - } - for(l=0;l<*n;l++) Dmu[l] = mujk[l]*Djk[l]; + for(l=0;l<*n;l++) Mstar[l][j] = MyMixMean[l][j][Zstar[l]]; + for(l=0;l<*n;l++) DMstar[l][j] = D[l]*Mstar[l][j]; + MDMstar[j] = MDM[j] - M[i][j]*DM[i][j] + Mstar[i][j]*DMstar[i][j]; + + // Create Rstar for(l=0;l<*n-1;l++) tempv[l] = v[l][j]; - - - maketri(tempv,*n,Djk,triupper,tridiag,trilower); - - trisolve (*n, trilower,tridiag, triupper, Dmu, VDmu); - - for(l=0;l<*n;l++) { - mujkstar[l] = MyMixMean[l][j][(int)Kstar[l]]; - - Djkstar[l] = MyMixPrec[l][(int)Kstar[l]]; + for(l=0;l<*n;l++) tempDM[l] = DMstar[l][j]; + for(l=0;l<*n;l++) tempR[l] = R[l][j]; + + maketri(tempv,*n,D,triupper,tridiag,trilower); + trisolve(*n,trilower,tridiag,triupper,tempDM,tempR); + + MDRstar[j] = 0; + for(l=0;l<*n;l++) MDRstar[j] += tempDM[l]*tempR[l]; + } + + logRz = log(MyMixPro[i][Zstar[i]]/MyMixPro[i][Z[i]]); + for(j=0;j<*m;j++) { + logRz += - 0.5*MDMstar[j] + 0.5*MDRstar[j] + 0.5*MDM[j] - 0.5*MDR[j]; + } + + accept = (int)UpdateMCMC(logRz,1,0,1); + if(accept==1) { + Z[i] = Zstar[i]; + for(j=0;j<*m;j++) { + for(l=0;l<*n;l++) M[l][j] = Mstar[l][j]; + for(l=0;l<*n;l++) DM[l][j] = DMstar[l][j]; + MDM[j] = MDMstar[j]; + MDR[j] = MDRstar[j]; + + for(l=0;l<*n-1;l++) tempv[l] = v[l][j]; + for(l=0;l<*n;l++) tempDM[l] = DMstar[l][j]; + for(l=0;l<*n;l++) tempR[l] = R[l][j]; + maketri(tempv,*n,D,triupper,tridiag,trilower); + trisolve(*n,trilower,tridiag,triupper,tempDM,tempR); + for(l=0;l<*n;l++) R[l][j] = tempR[l]; } - - for(l=0;l<*n;l++) Dmustar[l] = mujkstar[l]*Djkstar[l]; - for(l=0;l<*n-1;l++) tempv[l] = v[l][j]; - - maketri(tempv,*n,Djkstar,triupper,tridiag,trilower); - - - trisolve (*n, trilower,tridiag, triupper, Dmustar, VDmustar); - for(l=0;l<*n;l++) ones[l] = 0.0; + } + + // End of loop for z + } + + // Sample new phi1 and phi2 + for(j=0;j<*m;j++) { + phi1star = truncatedwalk(phi1[j], *phi1mhsd, 0, LARGE); + phi1starrat = truncatedrat(phi1[j], *phi1mhsd, 0, LARGE, phi1star); - ones[i] = 1.0; - trisolve (*n, trilower,tridiag, triupper, ones, Vone); + logRphi1 = dlnorm(phi1star,phi1dlmean[j],phi1dlsd[j],1) - dlnorm(phi1[j],phi1dlmean[j],phi1dlsd[j],1); + for(i=0;i<*n-1;i++) logRphi1 += dlinvgauss2(v[i][j],phi1star*diffchron[i],phi2[j]*diffchron[i]) - dlinvgauss2(v[i][j],phi1[j]*diffchron[i],phi2[j]*diffchron[i]); + //for(i=0;i<*n-1;i++) Rprintf("i=%i, v=%lf, first=%lf, second=%lf diffchron=%lf \n",i,v[i][j],dlinvgauss2(v[i][j],phi1star*diffchron[i],phi2[j]*diffchron[i]), dlinvgauss2(v[i][j],phi1[j]*diffchron[i],phi2[j]*diffchron[i]),diffchron[i]); + + //Rprintf("logRphi1=%lf j=%i, phi1star=%lf, phi1starrat=%lf, phi1[j]=%lf, dlnorm(phi1star,phi1dlmean[j],phi1dlsd[j],1)=%lf, dlnorm(phi1[j],phi1dlmean[j],phi1dlsd[j],1)=%lf, \n",logRphi1,j,phi1star,phi1starrat,phi1[j],dlnorm(phi1star,phi1dlmean[j],phi1dlsd[j],1),dlnorm(phi1[j],phi1dlmean[j],phi1dlsd[j],1)); + phi1[j] = UpdateMCMC(logRphi1,phi1star,phi1[j],phi1starrat); - oneVone = Vone[i]; - exp1bit = 0.0; - for(l=0;l<*n;l++) exp1bit = mujkstar[l]*Djkstar[l]*mujkstar[l]-mujk[l]*Djk[l]*mujk[l]; + phi2star = truncatedwalk(phi2[j], *phi2mhsd, 0, LARGE); + phi2starrat = truncatedrat(phi2[j], *phi2mhsd, 0, LARGE, phi2star); - exp2bit = 0.0; - for(l=0;l<*n;l++) exp2bit = mujk[l]*Djk[l]*VDmu[l]-mujkstar[l]*Djkstar[l]*VDmustar[l]; - logdetratio = -log(1+(Djkstar[i]-Djk[i])*oneVone); - - logRk += - 0.5*exp1bit-0.5*exp2bit+0.5*logdetratio + log(MyMixPro[i][(int)kstar])-log(MyMixPro[i][(int)K[i]])+0.5*log(Djkstar[i])-0.5*log(Djk[i]); + logRphi2 = dlnorm(phi2star,phi2dlmean[j],phi2dlsd[j],1) - dlnorm(phi2[j],phi2dlmean[j],phi2dlsd[j],1); + for(i=0;i<*n-1;i++) logRphi2 += dlinvgauss2(v[i][j],phi1[j]*diffchron[i],phi2star*diffchron[i]) - dlinvgauss2(v[i][j],phi1[j]*diffchron[i],phi2[j]*diffchron[i]); - - } // End of j loop updating k + //Rprintf("logRphi2=%lf \n",logRphi2); + phi2[j] = UpdateMCMC(logRphi2,phi2star,phi2[j],phi2starrat); - K[i] = UpdateMCMC(logRk,kstar,K[i],1.0); - - - } // End of i loop updating k - for(i=*n-10;i<*n;i++) { - for(l=*n-10;l<*n;l++) Rprintf("%4.2lf ",bigQ[i][l]); - //Rprintf("%4.2lf ",trilower[i]); - Rprintf("\n"); - } + // End of j loop updating phi1 and phi2 + } + + +// End of iterations loop +} -// for(l=0;l<*n;l++) Rprintf("diag=%lf, upp=%lf, choldiag=%lf, chollower=%lf \n",tridiag[l],triupper[l],choldiag[l],chollower[l]); - for(l=0;l<*n;l++) Rprintf("tempc = %lf \n",tempc[l]); - -*/ +// Change the RNG state +PutRNGstate(); + +Rprintf("\r"); +R_FlushConsole(); +Rprintf("Completed: 100.00 %%"); +Rprintf("\n"); +R_FlushConsole(); + +// End of function +} diff --git a/src/use.c b/src/use.c index 14ffb0a..6ed3415 100755 --- a/src/use.c +++ b/src/use.c @@ -366,6 +366,7 @@ free(B); return 2.0*sum; } + double dlognormal(double *x, double *mean, double **var,int len) { @@ -375,14 +376,12 @@ double dens,logdeter,**invvar,bigsum=0.0; int i,j; invvar = (double **)calloc(len, sizeof(double *)); -for (i=0;i