`1`. 회귀분석 (35점)

아래와 같은 수식을 고려하자.

$$y_i= 1.5+ 5\sin(t_i)+ 3 \cos(4 t_i)+\epsilon_i, \quad i=1,2,\dots, 1000$$

여기에서 $t_i=\frac{2\pi i}{1000}$ 이다. 그리고 $\epsilon_i \sim i.i.d~ N(0,1)$, 즉 서로 독립인 표준정규분포에서 추출된 샘플이다.

`(1)` `rnorm()`을 이용하여 $\{\epsilon_i \}_{i=1}^{1000}$ 을 만들고 이를 벡터로 만들어라. 즉 아래와 같은 벡터를 생성하라.

$\boldsymbol{\epsilon}=(\epsilon_1,\dots,\epsilon_{1000})$

epsilon = rnorm(1000, 0, 1)
head(epsilon)

`(2)` $\big\{\sin(t_i) \big\}_{i=1}^{1000}$와 $\big\{\cos(4t_i) \big\}_{i=1}^{1000}$를 각각 생성하고 이를 ${\bf x}_1$와 ${\bf x}_2$에 저장하라. 즉 아래와 같은 2개의 벡터를 생성하라.

${\bf x}_1=\big(\sin(\frac{2\pi}{1000}),\sin(\frac{4\pi}{1000}),\dots,\sin(2\pi)\big)$
${\bf x}_2=\big(\cos(\frac{8\pi}{1000}),\cos(\frac{16\pi}{1000}),\dots,\cos(8\pi)\big)$

t = c()
for (i in 1:1000) t[i] = (2*pi*i)/1000
head(t)

x1 = c()
for(i in 1:1000) x1[i] = sin(t[i])
head(x1)

x2 = c()
for(i in 1:1000) x2[i] = cos(4*t[i])
head(x2)

`(3)` (1),(2)로 부터 $y_i$를 계산하고, $(t_i,y_i)$를 scatter plot으로 시각화 하라. 색은 'gray60'을 사용하라. 시각화 예시는 아래와 같다.

y = c()
for(i in 1:1000) y[i] = 1.5+5*x1[i]+3*x2[i]+epsilon[i]
head(y)

plot(t,y)

`(4)` 두 벡터 ${\bf x}_1$, ${\bf x}_2$와 `cbind()`를 이용하여 아래와 같은 $1000 \times 3$ 매트릭스 ${\bf X}$를 만들어라.

$${\bf X}=\begin{bmatrix} 1 & \sin(\frac{2\pi}{1000}) & \cos(\frac{8\pi}{1000})\\ 1 & \sin(\frac{4\pi}{1000}) & \cos(\frac{16\pi}{1000})\\ \dots & \dots & \dots \\ 1 & \sin(2\pi) & \cos(8\pi)\\ \end{bmatrix} $$

x0 = rep(1, times = 1000)
head(x0)

X = cbind(x0, x1, x2)
head(X)

`(5)` $3 \times 1$ 매트릭스 $\boldsymbol{\beta}=\begin{bmatrix} 1.5 \\ 5 \\ 3 \end{bmatrix}$를 만들고 매트릭스 곱을 이용하여 ${\bf X} {\boldsymbol \beta}$를 계산하라. 계산결과를 벡터화하고 (3)의 $y$축에 붉은선으로 시각화 하라. 시각화 예시는 아래와 같다. (필요하다면 적절하게 선의 굵기를 조정할 것)

b = rbind(1.5, 5, 3)
b

rl = X %*% b

plot(t,y)
points(t, rl, type = 'l', lwd = 5, col = "red")

`(6)` 아래의 수식을 이용하여 $\boldsymbol{\hat\beta}$를 계산하고 계산결과를 $\boldsymbol{\beta}$와 비교하라.

$$ \boldsymbol{\hat\beta} = ({\bf X}^\top {\bf X})^{-1}{\bf X}^\top {\bf y}$$

여기에서 ${\bf X}$와 ${\bf y}$는 각각 아래와 같다.

${\bf X}=\begin{bmatrix} 1 & \sin(\frac{2\pi}{1000}) & \cos(\frac{4\pi}{1000})\\ 1 & \sin(\frac{8\pi}{1000}) & \cos(\frac{16\pi}{1000})\\ \dots & \dots & \dots \\ 1 & \sin(2\pi) & \cos(8\pi)\\ \end{bmatrix}$, $\quad {\bf y}=\begin{bmatrix} y_1\\ y_2\\ \dots \\ y_{1000}\\ \end{bmatrix}$

