#' Generate an unstructured signal
#'
#' @param m An integer value, the number of hypotheses
#' @param pi0 A numeric value in \eqn{[0,1]}, the proportion of true null signals
#' @param barmu A numeric value, the strength of the signal
#' @return A numeric vector of length \code{m}, the signal values for each hypothesis
#' @export
#' @examples
#' m <- 100
#' pi0 <- 0.7
#' barmu <- 4
#' mu <- gen.mu.noleaves(m = m, pi0 = pi0, barmu = barmu)
#' plot(mu)
gen.mu.noleaves <- function(m, pi0, barmu) {
mu <- numeric(m)
m0 <- floor(m * pi0)
m1 <- m - m0
mu[sample(seq(1, m), m1)] <- barmu
mu
}
#' Generate a block of deterministic signal of Gaussian shape
#'
#' @param barmu A numeric value, the strength of the signal
#' @param len An integer value, the length of the block
#' @return A numeric vector of length \code{len}, the signal values in the block
gauss_bloc <- function(barmu, len) {
# this choice of sigma2 is so that the values at the edge of the vector 'out' are = 1/200
sigma2 <- len^2/(8 * log(200))
out <- exp(-(seq_len(len) - floor(len/2))^2/(2 * sigma2))
return(barmu * out/mean(out))
}
#' Generate the signal in the leaves of the tree
#'
#' @param m An integer value, the number of hypotheses
#' @param K1 An integer value, the number of non-null leaves
#' @param d A numeric value in \eqn{[0,1]}, the proportion of non-null signals in each non-null leaf
#' @param grouped A boolean value, whether the non-null leaves are contiguous or not
#' @param setting A character value specifying the shape of the signal in each non-null leaf
#' @param barmu A numeric value, the strength of the signal
#' @param leaf_list A list of leaves as generated by \code{dyadic.from.Nnn} functions
#' @return A numeric vector of length \code{m}, the signal values for each hypothesis
#' @details If \code{setting == "const"} the non-null signal is constant equal to \code{barmu}. If \code{setting == "gauss"} and \eqn{d=1} the signal in an active leaf has the shape of a Gaussian bell with mean \code{barmu}, see the \code{\link{gauss_bloc}} function which generates it. If \eqn{d<1} it has the same shape but randomly pruned. If \code{setting == "poisson"} the non-null signal is randomly drawn according to a Poisson distribution of mean \code{barmu} (slightly modified to not yield 0). If \code{setting == "rgauss"} the non-null signal is randomly drawn according to a Normal distribution of mean \code{barmu} and variance \code{barmu}.
#' @seealso \code{\link{gauss_bloc}}
#' @export
#' @importFrom stats rpois rnorm
#' @examples
#' m <- 160
#' s <- 10
#' K1 <- floor(m/(s * 4))
#' d <- 1
#' barmu <- 4
#' dd <- dyadic.from.window.size(m, s, method = 2)
#' leaf_list <- dd$leaf_list
#' muC <- gen.mu.leaves(m = m, K1 = K1, d = d, grouped = FALSE,
#' setting = "const", barmu = barmu, leaf_list =leaf_list)
#' muC <- gen.mu.leaves(m = m, K1 = K1, d = d, grouped = TRUE,
#' setting = "const", barmu = barmu, leaf_list =leaf_list)
#' muG <- gen.mu.leaves(m = m, K1 = K1, d = d, grouped = FALSE,
#' setting = "gauss", barmu = barmu, leaf_list =leaf_list)
#' muP <- gen.mu.leaves(m = m, K1 = K1, d = d, grouped = FALSE,
#' setting = "poisson", barmu = barmu, leaf_list =leaf_list)
#' mu <- cbind(constant = muC, Gaussian = muG, Poisson = muP)
#' matplot(mu, t = 's')
gen.mu.leaves <- function(m, K1, d, grouped, setting, barmu, leaf_list) {
mu <- numeric(m)
K <- length(leaf_list)
if (K1 > K)
stop("K1>K,\nwe don't have so many leaves")
active_leaves <- numeric(K1)
if (grouped) {
active_leaves <- seq(1, K1) # + sample(seq(0, K - K1), 1)
} else {
active_leaves <- sample(seq(1, K), K1)
}
for (i in active_leaves) {
length_leaf <- length(leaf_list[[i]])
m1loc <- floor(length_leaf * d)
signal <- sample(length_leaf, m1loc)
mu[leaf_list[[i]][signal]] <- switch(setting, const = {
barmu
}, rgauss = {
rnorm(m1loc, mean = barmu, sd = sqrt(barmu))
}, gauss = {
gauss_bloc(barmu, m1loc)
# gauss_bloc(barmu, length_leaf)[signal] # old code, broken because
# that way the signal doesn't have mean barmu
# but lesser mean as soon as d<1
}, poisson = {
rpois(m1loc, 999 * barmu/1000) + barmu/1000
})
}
return(mu)
}
#' Generate one-sided p-values associated to a given signal with equi-correlated noise
#'
#' @param m An integer value, the number of hypotheses
#' @param mu A vector of \eqn{m} signal values
#' @param rho A numeric value in \eqn{[0,1]}, the level of equi-correlation between test statistics
#' @param alternative A character string specifying the alternative hypothesis.
#' Must be one of "two.sided" (default), "greater" or "less".
#' @return A vector of \eqn{m} one-sided \eqn{p}-values
#' @export
#' @importFrom stats pnorm
#' @examples
#' m <- 100
#' s <- 10
#' K1 <- floor(m/(s * 4))
#' d <- 1
#' barmu <- 4
#' dd <- dyadic.from.window.size(m, s, method = 2)
#' leaf_list <- dd$leaf_list
#' mu <- gen.mu.leaves(m = m, K1 = K1, d = d, grouped = FALSE,
#' setting = "const", barmu = barmu, leaf_list =leaf_list)
#' pvals <- gen.p.values(m = m, mu = mu, rho = 0)
#' plot(-log(pvals), t = 'b')
gen.p.values <- function(m, mu = rep(0, m), rho = 0, alternative = c("two.sided", "less", "greater")) {
alternative <- match.arg(alternative)
stopifnot(m == length(mu)) ## sanity check
Z <- rnorm(m + 1, 0, 1)
Y <- sqrt(1 - rho) * Z[seq_len(m)] + sqrt(rho) * Z[m + 1]
if (alternative == "two.sided"){
return(2 * pnorm(abs(Y + mu), lower.tail = FALSE))
}
else if (alternative == "less"){
return(pnorm(Y + mu, lower.tail = TRUE))
}
else if (alternative == "greater"){
return(pnorm(Y + mu, lower.tail = FALSE))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.