function(constraint.means, constraint.matrix, initial.probs = NA) { # MAxent function in R. Written by Bill Shipley # Bill.Shipley@USherbrooke.ca # Improved Iterative Scaling Algorithm from Della Pietra, S. # Della Pietra, V., Lafferty, J. (1997). Inducing features of random fields. # IEEE Transactions Pattern Analysis And Machine Intelligence # 19:1-13 #............................. # constraint.means is a vector of the community-aggregated trait values # constraint.matrix is a matrix with each trait defining a row and each # species defining a column. Thus, two traits and 6 species would produce # a constraint.matrix with two rows (the two traits) and 6 columns (the # six species) if(!is.vector(constraint.means)) return( "constraint.means must be a vector") if(is.vector(constraint.matrix)) { constraint.matrix <- matrix(constraint.matrix, nrow = 1, ncol = length( constraint.matrix)) } n.constraints <- length(constraint.means) dim.matrix <- dim(constraint.matrix) if(dim.matrix[2] == 1 & dim.matrix[1] > 1) { constraint.matrix <- t(constraint.matrix) dim.matrix <- dim(constraint.matrix) } if(n.constraints != dim.matrix[1]) return("Number of constraint means not equal to number in constraint.matrix" ) n.species <- dim.matrix[2] if(is.na(initial.probs)) { prob <- rep(1/n.species, n.species) } if(!is.na(initial.probs)) { prob <- initial.probs } # C.values is a vector containing the sum of each trait value over the species C.values <- apply(constraint.matrix, 1, sum) test1 <- test2 <- 1 while(test1 > 1e-008 & test2 > 1e-008) { denom <- constraint.matrix %*% prob gamma.values <- log(constraint.means/denom)/ C.values if(sum(is.na(gamma.values)) > 0) return("missing values in gamma.values" ) unstandardized <- exp(t(gamma.values) %*% constraint.matrix) * prob new.prob <- as.vector(unstandardized/sum( unstandardized)) if(sum(is.na(new.prob)) > 0) return("missing values in new.prob") test1 <- 100 * sum((prob - new.prob)^2) test2 <- sum(abs(gamma.values)) prob <- new.prob } list(probabilities = prob, final.moments = denom, constraint.means = constraint.means, constraint.matrix = constraint.matrix, entropy = -1 * sum(prob * log(prob))) }