또한 ${\bf X}^\top$는 ${\bf X}$의 전치행렬(transpose) 그리고 $({\bf X}^\top {\bf X})^{-1}$은 ${\bf X}^\top {\bf X}$의 역행렬을 의미한다.

X

y1 = cbind(y)
head(y1)

bhat = solve(t(X) %*% X) %*% t(X) %*% y1
bhat 

print(cbind(c(b),c(bhat)))

     [,1]     [,2]
[1,]  1.5 1.472656
[2,]  5.0 5.099967
[3,]  3.0 3.077521

`(7)` ${\bf X}\boldsymbol{\hat\beta}$를 계산하라. 계산결과를 벡터화 하고 (5)의 시각화 결과에 푸른점선으로 시각화 하라. 시각화 예시는 아래와 같다. (필요하다면 적절하게 선의 굵기를 조정할 것)

bl = X %*% bhat
bl

plot(t,y)
points(t, rl, type = 'l', lwd = 5, col = "red")
points(t, bl, type = 'l', lwd = 5, lty = 2, col = "blue")

`2`. 몬테카를로 적분 (10점)

`(1)` $\int_{-1.96}^{1.96}\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}x^2}dx$의 값을 `runif()`를 이용한 몬테카를로 적분으로 계산하라.

x=seq(-1.96, 1.96, 0.01)
y=(1/sqrt(2*pi))*exp((-1/2)*x^2)
min(y)
max(y)

plot(x,y,type = 'l')

xx = runif(10000, -2, 2)
yy = runif(10000, 0.05, 0.40)

plot(xx,yy)

plot(xx, yy)
lines(x,y,col = 'red', lwd = 3)

test = function(xx,yy){
    yy < (1/sqrt(2*pi))*exp((-1/2)*xx^2)
}

print(c(xx[1],yy[1]))
print((1/sqrt(2*pi))*exp((-1/2)*xx[1]^2))
test(xx[1],yy[1])

[1] -0.61375956  0.06942337
[1] 0.3304536

plot(xx,yy,col='gray')
lines(x,y,col='red',lwd=3)
points(xx[1],yy[1], col='blue')
points(xx[1],(1/sqrt(2*pi))*exp((-1/2)*xx[1]^2), col = 'red')

tst = c()
for (i in 1:10000) tst[i] = test(xx[i], yy[i])

head(tst)

plot(xx,yy,col='gray')
lines(x,y,col='red', lwd=3)
points(xx[tst],yy[tst],col='red')

sum(tst)

sum(tst)/10000 * 4*0.4

`(2)` `rnorm()`을 이용하여 표준정규분포에서 1000개의 확률변수를 생성하고 생성된 확률변수들 중 구간 $(-1.96 , 1.96)$사이에 있는 확률변수를 count하라.

X = rnorm(1000)
head(X)

Z = c()
for (i in 1:1000) Z[i] = -1.96 < X[i] & X[i] < 1.96 
Z

sum(Z)

`3`. 징검다리 (10점)

10명의 참가자가 (유리,강화유리)의 쌍으로 이루어진 징검다리를 건너는 게임을 시작한다. (유리,강화유리)의 쌍은 총 20개 있다. 따라서 유리는 총 40장이 있다. 참가자는 강화유리로 된 징검다리를 밟으면 살아남지만 유리로 된 징검다리를 밟으면 죽는다. 10명의 참가자의 특징은 아래와 같다.

1번~8번 참가자, 10번 참가자: 50퍼센트의 확률로 강화유리를 고를 수 있다.
9번 참가자: 95퍼센트의 확률로 강화유리를 고를 수 있다.

다리를 건너는 순서는 아래와 같이 번호의 오름차순으로 건너는 방법과 번호의 내림차순으로 건너는 2가지 방법이 있다고 하자.

Type A: 10번 참가자부터 다리를 건너고 그 뒤로 9번, 8번, ... 1번 순으로 다리를 건넌다.
Type B: 1번 참가자부터 다리를 건너고 그 뒤로 2번, 3번, ... 10번 순으로 다리를 건넌다.

8번 참가자는 Type A 와 Type B 중 어떠한 경우에 살아남을 확률이 높은가?

library(tidyverse)

