Fixed some bugs and changed the hyperparameter posterior sampler

no samples hyperparameters from joint distribution
eirikmn · Jul 5, 2023 · 6c4a284 · 6c4a284
1 parent 78e2018
commit 6c4a284
Show file tree

Hide file tree

Showing 192 changed files with 400,112 additions and 126 deletions.
diff --git a/.gitignore b/.gitignore
@@ -46,3 +46,4 @@ po/*~
 
 # DS_Store files
 **/.DS_Store
+inst/doc
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -38,9 +38,11 @@ Suggests:
     sn,
     sp (>= 1.4-5),
     rmarkdown,
-    graphics
+    graphics,
+    knitr
 Encoding: UTF-8
 LazyData: true
 URL: https://github.com/eirikmn/bremla
 BugReports: https://github.com/eirikmn/bremla/issues
 RoxygenNote: 7.2.3
+VignetteBuilder: knitr
diff --git a/R/bremla.R b/R/bremla.R
@@ -95,7 +95,7 @@
 #' depth = NGRIP_5cm$depth
 #' d18O = NGRIP_5cm$d18O
 #' proxy=d18O
-#' formula = dy~-1+depth2
+#' formula = dy~-1+depth2 + proxy
 #' depth2 = depth^2/depth[1]^2 #normalize for stability
 #'
 #' data = data.frame(age=age,dy=c(NA,diff(age)),depth=depth,depth2=depth2,proxy=proxy)

