Description Usage Arguments Value Author(s) See Also Examples
This function finds the two shape parameters for the Beta distribution of a random variable between 0 and 1. Note that the function uses a different method depending on the arguments supplied. The three methods are:
When the mean (mu
) and variance (sigma
) are supplied, the solution is found analytically.
When observed probabilities (probs
) at each quantile (quantiles
) are given, the
solution is found by minimizing the Sum of the Squared Errors (SSE) using the Nelder-Mead optimization
algorithm. Note that the fitting algorithm performs best when the five standard quantiles are supplied (Min, 25th, Median, 75th, Max).
When only observed probabilities (probs
) are supplied, the function uses Maximum Likelihood Estimation (MLE).
1 |
mu |
scalar giving the mean μ |
sigma |
scalar giving the variance σ^2 |
quantiles |
vector of quantiles for which proportions are observed. Expects: c('min', '25th', 'median', '75th', 'max'). |
probs |
vector of observed probabilities or proportions |
A list containing the two shape parameters of the Beta distribution
John Giles
Other beta_params:
get_beta_params_age()
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 | ### Analytic solution
mu <- 0.6
sigma <- 0.005
prm <- get_beta_params(mu=mu, sigma=sigma)
curve(dbeta(x, prm$shape1, prm$shape2), 0, 1, lwd=2,
xlab="Random variable", ylab='Probability density')
abline(v=mu, col='blue')
abline(v=c(mu-sqrt(sigma), mu+sqrt(sigma)), col='blue', lty=2)
### Numeric solution using SSE
quantiles <- c(0.001, 0.25, 0.5, 0.75, 0.999)
probs1 <- c(0.26, 0.57, 0.73, 0.81, 0.95) # WHO region AFR, <12mos, n=16
probs2 <- c(0.86, 0.86, 0.88, 0.92, 0.94) # WHO region AFR, >12mos, n=4
v <- rbind(probs1, probs2)
dimnames(v) <- list(Age=c('< 12mos', '> 12mos'),
Quantiles=c('min', '25th', 'median', '75th', 'max'))
par(mfrow=c(1,2))
prm <- get_beta_params(quantiles=quantiles, probs=probs1)
curve(dbeta(x, prm$shape1, prm$shape2), 0, 1, lwd=2,
xlab="Vaccine effectiveness", ylab='Probability density')
abline(v=probs1, col='blue', lty=2)
mtext(c('min', '25th', 'median', '75th', 'max'), 3, at=probs1, line=0.5, col='blue', cex=0.75)
mtext('<12mos', 3, at=0.05, line=-1.5)
prm <- get_beta_params(quantiles=quantiles, probs=probs2)
curve(dbeta(x, prm$shape1, prm$shape2), 0, 1, lwd=2,
xlab="Vaccine effectiveness", ylab='Probability density')
abline(v=probs2, col='blue', lty=2)
mtext('>12mos', 3, at=0.05, line=-1.5)
### Maximum Likelihood Estimation
get_beta_params(probs=rbeta(100, 1.5, 2))
# Measles vaccination coverage in Namibia
# WHO-UNICEF estiamtes from:
# http://apps.who.int/immunization_monitoring/globalsummary/estimates?c=NAM
mcv1 <- c(82, 80, 75, 85, 83, 82, 76, 74,
75, 76, 73, 69, 63, 73, 70, 70,
68, 58, 69, 65, 64, 59, 68)/100
mcv2 <- c(50, 32)/100
prm <- get_beta_params(probs=mcv1)
curve(dbeta(x, prm$shape1, prm$shape2), 0, 1)
abline(v=mcv1, lty=2, col='blue')
prm <- get_beta_params(probs=mcv2)
curve(dbeta(x, prm$shape1, prm$shape2), 0, 1)
abline(v=mcv2, lty=2, col='blue')
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.