R/simsamint.R

In MOM: Estimation and Testing of Hypothesis

#'Simulated confidence interval of a statistic
#'@description A function that returns a sampling interval for a statistic formed from random sample of certain probability distributions. The function generates the confidence interval using Monte Carlo simulations. The results might be unreliable if the resulting statistic has fat tailed distribution.
#'@param dist The parent population distribution
#'@param pop.par The value of the population parameters
#'@param FUN The statistic as a function of random data
#'@param side The type of the confidence interval (both sided, only lower bound or only upper bound)
#'@param range It controls the length of the interval in which the boundary points are searched for. One may increase the range in case the distribution of statistic is suspected to be fat tailed.
#'@param conf.coeff The confidence coefficient of the sampling interval
#'@param n sample size
#'@param sim.size simulation size, increasing it will gives more accuracy.
#'@importFrom stats runif rnorm rlnorm rgamma rcauchy rweibull rbeta rbinom rpois rnbinom rgeom rchisq rt rf sd
#'@importFrom actuar rpareto rlogarithmic
#'@importFrom utils head
#'@importFrom VGAM rrayleigh rlaplace
#'@details
#'The function asks the user to specify a distribution from which random sample is drawn and to specify a function of the random variables for which an approximate sampling Interval is to be provided. The function then uses Monte Carlo simulation technique to provide an approximate sampling interval of the statistic.
#'This function might be useful when the sampling distribution for a particular statistic is unknown, but that statistic might be useful in drawing meaningful inference. Although this function is inferior to other sophisticated techniques to deal with this problem, it might come handy for a beginner.
#'@return A list of class \code{"momint"} will be returned having the following components:
#'\describe{
#'  \item{Method}{The Method's Name}
#'  \item{Population.Distrbution}{The family of population distribution}
#'  \item{Paramater}{The parameter values of the population distribution}
#'  \item{Statistic}{The function of which the interval will be provided}
#'  \item{Sample.Size}{The sample size}
#'  \item{Confidence.Coefficient}{The confidence coefficient of the sampling interval}
#'  \item{Sampling.Interval}{The estimated sampling interval}
#'}
#'@examples
#'sim_sam_int(dist="normal",pop.par=c(0,1),FUN=mean,side="both")
#'sim_sam_int(dist="binomial",pop.par=c(5,0.5),FUN=sum,side="lower")
#'@export
sim_sam_int=function(dist=c("normal","lognormal", "gamma", "chisquare", "cauchy",
				"pareto", "weibull", "rayleigh", "laplace", "beta", "binomial",
				"poisson", "negativebinomial", "geometric","t","f","uniform"),
		     pop.par,FUN,side=c("lower","upper","both"),conf.coeff=0.95,range=1,n=100,sim.size=1000)
{
  dist=match.arg(dist)
  side=match.arg(side)
	p1=pop.par[1]
	p2=pop.par[2]
  FUN=match.fun(FUN)
  if(is.na(pop.par[2])) p2=0
  if(is.na(pop.par[3])) p3=0
 	sample.dist=function(n)
	{
		if(dist=="normal") s=rnorm(n,p1,p2)
		else if(dist=="lognormal") s=rlnorm(n,p1,p2)
		else if(dist=="gamma") s=rgamma(n,p1,p2)
		else if(dist=="chisquare") s=rchisq(n,p1)
		else if(dist=="cauchy") s=rcauchy(n,p1,p2)
		else if(dist=="pareto") s=rpareto(n,p1,p2)
		else if(dist=="weibull") s=rweibull(n,p1,p2)
		else if(dist=="rayleigh") s=rrayleigh(n,p1)
		else if(dist=="laplace") s=rlaplace(n,p1,p2)
		else if(dist=="beta") s=rbeta(n,p1,p2)
		else if(dist=="binomial") s=rbinom(n,p1,p2)
		else if(dist=="poisson") s=rpois(n,p1)
		else if(dist=="negativebinomial") s=rnbinom(n,p1,p2)
		else if(dist=="geometric") s=rgeom(n,p1)
		else if(dist=="t") s=rt(n,p1,p2)
		else if(dist=="f") s=rf(n,p1,p2,p3)
		else if(dist=="uniform") s=runif(n,min=p1,max=p2)
		return(s)
	}
	s=sample.dist(n*sim.size)
	samples=matrix(s,nrow=sim.size)
	theta1=apply(samples,1,FUN)
	theta= (theta1-mean(theta1))/sd(theta1)
	prob=function(c)
      {
		if(side=="upper") p= mean(theta<=c)
		else if(side=="lower") p=mean(c<=theta)
		return(p)
     	}
	prob1=function(c1,c2)
	{
		p=mean(c1<=theta&theta<=c2)
		return(p)
	}
	conf.int= function()
	{

	    c.pts=seq(min(theta)-2*range,max(theta)+2*range,length=500)
		  x1=NULL
		  x2=NULL
		  for(i in 1:500)
		  {
		    for(j in i:500)
		    {
		      p=prob1(c.pts[i],c.pts[j])
		      if(p<conf.coeff+0.01 & p>conf.coeff-0.01)
		      {
		        x1=c(x1,c.pts[i])
		        x2=c(x2,c.pts[j])
		      }
		    }
		  }
		  C=x2-x1
		  m=which.min(C)
      return(c(x1[m]*sd(theta1)+mean(theta1),x2[m]*sd(theta1)+mean(theta1)))

	}
	interval=function(n)
	{
		if(side=="upper") c.pts=seq(min(theta),max(theta)+4*range,length=500)
		else if(side=="lower")c.pts=seq(min(theta)-4*range,max(theta),length=500)
		if(side!="both")
		{
			c1=matrix(c.pts)
			x=apply(c1,1,prob)
			pivot=min(c.pts[x>conf.coeff-0.01 & x<conf.coeff+0.01])
		}
		if(side=="both") interval=conf.int()
		else if(side=="upper") interval= c(-Inf,pivot*sd(theta1)+mean(theta1))
		else interval= c(pivot*sd(theta1)+mean(theta1),Inf)
		return(interval)
	}
	c=interval(n)
	output=momint(dist=dist,pop.par,FUN,n,conf.coeff,c)
	return(output)
}