linrampfitter: 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

linrampfitter(object, control.linramp, print.progress = FALSE)

Arguments

object	List object which is the output of function bremla_chronology_simulation
control.linramp	Includes the data and specifications for the linear ramp model fitting procedure. Must include control.linramp\$proxy and control.linramp\$interval (preferably also control.linramp\$interval.unit), the rest can be set to default values. See control.linramp.default for details.
print.progress	Boolean. If TRUE progress will be printed to screen.

Value

Returns the same object list from the input, but appends another list of results from the linear ramp fit, including posterior marginal distributions of the hyperparameters, summary statistics, and simulations of linear ramp and transition end point (if any).

Author(s)

Eirik Myrvoll-Nilsen, [e-mail]

References

Erhardt et al. (2019)

See Also

bremla_chronology_simulation,events_depth_to_age

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")
synchronization=list(locations=depth[c(100,400,700)],method="gauss",
                           params=list(mean=c(age[c(100,400,700)]+c(30,-100,50)),
                                       sd=c(50,20,100)
                                       )
                       )
control.sim=list(synchronized=TRUE,
                 summary=list(compute=TRUE))

#simulate transition:
prox = rnorm(n,mean=c(rep(0,400),seq(0,4,length.out=200),rep(4,400)),sd=1)
window = 330:500
control.linramp = list(label="Simulated",proxy=prox,interval=window,interval.unit="index",
    depth.ref=depth[401])
object = bremla_prepare(formula,data,nsims=5000,reference.label="simulated timescale",
                        events = events,
                        synchronization=synchronization,
                        control.fit=control.fit,
                        control.sim=control.sim,
                        control.linramp=control.linramp)
object = bremla_modelfitter(object)
object = tiepointsimmer(object)
object = bremla_synchronized_simulation(object)
object = linrampfitter(object,print.progress=TRUE)
summary(object)
plot(object)
}

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