library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(nortest)

第一題

a小題

set.seed(122)
start<-(1000000*runif(1))%>% floor() ;start

## [1] 905020

k<-c()
for(i in 1:10000){
  k[i]<-start^2 %>% substr(.,2,7) %>% as.numeric()
  start<-start^2 %>% substr(.,2,7) %>% as.numeric()
}
nchar(k) %>% table()

## .
##    4    5    6 
##    4  657 9339

hist(k)

ks.test(k,"punif",0,999999)

## Warning in ks.test(k, "punif", 0, 999999): ties should not be present for
## the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  k
## D = 0.072443, p-value < 2.2e-16
## alternative hypothesis: two-sided

b小題

x<-c(1)
y<-c(1)
z<-c(1)
u<-c()
for(i in 2:10001){
  x[i]<-(x[i-1]*171)%%30269
  y[i]<-(y[i-1]*172)%%30307
  z[i]<-(z[i-1]*170)%%30323
  u[i]<-((x[i]/30269)+(y[i]/30307)+(z[i]/30323))%%1
}
u<-u[-1]

ks.test(u,"punif")

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  u
## D = 0.009332, p-value = 0.3486
## alternative hypothesis: two-sided

h<-hist(u)

h$counts

##  [1] 475 471 484 498 487 517 538 502 522 506 481 516 519 490 467 507 498
## [18] 543 483 496

chisq.test(h$counts)

## 
##  Chi-squared test for given probabilities
## 
## data:  h$counts
## X-squared = 17.052, df = 19, p-value = 0.5863

c小題*要解釋 >想把hist 一起比較發現a小提的數字有規律性 bc和normal很近比較p-value

uc<-c(1)
for(i in 2:10000){
  uc[i]<-((uc[i-1]+pi)^5)%%1
}
hist(uc)

ks.test(uc,"punif")

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  uc
## D = 0.0073868, p-value = 0.6462
## alternative hypothesis: two-sided

ks.test(uc,u)

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  uc and u
## D = 0.0143, p-value = 0.2582
## alternative hypothesis: two-sided

第二題

ALL

test<-function(k,j){
  N<-c(10,50,100)
  D<-c(10,20)
  a<-c();b<-c();c<-c();d<-c()
  for(i in 1:1000){
    tmp1<-rt(N[k],D[j]) %>% ks.test(.,"pnorm") %>% .[2] %>%`<`(0.05) %>% `*`(1) %>% as.numeric()
    tmp2<-rt(N[k],D[j]) %>% cvm.test() %>% .[2] %>%`<`(0.05) %>% `*`(1) %>% as.numeric()
    tmp3<-rt(N[k],D[j]) %>% ad.test()  %>% .[2] %>%`<`(0.05) %>% `*`(1) %>% as.numeric()
    tmp4<-rt(N[k],D[j]) %>% lillie.test() %>% .[2] %>% `<`(0.05) %>% `*`(1) %>% as.numeric()
    a<-sum(a,tmp1);b<-sum(b,tmp2);c<-sum(c,tmp3);d<-sum(d,tmp4)
  }
  paste(a,b,c,d,sep = "/")
}
norr<-function(k){
  N<-c(10,50,100)
  a<-c();b<-c();c<-c();d<-c()
  for(i in 1:1000){
    tmp1<-rnorm(N[k]) %>% ks.test(.,"pnorm") %>% .[2] %>%`<`(0.05) %>% `*`(1) %>% as.numeric()
    tmp2<-rnorm(N[k]) %>% cvm.test() %>% .[2] %>%`<`(0.05) %>% `*`(1) %>% as.numeric()
    tmp3<-rnorm(N[k]) %>% ad.test()  %>% .[2] %>%`<`(0.05) %>% `*`(1) %>% as.numeric()
    tmp4<-rnorm(N[k]) %>% lillie.test() %>% .[2] %>% `<`(0.05) %>% `*`(1) %>% as.numeric()
    a<-sum(a,tmp1);b<-sum(b,tmp2);c<-sum(c,tmp3);d<-sum(d,tmp4)
  }
  paste(a,b,c,d,sep = "/")
}


Alltest<-function(){
  N<-c(10,50,100)
  D<-c(10,20)
  m<-matrix(rep(0,9),nrow = 3)
  for(j in 1:2){  #D 自由度
    for(k in 1:3){ #n 樣本數
        m[j,k]<-test(k,j)
      }
  }
  for(k in 1:3){
    m[3,k]<-norr(k)
  }
  rownames(m)<-c("df = 10","df = 20","Normal");colnames(m)<-paste0("n = ",N);m
  }
Alltest()

##         n = 10        n = 50          n = 100         
## df = 10 "51/62/66/75" "56/119/128/85" "63/133/177/106"
## df = 20 "50/57/59/56" "47/71/78/50"   "53/76/99/60"   
## Normal  "53/50/52/43" "54/46/25/47"   "48/61/39/45"

