R의 확률함수, 정규분포

R은 확률 함수를 가지고 있어요.
우리가 보통(normal) 알고 있는 분포가 정규분포(normal distribution)이지요.

정규분포는 평균값을 중심으로 확률밀도가 대칭이에요.
평균과 분산만 알고 있다면 확률값을 계산할 수 있어 편리해요.

R은 분포와 관련해서 네 가지 종류의 머릿 글자를 가지고 있어요.

  1. d - probabilty Density: 확률을 계산
  2. p - cumulative Probability: 누적 확률을 계산
  3. q - probabilty Quantily: 누적 확률의 inverse
  4. r - probability Random: 주어진 정규분포로 랜덤 값 생성
위의 4번은 사실 분포 활용이고, 3번은 cdf의 역함수라 우리의 관심은 아니에요. 결국 1번 2번에 관심이 가지요. 조금 헷갈리나요? d가 확률이고 p가 누적확률이라는 것? 왜 cdf(cumulative probability density function)이니까 c로 시작해야 할 것 같은데 p로 시작하지요? 음... 이렇게 생각하지요. "확률"이라는 말은 너무 흔해서 그냥 안쓰기로 했어요. 그래서 probability density function(pdf)를 그냥 density function이라고 한 것이지요. 그래서 아! density function이니까 d로 시작해야겠다 이러고 d-로 시작하는 pdf를 정의했어요. 그럼 cdf는요? -____- 시작을 probability 없이 했으니 이것도 첫 단어는 지워버리고 시작해야겠다 생각을 한 것이지요.. 그래서 cumulative probability density function에서 cumulative를 지워버리고 p-로 시작했어요. 데헷.

그리하여...
정규분포를 나타내는 norm 에다가 d와 p를 붙여서...

d + norm -> dnorm : 정규분포의 PDF
p + norm -> pnorm : 정규분포의 CDF


이렇게 되었어요.

이제 정규분포의 PDF로 그래프를 그려 정말 정규분포인지 알아보지요.


