On 26-10-2012, at 12:50, Richard James wrote: > Dear Berend and Thomas, > > thank you for suggesting the lsei function. I found that the tlsce {BCE} > function also works very well: > > library("BCE") > tlsce(A=bmat,B=target) > > The limSolve package has an 'xsample' function for generating uncertainty > values via Monte-Carlo simulation, however it only works when specifying the > standard deviation on the target data (B). In my situation I have standard > deviations for the source areas (A) only. Therefore I need to generate the > uncertainty values manually. > > I've created a matrix of 1000 randonly distributed numbers for each of the > source area (A) properties: > > TSCa<-matrix(rnorm(1000, mean=0.03, sd=0.005),ncol=1) > TSMg<-matrix(rnorm(1000, mean=0.0073, sd=0.002),ncol=1) > CBCa<-matrix(rnorm(1000, mean=0.6, sd=0.1),ncol=1) > CBMg<-matrix(rnorm(1000, mean=0.0102, sd=0.005),ncol=1) > RCa<-matrix(rnorm(1000, mean=0.2, sd=0.05),ncol=1) > RMg<-matrix(rnorm(1000, mean=0.0141, sd=0.005),ncol=1) > DCa<-matrix(rnorm(1000, mean=0.35, sd=0.1),ncol=1) > DMg<-matrix(rnorm(1000, mean=0.012, sd=0.004),ncol=1) > > DistAll<-cbind(TSCa,TSMg,CBCa,CBMg,RCa,RMg,DCa,DMg) > colnames(DistAll)<-c("TopSoilCa","TopSoilMg","ChannelBankCa","ChannelBankMg","RoadCa","RoadMg","DrainCa","DrainMg") > > I now want to run the lsei model again: > > lsei(A=bmat,B=target,E=matrix(rep(1,4),ncol=4),F=1,G=matrix(rep(1,4),ncol=4),H=0, > fulloutput=TRUE)) >
Actually the G and H arguments are wrong. They should be G=diag(1,4) H=rep(0,4) since the weights should be >= 0 > but with A=bmat replaced by the appropriate random values in DistAll. I > assume I could then use the function "replicate" to then run the model 1000 > times to generate uncertainty values, e.g. > > replicate(n=1000,lsei(A=bmat,B=target,E=matrix(rep(1,4),ncol=4),F=1,G=matrix(rep(1,4),ncol=4),H=0, > fulloutput=TRUE)) > > Would anyone be able to help me write a function to replace bmat with new > values from DistAll each time the lsei model is run? > Maybe this will do what you need Nrep <- 10 # fors testing I think this is what you need ### start code kRow <- 1 bmat <- matrix(DistAll[kRow,],ncol=4,byrow=FALSE) bmat target <- c("Ca in river(%)"=0.33,"Mg in river (%)"=0.0114) target lsei(A=bmat,B=target,E=matrix(rep(1,4),ncol=4),F=1,G=diag(1,4),H=rep(0,4),fulloutput=TRUE) gen.single <- function(k,Distall,target) { bmat <- matrix(DistAll[k,],ncol=4,byrow=FALSE) z <- lsei(A=bmat,B=target,E=matrix(rep(1,4),ncol=4),F=1,G=diag(1,4),H=rep(0,4),fulloutput=TRUE) # insert tests of output of lsei to see if all is ok etc. z$X # weights } set.seed(2001) t(sapply(1:Nrep, FUN=gen.single,Distall=Distall,target=target)) ### end of code And now you can do whatever you wish with the columns of the output matrix Berend ______________________________________________ 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.