Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

14C calibration #2

Closed
MartinHinz opened this issue May 20, 2016 · 6 comments
Closed

14C calibration #2

MartinHinz opened this issue May 20, 2016 · 6 comments

Comments

@MartinHinz
Copy link
Member

Should 14C calibration be included at all?

If so, should we depend on eg. bchron, or shall we code a calibration by ourselves?

If so, including Bayesian Calibration? Maybe as future milestone?

@nevrome
Copy link
Member

nevrome commented May 20, 2016

Bchron is already on a very good way - nevertheless 14C-calibration and modelling is crucial. Maybe we can set up something in the future.

@nevrome
Copy link
Member

nevrome commented May 31, 2016

I realized that I need a very fast tool to do the calibration for neolithicR. My current solution with Bchron is much to slow.

The last time we talked about this, you described a pretty simple algorithm to perform an equally simple calibration. Do you think we could implement something like this as a Rcpp script quickly?

@MartinHinz
Copy link
Member Author

Sure. I am not sure if Rcpp is really necessary?

A short r-markdown piece below:

Quick Calibration

remark: step one and two are only necessary once, the calibration matrix could be stored

First: Get the cal curve, eg. from the internet (http://www.radiocarbon.org/IntCal13%20files/intcal13.14c)

con <- url("http://www.radiocarbon.org/IntCal13%20files/intcal13.14c")
calcurve<-read.csv(con,skip=11, header=F)
colnames(calcurve)<-c("CAL BP", "14C age","Error","Delta 14C","Sigma")

Second: make a matrix out of it

cal_bp <- calcurve[,1]
cal_matrix <- apply(calcurve,1,function(x){dnorm(cal_bp,mean = x[2], sd = x[3])})

Third: get your 14C date as probability vector

my_date <- c(4000,20)
my_date_prob<-dnorm(cal_bp,mean = my_date[1], sd = my_date[2])

Fourth: Multiply Vector by Matrix

res <- my_date_prob %*% cal_matrix

Fifth: Get the result in more meaningful shape

res.df <- cbind(cal_bp,res[1,])
res.df <- res.df[res.df[,2]>1e-05,]

Sixth: Plot

plot(res.df,type="l")

Seventh: Compare with Bcron:

library(Bchron)

res.bchron<-Bchron::BchronCalibrate(my_date[1],my_date[2],calCurves = "intcal13")
par(mfrow=c(2,1))
plot(res.df,type="l")
plot(res.bchron)
par(mfrow=c(1,1))

@MartinHinz
Copy link
Member Author

quickcal.pdf

@nevrome
Copy link
Member

nevrome commented May 31, 2016

Awesome - thanks! We should really consider to create an own calibration module. The simplicity of this code is impressive and list calibration is probably still the most commonly required task.

@MartinHinz
Copy link
Member Author

In the light of rcarbon and oxcaar, I think there is currently no need for an additional calibration function. So I will close this issue right now.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants