In purplelllll/StatComp21002: Several functions to estimate parameters and compute the distribution function

Overview

StatComp21002 is a simple R package developed to using admm method to estimate some parameters and estimate a parameter in a distribution and compute the distribution function of cauchy distribution. Three functions are considered, namely, df (compute the distribution function of cauchy) and est (estimate a parameter in a distribution and compute its bias and standard error) and admm(estimate some parameters in a regression model by the admm method). For each function,R versions are produced.

The R code for benchmarking est and df and admm is as follows.

est=function(sigma,n)
{
  ##generate random sample from the distribution using inverse transform
  u <- runif(n)
  x <- sqrt((-2)*((sigma)^2)*log(1-u))

  ##bootstrap
  B=10
  sigma_=numeric(B)
  for (b in 1:B)
  {
    i=sample(x,size=n,replace=TRUE)
    s1=sum(log(i))
    s2=sum(i^2)

    ##the function that needs to optimize
    f0=function(sigma)
      s1-2*n*log(sigma)-s2/(2*(sigma^2))
    ##the MLE
    sigma_[b]=optimize(f0,lower = 0,upper=sigma+4,maximum = TRUE)$maximum
  } 
  estimator=mean(sigma_)
  ##the bias and se of bootstrap
  bias=mean(sigma_)-sigma
  se=sd(sigma_)
  val=numeric(3)
  val=c(estimator,bias,se)
  return(val)
}

df=function(x,theta,eta)
{
  set.seed(34567)
  ## the density function
  f0=function(x,theta,eta)
  {
    a=((x-eta)/theta)^2
    b=theta*pi*(1+a)
    c=1/b
    return (c)
  }
  ## compute one part of the df
  s1=integrate(f0,lower=-Inf,upper=0,rel.tol = .Machine$double.eps^0.25,theta=theta,eta=eta)$value
  ## compute another part of the df
  u=runif(1e6)
  s2=mean(x/(theta*pi*(1+((u*x-eta)/theta)^2)))

  ##the df
  if (x>0)  
    s=s1+s2
  else
    s=integrate(f0,lower=-Inf,upper=x,rel.tol = .Machine$double.eps^0.25,theta=theta,eta=eta)$value
  return(c(s,pcauchy(x,eta,theta)))
}

admm=function(beta0,beta1,beta2,m,n)
{
  set.seed(78910)
  ##generate the sample
  X1=matrix(0,m,n)
  for (i in 1:m)
    for (j in 1:n)
      X1[i,j]=rnorm(1)

  X2=matrix(0,m,n)
  for (i in 1:m)
    for ( j in 1:n)
      X2[i,j]=rnorm(1)

  sigma=matrix(0,m,n)
  for (i in 1:m)
    for (j in 1:n)
      sigma[i,j]=rnorm(1,0,sd=0.8)

  #compute the yij
  y=matrix(0,m,n)
  for (i in 1:m)
    for (j in 1:n)
      y[i,j]=beta0+X1[i,j]*beta1+X2[i,j]*beta2+sigma[i,j]

  #set the initial value of the parameters
  beta0_=beta0
  beta1_=beta1
  beta2_=beta2

  theta_=matrix(0,m,n)
  for (i in 1:m)
    for ( j in 1:n)
      theta_[i,j]=rnorm(1)

  xi=1
  tau=0.5

  ##Aij Zij
  A=matrix(0,m,n)
  Z=matrix(0,m,n)

  N=1000
  tol <-1e-3
  r=0
  for (l in 1:N){
    for (i in 1:m)
    {
      for (j in 1:n)
      {
        A[i,j]=y[i,j]-(beta0_+X1[i,j]*beta1_+X2[i,j]*beta2_)
        if ((A[i,j]-theta_[i,j]/xi)>tau/xi)
          Z[i,j]=A[i,j]-tau/xi-theta_[i,j]/xi
        if ((A[i,j]-theta_[i,j]/xi)<=(tau-1)/xi)
          Z[i,j]=A[i,j]+(1-tau)/xi-theta_[i,j]/xi 
      }
    }

    ##iteration
    s1=sum(Z-y+X1*beta1_+X2*beta2_)
    s=sum(theta_)
    beta0_1=-(s1*xi+s)/(xi*m*n)
    e1=abs(beta0_1-beta0_)
    beta0_=beta0_1
    s2=sum((Z-y+beta0_+X2*beta2_)*X1)
    s20=sum(theta_*X1)
    beta1_1=-(s2*xi+s20)/(xi*sum(X1^2))
    e2=abs(beta1_1-beta1_)
    beta1_=beta1_1
    s3=sum((Z-y+beta0_+X1*beta1_)*X2)
    s30=sum(theta_*X2)
    beta2_1=-(s3*xi+s30)/(xi*sum(X2^2))
    e3=abs(beta2_1-beta2_)
    beta2_=beta2_1
    theta_1=theta_+xi*(Z-(y-(beta0_+X1*beta1_+X2*beta2_)))
    theta_=theta_1

    r=r+1
    if ((e1<tol)&&(e2<tol)&&(e3<tol))  break  
  }
  parameters=c(beta0,beta1,beta2)
  estimators=c(beta0_,beta1_,beta2_)
  r1=data.frame(parameters,estimators)

  return( list( (knitr::kable(t(r1),align="c")),r ) )
}

some examples for the functions

est(1,1e4)
df(1,1,0)
admm(1.5,1,0,20,5)

purplelllll/StatComp21002 documentation built on Dec. 24, 2021, 1:19 a.m.