mincrit <-function (obs, param, sumstats, obspar = NULL, abcmethod = abc,
crit = nn.ent, sumsubs = 1:ncol(sumstats), limit = length(sumsubs),
do.only = NULL, verbose = TRUE, do.crit = TRUE, do.err = FALSE,
final.dens = FALSE, errfn = rsse, ...)
{
argl <- list(...)
targind <- match(names(argl), "tol")
targind <- which(!is.na(targind))
margind <- match(names(argl), "method")
margind <- which(!is.na(margind))
if ((length(targind) == 0) & identical(abcmethod, abc)) {
argl$tol <- 0.01
}
if ((length(margind) == 0) & identical(abcmethod, abc)) {
argl$method <- "rejection"
}
if (!is.matrix(obs) | is.data.frame(obs)) {
obs <- matrix(obs, nrow = 1)
}
if (!is.matrix(param) | is.data.frame(param)) {
param <- as.matrix(param)
}
if (!is.matrix(sumstats) | is.data.frame(sumstats)) {
sumstats <- as.matrix(sumstats)
}
if (nrow(sumstats) != nrow(param)) {
stop("param and sumstats should be matrices/dataframes with the same number of rows.")
}
if (nrow(sumstats) <= 1 | nrow(param) <= 1) {
stop("Too few simulated datasets to perform ABC. The code will run with >=2, but a sensible number is at least thousands.")
}
if (!is.null(obspar) | is.data.frame(obspar)) {
if (!is.matrix(obspar)) {
obspar <- matrix(obspar, byrow = T, ncol = ncol(param))
}
if (nrow(obs) != nrow(obspar)) {
stop("Please supply observed statistics and observed parameter matrices with the same number of rows!\n")
}
}
if (any(is.na(sumstats)) | any(is.na(param)) | any(is.na(obs))) {
stop("Simulations and observations must not contain NAs")
}
argl$param=param
if (!length(colnames(param))) {
colnames(param) <- paste("P", 1:ncol(param), sep = "")
}
if (!length(colnames(sumstats))) {
colnames(sumstats) <- paste("C", 1:ncol(sumstats), sep = "")
}
sumstats <- sumstats[, sumsubs]
obs <- obs[, sumsubs]
data <- matrix(obs, nrow = 1)
npar <- ncol(param)
nr <- nrow(sumstats)
nstats <- length(sumsubs)
q <- (!is.null(obspar)) & (do.err)
if (!q) {
do.err <- FALSE
}
cm <- combmat(nstats, limit)
nc <- nrow(cm)
dom<-(0%in%do.only)
if(dom &!is.matrix(do.only)){
do.only<-matrix(do.only,nrow=1)
}
if (is.null(do.only) || max(do.only)>nc) {
do.only <- 1:nc
}
do.only<-as.matrix(do.only)
if(dom & (ncol(do.only) > nstats)){
stop("supplied subset matrix (do.only) does not match number of given summaries (ncol(sumsubs)) \n")
}
if(dom){
ldo<-nrow(do.only)
}
else{
ldo<-length(do.only)
}
vals <- err <- NULL
critvals <- matrix(0, 1, length(do.only))
best <- 1
# for (i in 1:length(do.only)) {
for (i in 1:ldo) {
if(dom){
I <- which(do.only[i, ] == 1)
}
else{
I <- which(cm[do.only[i], ] == 1)
}
if (verbose) {
cat("doing statistics: ", sumsubs[I], " (", i,
"/", ldo, ")\n")
}
argl$target=data[I]
argl$sumstat=sumstats[,I]
valsi <- do.call(abcmethod,argl)
# valsi <- eval(parse(text=paste("abcmethod(data[I], param, sumstats[, I],",addstr, "...)" )))
if (is.null(valsi$adj.values)) {
valsi <- valsi$unadj.values
}
else {
valsi <- valsi$adj.values
}
if (do.err) {
err[i] <- errfn(valsi, obspar, apply(param, 2, var))
}
if (do.crit) {
cat("doing criterion calculation...\n")
critvals[i] <- crit(valsi)
stick <- (critvals[best] <= critvals[i])
best <- ifelse(stick, best, i)
if (final.dens) {
if ((i == 1) | !stick) {
vals <- valsi
}
}
}
}
err <- err[best]
l <- list()
if (do.crit) {
l$critvals <- critvals
if(dom){
besti<-best
best <- matrix(which(do.only[besti, ]==1),nrow=1)
}
else{
besti <- do.only[best]
best <- matrix(which(cm[besti, ] == 1),nrow=1)
}
rownames(best)<-besti
l$best <- best
}
if (do.err) {
l$err <- err
}
if (final.dens) {
l$post.sample <- vals
}
l$posssubs <- do.only
l$sumsubs <- sumsubs
return(l)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.