1. round(0.125, 2)를 수행하여 값을 얻어라. round(0.1+0.025, 2)를 수행하여 값을 얻어라. 같은가?

round(0.125,2)
## [1] 0.12
round(0.1+0.025,2)
## [1] 0.12
round(0.135,2)
## [1] 0.14
round(0.145,2)
## [1] 0.14

10진법의 0.025는 2진법에서는 무한소수이다. 무한 자릿수를 컴퓨터에 넣을 수 없으니 뒷부분을 자르고 저장되면서 0.025는 0.025보다 아주 조금 작은 수로 입력이 된다. 하지만 10진법의 0.125는 1/8이 되어 2진법의 0.001 이 되어 유한소수이다. CPU 안에서는 80bit로 계산되고 CPU 밖으로 저장될 때는 64bit.

round(0.125,2) = round(0.1+0.025) = 0.12 의 동일한 결과를 얻었다. 10진법 수를 기준으로 반올림이 되었더라면 두 연산 모두 0.13이어야 한다. R round 함수의 설명을 참조해보니, 해당 함수는 round를 수행의 기준이 되는 자릿수(위의 경우에는 소수 세번째 자리수) 가 5인 경우 가까운 짝수 자릿수로 올림 혹은 내림하도록 설정되어 있기 때문인 것으로 확인되었다. 즉, 0.135를 소수 두번째 자리까지 round를 적용하는 경우 가까운 짝수인 0.14로, 0.145는 0.14가 되는 식이다.R help page의 원문은 아래 첨부하였다.

?round - Details: Note that for rounding off a 5, the IEC 60559 standard (see also ‘IEEE 754’) is expected to be used, ‘go to the even digit’. Therefore round(0.5) is 0 and round(-1.5) is -2. However, this is dependent on OS services and on representation error (since e.g. 0.15 is not represented exactly, the rounding rule applies to the represented number and not to the printed number, and so round(0.15, 1) could be either 0.1or 0.2).

2. 10^seq(100, 1000, 100)을 수행하여 R에서 어느 정도 큰 수까지 표시할 수 알아보아라. 앞으로 컴퓨터의 성능이 128비트로 높아지면 더 큰 수도 표시할 수 있을 것이다. 0이 아니면서 0에 가까운 수 또한 중요하다. 다음을 실행하여 machine \(\epsilon\) 값을 대충으로 찾아보아라.

10^seq(100,1000,100)
##  [1] 1e+100 1e+200 1e+300    Inf    Inf    Inf    Inf    Inf    Inf    Inf
10^seq(300,400,10)
##  [1] 1e+300    Inf    Inf    Inf    Inf    Inf    Inf    Inf    Inf    Inf
## [11]    Inf
10^seq(300,310,1)
##  [1] 1e+300 1e+301 1e+302 1e+303 1e+304 1e+305 1e+306 1e+307 1e+308    Inf
## [11]    Inf

위와 같이 순차적으로 코드를 시행한 결과,컴퓨터가 \(10^308\) 까지 표현할 수 있는 것으로 나타났다.

0.1^seq(100,1000,100)
##  [1] 1e-100 1e-200 1e-300  0e+00  0e+00  0e+00  0e+00  0e+00  0e+00  0e+00
0.1^seq(300,400,10)
##  [1] 1.000000e-300 1.000000e-310 9.999889e-321  0.000000e+00  0.000000e+00
##  [6]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [11]  0.000000e+00
0.1^seq(320,330,1)
##  [1] 9.999889e-321 9.980126e-322 9.881313e-323 9.881313e-324  0.000000e+00
##  [6]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [11]  0.000000e+00

