calculateConstantCVGenDirichletK.RdCalculate Generalized Dirichlet parameters such that the margins have constant Coefficient of Variation
calculateConstantCVGenDirichletK(p, most.likely.k)
| p | probabilities |
|---|---|
| most.likely.k | The k of the most likely branch or dirichlet marginal |
the k vector such that all margins have a constant coefficient of variation
calculateConstantCVGenDirichletK#> function(p, most.likely.k) #> { #> lenp <- length(p) #> stopifnot(lenp > 2) #> #> # if the p's are not probabilities, change them into probabilities #> if (sum(p) != 1) #> { #> warning("Renormalizing p in calculateConstantCVGenDiricletK") #> p <- p/sum(p) #> } #> #> # if the p's are not in order, order them #> ind <- order(p, decreasing = TRUE) #> p <- p[ind] #> #> # calculate CV of most likely branch #> varMostLikely <- p[1] * (1 - p[1]) / (most.likely.k + 1) #> CVMostLikely <- sqrt(varMostLikely) / p[1] #> #> k <- rep(1, lenp) #> k[1] <- most.likely.k #> theta <- p / (1 + p - cumsum(p)) #> tempVector <- (k * (1 - theta) + 1) / (k + 1) #> #> for (i in 2:(lenp - 1)) #> { #> # must re-calculate temp based on new k's #> tempVector <- (k * (1 - theta) + 1) / (k + 1) #> temp <- prod(tempVector[1:(i - 1)]) #> k[i] <- (p[i] * (CVMostLikely^2 + 1) - temp) / #> (temp * theta[i] - p[i] * (CVMostLikely^2 + 1)) #> if (k[i] < 0 || !is.finite(k[i])) #> { #> warning("Negative or Infinite k at index ", i, " k = ", k) #> # in this situation, set the k to the largest machine number #> k[i] <- .Machine$double.xmax #> } #> } #> k[lenp] <- k[lenp - 1] #> #> ktemp <- numeric(length(k)) #> ktemp[ind] <- k #> ktemp #> } #> <environment: namespace:dirichlet>