inst/doc/Introduction.R

In zhangmanustc/StatComp18024: Three functions for StatComp

## ----eval=FALSE----------------------------------------------------------
#  Antithetic.method<- function(sigma, n) {#construct funtions to generate samples from  Rayleigh distribution
#    X <- X_ <- numeric(n)
#    for (i in 1:n) {
#      U<- runif(n)#generate n observations from U~U(0,1)
#      V<- 1-U#then V~U(0,1)
#      X<- sigma*sqrt(-2*log(V))
#      X_<-sigma*sqrt(-2*log(U))
#      #according inverse transform from annalysis above,we obtain that $X=F_{X}^{-1}(U)=\sqrt{-2ln(1-U)}=\sqrt{-2ln(V)}$\qquad
#      var1<-var(X)#if X1 and X2 are independent,then var((X1+X2)/2)=var(X)
#      var2<-(var(X)+var(X_)+2*cov(X,X_))/4#compute the variance of (X+X')/2 generated by antithetic variables
#      reduction <-((var1-var2)/var1)#compute the reduction in variance by using antithetic variables
#    }
#    percent_reduction<-paste0(format(100*reduction, format = "fg", digits = 4), "%")#set format type by function format
#    return(percent_reduction)
#  }
#  set.seed(123)
#  sigma = 1#set the value of sigma-the scale parameter of Rayleigh distribution to 1
#  n <- 1000
#  Antithetic.method(sigma, n)

## ------------------------------------------------------------------------
Betacdf <- function( x, alpha=3, beta=3) {#construct cdf of Beta(3,3) by function()
m<-1e3
if ( any(x < 0) ) return (0)
stopifnot( x < 1 )#the range of x is between 0 and 1
t <- runif( m, min=0, max=x )#generate samples from U(0,x)
p<-(x-0)*(1/beta(alpha,beta)) * t^(alpha-1) * (1-t)^(beta-1)
cdf<-mean(p)#compute Monte Carlo estimate of cdf
return( min(1,cdf) )#ensure that cdf<=1
} #end function

set.seed(123)
for (i in 1:9) {
estimate<-Betacdf(i/10,3,3)#use function constructed to estimate F(x) when x=0.1,0.2,...,0.9
returnedvalues<-pbeta(i/10,3,3)#obtain values returned by the pbeta function
print( c(estimate,returnedvalues) )#compare the eatimates with returned values
} #end for loop



## ----eval=FALSE----------------------------------------------------------
#  #construct function to compute  the expected (theoretical) count
#  expected <- function(colsum, rowsum, total) {
#    (colsum / total) * (rowsum / total) * total
#  }
#  #construct function to compute the value of the test-statistic
#  chi_stat <- function(observed, expected) {
#    ((observed - expected) ^ 2) / expected
#  }
#  
#  
#  chisq_test2 <- function(x, y) {
#    total <- sum(x) + sum(y)
#    rowsum_x <- sum(x)
#    rowsum_y <- sum(y)
#    chistat <- 0
#  # computes the chi-square test statistic which is apparently different from chisq.test function
#    for (i in seq_along(x)) {
#      colsum <- x[i] + y[i]
#      expected_x <- expected(colsum, rowsum_x, total)
#      expected_y <- expected(colsum, rowsum_y, total)
#      chistat <- chistat + chi_stat(x[i], expected_x)
#      chistat <- chistat + chi_stat(y[i], expected_y)
#    }
#    chistat
#  }
#  o <- as.table(rbind(c(56, 67, 57,34), c(135, 125, 45,65)))
#  #compare the results of chi-square test
#  print(chisq_test2(c(56, 67, 57,34), c(135, 125, 45,65)))
#  print(chisq.test(as.table(rbind(c(56, 67, 57,34), c(135, 125, 45,65)))))
#  #compare the computing time
#  print(microbenchmark::microbenchmark(
#    chisq_test2(c(56, 67, 57,34), c(135, 125, 45,65)),
#    chisq.test(as.table(rbind(c(56, 67, 57,34), c(135, 125, 45,65))))
#  ))
#  
#