# NOT RUN {
##====================================
# example using bceInput data
# first try
X <- BCE(bceInput$Rat,bceInput$Dat,relsdRat=.2,relsdDat=.2,
iter=1000,outputlength=5000,jmpX=.01,jmpRat=.01)
## the number of accepted runs is too low;
## we play around with the jump lengths jmpX and jmpRat
X <- BCE(bceInput$Rat,bceInput$Dat,relsdRat=.2,relsdDat=.2,
iter=1000,outputlength=5000,jmpX=.02,jmpRat=.002)
## we inspect the output:
plot(X)
## For every element of X and Rat, we want to obtain a well-mixed,
## random trace. In this case, mixing is still a little poor.
## to optimize mixing in the ratio matrix, it is a good idea
## to make the jump length linear to the ratio matrix
## standard deviation (sdrat=.2*rat) :
X <- BCE(bceInput$Rat,bceInput$Dat,relsdRat=.2,relsdDat=.2,
iter=1000,outputlength=5000,jmpX=.02,
jmpRat=.2*(.2*bceInput$Rat))
plot(X)
## mixing improved a lot; we repeat the run with more iterations
## to improve the reliability of the results.
## the following run can take a few minutes - so it is toggled off
#X <- BCE(bceInput$Rat,bceInput$Dat,relsdRat=.2,relsdDat=.2,
# iter=100000,outputlength=5000,jmpX=.02,
# jmpRat=.2*(.2*bceInput$Rat))
#plot(X)
## you can see in the plots that traces for all elements of Rat and X
## are well-mixed. This run was saved in "bceOutput"
Sum <-summary(bceOutput)
# show results as mean with ranges
print(Sum$meanX)
# plot estimated means and ranges (lbX=lower, ubX=upper bound)
xlim <- range(c(Sum$lbX,Sum$ubX))
# first the mean
dotchart(x=t(Sum$meanX),xlim=xlim,
main="Taxonomic composition",
sub="using bce",pch=16)
# then ranges
nr <- nrow(Sum$meanX)
nc <- ncol(Sum$meanX)
for (i in 1:nr)
{ip <-(nr-i)*(nc+2)+1
cc <- ip : (ip+nc-1)
segments(t(Sum$lbX[i,]),cc,t(Sum$ubX[i,]),cc)
}
# show results as pairs plot
pairs(bceOutput,sample=3,main="Station 3")
# }
Run the code above in your browser using DataLab