andrewcparnell / Bchron

A Bayesian radiocarbon chronology model with R
http://andrewcparnell.github.io/Bchron/
36 stars 10 forks source link

extracting means and SD for calibrated dates? #4

Closed cmbarton closed 8 years ago

cmbarton commented 8 years ago

It is easy to get calibrated ranges from BchronCalibrate. But is there some way to extract both the mean and 1 or 2 SD for calibrated dates? It looks like the sum(agegrid * densities) = mean. But what about 1 or 2 SD. The ranges are the most accurate representation of course. But it would nice to be able to statistically summarize the density curve for some analyses.

Thanks Michaesl

andrewcparnell commented 8 years ago

Hi Michael,

Here are two different ways of doing it:

library(Bchron)

ages1 = BchronCalibrate(ages=11553,
                        ageSds=230,
                        calCurves='intcal13',
                        ids='Date-1')

# Get the mean and sd the maths-y way

# Mean
with(ages1[[1]], sum(ageGrid*densities))
# sd
with(ages1[[1]], sqrt(sum(densities*(ageGrid-mean)^2)))

# Do it the simulation way
samples = with(ages1[[1]], sample(ageGrid, size = 1e5, prob=densities, replace=TRUE))
mean(samples)
sd(samples)
cmbarton commented 8 years ago

Thanks much!

Michael Barton School of Human Evolution &Social Change Center for Social Dynamics & Complexity Arizona State University

...Sent from my iPad

On May 1, 2016, at 1:36 PM, Andrew Parnell notifications@github.com wrote:

Hi Michael,

Here are two different ways of doing it:

library(Bchron)

ages1 = BchronCalibrate(ages=11553, ageSds=230, calCurves='intcal13', ids='Date-1')

Get the mean and sd the maths-y way

Mean

with(ages1[[1]], sum(ageGrid*densities))

sd

with(ages1[[1]], sqrt(sum(densities*(ageGrid-mean)^2)))

Do it the simulation way

samples = with(ages1[[1]], sample(ageGrid, size = 1e5, prob=densities, replace=TRUE)) mean(samples) sd(samples) — You are receiving this because you authored the thread. Reply to this email directly or view it on GitHub

cmbarton commented 8 years ago

Looking at this again, my question is a little bit different from perhaps what it seemed. Using an example from the vignette:

ages2 = BchronCalibrate(ages=c(3445,11553,7456), ageSds=c(50,230,110), calCurves=c('intcal13','intcal13','shcal13'))

Is there a way to get the calibrated mean and SD (and 2SD) for each individual date, in addition to the calibrated ranges at CI=.5, .95, and .97 that are automatically created in the summary? Would it be creating some kind of apply function based on your response above? Thanks.

andrewcparnell commented 8 years ago

I've also created a new function called SampleAges to do this a bit more easily. See the example in the new vignette. Simply replace the quantile function with mean or sd if that's what you want. However, I wouldn't recommend mean +/- 2 sd for calibrated dates. You're better off either with a credible interval using quantile, or even better the highest density region as Bchron provides by default.

Andrew

cmbarton commented 8 years ago

Thanks. I'll test both of the new features soon.

Michael Barton

On Jun 5, 2016, at 10:47 AM, Andrew Parnell notifications@github.com wrote:

I've also created a new function called SampleAges to do this a bit more easily. See the example in the new vignette. Simply replace the quantile function with mean or sd if that's what you want. However, I wouldn't recommend mean +/- 2 sd for calibrated dates. You're better off either with a credible interval using quantile, or even better the highest density region as Bchron provides by default.

Andrew

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/andrewcparnell/Bchron/issues/4#issuecomment-223826728, or mute the thread https://github.com/notifications/unsubscribe/AIjVr-6h2ZGqoRZIJDr_o6-sjuXh_cvPks5qIwuUgaJpZM4IQZ5y.

cmbarton commented 8 years ago

Unfortunately, I can't test until you release it to CRAN as a binary. Your source code requires R 3.3 and I have not yet upgraded.

andrewcparnell commented 8 years ago

No worries. Just been told it will be on CRAN in the next 24 hours.

cmbarton commented 8 years ago

wonderful.

I'll check back. I haven't yet updated because IMHO, it takes 3-6 months for the packages I use to be properly updated too. So I usually just wait for the xxx.1 version.

Michael

On Jun 6, 2016, at 5:13 AM, Andrew Parnell notifications@github.com wrote:

No worries. Just been told it will be on CRAN in the next 24 hours.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/andrewcparnell/Bchron/issues/4#issuecomment-223941584, or mute the thread https://github.com/notifications/unsubscribe/AIjVr9eiFwAX5PXzzarJ4NSst3XIcK24ks5qJA7VgaJpZM4IQZ5y.

cmbarton commented 8 years ago

Hi Andrew,

I can't see where BChron can do this, but you keep showing me new functionality so I thought I'd ask.

Sometimes I have a situation where there are multiple C14 dates for a single level. Before running BChronDensity on a bunch of levels, I'd like to combine the multiple dates for a level into a single one.

Layer 1: date 1 Layer 2: date 2, date 3 Layer 3: date 4

I need to combine dates 2 and 3 into a single, uncalibrated date before running BChronDensity.

Is there a way to do this with BChron?

Thanks Michael


C. Michael Barton Director, Center for Social Dynamics & Complexity Professor of Anthropology, School of Human Evolution & Social Change Head, Graduate Faculty in Complex Adaptive Systems Science Arizona State University

voice: 480-965-6262 (SHESC), 480-965-8130/727-9746 (CSDC) fax: 480-965-7671 (SHESC), 480-727-0709 (CSDC) www: http://www.public.asu.edu/~cmbarton, http://csdc.asu.edu

andrewcparnell commented 8 years ago

Hi Michael,

The classic way would be to combine the dates (x1 +/- s1), (x2 +/- s2) as x3 +/- s3 where: x3 = (x1 + x2)/2 s3 = sqrt(s1^2 + s2^2)/2 You could do this by hand and then shove into BchronDensity.

A better way would be to fit a model where both dates contribute to a single calendar date. However, this isn't possible in Bchron (yet).

Andrew

cmbarton commented 8 years ago

Thanks. I didn't think I could get away with a simple average.

Michael Barton

On Oct 25, 2016, at 2:44 PM, Andrew Parnell notifications@github.com wrote:

Hi Michael,

The classic way would be to combine the dates (x1 +/- s1), (x2 +/- s2) as x3 +/- s3 where: x3 = (x1 + x2)/2 s3 = sqrt(s1^2 + s2^2)/2 You could do this by hand and then shove into BchronDensity.

A better way would be to fit a model where both dates contribute to a single calendar date. However, this isn't possible in Bchron (yet).

Andrew

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/andrewcparnell/Bchron/issues/4#issuecomment-256185119, or mute the thread https://github.com/notifications/unsubscribe-auth/AIjVr88IlekFUig8DCbhzPxZ6AIdsUOBks5q3ngygaJpZM4IQZ5y.