확률론2-1 : 코딩을 활용한 확률계산[in r]-1

R코딩으로 기본적인 확률계산을 해보겠습니다.

install.package("dplyr")
install.package("magrittr")
library("dplyr")
library("magrittr")
# 공정한 주사위 4개를 던지는 시행의 표본공간
omega1 <- expand.grid(dice1=1:6,dice2=1:6,
                      dice3=1:6,dice4=1:6)

# 네 주사위 눈의 합이 20 이상이 될 확률
#filter 사용.
(apply(omega1,1,sum)>=20) %>% sum # 1:TRUE, 0:FALSE임.
apply(omega1,1,sum)>=20
   [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
  [12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
  [23] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
  [34] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
  [45] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
  [56] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
  [67] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
  [78] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
  [89] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [100] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [111] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [122] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [144] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [155] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [166] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [177] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [188] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [199] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [210] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [221] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [232] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [243] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [254] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [265] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [276] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [287] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [298] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [309] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [320] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [331] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [342] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [353] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [364] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [375] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [386] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [397] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [408] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [419] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [430] FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [441] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [452] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [463] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [474] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [485] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [496] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [507] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [518] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [529] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [540] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [551] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [562] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [573] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [584] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [595] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [606] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE
 [617] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [628] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [639] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE
 [650] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [661] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [672] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [683] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [694] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [705] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [716] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [727] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [738] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [749] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [760] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [771] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [782] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
 [793] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [804] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [815] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE
 [826] FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [837] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [848] FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE
 [859] FALSE FALSE FALSE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE
 [870] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [881] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [892] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [903] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [914] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [925] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [936] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [947] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [958] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [969] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [980] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [991] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [ reached getOption("max.print") -- omitted 296 entries ]

 (apply(omega1,1,sum)>=20) %>% sum # 1:TRUE, 0:FALSE임.
[1] 70

라고 표현해도 되고 dplyr에서의 필터링을 통해서 필터해도 된다. subset도 가능. 각각 써보면 아래와 같다.

> A<-omega1 %>% filter(apply(.,1,sum)>=20) %>% nrow()
> A
[1] 70
> #subset 사용.
> subset(omega1,apply(omega1,1,sum)>=20) %>% nrow()
[1] 70
#눈 4개 합 20이상확률
>70/nrow(omega1)
[1] 0.05401235

> # "네 주사위 눈의 합"이 각각 나올 확률과 그래프
> ##=> dice4 옆에 sum이라는 변수를 만들고 합을 열로 합.
> sum<-apply(omega1,1,sum)
> omega1<-cbind(omega1,sum)
> omega1 %>% head()
  dice1 dice2 dice3 dice4 sum
1     1     1     1     1   4
2     2     1     1     1   5
3     3     1     1     1   6
4     4     1     1     1   7
5     5     1     1     1   8
6     6     1     1     1   9
#또는 다음과 같이 한번에 써도 된다.
> omega1 <- expand.grid(dice1=1:6,dice2=1:6,
+                       dice3=1:6,dice4=1:6)
> omega1$sum<-apply(omega1,1,sum) #바로 컬이름 적고 적용.
> omega1 %>% head()
  dice1 dice2 dice3 dice4 sum
1     1     1     1     1   4
2     2     1     1     1   5
3     3     1     1     1   6
4     4     1     1     1   7
5     5     1     1     1   8
6     6     1     1     1   9

이제 각 시행마다 눈이 얼마나 나온지 궁금하므로 table형태로 나타내면 빈도수를 알 수있다.

> (omega1$sum %>% table)
.
  4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21 
  1   4  10  20  35  56  80 104 125 140 146 140 125 104  80  56  35  20 
 22  23  24 
 10   4   1 
 #전체 개수로 나눠지면 그에 대한 확률이 나온다.
 > (omega1$sum %>% table)/nrow(omega1)
.
           4            5            6            7            8 
0.0007716049 0.0030864198 0.0077160494 0.0154320988 0.0270061728 
           9           10           11           12           13 
0.0432098765 0.0617283951 0.0802469136 0.0964506173 0.1080246914 
          14           15           16           17           18 
0.1126543210 0.1080246914 0.0964506173 0.0802469136 0.0617283951 
          19           20           21           22           23 
0.0432098765 0.0270061728 0.0154320988 0.0077160494 0.0030864198 
          24 
0.0007716049 
#이를 bar plot으로 그리면 다음과 같다.
> ((omega1$sum %>% table)/nrow(omega1)) %>% barplot()

4개 주사위 합 20이상일 확률의 bar plot.

# 2 ####
#동전을 열 번 던지는 시행에서 
# 앞면이 세 번 이상 나오는 사상을 A,
# 앞면이 뒷면보다 많이 나오는 사상을 B,
# 뒷면이 세 번 이상 나오는 사상을 C라 하자

> # 표본공간
> coin<-c("H","T")
> coin
[1] "H" "T"
> omega2<-expand.grid(c1=coin,c2=coin,c3=coin,c4=coin,
+                     c5=coin,c6=coin,c7=coin,c8=coin,
+                     c9=coin,c10=coin)

동전 열번 던졌을 때의 표본공간.(총 2^10=1024개)

# 사상 A, B, C
#앞면,뒷면 둘다 문자형이므로 그냥 셀 수는 없음.
#그러므로, H/T면 카운팅하는 함수를 만들자. 
#어떤 원소가 들어왔을 때 그것이 H/T인가?
#그리고 그것의 합은?(T/F는 1/0으로 취급.)
countH <- function(x) sum(x=="H")
countT <- function(x) sum(x=="T")

실제로 잘 돌아가는지 확인해 보자.

> omega2 %>% head
  c1 c2 c3 c4 c5 c6 c7 c8 c9 c10
1  H  H  H  H  H  H  H  H  H   H
2  T  H  H  H  H  H  H  H  H   H
3  H  T  H  H  H  H  H  H  H   H
4  T  T  H  H  H  H  H  H  H   H
5  H  H  T  H  H  H  H  H  H   H
6  T  H  T  H  H  H  H  H  H   H

> countH(omega2[1,])
[1] 10
> countH(omega2[2,])
[1] 9
> countT(omega2[1,])
[1] 0
> countT(omega2[2,])
[1] 1

카운팅 함수에는 문제 없다.

> #A:H가 3개이상
> A<-omega2 %>% filter(apply(omega2,1,countH)>=3)
#apply(data,가로or세로계산,무슨계산?)
> A %>% head()
  c1 c2 c3 c4 c5 c6 c7 c8 c9 c10
1  H  H  H  H  H  H  H  H  H   H
2  T  H  H  H  H  H  H  H  H   H
3  H  T  H  H  H  H  H  H  H   H
4  T  T  H  H  H  H  H  H  H   H
5  H  H  T  H  H  H  H  H  H   H
6  T  H  T  H  H  H  H  H  H   H

사상 A 은 앞면(H)가 3개이상인 사상이다. 그렇기에 표본공간에서 문자형을 카운팅 시킬 함수를 이용한 필터링을 통해 3개이상임을 찾아주면된다. 근데 이것이 진짜 3개이상인지는 아래와 같이 검산해보았다.

> #진짜로 3개 이상?
> (apply(A,1,countH)<3) %>% sum
[1] 0

내가 구한 사상A에서 3미만이 있는지 체크해두면 된다. 0이므로 내가 원하는 사상을 잘 구했다고 할 수 있겠다.

#B:앞면수가 뒷면수보다 많아야함.
>B<-omega2 %>% filter(apply(omega2,1,countH) >= apply(omega2,1,countT))
> B %>% head()
  c1 c2 c3 c4 c5 c6 c7 c8 c9 c10
1  H  H  H  H  H  H  H  H  H   H
2  T  H  H  H  H  H  H  H  H   H
3  H  T  H  H  H  H  H  H  H   H
4  T  T  H  H  H  H  H  H  H   H
5  H  H  T  H  H  H  H  H  H   H
6  T  H  T  H  H  H  H  H  H   H

B는 표본공간에서 앞면수가 뒷면수보다 많음을 체크해야하는데, 이는 표본공간의 데이터프레임 형태의 행의 합으로 앞서 구한 카운팅함수를 이용해 개수를 세어주면 된다.

#C:뒷면이 세번이상
C<-omega2 %>% filter(apply(omega2,1,countT)>=3)

> # Ac = Omega2 - A
> paste("Pr{A^c}=",1-(nrow(A)/nrow(omega2)))
[1] "Pr{A^c}= 0.0546875"
> # A∩Bc = A - B (순수A만 나타냄.)
> setdiff(A,B) %>% nrow()
[1] 330
> paste("Pr{A∩Bc} = A - B",(nrow(setdiff(A,B))/nrow(omega2)))
[1] "Pr{A∩Bc} = A - B 0.322265625"
> #Pr{A∩C} 
> paste("Pr{A∩C} = ",(nrow(intersect(A,C))/nrow(omega2)))
[1] "Pr{A∩C} =  0.890625"
> #Pr{B∩C}
> paste("Pr{B∩C} = ",(nrow(intersect(B,C))/nrow(omega2)))
[1] "Pr{B∩C} =  0.568359375"

위의 Ac는 A의 여사상(complement)를 의미한다. 그리고 차사상은 setdiff라는 함수를 이용해 바로 나타낼 수 있다.

> #Pr{A∩B∩C}
> paste("Pr{A∩B∩C} = ",(nrow(intersect(A,B,C))/nrow(omega2))) #error
Error in `intersect()`:
! `...` must be empty.
✖ Problematic argument:
• ..1 = C
ℹ Did you forget to name an argument?
Run `rlang::last_trace()` to see where the error occurred.
#한번에 교사상을 구하는 것은 불가능하다. A,B,C와 ∩ 하나로 A∩B∩C표현은 불가능하기에 ∩(intersect)를
#두번 사용해야한다.
> paste("Pr{A∩B∩C} = ",(nrow(intersect(intersect(A,B),C))/nrow(omega2)))
[1] "Pr{A∩B∩C} =  0.568359375"

 

-출처: AI소프트웨어학과 이두호교수님 강의파일