Two_Node_Process

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(BayesMassBal)

The function BMB is used with a two node process and simulated data.

yshift <- 1

rekt <- data.frame(matrix(NA, ncol = 4, nrow = 2))
names(rekt) <- c("xleft", "ybottom", "xright", "ytop")
rekt$xleft <- c(3,7)
rekt$ybottom <- 4
rekt$xright <- c(5,9)
rekt$ytop <- 6

aros <- data.frame(matrix(NA, ncol = 4, nrow = 5))
names(aros) <-c("x0","y0","x1","y1")
aros[1,] <- c(1,5,rekt$xleft[1],5)
aros[2,] <- c(rekt$xright[1],5,rekt$xleft[2],5)
aros[3,] <- c(rekt$xright[2],5,11,5)
aros[4,] <- c(rekt$xright[1] - 1, rekt$ybottom[1],rekt$xright[1] - 1, rekt$ybottom[1] - 2)
aros[5,] <- c(rekt$xright[2] - 1, rekt$ybottom[2], rekt$xright[2] - 1 , rekt$ybottom[2] - 2)

aros$y0 <- aros$y0
aros$y1 <- aros$y1

b.loc <- data.frame(matrix(NA, ncol = 2, nrow = 5))
names(b.loc) <- c("x","y")
b.loc[1,] <- c(0.5,aros$y0[1])
b.loc[2,] <- c(mean(c(aros$x0[2],aros$x1[2])),aros$y0[2] + 0.6)
b.loc[3,] <- c(aros$x1[3] + 0.5, aros$y1[3])
b.loc[4,] <- c(aros$x1[4], aros$y1[4] - 0.6)
b.loc[5,] <- c(aros$x1[5], aros$y1[5] - 0.6)

p.loc <- data.frame(matrix(NA, ncol = 2, nrow = 2))
names(p.loc) <- c("x","y")

p.loc$x <- rekt$xleft + 1
p.loc$y <- rekt$ybottom + 1

par(mar = c(0.1,0.1,0.1,0.1))
plot(1, type="n", xlab="", ylab="", xlim=c(0, 12), ylim=c(1, 6), axes = FALSE)
rect(xleft = rekt$xleft, ybottom = rekt$ybottom, xright = rekt$xright, ytop =rekt$ytop, col = "skyblue")
arrows(aros$x0,aros$y0, x1= aros$x1,y1 = aros$y1, code = 2, length = 0.1)
for(i in 1:5){
  text(b.loc$x[i], b.loc$y[i],labels = bquote(y[.(i)]), adj = c(0.5,0.5), cex = 1.2)
}
for(i in 1:2){
  text(p.loc$x[i], p.loc$y[i], labels = bquote(P[.(i)]), adj = c(0.5,0.5), cex = 1.2)
}
text(1.5,4.6,labels=  "F", adj = c(0.5,0.5), cex = 0.7)
text(c(5.5,9.5),c(4.6,4.6), labels = "C", adj= c(0.5,0.5), cex = 0.7)
text(c(4.4,8.4),c(3.5,3.5), labels = "T", adj= c(0.5,0.5), cex = 0.7)

The constraints around these process nodes are:

\begin{align} y_1 &= y_2 +y_4\ y_2 &= y_3 +y_5 \end{align}

Therefore the matrix of constraints, C is:

C <- matrix(c(1,-1,0,-1,0,0,1,-1,0,-1), nrow = 2, ncol = 5, byrow = TRUE)
C

The constrainProcess function in the BayesMassBal package is used to generate an X matrix based on C that will later be used with the BMB function.

X <- constrainProcess(C = C)
X

Constraints can also be imported from a .csv file. The path to a file, included in the BayesMassBal package, for this process can be found and constraints can be imported by specifying the location for the file argument for constrainProcess as shown below:

constraint_file_location <- system.file("extdata", "twonode_constraints.csv",package = "BayesMassBal")
X <- constrainProcess(file = constraint_file_location)

The previously simulated data is loaded from a .csv file using the importObservations() function. The local location of the the file imported below can be found by typing system.file("extdata", "twonode_example.csv",package = "BayesMassBal"). View the document in Excel to see how your data should be formatted for import. Note: it is not required that the entries into the *.csv file are separated by ";".

y <- importObservations(file = system.file("extdata", "twonode_example.csv",
                                  package = "BayesMassBal"),
                  header = TRUE, csv.params = list(sep = ";"))

Then, the BMB function is used to generate the distribution of constrained masses from the data with cov.structure = "indep".

indep.samples <- BMB(X = X, y = y, cov.structure = "indep", BTE = c(100,3000,1), lml = TRUE, verb = 0)

The output of BMB is a BayesMassBal object. Special instructions are designated when feeding a BayesMassBal object to the plot() function. Adding the argument layout = "dens" and indicating the mass balanced flow rate for CuFeS2 at $y_3$ should be plotted using a list supplied to sample.params, the desired distribution can be plotted with its 95% Highest Posterior Density Interval.

plot(indep.samples,sample.params = list(ybal = list(CuFeS2 = 3)),
    layout = "dens",hdi.params = c(1,0.95))

