R help - Feb 2013 - Bayesian mixing model

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]]