第三題

Gap test

gap.test=function(data,a,b){
  x<-ifelse(a<data & data<b,1,0)
   x1=which(x==1)
    y=x1[-1]-x1[-length(x1)]-1
     k<-unique(y) %>% sort() #看有哪些數字
      obs<-table(y) %>% as.numeric()# 實際值
       prob<-(b-a)*((1-b+a))^k #理論機率
      exp<-length(data)*prob
     stop<-which(exp<5) #從stop開始期望次數會小於5
    go<-!(1:length(exp)%in%stop)
   obs<-c(obs[go],sum(obs[stop]))
  prob<-c(prob[go],1-sum(prob[go]))
 chisq.test(obs,p=prob)
}
gap.test(runif(1000),0.1,0.4)

## Warning in chisq.test(obs, p = prob): Chi-squared approximation may be
## incorrect

## 
##  Chi-squared test for given probabilities
## 
## data:  obs
## X-squared = 11.735, df = 12, p-value = 0.4672

Permutation test

permute.test=function(data,k){
  y=rep(10,k)^c((k-1):0)
  x=matrix(data,ncol=k,byrow=T)
  x1=apply(x,1,rank)
  yy=apply(x1*y,2,sum)
  return(table(yy))
}

Run test

run<-function(r){
  r1 <- 1*(r[-1]>r[-1000] ) #變號位置
sum(r1[-1] != r1[-999]) + 1 # 共有幾個run
r2 <- which(r1[-1] != r1[-999]) #run的位置
r3 <- table(r2[-1]-r2[-length(r2)]) #看run＝n的個數
r3
Z <- (sum(r3)-(2*1000-1)/3)/((16*1000-29)/90)^0.5
1-pnorm(Z)
}
run(runif(1000))

## [1] 0.7342006

第四題

第一小題

k<-c()
for(i in 1:10000){
  k[i]<-runif(12) %>% sum() %>% `-`(6)
}
ks.test(k,"pnorm")

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  k
## D = 0.0066173, p-value = 0.7736
## alternative hypothesis: two-sided

nor<-rnorm(10000)
y<-cut(nor,breaks = c(-100,-2,-1,0,1,2,100)) %>% table() %>% `/`(10000)
x<-cut(k,breaks = c(-100,-2,-1,0,1,2,100)) %>% table() %>% `/`(10000)
p<-y %>% as.numeric()
chisq.test(x,p = p)

## Warning in chisq.test(x, p = p): Chi-squared approximation may be incorrect

## 
##  Chi-squared test for given probabilities
## 
## data:  x
## X-squared = 0.0023654, df = 5, p-value = 1

第二小題*未完成

第五題

第一小題

\(e-\frac{1}{3.5+x}\)

Bisection

fcn<-function(x){exp(1)-(1/(3.5+x))}
left<--3.5
right<- -2
i<-1
while(right-left>1e-6){
  m<-sum(left,right)/2
  ifelse(fcn(left)*fcn(m)<0,right<-m,left<-m)
  i<-i+1
}
paste0("The root is ",m,"\ and i = ",i)

## [1] "The root is -3.13212037086487 and i = 22"

False position

F.position<-function(){
 fcn<-function(x){exp(1)-(1/(3.5+x))}
  a<--3.4;b<--3
  i<-1
   while(abs(a-b)>1e-6|(fcn(a)*fcn(b)>0)){
    m1<-(fcn(a)-fcn(b))/(a-b)
     b1<-fcn(a)-m1*a
      tmp1<-(-(b1/m1))

      m2<-(fcn(tmp1)-fcn(a))/(tmp1-a)
     b2<-fcn(tmp1)-m2*tmp1
    tmp2<-(-(b2/m2))
   a<-tmp2;b<-tmp1
   i<-i+1
 }
  paste0("The root is ",tmp2,"\ and i = ",i)
}
F.position()

## [1] "The root is -3.13212055882856 and i = 6"

secant method

secant<-function(){
fcn<-function(x){exp(1)-(1/(3.5+x))}
a<--4;b<--3
i<-1
while(abs(a-b)>1e-6){
m1<-(fcn(a)-fcn(b))/(a-b)
b1<-fcn(a)-m1*a
tmp1<-(-(b1/m1))

m2<-(fcn(tmp1)-fcn(b))/(tmp1-b)
b2<-fcn(tmp1)-m2*tmp1
tmp2<-(-(b2/m2))
a<-tmp2;b<-tmp1
i<-i+1
}
paste0("The root is ",tmp2,"\ and i = ",i)
}
secant()

## [1] "The root is -3.13212055883988 and i = 6"

Uniroot

