R/tmp/estimate.R

In hosscine/phacm: The R Package for Analyzing Persistent Homology and Cycles of Manifold

#' Title
#'
#' @param X
#' @param maxdim
#' @param maxscale
#' @param const.band
#' @param maximum.thresh
#' @return
#' @export
#' @examples
bootstrap.homology <- function(X, maxdim, maxscale, const.band = 0, maximum.thresh = F) {
  # require(pracma)
  if (!("bootsSamples" %in% class(X)))
    stop("input must be bootsSamples")
  peak <- matrix(0, maxdim, length(X))
  # band <- ifelse(const.band > 0,const.band,hausdInterval(X, m=sample.size, B=times, alpha =
  # (1-confidence)))
  tseq <- seq(0, maxscale, length.out = 1000)
  diags <- lapply(X, function(x) calcPhom(x, maxdim, maxscale, ret = T, plot = F))
  print(sapply(diags, function(diag) calcDiagCentroid(diag)[3]))
  band <- ifelse(const.band == 0, max(sapply(diags, function(diag) calcDiagCentroid(diag)[3])),
    const.band)
  print(band)

  for (t in 1:length(X)) {
    land <- lapply(1:maxdim, function(d) TDA::landscape(diags[[t]][[1]], dimension = d, KK = 1,
      tseq = tseq))
    if (maximum.thresh)
      band <- max(sapply(land, max))/4
    for (d in 1:maxdim) {
      peak[d, t] <- countPeaksPL(X = land[[d]], thresh = (band/(2 * d)), tseq = tseq)
    }
  }

  dimnames(peak) <- list(paste0("dim", 1:maxdim), paste0("sample", 1:length(X)))
  bootstrap.summary <- list(peak = peak)
  bootstrap.summary <- append(bootstrap.summary, c(band = band, show.hole.density(peak)))
  class(bootstrap.summary) <- "smoothPhom"
  return(bootstrap.summary)
}

#' Title
#'
#' @param diag
#'
#' @return
#' @export
#'
#' @examples
calcDiagCentroid <- function(diag = diagram) {
  if (class(diag) == "list")
    diag <- diag[[1]]
  diag <- diag[-which(diag[, 1] == 0), ]
  centroid <- apply(diag[, 2:3], 2, mean)
  cpersistence <- (centroid[2] - centroid[1])
  ret <- c(centroid, cpersistence, cpersistence * 2)
  names(ret) <- c("birth", "death", "persistence", "noizes.thresh")
  return(ret)
}

#' Title
#'
#' @param x
#' @param ...
#'
#' @return
#' @export
#'
#' @examples
plot.smoothPhom <- function(x, ...) x <- show.hole.density(x$peak)

# ??<U+383C><U+3E34>次??<U+383C><U+3E33>??<U+383C><U+3E65>?<U+383C><U+3E31>?peak数??<U+383C><U+3E63>?<U+383C><U+3E38>???<U+393C><U+3E32>??<U+383C><U+3E36>度推定し<U+653C><U+3E66>??表示する??<U+383C><U+3E32>
# Internal FUnction
#' Title
#'
#' @param X
#'
#' @return
#' @export
#'
#' @examples
show.hole.density <- function(X) {
  dens <- apply(X, 1, stats::density)
  maxdim <- nrow(X)
  bootstrap.summary <- list()
  xlim = c(min(sapply(dens, function(den) min(den$x))), max(sapply(dens, function(den) max(den$x))))
  ylim = c(0, max(sapply(dens, function(den) max(den$y))))
  plot(1, type = "n", xlim = xlim, ylim = ylim, xlab = "betti number", ylab = "probability")
  for (d in 1:maxdim) {
    mhole <- mean(X[d, ])
    dhole <- dens[[d]][["x"]][which.max(dens[[d]][["y"]])]
    par(new = T)
    plot(dens[[d]], xlim = xlim, ylim = ylim, col = d + 1, ann = F, axes = F)
    print(paste0("dimension ", d, ", ", round(mhole, digits = 2), " mean hole, ", round(dhole),
      " density hole"))
    bootstrap.summary[[paste0("dim", d, "mhole")]] <- mhole
    bootstrap.summary[[paste0("dim", d, "dhole")]] <- dhole
  }
  legend("topright", legend = paste("dim", 1:maxdim), pch = 0, col = 2:(maxdim + 1))
  return(bootstrap.summary)
}

#' Title
#'
#' @param confidence
#' @param maxscale
#' @param haus.dist
#' @param digit
#'
#' @return
#' @export
#'
#' @examples
calc.bootstrap.confidence <- function(confidence, maxscale, haus.dist, digit = 3) {
  Lb <- function(t, haus.dist) sum(length(which(haus.dist > t)))/length(haus.dist)
  max.dist <- max(haus.dist)
  from <- 0
  to <- max.dist
  for (i in 1:digit) {
    tseq <- seq(from, to = to, length.out = 10)
    # ErrorText(from,to,haus.dist,seq(from,to = to,length.out = 10))
    prob <- sapply(tseq, Lb, haus.dist)
    # print(prob)
    from <- tseq[which(prob < confidence)[1] - 1]
    # print(from)
    to <- tseq[which(prob < confidence)[1]]
  }
  c <- 2 * from
  return(c)
}