bremla_biased_chronologies: Linear ramp fit

In eirikmn/bremla: Bayesian Regression Modeling of Layer-Counted Archives

Linear ramp fit

Description

Fits the linear ramp described by Erhardt et al. (2019) to proxy values using INLA.

Usage

bremla_biased_chronologies(object, control.bias, print.progress = FALSE)

Arguments

object	List object which is the output of function bremla_chronology_simulation
control.bias	List containing specifications for how the unknown random bias should be implemented. See control.bias.default for details.
print.progress	Boolean. If TRUE progress will be printed to the screen.

Value

Returns the same object list from the input, but appends another list of the summary statistics for each analysis (and samples if store.samples=TRUE), and inputs (settings). These are stored in object\$biases.

Author(s)

Eirik Myrvoll-Nilsen, [e-mail]

See Also

bremla_chronology_simulation

Examples

if(inlaloader()){
require(stats)
set.seed(1)
n <- 1000
phi <- 0.8
sigma <- 1.2
a_lintrend <- 0.3; a_proxy = 0.8
dy_noise <- as.numeric(arima.sim(model=list(ar=c(phi)),n=n,sd=sqrt(1-phi^2)))
lintrend <- seq(from=10,to=15,length.out=n)

proxy <- as.numeric(arima.sim(model=list(ar=c(0.9)),n=n,sd=sqrt(1-0.9^2)))
dy <- a_lintrend*lintrend + a_proxy*proxy + sigma*dy_noise

y0 = 11700;z0=1200
age = y0+cumsum(dy)
depth = 1200 + 1:n*0.05
depth2 = depth^2/depth[1]^2 #normalize for stability


formula = dy~-1+depth2 + proxy
data = data.frame(age=age,dy=dy,proxy=proxy,depth=depth,depth2=depth2)
data = rbind(c(y0,NA,NA,z0,NA),data) #First row is only used to extract y0 and z0.

events=list(locations=c(1210,1220,1240))
control.fit = list(ncores=2,noise="ar1")
control.sim=list(synchronized=2,
                 summary=list(compute=TRUE))

object = bremla_prepare(formula,data,nsims=5000,reference.label="simulated timescale",
                        events = events,
                        control.fit=control.fit,
                        control.sim=control.sim)
object = bremla_modelfitter(object)
object = bremla_chronology_simulation(object)
object = bremla_biased_chronologies(object, control.bias = list(bias.model="uniform",nsims=5000),
             print.progress=TRUE)
summary(object)
plot(object)
}

eirikmn/bremla documentation built on Jan. 25, 2025, 4:41 a.m.