fcn<-function(x){exp(1)-(1/(3.5+x))}
uniroot(fcn,c(-3.2,-3))$root

## [1] -3.132115

第二小題

Bisection

left<--10000
right<-10000
fcn2<-function(x){(exp(-x)/(sqrt((1+x^2))))-0.5}
i<-1
while(right-left>1e-10){
  m<-sum(left,right)/2
  ifelse(fcn2(left)*fcn2(m)<0,right<-m,left<-m)
  i<-i+1
}
paste0("The root is ",m,"\ and i = ",i)

## [1] "The root is 0.557728698922233 and i = 49"

False position

F.position<-function(){
fcn2<-function(x){(exp(-x)/(sqrt((1+x^2))))-0.5}
a<-1;b<-10
i<-1
while(abs(a-b)>1e-6|(fcn2(a)*fcn2(b)>0)){
m1<-(fcn2(a)-fcn2(b))/(a-b)
b1<-fcn2(a)-m1*a
tmp1<-(-(b1/m1))

m2<-(fcn2(tmp1)-fcn2(a))/(tmp1-a)
b2<-fcn2(tmp1)-m2*tmp1
tmp2<-(-(b2/m2))
a<-tmp2;b<-tmp1
i<-i+1
}
paste0("The root is ",tmp2,"\ and i = ",i)
}
F.position()

## [1] "The root is 0.557728698980331 and i = 7"

secant method

secant<-function(){
fcn2<-function(x){(exp(-x)/(sqrt((1+x^2))))-0.5}
a<-1;b<-4
i<-1
while(abs(a-b)>1e-6){
m1<-(fcn2(a)-fcn2(b))/(a-b)
b1<-fcn2(a)-m1*a
tmp1<-(-(b1/m1))

m2<-(fcn2(tmp1)-fcn2(b))/(tmp1-b)
b2<-fcn2(tmp1)-m2*tmp1
tmp2<-(-(b2/m2))
a<-tmp2;b<-tmp1
i<-i+1
}
paste0("The root is ",a,"\ and i = ",i)
}
secant()

## [1] "The root is 0.557728698980331 and i = 9"

Uniroot

fcn2<-function(x){(exp(-x)/(sqrt((1+x^2))))-0.5}
uniroot(fcn2,c(-3,3))$root

## [1] 0.5577375

第三小題

Bisection

left<-c(-0.6,-0.3,0.2)
right<-c(-0.3,0.2,14)
fcn3<- function(x){
((2-x)*(4-x)*(6-x)+3*4*5+4*3*5-4*(4-x)*4-5*5*(2-x)-(6-x)*3*3 )
}
test<-function(left,right){
  i<-1
  while(right-left>1e-10){
  m<-sum(left,right)/2
  ifelse(fcn3(left)*fcn3(m)<0,right<-m,left<-m)
  i<-i+1
  }
  return(c(left,i))
}
# fun(-0.6)
# fun(-0.3)
# fun(0.2)
# fun(14)
q<-test(left[1],right[1])
w<-test(left[2],right[2])
e<-test(left[3],right[3])

paste0("The root is ",q[1],"\ and i = ",q[2])

## [1] "The root is -0.480740698473528 and i = 33"

paste0("The root is ",w[1],"\ and i = ",w[2])

## [1] "The root is -1.16415210804632e-11 and i = 34"

paste0("The root is ",e[1],"\ and i = ",e[2])

## [1] "The root is 12.4807406983826 and i = 39"

False position

a<-c(-0.6,-0.3,0.2)
b<-c(-0.3,0.2,14)

test<-function(a,b){
  fcn3<- function(x){
((2-x)*(4-x)*(6-x)+3*4*5+4*3*5-4*(4-x)*4-5*5*(2-x)-(6-x)*3*3 )
}
  i<-1
while(abs(a-b)>1e-6|(fcn3(a)*fcn3(b)>0)){
m1<-(fcn3(a)-fcn3(b))/(a-b)
b1<-fcn3(a)-m1*a
tmp1<-(-(b1/m1))

m2<-(fcn3(tmp1)-fcn3(a))/(tmp1-a)
b2<-fcn3(tmp1)-m2*tmp1
tmp2<-(-(b2/m2))
a<-tmp2;b<-tmp1
i<-i+1
}
  return(c(a,i))
}

q<-test(a[1],b[1])
w<-test(a[2],b[2])
e<-test(a[3],b[3])

paste0("The root is ",q[1],"\ and i = ",q[2])

## [1] "The root is -0.480740698420936 and i = 4"

paste0("The root is ",w[1],"\ and i = ",w[2])##

## [1] "The root is 1.86992039053421e-15 and i = 7"

paste0("The root is ",e[1],"\ and i = ",e[2])##

## [1] "The root is 2.29797423893347e-15 and i = 5"

secant method