Warning message:
“package ‘tidyverse’ was built under R version 4.0.3”
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

✔ ggplot2 3.3.5     ✔ purrr   0.3.4
✔ tibble  3.1.6     ✔ dplyr   1.0.7
✔ tidyr   1.1.4     ✔ stringr 1.4.0
✔ readr   2.1.1     ✔ forcats 0.5.1

Warning message:
“package ‘ggplot2’ was built under R version 4.0.5”
Warning message:
“package ‘tibble’ was built under R version 4.0.5”
Warning message:
“package ‘tidyr’ was built under R version 4.0.5”
Warning message:
“package ‘readr’ was built under R version 4.0.5”
Warning message:
“package ‘purrr’ was built under R version 4.0.3”
Warning message:
“package ‘dplyr’ was built under R version 4.0.5”
Warning message:
“package ‘stringr’ was built under R version 4.0.5”
Warning message:
“package ‘forcats’ was built under R version 4.0.3”
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

rm(list=ls())

SIMULATE1 = function(){
### 변수들의 모음
ARR = c('N1','N2','N3','N4','N5','N6','N7','N8','N9','N10')
SURV = 10 
PLAYER = ARR[SURV]
STAGE = 0 
TOSSRSLT = NA

### Subfunction들의 모음 
toss = function(p) rbinom(n=1,size=1,prob=p) %>% as.logical
reset = function(){
    TOSSRSLT <<- NA
    SURV <<- 10
    STAGE <<- 0 
    PLAYER <<- ARR[SURV]
}
record = function(){
    list(PRE_TOSSRSLT=TOSSRSLT, SURV=SURV, STAGE=STAGE, PLAYER=PLAYER) 
}
go = function(){
    PROB <<- 0.5+ (PLAYER=='N9')*0.45
    TOSSRSLT <<- toss(PROB) 
    if (TOSSRSLT==FALSE) SURV <<- SURV - 1 
    STAGE <<- STAGE + 1 
    PLAYER <<- ARR[SURV]    
}
gogo = function() for(i in 1:10) go()

gogo_history = function(){
    rslt_ = as_tibble(record()) ## go를 시작하기전의 결과가 티블로 만들어짐 
    for(i in 1:10){
        go()
        rslt_ = rbind(rslt_, as_tibble(record())) 
    }
    print(rslt_)
}

simulate_once = function(){
    reset()
    gogo()
    return(record()$SURV )
}

### BODY
simrslt = c() 
for (i in 1:1000) simrslt[i] = simulate_once() > 7
sum(simrslt)
}

rslt1 = SIMULATE1()
rslt1/1000*100

SIMULATE2 = function(){
### 변수들의 모음
ARR = c('N10','N9','N8','N7','N6','N5','N4','N3','N2','N1')
SURV = 10 
PLAYER = ARR[SURV]
STAGE = 0 
TOSSRSLT = NA

### Subfunction들의 모음 
toss = function(p) rbinom(n=1,size=1,prob=p) %>% as.logical
reset = function(){
    TOSSRSLT <<- NA
    SURV <<- 10
    STAGE <<- 0 
    PLAYER <<- ARR[SURV]
}
record = function(){
    list(PRE_TOSSRSLT=TOSSRSLT, SURV=SURV, STAGE=STAGE, PLAYER=PLAYER) 
}
go = function(){
    PROB <<- 0.5+ (PLAYER=='N9')*0.45
    TOSSRSLT <<- toss(PROB) 
    if (TOSSRSLT==FALSE) SURV <<- SURV - 1 
    STAGE <<- STAGE + 1 
    PLAYER <<- ARR[SURV]    
}
gogo = function() for(i in 1:10) go()

gogo_history = function(){
    rslt_ = as_tibble(record()) ## go를 시작하기전의 결과가 티블로 만들어짐 
    for(i in 1:10){
        go()
        rslt_ = rbind(rslt_, as_tibble(record())) 
    }
    print(rslt_)
}

simulate_once = function(){
    reset()
    gogo()
    return(record()$SURV )
}

### BODY
simrslt = c() 
for (i in 1:1000) simrslt[i] = simulate_once() > 2
sum(simrslt)
}

rslt2 = SIMULATE2()
rslt2/1000*100

`4`. COVID19 (15점)

아래는 한국의 코로나확진자 발생수를 tibble로 불러오는 코드와 그 결과이다.

library(tidyverse)

