Fellow R users,
I'm using the BCE {BCE} function to run a Bayesian sediment mixing model. The
aim is to find the optimum % contribution from each of the 4 source areas that
can yield the target geochemistry.
I have geochemistry for 4 source areas called Rat:
Rat<-read.table(text="CaO MgO Na2O Al2O3
Topsoils 2.511250 0.7445500 0.7085500 14.10375
ChannelBanks 55.021500 0.8665000 0.3045000 10.19800
FieldDrains 17.958221 0.9415394 0.2979383 14.68013
RoadRunoff 9.177297 1.9034304 0.4618634 22.22206", header=TRUE)
...and geochemistry for 1 target sediment called Dat:
Dat<-read.table(text=" CaO MgO Na2O Al2O3
Targethf 15.741 1.346 0.368 19.872", header=TRUE)
I have absolute stdev on the sources called adssdRat:
abssdRat<-read.table(text=" CaO MgO Na2O Al2O3
Topsoils 1.8552515 0.1135672 0.06212094 1.491125
ChannelBanks 9.8400162 0.1401057 0.08599080 2.708710
FieldDrains 2.3896499 0.1961217 0.02545431 4.300644
RoadRunoff 0.7780579 0.1749869 0.02848264 1.116747", header =TRUE)
...and absolute stdev on the target called adssdDat:
abssdDat<-read.table(text=" CaO MgO Na2O Al2O3
1 0.877 0.531 0.264 2.439", header=TRUE)
Whilst the model results are OK, there is one problem. The model fails to take
into consideration covariance between the geochemistry in each source area
(Rat) and therefore the uncertainty in the output is likely to be overestimated
(e.g. Mg and Na are strongly positively correlated in each of the 4 source
areas)
My understanding is that this can be corrected for by changing the "userProb"
argument such that the multivariate normal distribution is called for 'Rat'
instead of the default gamma distribution. However I am having trouble in
writing a new "userProb" argument, so I'd be very grateful to anyone who could
help me out with this. Failing that, does anybody know of another package for
Bayesian optimisation that could account for covariance in the source area
geochemistry?
For anyone that can help, you may need the logprobability function that BCE
calls as the default when userProb is set to NULL:
##########
logProbabilityBCE1 <-
function (RAT, X, init.list)
{
with(init.list, {
if (!is.null(userProb))
logp <- log(userProb(RAT, X))
else {
y <- X %*% RAT
if (!is.null(positive)) {
dA.p <- dgamma(RAT[ind.pos], krat[ind.pos], lrat[ind.pos],
log = TRUE)
kdat <- y^2/sddat^2
ldat <- y/sddat^2
kdat[Dat == 0] <- 1
ldat[Dat == 0 & !is.na(Dat)] <- 1/y[Dat == 0 &
!is.na(Dat)]
dB.p <- dgamma(Dat[, positive], kdat[, positive],
ldat[, positive], log = TRUE)
}
else {
dA.p <- 0
dB.p <- 0
}
if (any(wholeranged)) {
dA.r <- dnorm(RAT[ind.r], Rat[ind.r], sdrat[ind.r],
log = TRUE)
dB.r <- dnorm(y[, wholeranged], Dat[, wholeranged],
sddat[, wholeranged], log = TRUE)
}
else {
dA.r <- 0
dB.r <- 0
}
dA.p[dA.p < -1e+08] <- -1e+08
dA.p[dA.p > +1e+08] <- +1e+08
dA.r[dA.r < -1e+08] <- -1e+08
dA.r[dA.r > +1e+08] <- +1e+08
drat <- sum(dA.p, dA.r, na.rm = TRUE)
if (unif)
drat <- 0
if (any(RAT < minRat | RAT > maxRat, na.rm = TRUE))
drat <- -Inf
dB.p[dB.p < -1e+08] <- -1e+08
dB.p[dB.p > +1e+08] <- +1e+08
dB.r[dB.r < -1e+08] <- -1e+08
dB.r[dB.r > +1e+08] <- +1e+08
ddat <- sum(dB.p, dB.r, na.rm = TRUE)/nst
logp <- drat + ddat
}
return(logp)
})}
##########
Regards
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
