Nothing
#' Row-wise Cumulative Sums
#' @description Row-wise Cumulative Sums in mpfr Array.
#'
#' @param mtrx mpfr two-dimensional array.
#' @return mpfr array with row-wise cumulative sums,
#' same dimension as the original array.
#' @examples
#' nill <- Rmpfr::mpfr(0, 120)
#' one <- Rmpfr::mpfr(1, 120)
#' two <- Rmpfr::mpfr(2, 120)
#' contents <- c(one,nill,nill, one,one,one, two,two,two)
#' mtrx3 <- Rmpfr::mpfr2array(contents, dim = c(3, 3))
#' print(mtrx3)
#' print(cumsumm(mtrx3))
#' @export
cumsumm <- function(mtrx) {
mtrs <- mtrx
nrw <- nrow(mtrx)
for (rw in 1:nrw) mtrs[rw, ] <- cumsum(mtrx[rw, ])
mtrs
}
#' Column-Wise Cumulative Sums
#' @description Column-wise cumulative sums in
#' mpfr array.
#'
#' @param mtrx mpfr two-dimensional array.
#' @return mpfr array with column-wise cumulative sums,
#' same dimension as the original array.
#' @examples
#' nill <- Rmpfr::mpfr(0, 120)
#' one <- Rmpfr::mpfr(1, 120)
#' two <- Rmpfr::mpfr(2, 120)
#' contents <- c(one,nill,nill, one,one,one, two,two,two)
#' mtrx3 <- Rmpfr::mpfr2array(contents, dim = c(3, 3))
#' print(mtrx3)
#' print(cumsummcol(mtrx3))
#' @export
cumsummcol <- function(mtrx) {
mtrs <- mtrx
ncl <- ncol(mtrx)
for (cl in 1:ncl) mtrs[, cl] <- cumsum(mtrx[, cl])
mtrs
}
#' Box Cumulative Sums
#' @description A box cumulative sum is defined as the
#' cumulative sum over a lower left rectangle. This function
#' is primarily for use when the components are point
#' probabilities for the number of crossings C and the longest
#' run L, then component (c,l) in the result is the
#' probability \eqn{P(C \ge c, L \le l)}.
#'
#' @param mtrx mpfr array
#' @return mpfr array
#' @examples
#' nill <- Rmpfr::mpfr(0, 120)
#' one <- Rmpfr::mpfr(1, 120)
#' two <- Rmpfr::mpfr(2, 120)
#' contents <- c(one,nill,nill, one,one,one, two,two,two)
#' mtrx3 <- Rmpfr::mpfr2array(contents, dim = c(3, 3))
#' print(mtrx3)
#' print(boxprobt(mtrx3))
#' @export
boxprobt <- function(mtrx) {
mtrs <- mtrx
rw <- nrow(mtrx)
cl <- ncol(mtrx)
for (row in 1:rw) mtrs[row, ] <- cumsum(mtrx[row, ])
for (col in 1:cl) mtrs[, col] <- rev(cumsum(rev(mtrs[, col])))
return(mtrs)
}
#' Exact Joint Probabilities for Low n
#' @description Exact joint probabilities, for low n,
#' of the number of crossings C and the longest run L
#' in n independent Bernoulli observations with success
#' probability p. Probabilites are multiplied by \eqn{2^{n-1}}.
#' @param n number, length of seqience, at most 6.
#' @param p success probability.
#' @param prec precision in mpfr calculations.
#' Default 120.
#' @return mpfr array
#' @examples
#' exactbin(n=6)
#' exactbin(n=5, p=0.6)
#' @export
exactbin <- function(n, p = 0.5, prec = 120) {
nill <- Rmpfr::mpfr(0, prec)
one <- Rmpfr::mpfr(1, prec)
two <- Rmpfr::mpfr(2, prec)
pm <- Rmpfr::mpfr(p, prec)
qm <- one - pm
res <- Rmpfr::mpfr2array(one, dim = c(1, 1))
if (n == 2)
res <- Rmpfr::mpfr2array(c(nill, two * pm * qm, pm^2 + qm^2, nill),
dim = c(2, 2))
if (n == 3)
res <- Rmpfr::mpfr2array(c(nill, nill, pm * qm, nill, two * pm * qm,
nill, pm^3 + qm^3, nill, nill), dim = c(3, 3))
if (n == 4)
res <- Rmpfr::mpfr2array(c(rep(nill, 3), 2 * pm^2 * qm^2, nill, 2 *
pm^2 * qm^2, two * pm * qm * (1 - pm * qm), nill, nill,
two * pm * qm * (pm^2 + qm^2), nill, nill, pm^4 + qm^4,
rep(nill, 3)), dim = c(4, 4))
if (n == 5)
res <- Rmpfr::mpfr2array(c(rep(nill, 4), pm^2 * qm^2, nill, nill, pm *
qm * (1 - pm * qm), 4 * pm^2 * qm^2, nill, nill, 2 * pm^2 *
qm^2, pm * qm * (2 * pm^3 + pm * qm + 2 * qm^3), nill, nill,
nill, 2 * pm * qm * (pm^3 + qm^3), rep(nill, 3), pm^5 +
qm^5, rep(nill, 4)), dim = c(5, 5))
if (n == 6)
res <- Rmpfr::mpfr2array(c(rep(nill, 5), 2 * pm^3 * qm^3, nill, nill,
pm^4 * qm^2 + pm^2 * qm^4, 2 * pm^4 * qm^2 + 8 * pm^3 *
qm^3 + 2 * pm^2 * qm^4, 3 * pm^4 * qm^2 + 4 * pm^3 *
qm^3 + 3 * pm^2 * qm^4, nill, nill, 2 * pm^3 * qm^3,
2 * pm^5 * qm + 2 * pm^4 * qm^2 + 4 * pm^3 * qm^3 + 2 *
pm^2 * qm^4 + 2 * pm * qm^5, 4 * pm^4 * qm^2 + 4 * pm^2 *
qm^4, nill, nill, nill, 2 * pm^4 * qm^2 + 2 * pm^2 *
qm^4, 2 * pm^5 * qm + pm^4 * qm^2 + pm^2 * qm^4 + 2 *
pm * qm^5, rep(nill, 3), nill, 2 * pm^5 * qm + 2 * pm *
qm^5, rep(nill, 4), pm^6 + qm^6, rep(nill, 5)), dim = c(6,
6))
rownames(res) <- c(0:(n - 1))
colnames(res) <- c(1:n)
return(res)
}
#' Number of Crossings and Longest Run
#' @description Auxiliary function for simclbin, computing
#' the number of crossings (type=0) or longest run (type=2)
#' in a sequence of independent normal observations. Crossings
#' and runs are related to whether the observations are above
#' a shift.
#' @param seri numeric; seri a sequence of random draws
#' @param shift numeric; shift for the observatoobs
#' @param type numeric; 0 number of crossings, 1 longest run
#' @return number of crossings or longest run, numeric
#' @export
clshift <- function(seri, shift = 0, type = 0) {
rle.sh <- rle(seri + shift > 0)$lengths
if (type == 0)
res <- length(rle.sh) - 1
if (type == 1)
res <- max(rle.sh)
return(res)
}
#' Simulation of Independent Bernoulli Observations
#' @description Simulation of a sequence of independent Bernoulli
#' Observations. To reduce the amount of random draws, each
#' simulation is based on a sequence of standard normal
#' variables, and whether each observation is above a shift
#' defined by the binomial probabilities assumed.
#' @param nser length of sequence simulated
#' @param nsim number of simulations
#' @param probs binomial probabilites
#' @return a data frame with the number of crossings and
#' longest run for each probability. For instance
#' the variables nc0.5 and lr0.5 are the number of
#' crossings and the longest run for success probability
#' 0.5. One row for each simulation.
#' @examples
#' cl30simbin <- simclbin(nser=30, nsim=100)
#' mean(cl30simbin$nc0.5) # mean number of crossings, p=0.5
#' mean(cl30simbin$lr0.9) # mean longest run, p=0.9
#' @export
simclbin <- function(nser = 100, nsim = 1e+05, probs = c(0.5, 0.6, 0.7,
0.8, 0.9)) {
nprob <- length(probs)
shifts <- stats::qnorm(probs)
series <- data.frame(matrix(stats::rnorm(nser * nsim), nrow = nsim))
res <- data.frame(matrix(rep(NA, 2 * nsim * nprob), nrow = nsim))
names(res) <- paste(rep(c("nc", "lr"), nprob), rep(probs, rep(2,
nprob)), sep = "")
for (shnr in 1:nprob) {
res[, 1 + 2 * (shnr - 1)] <- apply(series, 1, clshift, shift = shifts[shnr],
type = 0)
res[, 2 * shnr] <- apply(series, 1, clshift, shift = shifts[shnr],
type = 1)
}
return(res)
}
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.