plot(x=200:300,y=dnorm(200:300,mean=250,sd=8),type='l',xlab='Score,ylab='Probability')

평균이 250, 표준편차가 8인 정규분포에서 x값이 200과 300 사이인 확률 값들을 계산했어요. 250을 중심으로 좌우 대칭이네요. 다시 한번! pDf라서 d가 사용된 것이라는 점!

그럼 문제를 바꿔서 242와 266 사이에 x값이 있을 확률은 얼마나 될까요? 이때는 cPdf를 써야겠네요.

pnorm(266,250,8) - pnorm(242,250,8)
[1] 0.8185946

첫 파라미터가 x값, 두 번째는 평균, 세 번째는 표준편차에요.
pnorm() 함수로 cdf(혹은 cPdf)를 계산했어요.

약 82%의 확률로 우리가 242부터 266까지의 값을 가질 수가 있겠어요.

정규분포는 대단한 것이랍니다.
왜냐하면 우리가 알아야 할 대상이 10 만 개라고 해도 그것의 1%에 해당하는 1 천 개의 샘플을 잘 구하면 표본의 평균과 표본의 분산을 계산해서 정규분포를 활용할 수 있기 때문이에요. 즉, 1 천 개 표본의 평균이 250이라고 하고, 표본 분산이 8이라고 했을 때, 10 만 개에 해당하는 전체 대상의 점수가 242와 266 사이에 있을 확률을 0.82라고 계산할 수 있으니까요. 어림잡아 8 만 2 천 명이 저 점수를 받았다 짐작해요. 이것이 얼마나 믿을 만한 것인가 하는 이야기에는 조금 다른 생각이 필요하지요. ㅋㅋ

** Central Limit Theorem **

중심극한정리(CLT)는 샘플 수가 적당히 크면(n>=30 정도) 표본의 확률분포가 정규분포에 가깝다는 정리에요. 큰 수의 법칙(Law of Large Numbers)이 표본평균을 모평균처럼 사용해도 된다는 것이었다면 CLT는 이제 분포도 정규분포로 쓰면 된다는 식이니 참 좋죠.

R에서 평균을 구하는 함수 mean()이나 표준편차를 구하는 함수 sd()는 모두 "표본"에 대한 것이에요. 큰 수의 법칙 때문에 mean()은 그냥 모평균이겠거니 생각해도 되는데 sd()는 어떻게 하나요? 우리가 표준편차를 구하기 전에 분산을 구할 때 생각하면

Variance = sigma( (x-mu)^2 ) * 1/N

이렇게 하고 N은 전체 population의 숫자였어요. 그런데 표본분산을 모분산과 같다고 하려면, 표본의 크기를 n이라 할 때

Sample Variance = sigma( (x-mu)^2 ) * 1/(n-1)

이렇게 구해야 해요. R은 바로 Sample Variance의 제곱근을 구해서 돌려주는 sd() 함수를 가지고 있어요.

예를 들어 자료가 10,20,30이라면 평균은 20이고, 표본의 표준편차는

sqrt( ( (10-20)^2 + (20-20)^2 + (30-20)^2 ) * 1/(3-1) ) = 10

이고, R로 계산하면

sd( c(10,20,30) )
[1] 10

이렇답니다.

** 유의 수준 **

[문제-a] 수원대학교 축구부는 한 해 경기에서 평균 20번 이긴다. 또한 이것에 관한 모표준편차가 5라고 한다. 25번 경기를 한 뒤 표본평균이 21이었다. 한 해 평균 20번 이긴다고 정말 이야기할 수 있나?

이 분포의 샘플분산은 = 모분산 / n 로 계산하는 점에 주의!

대세가설(귀무가설)은 "평균이 20이다" 에요.

이것을 지지하는 근거를 값으로 나타낸 것이 p-value 에요.

pnorm(21, 20, 1, lower.tail = FALSE)
[1] 0.1586553

lower.tail=FALSE로 하면 1에서 cdf를 빼고 남은 값을 돌려줘요. 즉, 누적확률의 바깥 쪽 값을 계산해요. 평균이 20인데 x=21이라고 하면 이제 남은 확률이 15.8% 정도인 것이죠.

이제 p-value를 구해보지요. 먼저 표준화를 해서 Z값을 계산해야 해요. 사실을 말씀드리면 이미 우리가 구한 값이 표준화한 값이에요. R이 미리 그렇게 계산했어요.. 못 믿겠다면

Z = (x-mu) / sd 니까...

x = 21일 때, Z = (21 - 20)/1 = 1

pnorm(1,0,1,lower.tail=FALSE)
[1] 0.1586553

같아요.

p-value는

y <- pnorm(21,20,1,lower.tail=FALSE)
pvalue <- 2*y

이고, 이 때의 pvalue는  0.37 정도에요. 100명 중에서 37명 정도는 수원대학교 축구부가 한 해 평균 20번 정도 이긴다는 증거를 확실히 가지고 있다는 의미에요. 보통 우리가 100명 중에 5명만 증거를 확실히 가지고 있어도 괜찮다고 생각한다면 37명은 아주 큰 숫자에요. 따라서 표본평균이 21이지만 한 해 20번 이긴다는 주장을 반박하기 어려워요.

유의수준: 100명 중에 5명 정도가 확실하다는 것을 alpha=0.05 인 유의수준이라고 하고 100명 중에 1명 정도가 확실하다는 것을 alpha=0.01의 유의수준이라고 해요. 유의수준이란 넘기면 귀무가설 못 깬다는 기준이에요.

25번 경기 결과 평균이 20번 이긴다라고 나온다면

2*pnorm(20,20,1,lower.tail=F)
[1] 1

이라서 100명이면 100명 전부 지지한다는 것,

평균이 23으로 나온다면

2*pnorm(23,20,1,lower.tail=F)
[1] 0.0027

로 1만명 중에서 고작 27명만이 지지한다는 것.

p-value는 재판관 같은 것이에요.

이처럼 R에서 pnorm을 사용하면 p-value도 손쉽게 확인할 수 있어요.

이제 5% 지지자들의 신뢰구간까지 구해보지요.

qnorm() 함수를 써서 주어진 cdf의 quantile, 즉 x값을 구하지요.

5%지지자가 양쪽에 다 있으니 한 쪽을 기준으로 하면 2.5%에요.

K <- qnorm(0.025,lower.tail=F) * 5/(sqrt(25))

20 - K
[1] 18.35515

20 + K
[1] 21.64485

유의수준 5% 내에서는 19승, 20승, 21승까지는 20승이라고 인정해줄 만 하네요.

댓글

이 블로그의 인기 게시물

R에서 csv 파일 읽는 법

xlwings tutorial - 데이터 계산하여 붙여 넣기