source("fig_b.R") b1 <- function(i) bm[,sprintf("b[%i,1]", i)] plot(density(b1(1) - b1(3))) # log(NV) - log(EW) plot(density(b1(2) + b1(3) - 2 * b1(1))) # log(TA) + log(EW) - 2 * log(NV)