R/beta.est.R

In MOM: Estimation and Testing of Hypothesis

#'Method of Moments Estimation of Beta distribution
#'@description function to get the method of moment estimate(s) of beta distribution
#'@param data A numeric vector.
#'@param unknown A character string specifying which parameter is (are) unknown to the user.
#'@param shape1,shape2	Non-negative parameters of the Beta distribution.
#'@param plot logical which controls whether the histogram of the data along with the density curve of the theoretical beta distribution with the estimated parameters.
#'@param curvecol color of the theoretical density curve
#'@param ... additional plotting parameters
#'@importFrom stats var dbeta
#'@importFrom graphics hist curve
#'@return the estimated parameters by the method of moments of the data assuming the underlying distribution is beta distribution
#'@examples
#'beta_est(rbeta(1000,shape1=2,shape2=1),unknown="shape2",shape1=2)#shape1 is known
#'beta_est(rbeta(1000,shape1=2,shape2=1),unknown="shape1",shape2=1)#shape2 is known
#'beta_est(rbeta(1000,shape1=2,shape2=1),unknown="both")#both will be estimated
#'@export
beta_est=function(data,unknown=c("shape1","shape2","both"),shape1=NULL,shape2=NULL,plot=TRUE,curvecol="red",...)
{
  stopifnot(shape1>0&shape2>0)
  avg=mean(data)
  va.r=var(data)
  if(unknown=="shape1") shape1=((shape2*avg)/(1-avg))
  if(unknown=="shape2") shape2=((shape1*(1-avg))/avg)
  if(unknown=="both")
  {
    shape1=(avg^2*(1-avg)/va.r)-avg
    shape2=((shape1*(1-avg))/avg)
  }
  if(plot==TRUE)
  {
    hist(data,prob=TRUE,...)
    curve(dbeta(x,shape1=shape1,shape2=shape2),add=T,col=curvecol)
  }
  output=list(shape1,shape2)
  names(output)=c("shape1","shape2")
  return(output)
}