Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
osl_calibration
# 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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment