1.
## [1] 1
midSquareRand <- function(seed, length) {
  randvector <- NULL
  for(i in 1:length) {
    value <- seed * seed 
    seed <- (value %/% 1000) %% 1000000
    randvector <- c(randvector, seed)
  }   
  return(randvector)
}

random <- midSquareRand(123457, 10000)
ks.test(random , "punif")
## Warning in ks.test(random, "punif"): ties should not be present for the
## Kolmogorov-Smirnov test
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  random
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

x <- 3
y <- 67
z <- 4453
u <- NULL
for(i in 1:10000){
  x <- (x*171)%%30269
  y <- (y*172)%%30307
  z <- (z*170)%%30323
  u <- c(u,((x/30269)+(y/30307)+(z/30323))%%1)
}
hist(u)

ks.test(u ,"punif")
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  u
## D = 0.007003, p-value = 0.7107
## alternative hypothesis: two-sided
chisq.test( u , p =rep(1/length(u),length(u)))
## Warning in chisq.test(u, p = rep(1/length(u), length(u))): Chi-squared
## approximation may be incorrect
## 
##  Chi-squared test for given probabilities
## 
## data:  u
## X-squared = 1683.4, df = 9999, p-value = 1
U <- function(u){
  ((pi+u)^5)%%1
}

u <- 10
r <- NULL
for (i in 1:10000) {
  u <- U(u)
  r <- c(r,u)
}
hist(r)

ks.test(r ,"punif")
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  r
## D = 0.0075075, p-value = 0.6259
## alternative hypothesis: two-sided
library(nortest)
r <- NULL
r[[1]] <- rnorm(10)
r[[2]] <- rnorm(50)
r[[3]] <- rnorm(100)
r[[4]] <- rt(10,10)
r[[5]] <- rt(50,10)
r[[6]] <- rt(100,10)
r[[7]] <- rt(10,20)
r[[8]] <- rt(50,20)
r[[9]] <- rt(100,20)

ks <- NULL
ad <- NULL
cvm <- NULL
lillie <- NULL
pearson <- NULL
sf <- NULL
p_value <- NULL
for(i in 1:9){
  ks <- ks.test(r[[i]],"pnorm")$p.value
  ad <- ad.test(r[[i]])$p.value
  cvm <- cvm.test(r[[i]])$p.value
  lillie <- lillie.test(r[[i]])$p.value
  pearson <- pearson.test(r[[i]])$p.value
  sf <- sf.test(r[[i]])$p.value
  p_value <- rbind(p_value,c(ks,ad,cvm,lillie,pearson,sf))
}
colnames(p_value)[c(1:6)] <- c("ks.test","ad.test","cvm.test","lillie.test","pearson.test","sf.test")
rownames(p_value)[c(1:9)] <- c("n=10,N(0,1)","n=50,N(0,1)","n=100,N(0,1)","n=10,t(10)","n=50,t(10)","n=100,t(10)","n=10,t(20)","n=20,t(20)","n=100,t(20)")
p_value
##                 ks.test   ad.test  cvm.test lillie.test pearson.test
## n=10,N(0,1)  0.46956622 0.1331787 0.1265189   0.2140911   0.07855316
## n=50,N(0,1)  0.78978283 0.9318983 0.8928752   0.7663081   0.29864634
## n=100,N(0,1) 0.87204257 0.3148176 0.5997653   0.6232773   0.53780714
## n=10,t(10)   0.57312796 0.6160639 0.5318301   0.2951582   0.22138539
## n=50,t(10)   0.02571953 0.1719093 0.1691688   0.3167980   0.14758397
## n=100,t(10)  0.68956032 0.5627958 0.4783630   0.3332106   0.11682464
## n=10,t(20)   0.30295665 0.4739937 0.4903994   0.7284836   0.22138539
## n=20,t(20)   0.54868708 0.1499623 0.1735947   0.2822765   0.40835741
## n=100,t(20)  0.59582833 0.9020565 0.9174786   0.9745250   0.94734698
##                  sf.test
## n=10,N(0,1)  0.204674490
## n=50,N(0,1)  0.753345982
## n=100,N(0,1) 0.008058011
## n=10,t(10)   0.810119433
## n=50,t(10)   0.278466643
## n=100,t(10)  0.743163757
## n=10,t(20)   0.352663958
## n=20,t(20)   0.233070293
## n=100,t(20)  0.890714074
#power

