Nothing
#' calculate weights inverse to the distances from the specified observation.
#'
#' for calculation of a shape model by averaging the observations neighbouring
#' the configuration in question, it is necessary to calculate weights by
#' similarity.
#'
#' distances of zero will get a weight of 1e12 (this is scaled to all weights
#' summing to one), thus weights for observations further away are converging
#' to zero.
#'
#' @param data array containing landmark configurations
#' @param number integer: how many of the neighbours are to be involved.
#' @param ref integer: position in the array that is used as reference.
#' @param report logical: require report about name of the reference.
#' @param reg numeric: regularise mahalanobis distance by adding reg to the
#' diagonal of eigenvalues of the covariance matrix.
#' @param log logical: use the logarithm of the distances.
#' @param mahalanobis logical: use mahalanobis distance.
#' @param weightfun custom function that operates on a vector of distances (see examples) and generates weights accordingly.
#' @return
#' \item{data }{dataframe containing id, procrustes/mahalanobis distance
#' and weight according to the reference}
#' \item{reference }{returns observations' names if available}
#' \item{rho.all }{dataframe containing distances to references of all observations}
#' @examples
#'
#' if (require(shapes)) {
#' proc <- procSym(gorf.dat)
#' ##get weights for the the four specimen closest to the first observation.
#' weights <- proc.weight(proc$rotated,4,1)
#'
#' ##estimate the first specimen by weighted neighbour shapes.
#' estim <- proc$mshape*0;
#' for (i in 1:4)
#' {estim <-estim+proc$rotated[,,weights$data$nr[i]]*weights$data$weight[i]}
#'
#' ### visualise
#' plot(estim,asp=1)## show estimation
#' points(proc$rotated[,,1],col=3)##show original
#'
#' ## use a gaussian smoother to compute weights using a bandwidth of 0.05
#' gaussWeight <- function(r,sigma=0.05) {
#' sigma <- 2*sigma^2
#' return(exp(-r^2/ sigma))
#' }
#' weights <- proc.weight(proc$rotated,4,1,weightfun=gaussWeight)
#' }
#' @export
proc.weight <- function(data,number,ref,report=TRUE,reg=0,log=FALSE,mahalanobis=FALSE,weightfun=NULL) {
lmdat <- FALSE
col <- 3
rho <- 0
if (length(dim(data))==3) {
data <- vecx(data)
lmdat <- TRUE
}
l <- dim(data)[1]
obsnames <- dimnames(data)[[1]]
if (lmdat && !mahalanobis) {
for(i in 1:l)
rho[i] <- angle.calc(data[ref,],data[i,])
if (is.null(obsnames))
id <- as.character(c(1:l))
else
id <- obsnames
} else {
if (mahalanobis) {
covtmp <- cov(data)
if (reg != 0) {
eig <- eigen(covtmp,symmetric=TRUE)
covtmp <- t(eig$vectors)%*%diag(eig$values+reg)%*%eig$vectors
}
checksing <- try(covtmp <- solve(covtmp),silent = TRUE)
if (inherits(checksing,"try-error")) {
covtmp <- armaGinv(covtmp)
cat("singular covariance matrix: using general inverse\n")
}
}
else
covtmp <- diag(ncol(data))
rho <- sqrt(mahalanobis(data,cov=covtmp,center=data[ref,],inverted = TRUE))
}
if (is.null(obsnames))
id <- as.character(c(1:l))
else
id <- obsnames
if (log)
rho <- log(rho)
nr <- c(1:l)
data <- data.frame(nr,id,rho)
dat.sort.i <- data[order(data[,col]),]
dat.which <- dat.sort.i[2:(number+1),]
if (is.null(weightfun)) {
weightfun <- function(x) {
if (0 %in% x)
x[which(x == 0)] <- 1/1e12
weight <- 1/x
weight <- weight/sum(weight)
return(weight)
}
}
weight <- weightfun(dat.which[,col])
weight <- weight/sum(weight)
out <- data.frame(dat.which,weight)
if (report)
cat(paste(" reference was",id[ref],"\n"))
return(list(data=out,reference=id[ref],rho.all=data))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.