osl_calibration

osl_calibration.r

# data ingest ---------------------------------------

grain_1 <- read.csv("CMC1.csv", header = FALSE)
grain_2 <- read.csv("CMC2.csv", header = FALSE)
grain_3 <- read.csv("CMC3.csv", header = FALSE)

grains <- list(grain_1 = grain_1, 
               grain_2 = grain_2, 
               grain_3 = grain_3)

dose_rate <- read.csv("dr.csv", header = FALSE)$V1
random_error <- read.csv("drs.csv", header = FALSE)$V1

# data structure for stan model  ----------------------

N <-  length(grains) # number of grains we are analysing 
J <-  unname(sapply(grains, nrow)) # number of analyses per grain 
y <-  stack(sapply(grains, `[[`, 1))$values # stack of 1st col of each list item
sigma <- stack(sapply(grains, `[[`, 2))$values # stack of 2nd col of each list item
dr <- dose_rate
drs <- random_error

# starting values for the parameters ------------------

mu <- unname(sapply(grains, function(i) median(i$V1)))
theta <- rnorm(sum(J)) * (0.03 * y) + y
theta <- ifelse(theta >= 0.01, theta, 0.01)
tau1 <- 0.20
tau2 <- 0.50
lamda <- rnorm(sum(J)) * 0.05 + 0.5
delta <- rnorm(N) * drs + dr

# stan... ------------------

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
model <- stan_demo()

# test to see if my objects are similar to his
df <- list(data.frame(mu = mu,
                      al_mu = al_mu),
           data.frame(theta = theta,
                      al_theta = al_theta),
           data.frame(lamda = lamda,
                      al_lamda = al_lamda),
           data.frame(delta = delta,
                      al_delta = al_delta))

library(ggplot2)
# theta
ggplot(stack(df[[2]]), aes(colour = ind,  values)) +
  geom_density()

# lamda
ggplot(stack(df[[3]]), aes(colour = ind,  values)) +
  geom_density()

# delta
ggplot(stack(df[[4]]), aes(colour = ind,  values)) +
  geom_density()

# grains
grains_df <- do.call(rbind, grains)
grains_df$sample <- gsub("\\..*$", "", rownames(grains_df))
ggplot(grains_df, aes(V1, fill = sample)) +
  geom_histogram()

# explore a negative binomial fit...

N <- 100
x <- rnbinom(N, 10, .25)
hist(x, 
     xlim=c(min(x),max(x)), probability=T, nclass=max(x)-min(x)+1, 
     col='lightblue',
     main='Negative binomial distribution, n=10, p=.25')
lines(density(x,bw=1), col='red', lwd=3)

# load library
library(fitdistrplus)

# fit the negative binomial distribution
samp <- grain_2$V1
fit <- fitdist(samp, "nbinom",  method = "mme")

# get the fitted densities. mu and size from fit.
fitD <- dnbinom(0:(length(samp)*2), size = unname(fit$estimate[1]) , mu = unname(fit$estimate[2]))

# add fitted line (blue) to histogram
plot(density(samp), col='red', lwd=3)
lines(fitD, lwd="3", col="blue")