df=read_csv('https://raw.githubusercontent.com/guebin/2021IR/master/_notebooks/covid19.csv')
head(df)

Rows: 12294 Columns: 5

── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): prov
dbl (4): year, month, day, cases


ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

`(1)` 2020년의 확진자(cases)총합과 2021년의 확진자(cases)총합을 구하라.

year1 = df %>% filter(year == 2020 | year == 2021)

total_cases = year1 %>% summarize(sum_cases = sum(cases))
total_cases

`(2)` 2020년 2월1일 ~ 2020년 2월15일까지의 기간동안 확진자의 합을 지역별로 나누어 구하라 (아래예시 참고). 가장 많은 확진자가 발견된 지역은 어디인가?

즉 아래를 각각 계산하고 2020년 2월1일 ~ 2020년 2월15일 기간동안 발생한 확진자의 합이 가장 큰 지역을 찾으면 된다.

2020년2월1일 ~ 2020년2월15일 기간 서울에서 발생한 확진자의 합 = 2020년2월1일 서울에서 발생한 확진자수 + ... + 2020년2월15일 서울에서 발생한 확진자수
...
2020년2월1일 ~ 2020년2월15일 기간 제주에서 발생한 확진자의 합 = 2020년2월1일 제주에서 발생한 확진자수 + ... + 2020년2월15일 제주에서 발생한 확진자수

year2 = df %>% filter(year == 2020 & month == 2 & day > 0 & day < 16 )

prov_cases1 = year2 %>% group_by(prov) %>% summarize(sum_cases = sum(cases)) %>% arrange(desc(sum_cases)) %>% head(1)
prov_cases1

`(3)` 2020년 2월16일 ~ 2020년 2월29일까지의 기간동안 확진자의 총합을 지역별로 나누어 구하라 (아래예시 참고). 가장 많은 확진자가 발견된 지역은 어디인가?

즉 아래를 각각 계산하고 2020년 2월16일 ~ 2020년 2월29일 기간동안 발생한 확진자의 합이 가장 큰 지역을 찾으면 된다.

2020년2월16일 ~ 2020년2월29일 기간 서울에서 발생한 확진자의 합 = 2020년2월16일 서울에서 발생한 확진자수 + ... + 2020년2월29일 서울에서 발생한 확진자수
...
2020년2월16일 ~ 2020년2월29일 기간 제주에서 발생한 확진자의 합 = 2020년2월16일 제주에서 발생한 확진자수 + ... + 2020년2월29일 제주에서 발생한 확진자수

year3 = df %>% filter(year == 2020 & month == 2 & day > 15 & day < 30 )

prov_cases2 = year3 %>% group_by(prov) %>% summarize(sum_cases = sum(cases)) %>% arrange(desc(sum_cases)) %>% head(1)
prov_cases2

x0	x1	x2
1	0.006283144	0.9996842
1	0.012566040	0.9987370
1	0.018848440	0.9971589
1	0.025130095	0.9949510
1	0.031410759	0.9921147
1	0.037690183	0.9886517

x0	x1	x2
1	0.006283144	0.9996842
1	0.012566040	0.9987370
1	0.018848440	0.9971589
1	0.025130095	0.9949510
1	0.031410759	0.9921147
1	0.037690183	0.9886517
1	0.043968118	0.9845643
1	0.050244318	0.9798551
1	0.056518534	0.9745269
1	0.062790520	0.9685832
1	0.069060026	0.9620277
1	0.075326806	0.9548645
1	0.081590612	0.9470983
1	0.087851197	0.9387339
1	0.094108313	0.9297765
1	0.100361715	0.9202318
1	0.106611154	0.9101060
1	0.112856385	0.8994053
1	0.119097160	0.8881364
1	0.125333234	0.8763067
1	0.131564359	0.8639234
1	0.137790291	0.8509945
1	0.144010783	0.8375280
1	0.150225589	0.8235326
1	0.156434465	0.8090170
1	0.162637165	0.7939904
1	0.168833445	0.7784623
1	0.175023059	0.7624425
1	0.181205764	0.7459411
1	0.187381315	0.7289686
⋮	⋮	⋮
1	-1.812058e-01	0.7459411
1	-1.750231e-01	0.7624425
1	-1.688334e-01	0.7784623
1	-1.626372e-01	0.7939904
1	-1.564345e-01	0.8090170
1	-1.502256e-01	0.8235326
1	-1.440108e-01	0.8375280
1	-1.377903e-01	0.8509945
1	-1.315644e-01	0.8639234
1	-1.253332e-01	0.8763067
1	-1.190972e-01	0.8881364
1	-1.128564e-01	0.8994053
1	-1.066112e-01	0.9101060
1	-1.003617e-01	0.9202318
1	-9.410831e-02	0.9297765
1	-8.785120e-02	0.9387339
1	-8.159061e-02	0.9470983
1	-7.532681e-02	0.9548645
1	-6.906003e-02	0.9620277
1	-6.279052e-02	0.9685832
1	-5.651853e-02	0.9745269
1	-5.024432e-02	0.9798551
1	-4.396812e-02	0.9845643
1	-3.769018e-02	0.9886517
1	-3.141076e-02	0.9921147
1	-2.513010e-02	0.9949510
1	-1.884844e-02	0.9971589
1	-1.256604e-02	0.9987370
1	-6.283144e-03	0.9996842
1	-2.449294e-16	1.0000000