diff --git a/R/bremla_chronology_simulation.R b/R/bremla_chronology_simulation.R
@@ -86,45 +86,6 @@ bremla_chronology_simulation = function(object, control.sim,print.progress=FALSE
 
   time.start = Sys.time()
 
-  if(tolower(method) == "inla"){ #if INLA is used
-    reg.model = object$.internal$lat.selection
-
-    if(print.progress) cat("Simulating ",nsims, " hyperparameters from INLA posterior...",sep="")
-
-    if(tolower(noise) %in% c("rgeneric","custom")){
-      hypersamples = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
-      param.names = object$.args$control.fit$rgeneric$param.names
-      for(i in 1:length(object$.args$control.fit$rgeneric$from.theta)){
-        paramsamp = object$.args$control.fit$rgeneric$from.theta[[i]](hypersamples[,i])
-        if(is.null(param.names[i]) || is.na(param.names[i])){
-          tempname = paste0("hyperparameter",i)
-          object$simulation[[tempname]] = paramsamp
-        }else{
-          object$simulation[[param.names[i] ]] = paramsamp
-
-        }
-      }
-
-    }else if(tolower(noise) %in% c(0,"ar(0)","ar0","iid","independent")){
-      hypersamples = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
-      object$simulation = list(sigma = 1/sqrt(hypersamples[,1]))
-    }else if (tolower(noise) %in% c(1,"ar1","ar(1)")){
-      hypersamples = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
-      object$simulation = list(sigma = 1/sqrt(hypersamples[,1]), phi=hypersamples[,2])
-
-    }else if (tolower(noise) %in% c(2,"ar2","ar(2)")){
-      hypersamples = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
-      p=2
-      hypersamplesar2 = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
-      phii = hypersamplesar2[, 2L:(2L+(p-1L))]
-      phis = apply(phii, 1L, INLA::inla.ar.pacf2phi)
-      object$simulation = list(sigma = 1/sqrt(hypersamples[,1]),phi1=phis[1,],phi2=phis[2,])
-    }
-
-
-    if(print.progress) cat(" completed!\n",sep="")
-  }
-
 
   ## sample mean ("fixed") vector and ("stochastic") noise component
 
@@ -134,23 +95,22 @@ bremla_chronology_simulation = function(object, control.sim,print.progress=FALSE
 
     ##sample fixed parameters first
     latentselection = object$.internal$lat.selection
-    # latentselection = list()
-    # if(reg.model$const) latentselection$`(Intercept)`=1
-    # if(reg.model$depth1) latentselection$z=1
-    # if(reg.model$depth2) latentselection$z2=1
-    # if(reg.model$proxy) latentselection$x = 1
-
-    # for(i in 2:object$.args$nevents){
-    #   if(reg.model$psi0) latentselection[[paste0("a",i-1)]] = 1
-    #   if(reg.model$psi1)latentselection[[paste0("c",i-1)]] = 1
-    # }
-    #latentsamples = inla.posterior.sample(nsims,object$fitting$fit,selection=latentselection,verbose=FALSE,add.names=FALSE)
+    reg.model = object$.internal$lat.selection
+
     ncores_postsamp = max(1,object$.args$control.sim$ncores)
+
+
+
     latentsamples = INLA::inla.posterior.sample(nsims,object$fitting$inla$fit,
                                           selection=latentselection,verbose=FALSE,
                                           add.names=FALSE,num.threads = ncores_postsamp)
 
+    int_hyper = data.frame(matrix(NA, nrow=nsims,ncol=length(latentsamples[[1]]$hyperpar))) #get hyperparameters from joint distribution
+
+    colnames(int_hyper) = names(latentsamples[[1]]$hyperpar)
+
     n=nrow(object$data)
+
     if(object$.args$control.sim$store.everything) object$simulation$dmean = matrix(NA,nrow=n,ncol=nsims)
     time.endmean = Sys.time()
     if(print.progress) cat(" completed in ",difftime(time.endmean,time.startmean,units="secs")[[1]]," seconds!\n",sep="")
@@ -162,12 +122,15 @@ bremla_chronology_simulation = function(object, control.sim,print.progress=FALSE
 
 
 
+
+
     for(i in 1:nsims){
       if(i %% 1000 == 0 && print.progress){
         cat("Age simulation ",i,"/",nsims,". Elapsed time: ",difftime(Sys.time(),time.startage,units="secs")[[1]]," seconds...","\n",sep="")
       }
       ## from fixed parameters compute fixed model component using 'meanmaker' function
       coefs = latentsamples[[i]]$latent
+      int_hyper = latentsamples[[i]]$hyperpar
 
       dmeansim = meanmaker( coefs, reg.model, data = object$data )
 
@@ -177,27 +140,67 @@ bremla_chronology_simulation = function(object, control.sim,print.progress=FALSE
 
       ##sample noise component
       if(tolower(noise) %in% c("rgeneric","custom")){
-        theta = hypersamples[i,]
+
+        theta = int_hyper
+        param.names = object$.args$control.fit$rgeneric$param.names
+        hypersamples = int_hyper
+        for(i in 1:length(object$.args$control.fit$rgeneric$from.theta)){
+          paramsamp = object$.args$control.fit$rgeneric$from.theta[[i]](hypersamples[i])
+          if(is.null(param.names[i]) || is.na(param.names[i])){
+            tempname = paste0("hyperparameter",i)
+            object$simulation$params[[tempname]] = paramsamp
+          }else{
+            object$simulation$params[[param.names[i] ]] = paramsamp
+          }
+        }
+
 
         Q = object$.args$control.fit$rgeneric$Q(theta,n,ntheta=length(theta))
         muvek = numeric(n)
         noisesim = Qsimmer(1, Q, muvek)
+        #hypersamples = int_hyper
+      }else{
+        hypersamples = int_hyper
+
+        if(tolower(noise) %in% c(0,"iid","independent","ar0","ar(0)")){
+          object$simulation$params$sigma[i] = 1/sqrt(hypersamples[1])
+          #object$simulation$params = list(sigma = 1/sqrt(hypersamples[,1]))
+
+          noisesim = rnorm(n,mean=0,sd=object$simulation$params$sigma[i])
+
+
+        }else if(tolower(noise) %in% c(1,"ar1","ar(1)")){
+          hypersamples = int_hyper
+          object$simulation$params$sigma[i] = 1/sqrt(hypersamples[1])
+          object$simulation$params$phi[i] = hypersamples[2]
+
+          noisesim = arima.sim(n=n,list(ar=c(object$simulation$params$phi[i])),
+                               sd = object$simulation$params$sigma[i]*sqrt(1-object$simulation$params$phi[i]^2))
+
+
+        }else if(tolower(noise) %in% c(2,"ar2","ar(2)")){
+          p=2
+          hypersamplesar2 = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
+          phii = hypersamplesar2[, 2L:(2L+(p-1L))]
+          phis = apply(phii, 1L, INLA::inla.ar.pacf2phi)
+          object$simulation$params$sigma[i] = 1/sqrt(hypersamples[1])
+          object$simulation$params$phi1[i] = phis[1]
+          object$simulation$params$phi2[i] = phis[2]
 
-      }else if(tolower(noise) %in% c(0,"iid","independent","ar0","ar(0)")){
-        noisesim = rnorm(n,mean=0,sd=object$simulation$sigma[i])
 
-      }else if(tolower(noise) %in% c(1,"ar1","ar(1)")){
-        noisesim = arima.sim(n=n,list(ar=c(object$simulation$phi[i])),
-                             sd = object$simulation$sigma[i]*sqrt(1-object$simulation$phi[i]^2))
+          gamma0 = (1-phis[2,i])/((1+phis[2,i])*(1-phis[1,i]-phis[2,i])*(1+phis[1,i]-phis[2,i]))
 
-      }else if(tolower(noise) %in% c(2,"ar2","ar(2)")){
-        gamma0 = (1-phis[2,i])/((1+phis[2,i])*(1-phis[1,i]-phis[2,i])*(1+phis[1,i]-phis[2,i]))
+          noisesim = arima.sim(n = n, list(ar = c( phis[1,i],phis[2,i])),
+                               sd = object$simulation$params$sigma[i]*sqrt(1/gamma0))
+
+
+          #object$simulation$params = list(sigma = 1/sqrt(hypersamples[,1]),phi1=phis[1,],phi2=phis[2,])
+        }
 
-        noisesim = arima.sim(n = n, list(ar = c( phis[1,i],phis[2,i])),
-                             sd = object$simulation$sigma[i]*sqrt(1/gamma0))
       }
 
 
+
       ## Take cumulatives. If log transformation is used, transform back first
       if(object$.args$control.fit$transform == "log"){
         object$simulation$age[,i] = object$initials$age + cumsum(exp(dmeansim+noisesim))
@@ -206,6 +209,30 @@ bremla_chronology_simulation = function(object, control.sim,print.progress=FALSE
       }
 
 
+
+    #if(print.progress) cat("Simulating ",nsims, " hyperparameters from INLA posterior...",sep="")
+
+    # if(tolower(noise) %in% c("rgeneric","custom")){
+      #hypersamples = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
+
+
+
+
+    # }else if(tolower(noise) %in% c(0,"ar(0)","ar0","iid","independent")){
+    #   hypersamples = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
+    #   object$simulation = list(sigma = 1/sqrt(hypersamples[,1]))
+    # }else if (tolower(noise) %in% c(1,"ar1","ar(1)")){
+    #   hypersamples = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
+    #   object$simulation = list(sigma = 1/sqrt(hypersamples[,1]), phi=hypersamples[,2])
+    #
+    # }else if (tolower(noise) %in% c(2,"ar2","ar(2)")){
+    #   hypersamples = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
+    #   p=2
+    #   hypersamplesar2 = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
+    #   phii = hypersamplesar2[, 2L:(2L+(p-1L))]
+    #   phis = apply(phii, 1L, INLA::inla.ar.pacf2phi)
+    #   object$simulation = list(sigma = 1/sqrt(hypersamples[,1]),phi1=phis[1,],phi2=phis[2,])
+    # }
     }
   }
 

diff --git a/R/bremla_synchronized_simulation.R b/R/bremla_synchronized_simulation.R
@@ -125,83 +125,66 @@ time.start = Sys.time()
   free_n = n-m
   free_indexes = (1:n)[-tie_indexes]
 
-  if(tolower(object$.args$control.fit$method) == "inla" && is.null(object$simulation)){ #if INLA is used
+  if(tolower(object$.args$control.fit$method) == "inla"){ # && is.null(object$simulation)){ #if INLA is used
+
     reg.model = object$.internal$lat.selection
+    latentselection = object$.internal$lat.selection
 
-    if(print.progress) cat("Simulating ",nsims, " hyperparameters from INLA posterior...",sep="")
 
-    if(tolower(noise) %in% c("rgeneric","custom")){
-      hypersamples = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
-      param.names = object$.args$control.fit$rgeneric$param.names
-      for(i in 1:length(object$.args$control.fit$rgeneric$from.theta)){
-        paramsamp = object$.args$control.fit$rgeneric$from.theta[[i]](hypersamples[,i])
-        if(is.null(param.names[i]) || is.na(param.names[i])){
-          tempname = paste0("hyperparameter",i)
-          object$simulation[[tempname]] = paramsamp
-        }else{
-          object$simulation[[param.names[i] ]] = paramsamp
+    timecoef.start = Sys.time()
+    ncores_postsamp = max(1,object$.args$control.sim$ncores)
+    latentsamples = INLA::inla.posterior.sample(nsims,object$fitting$inla$fit,
+                                                selection=latentselection,verbose=FALSE,
+                                                add.names=FALSE,num.threads = ncores_postsamp)
+    #int_hyper = latentsamples[[i]]
+
+    #int_hyper = data.frame(matrix(NA, nrow=nsims,ncol=length(latentsamples[[1]]$hyperpar))) #get hyperparameters from joint distribution
+
+    #colnames(int_hyper) = names(latentsamples[[1]]$hyperpar)
+    timecoef.end = Sys.time()
+    if(print.progress) cat(" completed in ",difftime(timecoef.end,timecoef.start,units="secs")[[1]]," seconds.\n",sep="")
+    samples = matrix(NA,nrow=n,ncol=nsims)
 
-        }
-      }
-    }else if(tolower(noise) %in% c(0,"ar(0)","ar0","iid","independent")){
-      hypersamples = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
-      object$simulation = list(sigma = 1/sqrt(hypersamples[,1]))
-    }else if (tolower(noise) %in% c(1,"ar1","ar(1)")){
-      hypersamples = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
-      object$simulation = list(sigma = 1/sqrt(hypersamples[,1]), phi=hypersamples[,2])
-
-    }else if (tolower(noise) %in% c(2,"ar2","ar(2)")){
-      hypersamples = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
-      p=2
-      hypersamplesar2 = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
-      phii = hypersamplesar2[, 2L:(2L+(p-1L))]
-      phis = apply(phii, 1L, INLA::inla.ar.pacf2phi)
-      object$simulation = list(sigma = 1/sqrt(hypersamples[,1]),phi1=phis[1,],phi2=phis[2,])
-    }
 
 
-    if(print.progress) cat(" completed!\n",sep="")
-  }
 
   if(print.progress) cat("Sampling fixed coefficients...",sep="")
-  latentselection = object$.internal$lat.selection
-  # latentselection = list()
-  # reg.model = object$.args$reg.model
-  # if(reg.model$const) latentselection$`(Intercept)`=1
-  # if(reg.model$depth1) latentselection$z=1
-  # if(reg.model$depth2) latentselection$z2=1
-  # if(reg.model$proxy) latentselection$x = 1
-  #
-  # for(i in 2:object$.args$nevents){
-  #   if(reg.model$psi0) latentselection[[paste0("a",i-1)]] = 1
-  #   if(reg.model$psi1)latentselection[[paste0("c",i-1)]] = 1
-  # }
-  timecoef.start = Sys.time()
-  ncores_postsamp = max(1,object$.args$control.sim$ncores)
-  latentsamples = INLA::inla.posterior.sample(nsims,object$fitting$inla$fit,
-                                        selection=latentselection,verbose=FALSE,
-                                        add.names=FALSE,num.threads = ncores_postsamp)
-
-  timecoef.end = Sys.time()
-  if(print.progress) cat(" completed in ",difftime(timecoef.end,timecoef.start,units="secs")[[1]]," seconds.\n",sep="")
-  samples = matrix(NA,nrow=n,ncol=nsims)
-
-
-  if(print.progress) cat("Simulating synchronized chronologies...\n",sep="")
+
+
+  if(print.progress) cat("\nSimulating synchronized chronologies...\n",sep="")
   timeage.start = Sys.time()
   for(r in 1:nsims){
+    hypersamples = latentsamples[[r]]$hyperpar
+    int_hyper = hypersamples
 
     if(tolower(noise) %in% c("rgeneric","custom")){
 
-      theta = hypersamples[r,]
+      param.names = object$.args$control.fit$rgeneric$param.names
+      for(i in 1:length(object$.args$control.fit$rgeneric$from.theta)){
+        paramsamp = object$.args$control.fit$rgeneric$from.theta[[i]](hypersamples[i])
+        if(is.null(param.names[i]) || is.na(param.names[i])){
+          tempname = paste0("hyperparameter",i)
+          object$simulation$params[[tempname]][i] = paramsamp
+        }else{
+          object$simulation$params[[param.names[i] ]][i] = paramsamp
+        }
+      }
+
+      theta = hypersamples
 
       Qx = object$.args$control.fit$rgeneric$Q(theta,n,ntheta=length(theta))
 
       Qfull = Qymaker(Qx)
 
-    }else if(tolower(noise) %in% c(1,"ar1","ar(1)")){
-      sigma_sample = object$simulation$sigma[r]
-      phi_sample = object$simulation$phi[r]
+    }else if( tolower(noise) %in% c(1,"ar1","ar(1)") ){
+
+      hypersamples = int_hyper
+      #hypersamples = INLA::inla.hyperpar.sample(nsims,object$fitting$inla$fit)
+      object$simulation$params$sigma[r] = 1/sqrt(hypersamples[1])
+      object$simulation$params$phi[r] = hypersamples[2]
+
+      sigma_sample = object$simulation$params$sigma[r]
+      phi_sample = object$simulation$params$phi[r]
 
       if(object$.args$control.fit$noise=="ar1"){
         Qfull = Qmaker_ar1cum(n,sigma_sample,phi_sample)
@@ -214,6 +197,7 @@ time.start = Sys.time()
     }else{
       stop("Currently, only the 'ar1' noise model is implemented. Other models can be specified via the 'rgeneric' model, see '?rgeneric.fitting' for details.")
     }
+
     Qa = Qfull[-tie_indexes,-tie_indexes]
     Qab = Qfull[-tie_indexes,tie_indexes]
     if(m==1){
@@ -244,6 +228,7 @@ time.start = Sys.time()
       cat("Synchronous age simulation ",r,"/",nsims,". Elapsed time: ",difftime(Sys.time(),timeage.start,units="secs")[[1]]," seconds...","\n",sep="")
     }
   }
+  }
 
   timeage.end = Sys.time()
 

diff --git a/R/default_arguments.R b/R/default_arguments.R
@@ -87,9 +87,8 @@ control.fit.default <- function(){
   return(list(
     noise="ar1",
     method="inla",
-    hyperprior = NULL,
     verbose=FALSE,
-    log.theta.prior=NULL,
+    hyperprior = NULL,
     improve.fixed=FALSE,
     ncores=1,
     transform="identity",

diff --git a/R/helpfulfunctions.R b/R/helpfulfunctions.R
@@ -552,10 +552,11 @@ tiepointsimmer = function(object, synchronization,print.progress=FALSE,...){
       samples = adolphi_tiepoint_simmer(nsims=synchronization$nsims,
                                         tieshifts=adolphilocs,...)
       samples = samples[,!is.na(tieinclude)] #only use tie-points where location is not NA
-      adolphilocs = adolphilocs[!is.na(tieinclude)]
+      #adolphilocs = adolphilocs[!is.na(tieinclude)]
+      adolphilocs = adolphilocs[tieinclude]
       synchronization$locations_unit="age"
       locations_indexes = which.index(adolphilocs,object$data$age)
-      samples = samples[,!is.na(locations_indexes)]
+      samples = samples[,tieinclude]
       locations_indexes = locations_indexes[!is.na(locations_indexes)]
       synchronization$locations = object$data$age[locations_indexes]
       #synchronization$locations = synchronization$locations[!is.na(locations_indexes)]

diff --git a/inst/reproduce_results/corricdata.xlsx b/inst/reproduce_results/corricdata.xlsx
diff --git a/inst/reproduce_results/corricktable.pdf b/inst/reproduce_results/corricktable.pdf