R/mean.mu.R

In ctmm-initiative/ctmm: Continuous-Time Movement Modeling

# match COV.mu structure to flattened mean c(mu) # [axes*m]
flatten.cov.mu <- function(COV.mu)
{
  if(length(dim(COV.mu))==4) # [axes,m,m,axes]
  {
    COV.mu <- aperm(COV.mu,c(2,1,3,4)) # [m,axes,m,axes]
    DIM <- dim(COV.mu)
    DIM <- prod(DIM[1:2])
    dim(COV.mu) <- c(DIM,DIM) # [axes*m,axes*m]
  }

  return(COV.mu)
}

##########
mean.mu <- function(x,debias=TRUE,weights=NULL,trace=FALSE,IC="AICc",...)
{
  if(is.null(weights))
  { weights <- rep(1,length(x)) }

  axes <- x[[1]]$axes
  AXES <- length(axes)
  N <- length(x)

  mean <- sapply(x,function(y){y$mean})
  mean <- unique(mean)

  # dimensions of mean functions
  m <- sapply(x,function(y){nrow(y$mu)})
  M <- max(m)

  # Gaussian-Gaussian in all cases
  MU <- array(0,c(N,AXES*M))
  SIGMA <- array(0,c(N,AXES*M,AXES*M))

  if(M==1)
  {
    colnames(MU) <- axes
    dimnames(SIGMA) <- list(NULL,axes,axes)
  }

  for(i in 1:N)
  {
    # fill in for missing variances
    diag(SIGMA[i,,]) <- Inf

    # existing indices
    SUB <- logical(M*AXES)
    for(j in 1:AXES) { SUB[(j-1)*M+1:m[i]] <- TRUE }
    SUB <- which(SUB)

    MU[i,SUB] <- c(x[[i]]$mu)
    SIGMA[i,SUB,SUB] <- flatten.cov.mu(x[[i]]$COV.mu)
  }

  # list of candidate models
  MM <- list()

  # no variance in mean location
  S <- "Dirac-\u03B4(\u03BC)"
  if(trace) { message("Fitting location-mean model ",S) }
  MM[[S]] <- meta.normal(MU,SIGMA,VARS=FALSE,isotropic=TRUE,debias=debias,weights=weights,WARN=FALSE,...)
  GUESS <- MM[[S]]

  # symmetric mean distribution
  if(AXES*M*N >= AXES*M+1)
  {
    S <- "isotropic-\u03BC"
    if(trace) { message("Fitting location-mean model ",S) }
    MM[[S]] <- meta.normal(MU,SIGMA,isotropic=TRUE,debias=debias,weights=weights,GUESS=GUESS,WARN=FALSE,...)
    GUESS <- MM[[S]]

    # general distribution
    if(AXES*M*N >= AXES*M+(AXES*M*(AXES*M+1))/2)
    {
      S <- "anisotropic-\u03BC"
      if(trace) { message("Fitting location-mean model ",S) }
      MM[[S]] <- meta.normal(MU,SIGMA,debias=debias,weights=weights,GUESS=GUESS,WARN=FALSE,...)
    }
  }

  ICS <- sapply(MM,function(m){m[[IC]]})
  names(ICS) <- names(MM)
  i <- which.min(ICS)
  i <- names(MM)[i]
  MM <- MM[[i]]
  MM$name <- i

  # report model selection
  if(trace)
  {
    i <- sort(ICS,index.return=TRUE)$ix
    ICS <- ICS[i] # sorted
    ICS <- ICS - ICS[1] # start at zero
    ICS <- data.frame(ICS)
    colnames(ICS) <- paste0("\u0394",IC)
    message("* Model selection for location-mean \u03BC distribution.")
    print(ICS)
  }

  # R$mu # population mean of mean locations
  # R$COV.mu # uncertainty in mean of means estimate
  names(MM)[ which(names(MM)=="sigma") ] <- "POV.mu" # dispersion of means
  names(MM)[ which(names(MM)=="COV.sigma") ] <- "COV.POV.mu" # uncertainty in dispersion of means

  if(M>1)
  {
    MM$mu <- array(MM$mu,c(M,AXES))

    dim(MM$COV.mu) <- c(M,AXES,M,AXES)
    MM$COV.mu <- aperm(MM$COV.mu,c(2,1,3,4))

    dim(MM$POV.mu) <- c(M,AXES,M,AXES)
    MM$POV.mu <- aperm(MM$POV.mu,c(2,1,3,4))
  }

  MM$mean <- mean

  return(MM)
}