위의 코드를 수행한 결과 `machine \(\epsilon\)은 약\(10^{-320} < \epsilon < 10^{-310}\) 수준인 것으로 나타났다.

3. 초기값은 본인 학번을 사용하고 congruential (A=7^5, c=0, M=2^63) 방법으로 0과 1사이의 난수를 만들어라.

unif.cong = function(n){
  x = vector(mode="numeric",length = n)
  x[1]=2014122060
  for(i in 2:n){
    x[i]=(((7^5)*x[i-1])%%(2^63))
  }
  return(x/(2^63))
}
  1. 주사위 변수를 5개 발생시켜 평균을 구하는 것을 100개 하여 평균의 히스토그램을 그려라.
  2. 주사위 변수를 30개 발생시켜 평균을 구하는 것을 100개 하여 평균의 히스토그램을 그려라.
dice_sample = function(s,n){
  random = as.integer(unif.cong(n*s)*6)+1
  x = colMeans(matrix(random,s))
  return(x)
}


set.seed(2014122060)
x = dice_sample(5,100)
set.seed(2014122060)
y = dice_sample(30,100)
par(mfrow=c(1,2))
hist(x, main = "100 Sample Means with Sample Size 5",xlim=c(1,6))
abline(v=3.5,col="red",lty=3,lwd=3)

hist(y, main = "100 Sample Means with Sample Size 30",xlim=c(1,6))
abline(v=3.5,col="red",lty=3,lwd=3)

  1. 위의 두 히스토그램을 중심극한정리에 근거하여 비교하여라.
summary(x)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.200   3.000   3.600   3.502   4.000   5.200
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.433   3.333   3.500   3.468   3.675   4.100
summary = matrix(c(sum(1:6)/6,25/12,mean(x),var(x),mean(y),var(y)),
                 ncol=3,nrow=2,byrow=FALSE,
                 dimnames=list(c("mean","var"),c("theory","n=5","n=30")))
summary
##        theory       n=5       n=30
## mean 3.500000 3.5020000 3.46766667
## var  2.083333 0.6266626 0.09426386
boxplot(x,y,names=c("small sample","large sample"))

위의 상자그림을 살펴보면, 표본의 갯수가 5개인 소표본 자료의 그래프보다, 표본의 갯수가 30개인 대표본 자료의 그래프가 더욱 중심에 모여있는 형태를 보이고 있다. 또한 소표본 자료 그래프의 경우 \(median = 3.6>mean = 3.502\)으로 중위값이 평균보다 더 크므로 음의 왜도를 보인다고 할 수 있다. 상자그림 역시 중위값을 기준으로 상자의 하단부가 더 길게 나타나고 있다. 대표본 자료의 그래프 역시 중위값이 평균보다 약간 큰 값을 보이지만, 소표본의 경우 그 차이가 0.1(=3.60-3.502)이었던 반면, 대표본의 경우 0.032(=3.50-3.468) 수준으로 그 차이가 60% 가량 줄어들었다. 대표본 자료의 상자그림 역시 중위값을 기준으로 상자 및 수염의 길이가 대칭적인 형태를 보임을 확인할 수 있다.

히스토그램과 상자그림,summary 를 종합적으로 살펴볼 때, 표본의 갯수가 늘어남에 따라 모집단의 분포(해당 문제의 경우 uniform distribution)와 무관하게 표본 평균의 분포는 평균이 \(\mu\) , 분산이 \(\sigma^2/n\)인 정규분포에 근사한다는 것(중심극한정리)을 확인할 수 있다. 히스토그램으로 살펴보면, 빨간색 점선으로 표시된 것이 모평균 \(\mu\) =3.5 이다. 평균은 n=5,n=30의 경우 모두 모평균 \(\mu=3.5\)와 같거나 매우 유사한 값을 보인다. 분산의 경우 역시 중심극한 정리를 적용한 경우의 값과 유사한 값을 보인다.. 소표본의 경우 \(0.6266 \approx \frac{2.083}{5}\), 대표본의 경우 \(0.0942 \approx \frac{2.083}{30}\)이다.

par(mfrow=c(1,2))
qqnorm(x)
qqline(x,col="red")
qqnorm(y)
qqline(y,col="red")

shapiro.test(x)
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.98395, p-value = 0.2661
shapiro.test(y)
## 
##  Shapiro-Wilk normality test
## 
## data:  y
## W = 0.97422, p-value = 0.0469

엄밀하게 정규분포를 따르는지 확인하기 위해 Normal QQ plot과 정규성 검정을 수행하였다. Normal QQ plot을 살펴보면, 소표본과 대표본 경우 모두 직선에 가까운 형태를 보이고 있는 것을 확인할 수 있다. shapiro.test() 의 검정 결과를 살펴보면, 99% 유의 수준 하에 대표본과 소표본 모두 정규분포를 따른다. Shapiro-Wilk 검정의 귀무가설은 “해당 자료가 정규분포를 따른다.”이고, 대표본과 소표본의 검정 결과 p-value가 모두 유의수준 0.01보다 크기 때문에 두 경우 모두 정규분포를 따른다는 주장은 통계적으로 유의하다. 중심 극한 정리는 sample size(n)이 무한대에 가까워질 때의 limiting distribution이 정규분포를 따른다는 정리이지만, 해당 simulation의 결과를 통해 유한개의 sample size를 가지고도 통계적으로 유의미한 수준에서 정규 분포 근사가 가능함을 확인했다. 특히, 주사위 던지기의 경우에는 표본의 크기(sample size)가 5개만 되어도 평균, 분산의 형태 뿐만 아니라 통계적으로도 유의미한 수준에서 정규분포를 따른다는 것(CLT)을 확인할 수 있다.

4. A deck of 100 cards—numbered 1, 2, . . . , 100—is shuffled and then turned over one card at a time. Say that a “hit” occurs whenever card i is the i th card to be turned over, i = 1, . . . , 100.

  1. Write a simulation program to estimate the expectation and variance of the total number of hits. Run the program.
rand.perm = function(x){
  for (i in length(x):1){
    I = as.integer(runif(1)*i)+1
    x[c(i,I)]=x[c(I,i)]
  }
  return(x)
}

hit = function(x){
  n = 0
  for(i in 1:length(x)){
    if(i == x[i]) {
      n = n + 1 }
  }
  print(n)
}


hit_dist = function(data,n){
  dist = vector(mode="numeric",length =n)
  for(i in 1:n){
    set = rand.perm(data)
    dist[i] = hit(set)
  }
  hit_dist = dist
}
set.seed(2014122060)
hit_1000 = hit_dist(1:100,1000)
summary = data.frame(mean = mean(hit_1000),
                     var = var(hit_1000))
summary

simulation을 이용하여 구한 hit에 대한 1000개의 확률변수를 히스토그램으로 그리면 아래와 같다. 위에 붉은 점선으로 덧대어진 곡선은 Hit에 대한 이론적 확률 분포(Poisson Distribution with $ = 1$) 를 이용하여 구한 것이다.simaultion을 통해 구현한 확률변수의 평균은 0.977 \(\approx 1\), 분산은 0.9394104 \(\approx 1\)으로 이론적인 평균(1)과 분산(1)에 근사하는 값을 보이고 있다.

hist(hit_1000,main = "Simulated # of Hits(n=1000)",
    xlab="X - Values", ylim = c(0,400),
    breaks=c(seq(-0.5,5.5,1)),xlim=c(0,6))

curve(1000*exp(-1)/factorial(x),0,6,col="red",
      lwd=2, lty=2, add=T)

  1. Find the exact answers and compare them with your estimates.

<손으로 계산해 제출할 예정>

5. Give a method for simulating X, having the probability mass function pj , j =5, 6, . . . , 14, where pj = 0.11 when j is odd and 5 ≤ j ≤ 13, pj =0.09 when j is even and 6 ≤ j ≤ 14. Use the text’s random number sequence to generate X.

  1. Inverse Transformation Method
sim.xx = function(n){
  unif = runif(n)
  cdf = matrix(c(5:14,0.11,0.2,0.31,0.4,0.51,
           0.6,0.71,0.8,0.91,1),10,2)
  index = vector(mode="numeric", length = n)
  result = vector(mode="numeric", length = n)
  for (i in 1:n){
   index[i] = sum(cdf[,2]<unif[i])+1
   result[i] = cdf[index[i],1]
  }
  return(result)
}

set.seed(2014122060)
inverse.result = sim.xx(1000)
  1. Rejection Method
sim.yy = function(n){
  result = vector(mode ="numeric",length = n)
  data = data.frame(X =5:14, prob = rep(c(0.11,0.09),5))
  for(j in 1:n){
    y = as.integer(runif(n)*10)+5
    u2 = runif(n)
      for (i in 1:n){
        if (u2[i] > (data[data$X == y[i],2]/0.11)) next
        result[j] = y[i]
        break }
  }
  return(result)
}

set.seed(2014122060)
reject.result = sim.yy(1000)
  1. 위의 두 가지 방식으로 생성한 X 확률변수의 히스토그램에 이론값을 덧대어 표시하고 일치하는지 검토하라.
par(mfrow=c(1,2))

theor.sample = c(rep(seq(5,13,2),110),rep(seq(6,14,2),90))

hist(inverse.result,breaks=c(4:14), xlab = "Simulated X-values",
    main = "Inverse Trans Method",col= "mistyrose",
    ylim = c(0,140))

hist(theor.sample,breaks =c(4:14),col=rgb(1,1,0,0.2),
    ylim= c(0,140), add=T)

hist(reject.result,breaks=c(4:14), xlab = "Simulated X-values",
    main = "Rejection Method", col= "lightskyblue1",
    ylim=c(0,140))

hist(theor.sample,breaks =c(4:14),col=rgb(1,1,0,0.2),
    ylim= c(0,140),add=T)

set.seed(2014122060)
system.time(sim.xx(1000))
##    user  system elapsed 
##       0       0       0
set.seed(2014122060)
system.time(sim.yy(1000))
##    user  system elapsed 
##    0.14    0.00    0.14
diff = c(sum((sort(theor.sample)-sort(inverse.result))^2)/1000,
         sum((sort(theor.sample)-sort(reject.result))^2)/1000)
diff
## [1] 0.078 0.082

위 그래프에서 분홍색으로 표시된 것은 Inverse Transformation Method로 생성한 확률변수, 노란색은 이론적인 확률 질량 함수로부터 구현한 확률변수, 하늘색으로 표시한 것은 Rejection Method를 이용하여 생성한 확률변수를 뜻한다. 두 그래프 모두 크게는 이론적 확률변수의 히스토그램의 경향성을 따라가고 있다. 짝수에서 낮은 빈도수, 홀수에서 높은 빈도수를 나타낸다. 서로 다른 랜덤 변수 생성 방법에 따라 이론적인 분포와 얼마나 차이가 나는지를 확인하기 위해 Mean Square Difference를 계산해보았다. Inverse Transformation Method의 경우 그 값이 0.078, Rejection Method의 경우 그 값이 0.082로 나타났다. 이를 통해 Inverse Transformation Method가 이 경우에는 이론적 확률분포와 더욱 유사한 확률 변수를 생성했을 것임을 기대할 수 있고, 그래프를 통한 시각적 확인이 가능하다. 또한 각 방법에 따라 확률 변수 생성에 걸린 시간을 확인해 볼 수 있다. 1000개의 랜덤변수를 생성하는데 Inverse Method의 경우 0.02초, Rejection Method의 경우 0.15초가 소요됬다. 이 경우 Inverse Transformation Method로 생성한 결과가 이론적 분포와 더욱 유사했으며, 그 알고리즘의 계산속도도 우월하게 나타났다.

par(mfrow=c(1,3))

hist(inverse.result,breaks=c(4:14), xlab = "Simulated X-values",
    main = "Inverse Trans Method",col= "mistyrose",
    ylim = c(0,140))

hist(theor.sample,breaks =c(4:14),col="lightyellow",
    ylim= c(0,140), xlab = "All X-Values",
    main = "Theoretical Dist'n")

hist(reject.result,breaks=c(4:14), xlab = "Simulated X-values",
    main = "Rejection Method", col= "lightskyblue1",
    ylim=c(0,140))