nfun <- function(n){
  c(ks.test(rnorm(n),"pnorm")$p.value,
    ad.test(rnorm(n))$p.value,
    cvm.test(rnorm(n))$p.value,
    lillie.test(rnorm(n))$p.value,
    pearson.test(rnorm(n))$p.value,
    sf.test(rnorm(n))$p.value)
}
tfun <- function(n,df){
  c(ks.test(rt(n,df),"pnorm")$p.value,
    ad.test(rt(n,df))$p.value,
    cvm.test(rt(n,df))$p.value,
    lillie.test(rt(n,df))$p.value,
    pearson.test(rt(n,df))$p.value,
    sf.test(rt(n,df))$p.value)
}
q01 <- c(0,0,0,0,0,0);q02 <- c(0,0,0,0,0,0);q03 <- c(0,0,0,0,0,0);q11 <- c(0,0,0,0,0,0); q12 <- c(0,0,0,0,0,0) ;q13 <- c(0,0,0,0,0,0) ; q21 <- c(0,0,0,0,0,0) ; q22 <- c(0,0,0,0,0,0) ; q23 <-c(0,0,0,0,0,0)
for(i in 1:1000){
  q01 <- q01+(1*nfun(10)<0.05)
  q02 <- q02+(1*nfun(50)<0.05)
  q03 <- q03+(1*nfun(100)<0.05)
  q11 <- q11+(1*(tfun(10,10)<0.05))
  q12 <- q12+(1*(tfun(50,10)<0.05))
  q13 <- q13+(1*(tfun(100,10)<0.05))
  q21 <- q21+(1*(tfun(10,20)<0.05))
  q22 <- q22+(1*(tfun(50,20)<0.05))
  q23 <- q23+(1*(tfun(100,20)<0.05))
}
power <- rbind(q01/1000,q02/1000,q03/1000,q11/1000,q12/1000,q13/10000,q21/1000,q22/1000,q23/1000) 
colnames(power)[c(1:6)] <- c("ks.test","ad.test","cvm.test","lillie.test","pearson.test","sf.test")
rownames(power)[c(1:9)] <- c("n=10,N(0,1)","n=50,N(0,1)","n=100,N(0,1)","n=10,t(10)","n=50,t(10)","n=100,t(10)","n=10,t(20)","n=20,t(20)","n=100,t(20)")
power
##              ks.test ad.test cvm.test lillie.test pearson.test sf.test
## n=10,N(0,1)   0.0480  0.0560   0.0530      0.0550       0.0490  0.0490
## n=50,N(0,1)   0.0560  0.0480   0.0480      0.0430       0.0510  0.0480
## n=100,N(0,1)  0.0440  0.0440   0.0400      0.0480       0.0560  0.0420
## n=10,t(10)    0.0500  0.0720   0.0700      0.0830       0.0800  0.0810
## n=50,t(10)    0.0530  0.1080   0.1130      0.0790       0.0670  0.1790
## n=100,t(10)   0.0046  0.0167   0.0138      0.0123       0.0081  0.0302
## n=10,t(20)    0.0550  0.0580   0.0570      0.0540       0.0730  0.0660
## n=20,t(20)    0.0460  0.0660   0.0770      0.0550       0.0470  0.1240
## n=100,t(20)   0.0530  0.0690   0.0540      0.0740       0.0580  0.1500
############ 3.
#gap
set.seed(10)
x <- runif(1000)
r <- runif(2)
x1 <- which( x >= min(r) & x < max(r))
length(x1)
## [1] 96
x2 <- x1[-1]-x1[-length(x1)]
k <- table(x2) %>% as.data.frame() %>% .[,2]
kk <- k[c(1:9)] 
kk[10] <- sum(k[c(10:length(k))])
kkk <- kk
thm <- function(k){
  (max(r)-min(r))*(1-(max(r)-min(r)))^k
}
gap <- NULL
for(i in 1:length(k)){
  gap <- c(gap,thm(i))
}
gap1 <- gap[c(1:9)]
gap1[10] <- 1-sum(gap1)
sum(((kkk-gap1*sum(kkk))^2)/(gap1*sum(kkk)))
## [1] 15.21951
qchisq(0.95,9)
## [1] 16.91898
#run

set.seed(9)
r <- runif(1000)
r1 <- 1*(r[-1]>r[-1000] ) #變號位置
sum(r1[-1] != r1[-999]) + 1 # 共有幾個run
## [1] 664
r2 <- which(r1[-1] != r1[-999]) #run的位置
r3 <- table(r2[-1]-r2[-length(r2)]) #看run＝n的個數
r3
## 
##   1   2   3   4   5 
## 420 171  54  13   4
Z <- (sum(r3)-(2*1000-1)/3)/((16*1000-29)/90)^0.5
#permutation