It is also possible to generate trace plots to inspect convergence of the Gibbs sampler. Here are trace plots for $\beta$

plot(indep.samples,sample.params = list(beta = list(CuFeS2 = 1:3, gangue = 1:3)),layout = "trace",hdi.params = c(1,0.95))

A quantitative diagnostics for convergence and autocorrelation are available as part of the output from BMB:

indep.samples$diagnostics

The model with independent variances may not be the best fitting model. Models specifying covariance between sample locations for a single component, and covariance between components at a single location are fit.

component.samples <- BMB(X = X, y = y, cov.structure = "component", BTE = c(100,3000,1), lml = TRUE, verb = 0)
location.samples <- BMB(X = X, y = y, cov.structure = "location", BTE = c(100,3000,1), lml = TRUE, verb = 0)

Computing $\log(\mathrm{Bayes Factor})$ for $BF = p(y|\texttt{indep})/p(y|\texttt{component})$:

indep.samples$lml - component.samples$lml

Then comparing $p(y|\texttt{component})$ to $p(y|\texttt{location})$

component.samples$lml - location.samples$lml

Shows there is little difference between the models where cov.structure = "location" and cov.structure = "component", but both of these models better explain the data than cov.structure = "indep".

We can view a summary of the favored model by passing a BayesMassBal object to the summary function. While not done in this case, the summary table can be saved by passing the desired name of a *.csv file to the export argument.

summary(component.samples, export = NA)

The main effect of a variable independent of the process can be calculated by supplying a function, fn that takes the arguments of mass balanced flow rates ybal, and the random independent and uniformly distributed variables x. Information can be gained on the main effect of a particular element of x, xj, on fn using the mainEff function. Output from mainEff includes information on the distribution of $E_x\lbrack f(x,y_{\mathrm{bal}})|x_j \rbrack$.

fn_example <- function(X,ybal){
      cu.frac <- 63.546/183.5
      feed.mass <- ybal$CuFeS2[1] + ybal$gangue[1]
      # Concentrate mass per ton feed
      con.mass <- (ybal$CuFeS2[3] + ybal$gangue[3])/feed.mass
      # Copper mass per ton feed
      cu.mass <- (ybal$CuFeS2[3]*cu.frac)/feed.mass
      gam <- c(-1,-1/feed.mass,cu.mass,-con.mass,-cu.mass,-con.mass)
      f <- X %*% gam
      return(f)
      }

rangex <- matrix(c(4.00 ,6.25,1125,1875,3880,9080,20,60,96,208,20.0,62.5),
                   ncol = 6, nrow = 2)
mE_example <- mainEff(indep.samples, fn = "fn_example",rangex =  rangex,xj = 3, N = 25, res = 25)

A plot of the output can be made. To get lines that are better connected, change increase N in the mainEff function.

m.sens<- mE_example$fn.out[2,]
hpd.sens <- mE_example$fn.out[c(1,3),]
row.names(hpd.sens) <- c("upper", "lower")
g.plot <- mE_example$g/2000

y.lim <- range(hpd.sens)

lzero.bound <- apply(hpd.sens,1,function(X){which(X <= 0)})
lzero.mean <- which(m.sens <= 0)

main.grid <- pretty(g.plot)
minor.grid <- pretty(g.plot,25)
minor.grid <- minor.grid[-which(minor.grid %in% main.grid)]

y.main <- pretty(hpd.sens)

opar <- par(no.readonly =TRUE) 
par(mar = c(4.2,4,1,1))
plot(g.plot,m.sens, type = "n", xlim = range(g.plot), ylim = y.lim, ylab = "Net Revenue ($/ton Feed)", xlab=  "Cu Price ($/lb)")

abline(v = main.grid, lty = 6, col = "grey", lwd = 1)
abline(v = minor.grid, lty =3, col = "grey", lwd = 0.75)

abline(h = 0, col = "red", lwd = 1, lty = 6)

lines(g.plot[lzero.mean],m.sens[lzero.mean],col = "red", lwd =2)
lines(g.plot[-lzero.mean[-length(lzero.mean)]],m.sens[-lzero.mean[-length(lzero.mean)]],col = "darkgreen", lwd =2)

lines(g.plot[lzero.bound$lower],hpd.sens[2,][lzero.bound$lower], lty = 5, lwd = 2, col = "red")
lines(g.plot[-lzero.bound$lower],hpd.sens[2,][-lzero.bound$lower], lty = 5, lwd = 2, col = "darkgreen")

lines(g.plot[lzero.bound$upper],hpd.sens[1,][lzero.bound$upper], lty = 5, lwd = 2, col = "red")
lines(g.plot[-lzero.bound$upper],hpd.sens[1,][-lzero.bound$upper], lty = 5, lwd = 2, col= "darkgreen")

legend("topleft", legend = c("Expected Main Effect", "95% Bounds", "Net Revenue < $0", "Net Revenue > $0"), col = c("black","black","red", "darkgreen"), lty = c(1,6,1,1), lwd = c(2,2,2,2), bg = "white")

par(opar)


Try the BayesMassBal package in your browser

Any scripts or data that you put into this service are public.

BayesMassBal documentation built on June 18, 2022, 1:08 a.m.