prov	sum_cases
<chr>	<dbl>
경기	9

prov	sum_cases
<chr>	<dbl>
대구	2055

A matrix: 6 × 1 of type dbl
y
3.598227
5.038326
5.436468
6.595559
5.668046
3.503458

A matrix: 1000 × 1 of type dbl
y
4.581249
4.610377
4.637560
4.662801
4.686104
4.707471
4.726910
4.744425
4.760026
4.773721
4.785520
4.795436
4.803480
4.809667
4.814012
4.816530
4.817240
4.816159
4.813306
4.808704
4.802372
4.794335
4.784616
4.773241
4.760233
4.745622
4.729435
4.711701
4.692449
4.671711
⋮
2.844162
2.926477
3.007345
3.086734
3.164612
3.240950
3.315716
3.388884
3.460425

1. 회귀분석 (35점)

(1) rnorm()을 이용하여 $\{\epsilon_i \}_{i=1}^{1000}$ 을 만들고 이를 벡터로 만들어라. 즉 아래와 같은 벡터를 생성하라.

(2) $\big\{\sin(t_i) \big\}_{i=1}^{1000}$와 $\big\{\cos(4t_i) \big\}_{i=1}^{1000}$를 각각 생성하고 이를 ${\bf x}_1$와 ${\bf x}_2$에 저장하라. 즉 아래와 같은 2개의 벡터를 생성하라.

(3) (1),(2)로 부터 $y_i$를 계산하고, $(t_i,y_i)$를 scatter plot으로 시각화 하라. 색은 'gray60'을 사용하라. 시각화 예시는 아래와 같다.

(4) 두 벡터 ${\bf x}_1$, ${\bf x}_2$와 cbind()를 이용하여 아래와 같은 $1000 \times 3$ 매트릭스 ${\bf X}$를 만들어라.

(6) 아래의 수식을 이용하여 $\boldsymbol{\hat\beta}$를 계산하고 계산결과를 $\boldsymbol{\beta}$와 비교하라.

(7) ${\bf X}\boldsymbol{\hat\beta}$를 계산하라. 계산결과를 벡터화 하고 (5)의 시각화 결과에 푸른점선으로 시각화 하라. 시각화 예시는 아래와 같다. (필요하다면 적절하게 선의 굵기를 조정할 것)

2. 몬테카를로 적분 (10점)

(1) $\int_{-1.96}^{1.96}\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}x^2}dx$의 값을 runif()를 이용한 몬테카를로 적분으로 계산하라.

(2) rnorm()을 이용하여 표준정규분포에서 1000개의 확률변수를 생성하고 생성된 확률변수들 중 구간 $(-1.96 , 1.96)$사이에 있는 확률변수를 count하라.

3. 징검다리 (10점)

4. COVID19 (15점)

(1) 2020년의 확진자(cases)총합과 2021년의 확진자(cases)총합을 구하라.

(2) 2020년 2월1일 ~ 2020년 2월15일까지의 기간동안 확진자의 합을 지역별로 나누어 구하라 (아래예시 참고). 가장 많은 확진자가 발견된 지역은 어디인가?

(3) 2020년 2월16일 ~ 2020년 2월29일까지의 기간동안 확진자의 총합을 지역별로 나누어 구하라 (아래예시 참고). 가장 많은 확진자가 발견된 지역은 어디인가?