p <- matrix(runif(1600),4,400)
p1 <- apply(p, 2, rank)
per <- p1[1,]*1000+p1[2,]*100+p1[3,]*10+p1[4,]
x <- table(per) %>% as.data.frame() %>% .[,2]
Q <- sum(((x-(400/24))^2)/400/24)
qchisq(0.95,23)
## [1] 35.17246
#5
#bisection
f <- function(x){
  exp(1)-(1/(3.5+x))
}
g <- function(x){
  (exp(-x)/sqrt(1+x^2))-0.5
}
Ei=function(x){
  (2-x)*(4-x)*(6-x)+60+60-16*(4-x)-9*(6-x)-25*(2-x)
}
bis <- function(f,a,b){
  z <- abs(b-a)
  i <- 0
  print(z)
  while (z > 0.000000001) {
    c <- (a+b)/2
    ifelse(f(c)<0, a <- c , b <- c)
    z <- abs(b-a)
    i <- i+1
    
  }
  print(c(c,i))
}
curve(g,0,1);abline(h=0)

bis(f,-3.2,-3)
## [1] 0.2
## [1] -3.132121 28.000000
bis(g,0.6,0.4)
## [1] 0.2
## [1]  0.5577287 28.0000000
bis(Ei,-0.5,-0.4)
## [1] 0.1
## [1] -0.4 27.0
bis(Ei,-0.4,0)
## [1] 0.4
## [1] -7.450581e-10  2.900000e+01
bis(Ei,13,12)
## [1] 1
## [1] 12.48074 30.00000
#false positions
f <- function(x){
  exp(1)-(1/(3.5+x))
}
g <- function(x){
  (exp(-x)/sqrt(1+x^2))-0.5
}
Ei=function(x){
  (2-x)*(4-x)*(6-x)+60+60-16*(4-x)-9*(6-x)-25*(2-x)
}
fal <- function(f,a,b){
  z <- abs(b-a)
  i <- 0
  while(z >0.000001){
    c <- b-(f(b)*((b-a)/(f(b)-f(a))))
    if(f(b)*f(c)>0) b <- c
    z <- f(b)
    i <- i+1
    
  }
  print(c(c,i))
}
curve(f,-5,0)
abline(h=0)

fal(f,-3.4,-3)
## [1] -3.13212 44.00000
fal(g,0.4,0.6)
## [1] 0.5612611 1.0000000
fal(Ei,-0.5,-0.4)
## [1] -0.4768946  1.0000000
fal(Ei,-0.4,0)
## [1] 0 1
fal(Ei,13,12)
## [1] 12.48074  8.00000
curve(Ei,-1,0.01);abline(h=0)

#
#secant
f <- function(x){
  exp(1)-(1/(3.5+x))
}
g <- function(x){
  (exp(-x)/sqrt(1+x^2))-0.5
}
Ei=function(x){
  (2-x)*(4-x)*(6-x)+60+60-16*(4-x)-9*(6-x)-25*(2-x)
}
sec <- function(f,a,b){
  z <- abs(f(b)-f(a))
  i <- 0
  while (z > 0.000001) {
    c <- b-(f(b)*((b-a)/(f(b)-f(a))))
    a <- b
    b <- c
    z <- abs(f(b)-f(a))
    i <- i+1
    
  }
  print(c(c,i))
}
sec(f,-3.2,-3)
## [1] -3.132121  6.000000
sec(g,0.4,0.6)
## [1] 0.5577287 4.0000000
sec(Ei,-0.5,-0.4)
## [1] -0.4807407  5.0000000
sec(Ei,-0.4,0)
## [1] 0 1
sec(Ei,13,12)
## [1] 12.48074  6.00000
#6
add=function(a){
  (1997/(2+a))-(906/(1-a))-(904/(1-a))+(32/(a))
}

add2=function(a){
  (-1997/((2+a)^2))-(1810/((1-a)^2))-(32/(a^2))}

newton=function(f,fp,start){
  i=0
  new=start
  r=c(i,new,f(new))
  while(abs(f(new)>10^(-7))){
    i=i+1
    new=new-(f(new)/fp(new))
    r=rbind(r,c(i,new,f(new)))}
  r}
newton(add,add2,0.5)
## [1]     0.0     0.5 -2757.2