a<-c(-0.6,-0.3,0.2)
b<-c(-0.3,0.2,14)

secant<-function(a,b){
fcn3<-function(x){((2-x)*(4-x)*(6-x)+3*4*5+4*3*5-4*(4-x)*4-5*5*(2-x)-(6-x)*3*3 )}
i<-1
while(abs(a-b)>1e-6){
m1<-(fcn3(a)-fcn3(b))/(a-b)
b1<-fcn3(a)-m1*a
tmp1<-(-(b1/m1))

m2<-(fcn3(tmp1)-fcn3(b))/(tmp1-b)
b2<-fcn3(tmp1)-m2*tmp1
tmp2<-(-(b2/m2))
a<-tmp2;b<-tmp1
i<-i+1
}
paste0("The root is ",a,"\ and i = ",i)
}

secant(a[1],b[1])

## [1] "The root is -0.480740698407861 and i = 5"

secant(a[2],b[2])##

## [1] "The root is -9.50859926297007e-11 and i = 5"

secant(a[3],b[3])##

## [1] "The root is 1.67215864063299e-15 and i = 6"

Uniroot

a<-c(-0.6,-0.3,0.2)
b<-c(-0.3,0.2,14)
fcn3<-function(x){((2-x)*(4-x)*(6-x)+3*4*5+4*3*5-4*(4-x)*4-5*5*(2-x)-(6-x)*3*3 )}

uniroot(fcn3,c(a[1],b[1]))$root

## [1] -0.480741

uniroot(fcn3,c(a[2],b[2]))$root

## [1] 7.944986e-07

uniroot(fcn3,c(a[3],b[3]))$root

## [1] 12.48074

第六題

Secant method

fcn<-function(x){#這部分偷懶直接算
  1997/(2+x)-1810/(1-x)+32/x
}

secant<-function(){
a<-runif(1);b<-runif(1)
i<-1
while(abs(a-b)>1e-10){
m1<-(fcn(a)-fcn(b))/(a-b)
b1<-fcn(a)-m1*a
tmp1<-(-(b1/m1))

m2<-(fcn(tmp1)-fcn(b))/(tmp1-b)
b2<-fcn(tmp1)-m2*tmp1
tmp2<-(-(b2/m2))
a<-tmp2;b<-tmp1
a<-ifelse((a-b)=="NaN",0,a)
b<-ifelse((a-b)=="NaN",0,b)
i<-i+1
}
return(c(a,i))
# paste0("The root is ",a,"\ and i = ",i)
}

k<-secant()
while(k[1]<=0){
 k<-secant() 
}
paste0("The root is ",k[1],"\ and i = ",k[2])

## [1] "The root is 0.0357123022406281 and i = 6"

(Another method )

牛頓法

fcn<-function(x){#這部分直接算
  1997/(2+x)-1810/(1-x)+32/x
}

intgo<-function(){
  r<-c(1,2)
  k<-2
  int<-runif(1)
  plus<-0.0000001
  while(abs(r[k]-r[k-1])>1e-6){#T 繼續做 F停止
    x1<-c(int,fcn(int))
    x2<-c((int+plus),fcn(int+plus))
    m<-(x2[2]-x1[2])/plus
    b<-x2[2]-m*x2[1]
    int<-(-b/m)
    r[k+1]<-int
    k<-k+1
  }
  return(c(int,length(r)))
}

k<-intgo()
while(k[1]>1|k[1]<0){
  k<-intgo()
}
paste0("The root is ",k[1]," and i =",k[2])

## [1] "The root is 0.0357123022397858 and i =12"

統計模擬大全

Yu & Hsu

2018/3/17

第一題

a小題

b小題

c小題*要解釋 >想把hist 一起比較發現a小提的數字有規律性 bc和normal很近比較p-value

第二題

ALL

第三題

Gap test

Permutation test

Run test

第四題

第一小題

第二小題*未完成

第五題

第一小題

Bisection

False position

secant method

Uniroot

第二小題

Bisection

False position

secant method

Uniroot

第三小題

Bisection

False position

secant method

Uniroot

第六題

Secant method

(Another method )

牛頓法

統計模擬大全

Yu & Hsu

2018/3/17

第一題

a小題

b小題

c小題*要解釋 >想把hist 一起比較 發現a小提的數字有規律性 bc和normal很近 比較p-value

第二題

ALL

第三題

Gap test

Permutation test

Run test

第四題

第一小題

第二小題*未完成

第五題

第一小題

Bisection

False position

secant method

Uniroot

第二小題

Bisection

False position

secant method

Uniroot

第三小題

Bisection

False position

secant method

Uniroot

第六題

Secant method

(Another method )

牛頓法

c小題*要解釋 >想把hist 一起比較發現a小提的數字有規律性 bc和normal很近比較p-value