R/simmptest.R

In MOM: Estimation and Testing of Hypothesis

#'Most Powerful Test by Neyman-Pearson Lemma
#'@description It can be used to check whether a data comes from null distribution or from the alternative distribution
#'@param data A numeric vector
#'@param null.dist The family of null distribution
#'@param null.par The parameter values of the null distribution
#'@param alter.dist The family of alternative distribution
#'@param alter.par The parameter values of the alternative distribution
#'@param test.level The level of significance of the test
#'@param sim.size simulation size, increasing it will gives more accuracy.
#'@param power A loogical vector, whether power of the test will be calculated.
#'@importFrom stats quantile runif rnorm rlnorm rgamma rcauchy rweibull rbeta rbinom rpois rnbinom rgeom rchisq rt rf dunif dnorm dlnorm dgamma dcauchy dweibull dbeta dbinom dpois dnbinom dgeom dchisq dt df sd
#'@importFrom actuar rpareto dpareto rlogarithmic dlogarithmic
#'@importFrom utils head
#'@importFrom VGAM rrayleigh rlaplace drayleigh dlaplace
#'@details
#'This function mainly test whether data comes from the null distribution or alternative distribution.
#'It uses the theory of the Most Powerful (MP) test. It basically uses simulation to get the p value and make the decision.
#'Increasing sim.size give more accuracy as well as test can be failed if you increase it heavily.
#'@return A list of class \code{"momtest"} will be returned having the following components:
#'\describe{
#'  \item{Method}{The Method's Name}
#'  \item{Data}{The first 6 elements of input data}
#'  \item{Null.Distrbution}{The family of null distribution}
#'  \item{Null.Parameter}{The parameter values of the null distribution}
#'  \item{Alternative.Distrbution}{The family of alternative distribution}
#'  \item{Alternative.Parameter}{The parameter values of the alternative distribution}
#'  \item{Sample.Size}{The sample size}
#'  \item{Significance.level}{The level of the significance of the test}
#'  \item{Decision}{The Test Result, wheter the null hypotheis is rejected or not}
#'  \item{Power}{Power of the Test}
#'}
#'@examples
#'sim_mp_test(rnorm(100),null.dist="normal",null.par=c(0,1),alter.dist="cauchy",alter.par=c(0,1))
#'sim_mp_test(rnorm(100),null.dist="nor",null.par=c(2,1),alter.dist="nor",alter.par=c(0,1))
#'@export
sim_mp_test=function(data, null.dist=c("uniform","normal","lognormal", "gamma", "cauchy",
                               "pareto", "weibull", "rayleigh", "laplace", "beta", "binomial",
                               "poisson", "negativebinomial","geometric","t","f","logarithmic"),
                   null.par,
                   alter.dist=c("uniform","normal","lognormal", "gamma", "cauchy",
                               "pareto", "weibull", "rayleigh", "laplace", "beta", "binomial",
                               "poisson", "negativebinomial","geometric","t","f","logarithmic"),
                   alter.par,test.level=0.95,sim.size=1,power=TRUE)
{
  null.dist=match.arg(null.dist)
  alter.dist=match.arg(alter.dist)
  np1=null.par[1]
  np2=null.par[2]
  np3=null.par[3]
  ap1=alter.par[1]
  ap2=alter.par[2]
  ap3=alter.par[3]
  if(is.na(null.par[2])) np2=0
  if(is.na(null.par[3])) np3=0
  if(is.na(alter.par[2])) ap2=0
  if(is.na(alter.par[3])) ap3=0
  n=length(data)
  if(null.dist=="uniform")
  {
    null=prod(dunif(data,min=np1,max=np2))
    sim=matrix(runif(n*sim.size,min=np1,max=np2),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(dunif(x,min=np1,max=np2)))
  }
  else if(null.dist=="normal")
    {
    null=prod(dnorm(data,mean=np1,sd=np2))
    sim=matrix(rnorm(n*sim.size,mean=np1,sd=np2),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(dnorm(x,mean=np1,sd=np2)))
    }
  else if(null.dist=="lognormal")
    {
    null=prod(dlnorm(data,meanlog=np1,sdlog=np2))
    sim=matrix(rlnorm(n*sim.size,meanlog=np1,sdlog=np2),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(dlnorm(x,meanlog=np1,sdlog=np2)))
    }
  else if(null.dist=="gamma")
    {
    null=prod(dgamma(data,shape=np1,rate=np2))
    sim=matrix(rgamma(n*sim.size,shape=np1,rate=np2),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(dgamma(x,shape=np1,rate=np2)))
    }
  else if(null.dist=="cauchy")
    {
    null=prod(dcauchy(data,location=np1,scale=np2))
    sim=matrix(rcauchy(n*sim.size,location=np1,scale=np2),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(dcauchy(x,location=np1,scale=np2)))
    }
  else if(null.dist=="pareto")
    {
    null=prod(dpareto(data,shape=np1,scale=np2))
    sim=matrix(rpareto(n*sim.size,shape=np1,scale=np2),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(dpareto(x,shape=np1,scale=np2)))
    }
  else if(null.dist=="weibull")
    {
    null=prod(dweibull(data,shape=np1,scale=np2))
    sim=matrix(rweibull(n*sim.size,shape=np1,scale=np2),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(dweibull(x,shape=np1,scale=np2)))
    }
  else if(null.dist=="rayleigh")
    {
    null=prod(drayleigh(data,scale=np1))
    sim=matrix(rrayleigh(n*sim.size,scale=np1),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(drayleigh(x,scale=np1)))
    }
  else if(null.dist=="laplace")
    {
    null=prod(dlaplace(data,location=np1,scale=np2))
    sim=matrix(rlaplace(n*sim.size,location=np1,scale=np2),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(dlaplace(x,location=np1,scale=np2)))
    }
  else if(null.dist=="beta")
    {
    null=prod(dbeta(data,shape1=np1,shape2=np2,ncp=np3))
    sim=matrix(rbeta(n*sim.size,shape1=np1,shape2=np2,ncp=np3),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(dbeta(x,shape1=np1,shape2=np2,ncp=np3)))
    }
  else if(null.dist=="t")
    {
    null=prod(dt(data,df=np1,ncp=np2))
    sim=matrix(rt(n*sim.size,df=np1,ncp=np2),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(dt(x,df=np1,ncp=np2)))
    }
  else if(null.dist=="f")
    {
    null=prod(df(data,df1=np1,df2=np2,ncp=np3))
    sim=matrix(rf(n*sim.size,df1=np1,df2=np2,ncp=np3),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(df(x,df1=np1,df2=np2,ncp=np3)))
    }
  else if(null.dist=="binomial")
    {
    null=prod(dbinom(data,size=np1,prob=np2))
    sim=matrix(rbinom(n*sim.size,size=np1,prob=np2),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(dbinom(x,size=np1,prob=np2)))
    }
  else if(null.dist=="poisson")
    {
    null=prod(dpois(data,lambda=np1))
    sim=matrix(rpois(n*sim.size,lambda=np1),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(dpois(x,lambda=np1)))
    }
  else if(null.dist=="geometric")
    {
    null=prod(dgeom(data,prob=np1))
    sim=matrix(rgeom(n*sim.size,prob=np1),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(dgeom(x,prob=np1)))
    }
  else if(null.dist=="negativebinomial")
    {
    null=prod(dnbinom(data,size=np1,prob=np2))
    sim=matrix(rnbinom(n*sim.size,size=np1,prob=np2),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(dnbinom(x,size=np1,prob=np2)))
    }
  else if(null.dist=="logarithmic")
    {
    null=prod(dlogarithmic(data,prob=np1))
    sim=matrix(rlogarithmic(n*sim.size,prob=np1),nrow=n,byrow=T)
    crvald=apply(sim,MARGIN=2,FUN=function(x) prod(dlogarithmic(x,prob=np1)))
  }
  if(alter.dist=="uniform")
  {
    alter=prod(dunif(data,min=ap1,max=ap2))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(dunif(x,min=ap1,max=ap2)))
    if(power==TRUE)
    {
      simpow=matrix(runif(n*sim.size,min=ap1,max=ap2),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(dunif(x,min=ap1,max=ap2)))
    }
  }
  else if(alter.dist=="normal")
    {
    alter=prod(dnorm(data,mean=ap1,sd=ap2))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(dnorm(x,mean=ap1,sd=ap2)))
    if(power==TRUE)
      {
      simpow=matrix(rnorm(n*sim.size,mean=ap1,sd=ap2),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(dnorm(x,mean=ap1,sd=ap2)))
      }
    }
  else if(alter.dist=="lognormal")
    {
    alter=prod(dlnorm(data,meanlog=ap1,sdlog=ap2))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(dlnorm(x,meanlog=ap1,sdlog=ap2)))
    if(power==TRUE)
      {
      simpow=matrix(rlnorm(n*sim.size,meanlog=ap1,sdlog=ap2),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(dlnorm(x,meanlog=ap1,sdlog=ap2)))
      }
    }
  else if(alter.dist=="gamma")
    {
    alter=prod(dgamma(data,shape=ap1,rate=ap2))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(dgamma(x,shape=ap1,rate=ap2)))
    if(power==TRUE)
    {
      simpow=matrix(rgamma(n*sim.size,shape=ap1,rate=ap2),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(dgamma(x,shape=ap1,rate=ap2)))
    }
    }
  else if(alter.dist=="cauchy")
    {
    alter=prod(dcauchy(data,location=ap1,scale=ap2))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(dcauchy(x,location=ap1,scale=ap2)))
    if(power==TRUE)
    {
      simpow=matrix(rcauchy(n*sim.size,location=ap1,scale=ap2),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(dcauchy(x,location=ap1,scale=ap2)))
    }
    }
  else if(alter.dist=="pareto")
    {
    alter=prod(dpareto(data,shape=ap1,scale=ap2))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(dpareto(x,shape=ap1,scale=ap2)))
    if(power==TRUE)
    {
      simpow=matrix(rpareto(n*sim.size,shape=ap1,scale=ap2),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(dpareto(x,shape=ap1,scale=ap2)))
    }
    }
  else if(alter.dist=="weibull")
    {
    alter=prod(dweibull(data,shape=ap1,scale=ap2))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(dweibull(x,shape=ap1,scale=ap2)))
    if(power==TRUE)
    {
      simpow=matrix(rweibull(n*sim.size,shape=ap1,scale=ap2),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(dweibull(x,shape=ap1,scale=ap2)))
    }
    }
  else if(alter.dist=="rayleigh")
    {
    alter=prod(drayleigh(data,scale=ap1))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(drayleigh(x,scale=ap1)))
    if(power==TRUE)
    {
      simpow=matrix(rrayleigh(n*sim.size,scale=ap2),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(drayleigh(x,scale=ap2)))
    }
    }
  else if(alter.dist=="laplace")
    {
    alter=prod(dlaplace(data,location=ap1,scale=ap2))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(dlaplace(x,location=ap1,scale=ap2)))
    if(power==TRUE)
    {
      simpow=matrix(rlaplace(n*sim.size,location=ap1,scale=ap2),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(dlaplace(x,location=ap1,scale=ap2)))
    }
    }
  else if(alter.dist=="beta")
    {
    alter=prod(dbeta(data,shape1=ap1,shape2=ap2,ncp=ap3))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(dbeta(x,shape1=ap1,shape2=ap2,ncp=ap3)))
    if(power==TRUE)
    {
      simpow=matrix(rbeta(n*sim.size,shape1=ap1,shape2=ap2,ncp=ap3),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(dbeta(x,shape1=ap1,shape2=ap2,ncp=ap3)))
    }
    }
  else if(alter.dist=="t")
    {
    alter=prod(dt(data,df=ap1,ncp=ap2))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(dt(x,df=ap1,ncp=ap2)))
    if(power==TRUE)
    {
      simpow=matrix(rt(n*sim.size,df=ap1,ncp=ap2),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(dt(x,df=ap1,ncp=ap2)))
    }
    }
  else if(alter.dist=="f")
    {
    alter=prod(df(data,df1=ap1,df2=ap2,ncp=ap3))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(df(x,df1=ap1,df2=ap2,ncp=ap3)))
    if(power==TRUE)
    {
      simpow=matrix(rf(n*sim.size,df1=ap1,df2=ap2,ncp=ap3),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(df(x,df1=ap1,df2=ap2,ncp=ap3)))
    }
    }
  else if(alter.dist=="binomial")
    {
    alter=prod(dbinom(data,size=ap1,prob=ap2))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(dbinom(x,size=ap1,prob=ap2)))
    if(power==TRUE)
    {
      simpow=matrix(rbinom(n*sim.size,size=ap1,prob=ap2),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(dbinom(x,size=ap1,prob=ap2)))
    }
    }
  else if(alter.dist=="poisson")
    {
    alter=prod(dpois(data,lambda=ap1))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(dpois(x,lambda=ap1)))
    if(power==TRUE)
    {
      simpow=matrix(rpois(n*sim.size,lambda=ap1),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(dpois(x,lambda=ap1)))
    }
    }
  else if(alter.dist=="geometric")
    {
    alter=prod(dgeom(data,prob=ap1))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(dgeom(x,prob=ap1)))
    if(power==TRUE)
    {
      simpow=matrix(rgeom(n*sim.size,prob=ap1),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(dgeom(x,prob=ap1)))
    }
    }
  else if(alter.dist=="negativebinomial")
    {
    alter=prod(dnbinom(data,size=ap1,prob=ap2))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(dnbinom(x,size=ap1,prob=ap2)))
    if(power==TRUE)
    {
      simpow=matrix(rnbinom(n*sim.size,size=ap1,prob=ap2),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(dnbinom(x,size=ap1,prob=ap2)))
    }
    }
  else if(alter.dist=="logarithmic")
    {
    alter=prod(dlogarithmic(data,prob=ap1))
    crvaln=apply(sim,MARGIN=2,FUN=function(x) prod(dlogarithmic(x,prob=ap1)))
    if(power==TRUE)
    {
      simpow=matrix(rlogarithmic(n*sim.size,prob=ap1),nrow=n,byrow=T)
      pown=apply(simpow,MARGIN=2,FUN=function(x) prod(dlogarithmic(x,prob=ap1)))
    }
  }
  if(power==TRUE)
  {
    if(null.dist=="uniform") powd=apply(simpow,MARGIN=2,FUN=function(x) prod(dunif(x,min=np1,max=np2)))
    else if(null.dist=="normal") powd=apply(simpow,MARGIN=2,FUN=function(x) prod(dnorm(x,mean=np1,sd=np2)))
    else if(null.dist=="lognormal") powd=apply(simpow,MARGIN=2,FUN=function(x) prod(dlnorm(x,meanlog=np1,sdlog=np2)))
    else if(null.dist=="gamma") powd=apply(simpow,MARGIN=2,FUN=function(x) prod(dgamma(x,shape=ap1,rate=ap2)))
    else if(null.dist=="cauchy")powd=apply(simpow,MARGIN=2,FUN=function(x) prod(dcauchy(x,location=ap1,scale=ap2)))
    else if(null.dist=="logarithmic")powd=apply(simpow,MARGIN=2,FUN=function(x) prod(dlogarithmic(x,prob=ap1)))
    else if(null.dist=="negativebinomial")powd=apply(simpow,MARGIN=2,FUN=function(x) prod(dnbinom(x,size=ap1,prob=ap2)))
    else if(null.dist=="geometric")powd=apply(simpow,MARGIN=2,FUN=function(x) prod(dgeom(x,prob=ap1)))
    else if(null.dist=="poisson")powd=apply(simpow,MARGIN=2,FUN=function(x) prod(dpois(x,lambda=ap1)))
    else if(null.dist=="binomial")powd=apply(simpow,MARGIN=2,FUN=function(x) prod(dbinom(x,size=ap1,prob=ap2)))
    else if(null.dist=="f")powd=apply(simpow,MARGIN=2,FUN=function(x) prod(df(x,df1=ap1,df2=ap2,ncp=ap3)))
    else if(null.dist=="t")powd=apply(simpow,MARGIN=2,FUN=function(x) prod(dt(x,df=ap1,ncp=ap2)))
    else if(null.dist=="beta")powd=apply(simpow,MARGIN=2,FUN=function(x) prod(dbeta(x,shape1=ap1,shape2=ap2,ncp=ap3)))
    else if(null.dist=="laplace")powd=apply(simpow,MARGIN=2,FUN=function(x) prod(dlaplace(x,location=ap1,scale=ap2)))
    else if(null.dist=="rayleigh")powd=apply(simpow,MARGIN=2,FUN=function(x) prod(drayleigh(x,scale=ap2)))
    else if(null.dist=="weibull")powd=apply(simpow,MARGIN=2,FUN=function(x) prod(dweibull(x,shape=ap1,scale=ap2)))
    else if(null.dist=="pareto")powd=apply(simpow,MARGIN=2,FUN=function(x) prod(dpareto(x,shape=ap1,scale=ap2)))
  }

  crdata=alter/null
  crval=crvaln/crvald
  cutoff=quantile(crval,prob=test.level)
  if(crdata>cutoff)
    {
      d="Reject null hypothesis"
  }else
  d= "Do not reject null hypothesis"
  power.val="NA"
  if(power==TRUE)
  {
    powsim=pown/powd
    count=sum(powsim>cutoff)
    power.val=count/length(powsim)
  }
  output=momtest(data,null.dist,null.par,alter.dist,alter.par,n,test.level,d,power.val)
  return(output)
}