통계적 검정 (statistical testing) 은 모집단의 모수 또는 분포 형태에 대한 추정에 대해 그것이 옳은지 그른지를 임의로 추출한 표본으로부터 통계량을 측정하여 판단하는 통계적 절차를 말합니다.

 

단일 모집단에 대한 통계적 추론 (추정과 검정) 과 관련하여

 

[표본이 크고 정규성 충족 시]

- 단일 모집단의 모평균에 대한 신뢰구간 추정과 검정

   : t.test()

 

- 단일 모집단의 모분산에 대한 신뢰구간 추정과 검정

   : chi-square test

 

- 단일 모집단의 모비율에 대한 신뢰구간 추정과 검정

   : prop.test()

 

 

[정규성 미충족 시]

- 단일 모집단 중심에 대한 비모수 검정

  (Nonparametric test on one sample with median)

  : wilcox.test()

 

 

[정규성 여부 검정]

- 단일 모집단 분포의 정규성 검정

  : shapiro.test(), qqnorm(), qqline()

 

을 차례로 살펴보겠습니다.

 

지난번 포스팅의 정규분포 형태를 띠는 '단일 모집단의 모평균, 모분산, 모비율에 대한 신뢰구간 추정과 검정'에 이어서, 이번 포스팅에서는 정규분포 형태를 띠지 않는 단일 모집단의 비모수 검정에 대해 알아보고, R의 wilcox.test() 함수를 사용해 예를 들어보겠습니다.

 

표본의 개수가 많은 경우에는 중심극한의 정리에 의거해서 정규분포로 근사하게 되어 문제가 안되는데요, 표본의 개수가 작은 경우에는 정규성을 충족시키지 못하는 문제가 발생할 수 있어 t-검정을 사용할 수 없게 되며 이럴 때 사용하는 것이 비모수 검정이 되겠습니다. 혹은 모집단이 정규분포를 띠는지 아닌지를 알지 못하는 경우에도 비모수 검정 (nonparametric test)를 사용합니다.

 

모수 검정과 비모수 검정법은 아래 비교표처럼 서로 짝을 지어서 보면 좀더 쉽게 이해할 수 있을 것입니다.

 

 

[ 모수 검정과 비모수 검정 비교표 ]

 

 구분

 모수 검정

(Parametric Test)

비모수 검정

(Nonparametric Test) 

 When to use

 정규분포 가정 만족 시

 (normal distribution)

 정규분포 가정 불충족 시,

 혹은 모집단이 어떤 분포를 따르는지 모를 때

 (non-normal distribution, or un-know distribution,

  or very small sample size, or rankded data)

Statisctic

 평균 (mean)

 중앙값 (median) 

 1 sample

 1 sample t-test 

 1 sample Wilcoxon signed rank test

 2 samples

 2 sample t-test

 Mann-Whitney test

 paired 2-sample t-test

 Wilcoxon signed rank test

 more than

2 samples

 one-way ANOVA

 Kruskal-Wallis test

 

 

단일 모집단에 중심에 대한 비모수 검정에는 Wilcoxon signed rank 검정을 실시하면 되며, 아래에 R의 wilcox.test() 함수를 사용한 예를 소개하였습니다.

 

 

 

 

 


 

 

문제 ) 2010년도에 서울 지역의 고등학교 1학년 남학생의 몸무게를 전수 조사하였더니 평균이 63.0 kg, 표준편차가 8.0 kg이 나왔다.  2015년에 15명의 남학생을 무작위로 표본을 추출하여 몸무게를 재어보니 {70.2, 54.9, 67.0, 60.5, 63.4, 61.9, 71.8, 66.1, 72.6, 73.0, 68.7, 70.3, 66.2, 55.6, 65.9} 와 같이 나왔을 때 서울 지역 고등학교 1학년 남학생의 몸무게가 5년 전보다 더욱 증가했다고 볼 수 있는가?

 

(대립가설 1, 단측 검정 "greater") 서울 지역 고등학교 1학년 남학생의 몸무게가 5년 전보다 더욱 증가하였다는 대립가설을 검정하시오.

 

H0 : 서울 지역 고등학교 1학년 남학생의 몸무게는 5년 전보다 증가하지 않음

H1 : 서울 지역 고등학교 1학년 남학생의 몸무게가 5년 전보다 더욱 증가

 

 

> ##----------------------------------------------------------
> ## 단일 모집단의 모평균에 대한 비모수 검정: wilcox.test()
> ##----------------------------------------------------------
> 
> # x : random sample of 15 students's weight in high school
> x <- c(70.2, 54.9, 67.0, 60.5, 63.4, 61.9, 71.8, 66.1, 72.6, 73.0, 68.7, 70.3, 66.2, 55.6, 65.9)
> mean(x)
[1] 65.87333
> var(x); sd(x)
[1] 32.54495
[1] 5.704818
> 
> stem(x) # stem-and-leaf plot

  The decimal point is 1 digit(s) to the right of the |

  5 | 
  5 | 56
  6 | 123
  6 | 66679
  7 | 00233

> 
> 
> # wilcox signed rank test, one-sided test : "greater"
> wilcox.test(x, # weight vector for wilcox-test
+             alternative = c("greater"), #  alternative = c("less", "greater", "two-sided")
+             mu = 63.0, # mu of population
+             conf.int = TRUE) # 95% pseudomedian confidence interval

	Wilcoxon signed rank test

data:  x
V = 92, p-value = 0.0365
alternative hypothesis: true location is greater than 63
95 percent confidence interval:
 63.35   Inf
sample estimates:
(pseudo)median 
        66.175

 

 

위의 비모수 검정 결과 P-value가 0.0365 로서 5% 유의수준 하에 귀무가설을 기각하게 되고 대립가설을 채택하게 되어 5년 전보다 고등학교 남학생의 몸무게가 증가했다고 판단할 수 있겠습니다.

 

 

 

t-검정(t-Test)와 비교해보면 아래와 같습니다. 결과는 크게 다르지 않게 나왔습니다.

 

> # t-test, one-sided test : "greater"
> t.test(x, # weight vector for t-test
+        alternative = c("greater"), #  alternative = c("less", "greater", "two-sided")
+        mu = 63.0, # mu of population
+        conf.level = 0.95) # confidence level or confidence coefficient (1-α)

	One Sample t-test

data:  x
t = 1.9507, df = 14, p-value = 0.0357
alternative hypothesis: true mean is greater than 63
95 percent confidence interval:
 63.27896      Inf
sample estimates:
mean of x 
 65.87333

 

 

 

 

Wilcoxon signed rank test 의 lower confidence limit, upper confidence limit 을 indexing 해오는 방법은 아래와 같습니다.

 

> # indexing of 95% confidence level value > wilcox.test_confi_95 <- wilcox.test(x, alternative = c("greater"), mu = 63.0, conf.int = TRUE) > > names(wilcox.test_confi_95) # statistics [1] "statistic" "parameter" "p.value" "null.value" "alternative" "method" "data.name" [8] "conf.int" "estimate" > > wilcox.test_confi_95$conf.int # confidence interval at 95% confidence level [1] 63.35 Inf attr(,"conf.level") [1] 0.95 > > wilcox.test_confi_95$conf.int[1] # lower confidence limit [1] 63.35 > wilcox.test_confi_95$conf.int[2] # upper confidence limit [1] Inf

 

 

많은 도움 되었기를 바랍니다.

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

Posted by R Friend Rfriend

댓글을 달아 주세요

통계적 검정 (statistical testing) 은 모집단의 모수 또는 분포 형태에 대한 추정에 대해 그것이 옳은지 그른지를 임의로 추출한 표본으로부터 통계량을 측정하여 판단하는 통계적 절차를 말합니다.

 

단일 모집단에 대한 통계적 추론 (추정과 검정) 과 관련하여

 

[표본이 크고 정규성 충족 시]

- 단일 모집단의 모평균에 대한 신뢰구간 추정과 검정

   : t.test()

 

- 단일 모집단의 모분산에 대한 신뢰구간 추정과 검정

   : chi-square test

 

- 단일 모집단의 모비율에 대한 신뢰구간 추정과 검정

   : prop.test()

 

 

[정규성 미충족 시]

- 단일 모집단 중심에 대한 비모수 검정

  : wilcox.test()

 

 

[정규성 여부 검정]

- 단일 모집단 분포의 정규성 검정

  : shapiro.test(), qqnorm(), qqline()

 

을 차례로 살펴보겠습니다.

 

지난번 포스팅의 '단일 모집단의 모평균, 모분산에 대한 신뢰구간 추정과 검정'에 이어서, 이번 포스팅에서는 정규분포 형태를 띠는 단일 모집단에서 충분히 큰 규모로 임의로 표본을 추출하여 표본비율(one sample proportion)을 분석하여 미지의 모수인 모비율(one population proportion)에 대한 95% 신뢰계수의 신뢰구간 추정과 검정을 R의 prop.test() 함수를 사용해보겠습니다.

 

이번 포스팅은 선거철 특정 정당이나 후보에 대한 "지지율"(지지 여부 yes, no)이라든지, 특정 제품의 "불량률"(양품 여부 yes, no)과 같은 모집단 내 개체들의 특정 속성에 대한 구성 비율에 추정과 검정에 사용합니다.

 

 

[ 단일 모집단의 모비율에 대한 검정통계량 및 대립가설 형태별 P-value ]

 

 

확률 변수 X가 모수 n (시행 횟수), p (성공 확률) 인 이항분포 (binomial distribution)을 따를 때

모비율 p의 추정량으로는 표본비율 p^ = X / n 을 사용함.

 

표본비율 P^ 은 중심극한의 정리에 의거하여 n이 충분히 크면 평균이 p, 분산이 p(1-p)/n 인

정규분포로 근사하게 됨.

 

따라서, 모비율 p에 대한 검정을 위해 사용하는 통계량 및 대립가설 형태별 P-value는 다음과 같음

 

 

 

 

문제) 100원짜리 동전 던지기를 1000번 했는데 앞면이 485번 나왔다.  그렇다면 이 동전은 앞면과 뒷면이 균일한 동전이라고 말할 수 있는지를 유의수준 5%로 검정하고, 신뢰계수 95%의 신뢰구간을 구하여라.

 

> ##------------------------------------------------
> ## one population proportion test : prop.test()
> ##------------------------------------------------
> 
> prop.test(x = 485, # number of success
+           n = 1000, # sample size
+           p=0.50, # proportion of success
+           
+           alternative = c("two.sided"), # two-sided test
+           # alternative = c("greater"), # right-sided test
+           # alternative = c("less"),# left-sided test
+           
+           conf.level = 0.95) # confidence level

	1-sample proportions test with continuity correction

data:  485 out of 1000, null probability 0.5
X-squared = 0.841, df = 1, p-value = 0.3591
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.4536437 0.5164729
sample estimates:
    p 
0.485 

> 
> 
> prop.test_confi_95 <- prop.test(x = 485, n = 1000, p=0.50, alternative = c("two.sided"), conf.level = 0.95)
> 
> names(prop.test_confi_95)  # statistics
[1] "statistic"   "parameter"   "p.value"     "estimate"    "null.value"  "conf.int"    "alternative"
[8] "method"      "data.name"  
> 
> prop.test_confi_95$conf.int  # confidence interval at 95% confidence level
[1] 0.4536437 0.5164729
attr(,"conf.level")
[1] 0.95
> 
> prop.test_confi_95$conf.int[1] # lower confidence limit
[1] 0.4536437
> prop.test_confi_95$conf.int[2] # upper confidence limit
[1] 0.5164729

 

 

P-value가 0.3591 이므로 유의수준(significance level) 5%에서 귀무가설을 기각하지 못하고 채택하게 됩니다. 즉, 동전은 앞, 뒤가 균일하다고 말할 수 있겠습니다.

 

그리고 신뢰계수 95%의 신뢰구간은 0.453 ~ 0.516 가 되겠습니다. 즉, 1000번 동전던지기를 했는데, 453번 보다 앞면이 덜 나오거나, 혹은 516번 보다 앞면이 더 나오면 동전의 앞, 뒷가 균일하지 않고, 누군가가 사기의 목적을 가지고 조작을 했다고 유의수준 5% 기준 하에 의심해 볼 수 있겠습니다.

 

prop.test() 함수 말고 분석가가 직접 위의 표본비율 p^을 표준화한 통계량 Z를 직접 사용자 정의함수로 짜서 검정과 95% 신뢰구간 추정을 해도 동일한 결과가 나올겁니다.

 

많은 도움 되었기를 바랍니다.

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

Posted by R Friend Rfriend

댓글을 달아 주세요

통계적 검정 (statistical testing) 은 모집단의 모수 또는 분포 형태에 대한 추정에 대해 그것이 옳은지 그른지를 임의로 추출한 표본으로부터 통계량을 측정하여 판단하는 통계적 절차를 말합니다.

 

단일 모집단에 대한 통계적 추론 (추정과 검정) 과 관련하여

 

[표본이 크고 정규성 충족 시]

- 단일 모집단의 모평균에 대한 신뢰구간 추정과 검정

   : t.test()

 

- 단일 모집단의 모분산에 대한 신뢰구간 추정과 검정

   : chi-square test

 

- 단일 모집단의 모비율에 대한 신뢰구간 추정과 검정

   : prop.test()

 

 

[정규성 미충족 시]

- 단일 모집단 중심에 대한 비모수 검정

  : wilcox.test()

 

 

[정규성 여부 검정]

- 단일 모집단 분포의 정규성 검정

  : shapiro.test(), qqnorm(), qqline()

 

을 차례로 살펴보겠습니다.

 

지난번 포스팅의 '단일 모집단의 모평균에 대한 신뢰구간 추정과 검정'에 이어서, 이번 포스팅에서는 정규분포 형태를 띠는 단일 모집단에서 임의로 표본을 추출하여 미지의 모수인 모분산(one population variance)에 대한 95% 신뢰계수의 신뢰구간 추정과 검정을 R의 사용자 정의 함수를 적용해서 해보겠습니다.

 

제조업에서 품질을 평가하거나 금융업에서 Risk를 평가할 때 많이 사용하는 것이 퍼짐 정도를 나타내는 분산이므로 기준치, 혹은 허용치를 벗어나는지를 판단(모분산과 표본분산이 동일한가 여부)하는데 이번 포스팅을 활용할 수 있겠습니다.

 

연속형 확률분포 중에서 이전에 '카이제곱 분포'를 포스팅한적이 있는데요, 단일 모집단의 모비율 추정과 검정에 '카이제곱 분포 통계량'을 사용합니다.

 

 

정규분포 (μ, σ^2)을 따르는 모집단으로부터 표본크기 n인 표본을 랜덤하게 추출하면

 

  V = (n-1)* (S/σ)^2 인 통계량은

 

자유도가 n-1 인 카이제곱 분포(chisquare distribution)를 따름

 

 

위의 V 통계량에서 (S/σ)^2  부분에 주목해서 보면, 표본분산을 모분산으로 나눈 값을 제곱하였으므로 모분산 대비 표본분산이 크거나 작게 되면 1로부터 멀어지게끔 되어있습니다.

 

그리고, 참고로 카이제곱 분포를 따른다는 것은 정규분포나 t-분포처럼 좌우 대칭 형태의 분포가 아니라, 왼쪽으로 많이 몰려있고, 오른쪽으로 꼬리가 긴 분포를 따른다는 의미입니다.  따라서 양측검정할 때는 정규분포의 z통계량이나 t-분포의 T통계량을 사용할 때 처럼 '±통계량'의 형태가 아니구요, 하한 신뢰한계(lower confidence limit)와 상한 신뢰한계(upper confidence limit)을 따로 따로 계산해야만 합니다. 

 

또한 t-검정의 경우 모집단으로 부터 추출한 샘플들이 정확하게 정규분포를 따르지 않더라도 좌우 대칭이며 특이값만 없다면 검정 결과는 덜 영향을 받는 반면에, 카이제곱 검정은 정규성 가정에 매우 민감하게 영향을 받으므로 정규분포를 따르는지에 대해서 반드시 확인 후에 조심해서 사용해야 하겠습니다.

 

 

 

[ 카이제곱 검정 (Chi-square test) 의 검정 통계량 및 기각역 ]

 

 

 

 

 

문제) 빛나리 공장에서 생산하는 전구 수명은 정규분포를 따르는 것으로 알려져 있다.  빛나리 공장 전구의 분산이 30시간 이내이면 품질 기준을 만족한다.  완성품 전구에서 100개를 무작위로 추출하여 전구의 수명을 측정한 결과 분산이 38시간 이었다.  그러면, 빛나리 공장에서 생산한 전구의 수명의 분산이 기준치 30시간을 초과하는지 여부를 유의수준(significance level) 5%로 검정하여라.

 

 

위의 문제는 분산에 대해 우측검정(one-sided test, greater) 을 하는 문제이므로 카이제곱 검정을 하는 R 사용자 정의함수를 만들어보았습니다.  아래 R 사용자 정의 함수에 샘플 개수(n), 모분산(Sigma0Squared), 표본분산(SSquared), 유의수준(alpha) 부분을 다른 문제의 값으로 바꿔치기 해서 사용하시면 되겠습니다.

 

 

> # function for right-sided chi-square test with significance level 'alpha'
> var_test_TwoSided <- function(n, Sigma0Squared, SSquared, alpha)
+ {
+   df <- n - 1
+   v <- df*(SSquared)/Sigma0Squared
+   
+   upper.critical <- qchisq((1-alpha), df, lower.tail = FALSE) # right-sided test
+ 
+   # lower.critical <- qchisq((alpha), df, lower.tail = FALSE) # left-sided test
+   
+   # upper.critical <- qchisq((1-alpha)/2, df, lower.tail = FALSE) # two-sided test
+   # lower.critical <- qchisq(alpha/2, df, lower.tail = FALSE) # two-sided test
+   
+   print(paste("degrees of freedom = ", df), quote = FALSE)
+   print(paste("population variance = ", Sigma0Squared), quote = FALSE)
+   print(paste("sample variance = ", SSquared), quote = FALSE) 
+   print(paste("significance level = ", alpha), quote = FALSE)
+   print(paste("confidence level = ", 1-alpha), quote = FALSE)
+   print("  ")
+   print(paste((1-alpha)*100, "% confidence interval for variance"), quote=FALSE)
+   print(paste("   upper critical limit = ", round(upper.critical, 2)), quote = FALSE)
+   print(paste("   chisq statistic      = ", round(v, 2)), quote = FALSE)
+   print("  ")
+   print(paste("P-value =", round(pchisq(v, df, lower.tail=FALSE),4)), quote = FALSE)
+ }
> 
> 
> var_test_TwoSided(n=100, Sigma0Squared=30, SSquared=38, alpha=0.05)
[1] degrees of freedom =  99
[1] population variance =  30
[1] sample variance =  38
[1] significance level =  0.05
[1] confidence level =  0.95
[1] "  "
[1] 95 % confidence interval for variance
[1]    upper critical limit =  77.05
[1]    chisq statistic      =  125.4
[1] "  "
[1] P-value = 0.0377

 

 

P-value 가 0.0377 이므로 유의수준(significance level) 0.05에서 귀무가설은 기각(H0 is rejected), 대립가설 채택(H1 is accepted) 되었습니다.  따라서 빛나리 전구회사의 전구는 분산이 기준치 30시간을 유의수준 5%에서 초과했으므로 시장에 내다팔 수가 없게 되었네요.

 

많은 도움 되었기를 바랍니다.

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

 

 

Posted by R Friend Rfriend

댓글을 달아 주세요

통계적 검정 (statistical testing) 은 모집단의 모수 또는 분포 형태에 대한 추정에 대해 그것이 옳은지 그른지를 임의로 추출한 표본으로부터 통계량을 측정하여 판단하는 통계적 절차를 말합니다.

 

단일 모집단에 대한 통계적 추론 (추정과 검정) 과 관련하여

 

[표본이 크고 정규성 충족 시]

- 단일 모집단의 모평균에 대한 신뢰구간 추정과 검정

   : t.test()

 

- 단일 모집단의 모분산에 대한 신뢰구간 추정과 검정

 

   : chisq test

 

- 단일 모집단의 모비율에 대한 신뢰구간 추정과 검정

   : prop.test()

 

 

[정규성 미충족 시]

- 단일 모집단 중심에 대한 비모수 검정

  : wilcox.test()

 

 

[정규성 여부 검정]

- 단일 모집단 분포의 정규성 검정

  : shapiro.test(), qqnorm(), qqline()

 

을 차례로 살펴보겠습니다.

 

이번 포스팅에서는 정규분포 형태를 띠는 단일 모집단에서 임의로 많은 갯수의 표본을 추출하여 미지의 모수인 모평균 μ 에 대한 95% 신뢰계수의 신뢰구간 추정과 검정을 R의 t.test() 함수를 활용해서 해보겠습니다.

(표본의 크기가 충분히 크면 중심 극한의 정리에 의거하여 정규분포를 가정)

 

 

먼저, 기본 용어부터 살펴보자면

 

[ 귀무가설 (null hypothesis: H0), 대립가설 (alternative hypothesis : H1) ]

 

 

  • 귀무가설 (null hypothesis: H0) : 현재까지 참이라고 알려진 사실

  • 대립가설 (alternative hypothesis : H1) : 귀무가설에 반대되는 가설로, 연구자가 표본자료로 부터 입증하고자 하는 주장, 추정, 가설

 

 

 

 

[ 대립가설의 형태에 따른 검정의 분류 (Types of alternative hypothesis) ]

 

 

  • 단측 검정 (One-sided test) : H1 : population mean > (greater) or < (less) than specific value
    (greater : 우측 검정, less : 좌측 검정)
  • 양측 검정 (Tow-sided test) : H1 : population mean =! (not equal) with specific value

 

 

 

 

[ 검정 통계량 (test statistic)과 P-value]

  • 검정 통계량 (test statistic) : 귀무가설의 채택 또는 기각 여부를 판단하는데 사용하는 통계량
  • 표본자료로부터 관측된 모수의 추정값이 귀무가설 하에서 일어나기 힘든 값, 매우 희귀한 값이면 귀무가설을 기각하고 대립가설을 채택하는 원리
  • 임의적 판단이 아닌, 검정 통계량의 확률분포와 유의 수준(significance level)를 이용해 판단
  • P-value : 귀무가설이 옳다는 가정 하에서 표본으로부터 계산된 검정 통계량의 값보다 더 극단적인 결과가 나올 확률로서, P-value가 유의 수준보다 작으면 귀무가설을 기각하고 대립가설을 채택

 

---------------------------------------------------------------------

 

문제 하나를 가지고 예를 들어서 R t.test() 함수를 활용해 추정과 검정 해보겠습니다.

 

 

문제 ) 2010년도에 서울 지역의 고등학교 1학년 남학생의 몸무게를 전수 조사하였더니 평균이 63.0 kg, 표준편차가 8.0 kg이 나왔다.  2015년에 15명의 남학생을 무작위로 표본을 추출하여 몸무게를 재어보니 {70.2, 54.9, 67.0, 60.5, 63.4, 61.9, 71.8, 66.1, 72.6, 73.0, 68.7, 70.3, 66.2, 55.6, 65.9} 와 같이 나왔을 때 서울 지역 고등학교 1학년 남학생의 몸무게가 5년 전보다 더욱 증가했다고 볼 수 있는가?

 

(대립가설 1, 단측 검정 "greater") '서울 지역 고등학교 1학년 남학생의 몸무게가 5년 전과 동일하다'는 귀무가설(H0), '서울 지역 고등학교 1학년 남학생의 몸무게가 5년 전보다 더욱 증가하였다'는 대립가설(H1)을 검정하시오.

 

  • H0 : 서울 지역 고등학교 1학년 남학생의 몸무게는 5년 전보다 증가하지 않음
  • H1 : 서울 지역 고등학교 1학년 남학생의 몸무게가 5년 전보다 더욱 증가

 

 

# Normal distribution plot, X~N(63, 8^2), right-sided test
x1 <- c(33:93)
plot(x1, dnorm(x1, mean=63, sd=8), type='l',
     main="Normal distribution, X~N(63,8^2), right-sided test")

abline(v=63, col="blue", lty=3)
abline(v=63 + 1.96*8, col="red", lty=2)
text(82, 0.003, labels = "------->")

 

 

 

 

(1) 단측 검정 (one-sided test) 

 

 
> ##----------------------------------------------------------
> ## 단일 모집단의 모평균에 대한 신뢰구간 추정과 검정: t.test()
> ##----------------------------------------------------------
> 
> # x : random sample of 15 students's weight in high school
> x <- c(70.2, 54.9, 67.0, 60.5, 63.4, 61.9, 71.8, 66.1, 72.6, 73.0, 68.7, 70.3, 66.2, 55.6, 65.9)
> mean(x)
[1] 65.87333
> var(x); sd(x)
[1] 32.54495
[1] 5.704818
> 

> stem(x) # stem-and-leaf plot The decimal point is 1 digit(s) to the right of the | 5 | 5 | 56 6 | 123 6 | 66679 7 | 00233

> 
> 
> # one-sided test : "greater"
> t.test(x, # weight vector for t-test
+        alternative = c("greater"), #  alternative = c("less", "greater", "two-sided")
+        mu = 63.0, # mu of population
+        conf.level = 0.95) # confidence level or confidence coefficient (1-α)

	One Sample t-test

data:  x
t = 1.9507, df = 14, p-value = 0.0357
alternative hypothesis: true mean is greater than 63
95 percent confidence interval:
 63.27896      Inf
sample estimates:
mean of x 
 65.87333

 

 

(해석) 위 t.test()를 활용한 분석 결과 t 통계량은 1.9507, P-value 는 0.0357로서 유의수준 (significance level) 5% 하에서는 귀무가설이 기각, 대립가설이 채택되어 "5년 전보다 남학생들 몸무게가 증가"하였다고 판단할 수 있습니다.

 

P-value가 0.0357 이므로 15개의 임의 표본을 추출해서 100 번 반복 조사를 해보면 약 3번~4번 밖에 위와 같은 경우가 발생한다는 뜻이므로, 귀무가설 하에서는 일어나기에 쉽지 않은 경우라고 하겠지요?

 

95% 신뢰구간은 아래 값처럼 (63.27896 ~ 무한대)가 되겠네요.  단측검정으로 우측검정을 했기 때문에 이런 신뢰구간 값이 나온겁니다. 

 

 95 percent confidence interval:
63.27896      Inf

 

5년 전의 모집단의 몸무게 평균이 63.0 kg 이라고 했는데요, 5년 후에 조사한 표본집단의 몸무게 평균의 95% 신뢰구간이 (63.27896~무한대) 이므로 서로 겹치지 않습니다.  그만큼 귀무가설에서는 일어나기 힘든 일이 벌어진 것이므로, 대립가설("greater")을 채택하게 됩니다.

 

 

---------------------------------------------------------------------

 

(대립가설 2, 양측 검정 "two.sided))  서울 지역 고등학교 1학년 남학생의 몸무게는 5년 전과 동일하다는 대립가설을 검정하시오.

 

H0 : 서울 지역 고등학교 1학년 남학생의 몸무게는 5년 전과 동일

H1 : 서울 지역 고등학교 1학년 남학생의 몸무게가 5년 전과 다름 (증가 혹은 감소)

 

 

# Normal distribution plot, X~N(63, 8^2), two-sided test
x1 <- c(33:93)
plot(x1, dnorm(x1, mean=63, sd=8), type='l',
     main="Normal distribution, X~N(63,8^2), two-sided test")

abline(v=63, col="blue", lty=3)
abline(v=63 - 1.96*8, col="red", lty=2)
abline(v=63 + 1.96*8, col="red", lty=2)

text(82, 0.003, labels = "------->")
text(44.5, 0.003, labels = "<-------")

 

 

 

 

 

 

 

(2) 양측 검정 (two-sided test) 

 

 
> # two-sided test : "two.sided"
> t.test(x, # weight vector for t-test
+        alternative = c("two.sided"), #  alternative = c("less", "greater", "two-sided")
+        mu = 63.0, # mu of population
+        conf.level = 0.95) # confidence level or confidence coefficient (1-α)

	One Sample t-test

data:  x
t = 1.9507, df = 14, p-value = 0.0714
alternative hypothesis: true mean is not equal to 63
95 percent confidence interval:
 62.71411 69.03256
sample estimates:
mean of x 
 65.87333

 

 

위 양측검정 (two-sided testing) 결과 P-value 가 0.0714 이므로 유의수준 (significance level) 5% 하에서는 귀무가설을 채택(대립가설 기각, 즉 5년 전과 몸무게 평균 같다)하고, 유의수준 10% 하에서는 귀무가설을 기각(대립가설 채택, 즉 5년 전과 몸무게 평균 다르다)하게 됩니다.

 

 

 

---------------------------------------------------------------------

 

(3) 95% 신뢰구간 값 구하기 (indexing of 95% confidence level value)

 

 

> # indexing of 95% confidence level value
> t.test_confi_95 <- t.test(x, alternative = c("two.sided"), mu = 63.0, conf.level = 0.95) 
> 
> 
> names(t.test_confi_95) # statistics
[1] "statistic"   "parameter"   "p.value"     "conf.int"    "estimate"    "null.value"  "alternative"
[8] "method"      "data.name"  
> 
> t.test_confi_95$conf.int # confidence interval at 95% confidence level
[1] 62.71411 69.03256
attr(,"conf.level")
[1] 0.95
> 
> t.test_confi_95$conf.int[1] # lower limit
[1] 62.71411
> t.test_confi_95$conf.int[2] # upper limit
[1] 69.03256

 

 

 

많은 도움 되었기를 바랍니다.

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 


Posted by R Friend Rfriend

댓글을 달아 주세요

  1. ㅇㅇㄹ 2017.10.19 23:33  댓글주소  수정/삭제  댓글쓰기

    양측검정이니까 유의수준 10%하에서도 p값이 0.05보다 작아야 하므로 귀무가설 기각되지 않나요?

  2. 감사합니다! 2017.12.12 17:09  댓글주소  수정/삭제  댓글쓰기

    정말로 차근차근 잘 설명해주셔서 덕분에 잘 공부하고 있습니다! 좋은 자료 제공해주셔서 진심으로 감사드려요! ㅎㅎㅎ

  3. yes89929@gmail.com 2018.01.08 21:39  댓글주소  수정/삭제  댓글쓰기

    # Normal distribution plot, X~N(63, 8^2), two-sided test
    위 내용을 보면
    abline(v=63 - 1.96*8, col="red", lty=2)에서
    Z(시그마 / 루트 n)을 안하셔서 제가 알고있는 신뢰구간 식과는 다른데
    그림을 제가 원하는 형태가 나왔어요.

    왜 루트n을 계산하지 않았는지, 그리고 루트n을 계한하지 않아도 되는지 알려주시면 감사하겠습니다.

    추가로 95% t 검정을 했을때는 62~69가 나오는데 위에 값가고 많이 달라서 왜그런지도 알려주시면 감사하겠습니다.

  4. aa 2018.11.05 23:04  댓글주소  수정/삭제  댓글쓰기

    P-value : 귀무가설이 옳다는 가정 하에서 표본으로부터 계산된 검정 통계량의 값을 기준으로 대립가설 방향으로의 확률로서, P-value가 유의 수준보다 작으면 귀무가설을 기각하고 대립가설을 채택함

    라고 되어있는데 '대립 가설 방향으로의 확률' 보다는 '극단값이 나올 확률' 이 아닌가요??

  5. 세진컴퓨터랜드 2018.12.18 18:08  댓글주소  수정/삭제  댓글쓰기

    와 정말 좋은 포스팅 잘 읽고 갑니다 !!
    감사합니다 !!

  6. 배우는 자 2019.12.27 11:24  댓글주소  수정/삭제  댓글쓰기

    (1) 대립가설 2번에 H0, H1 설정이 반대로 되었어요~
    - 5년 전과 동일하다는 것이 대립가설이라고 했으므로.


    (2) 아래 본문 내용에 이해하기 쉽게 수정이 필요할것 같습니다~

    "위 양측검정 (two-sided testing) 결과 P-value 가 0.0714 이므로 유의수준 (significance level) 5% 하에서는 귀무가설을 채택(대립가설 기각, 즉 5년 전과 몸무게 평균 같다)하고, 유의수준 10% 하에서는 귀무가설을 채택(대립가설 기각, 즉 5년 전과 몸무게 평균 다르다)하게 됩니다."

    ->

    "위 양측검정 (two-sided testing) 결과 P-value 가 0.0714 이므로 유의수준 (significance level) 10% 하에서는 귀무가설을 채택(대립가설 기각, 즉 5년 전과 몸무게 평균 다르다)하게 됩니다."

    • R Friend Rfriend 2019.12.27 22:44 신고  댓글주소  수정/삭제

      안녕하세요.

      귀무가설(H0)는 5년 전과 동일,
      대립가설(H1)은 5년 전보다 증가

      입니다. 본문 내용 중에 혼선이 있을 수 있는 부분을 보다 명확하게 수정하였습니다.

  7. HG 2021.01.10 21:37  댓글주소  수정/삭제  댓글쓰기

    stem(x)가 무엇을 뜻하는건지 모르겠습니다.
    stem(x,scale=2)로 해야 되는거 아닌가요?
    54 | 96
    56 |
    58 |
    60 | 59
    62 | 4
    64 | 9
    66 | 120
    68 | 7
    70 | 238
    72 | 60

    • R Friend Rfriend 2021.01.11 00:05 신고  댓글주소  수정/삭제

      안녕하세요.

      stem() 함수는 줄기 잎 그림 (Stem-and-leaf Plot) 을 그릴때 사용합니다.

      결과에서
      The decimal point is 1 digit(s) to the right of the |
      5 |
      5 | 56 ==> 55, 56
      6 | 123 ==> 61, 62, 63
      6 | 66679 ==> 66, 66, 66, 67, 69
      7 | 00233 ==> 70, 70, 72, 73, 73

      을 의미합니다. 데이터의 분포 형태를 가늠해볼 수 있습니다.

이전 포스팅들에서 이산형 확률분포, 연속형 확률분포에 대해서 R을 가지고 그래프도 그려보고 확률, 분위수 계산, 난수 생성하는 방법을 모두 살펴보았습니다.

 

확률분포, 확률, 추정과 검정을 좀더 들어가게 되면 신뢰구간(confidence interval), 신뢰계수(confidence coefficient), 유의수준(significance level), 신뢰한계(conficence limits) 등과 맞부딪치게 되는데요, 비슷비슷한 용어들을 알듯 하면서도 정확히 모르고 있는 분들이 많은 통계용어가 아닌가 합니다.

 

통계의 (1) 큰 줄기 하나가 모집단에서 표본을 추출해서 주요 통계량을 가지고 표본의 특성(중심화 경향, 퍼짐 정도, 분포모양, 상관성 등)을 알아보는 기술통계(descriptive statistics)가 있고, (2) 또 하나의 큰 줄기가 이렇게 파악한 표본집단의 기술통계량을 가지고 모집단의 모수(parameter)를 추정(estimation)하고 가설을 검정(test)하는 추론통계(inferential staticstics) 입니다. 

 

추론통계의 모집단 모수 추정 시, 가령 모집단의 평균과 분산 추정 시에 특정 값 하나를 명시하는 점 추정 (point estimation)을 하게 되면 틀릴 확률이 아주 아주 높게 됩니다.  그래서 대부분은 특정 값 하나 대신에 어느 값부터 어느 값 사이에 어느 수준의 확률로 모집단의 모수가 들어가 있을 거라고 추정하는 구간 추정 (interval estimation)을 하게 됩니다. 

 

그러면, 아래의 그림을 참고해가면서 개념을 하나씩 설명드리도록 하겠습니다.

 

 

 

1) 신뢰구간 (confidence interval) : 주어진 확률 1-α (신뢰계수)에 대하여 표본분포의 통계량이 모집단 모수에 포함되는 구간 (θ1 ~ θ2)

 

2) 신뢰계수 (confidence coefficient) : n개의 확률표본을 추출하여 신뢰구간을 구하는 잡업을 N번 반복하여 얻은 N개의 신뢰구간 중 (1-α)%에 미지의 모수 μ 가 포함되어 있을 확률

 

3) 유의수준 (significance level) : n개의 확률표본을 추출하여 신뢰구간을 구하는 잡업을 N번 반복하여 얻은 N개의 신뢰구간 중 미지의 모수 μ 가 포함되어 있지 않을 확률 (α)

 

4) 신뢰한계 (confidence limits) : 신뢰구간의 하한값(θ1)과 상한값(θ2)

 

 

위의 개념을 좀더 알기 쉽도록 하기 위해서, X ~ N(mean=10, sd=3)인 정규분포에서 20개(n)의 확률표본을 추출하여 신뢰구간을 구하는 작업을 100번(N) 번 반복 수행했을 때 얻은 100개(N)의 신뢰구간(confidence interval) 중에서 95% 의 신뢰계수(confidence coefficient, 1-α)만큼은 미지의 모수 μ (모집단의 평균) 가 포함되어 있고, 5%의 유의수준(significance level, α) 만큼은 미지의 모수 μ(모집단의 평균)가 포함되지 않는 경우를, R을 이용해서 simulation 해보도록 하겠습니다. (신뢰구간을 벗어난 case는 색깔, 선 형태가 다르게 그려지도록 프로그래밍 해놨습니다)

 

 

 

> ##----------------------------------------------------------
> ## 신뢰구간(confidence interval), 신뢰계수 (confidence coefficient)
> ##----------------------------------------------------------
> 
> 
> # confidence interval simulation, parameter setting
> alpha <- 0.05 # significance level
> N <- 100 # total simulation frequency
> n <- 20 # sampling number
> mu <- 10 # mean of population
> sigma <- 3 # standard deviation of population
> 
> # plot
> # 면 분할 
> layout(matrix(1:2), heights=c(1,2))
> 
> 
> # Normal distribution plot, X~N(10,3^2)
> x <- seq(0, 20, length=100)
> plot(x, dnorm(x, mean=10, sd=3), type='l', main="Normal distribution, X~N(10,3^2)")
> 
> abline(v=mu, col="red", lty=2)
> abline(v=mu-1.96*3, col="red", lty=2)
> abline(v=mu+1.96*3, col="red", lty=2)
> 
> 
> plot(c(0, 20), c(1, N), 
+      xlab = "confidence interval", 
+      ylab = "simulation freqneucy", 
+      type = 'n')
> 
> abline(v=mu, col="blue", lty=2)
> 
> 
> within_yn <- 0
> 
> for(i in 1 : N) {
+   x = rnorm(n, mu, sigma)
+   conf_limit_lower = mean(x) - qnorm(1-alpha/2)*sd(x)/sqrt(n) # lower confidence limit
+   conf_limit_upper = mean(x) + qnorm(1-alpha/2)*sd(x)/sqrt(n) # upper confidence limit
+   within_yn_eval = mu <= conf_limit_upper & mu >= conf_limit_lower
+   if(within_yn_eval) within_yn = within_yn + 1
+   segments(conf_limit_lower, i, conf_limit_upper, i, col=(! within_yn_eval)+1, lty=(! within_yn_eval)+1)
+   # segments(x1, y1, x2, y2) ; line from (x1, y1) to (x2, y2)
+ }
> 
> mtext(paste("confidence coefficient(1-alpha) = ", within_yn/N), side=3, line = 0, outer = FALSE, cex = 1.3)

 

 

 

 

 

 

 

실생활 예시로서, 만약 아래와 같은 선거 지지율 설문조사가 있다고 했을 때, 이를 해석해보겠습니다.

 

“A조사기관에 따르면 이달 1~3일 전국 19세 이상 성인 남녀 1000명을 대상으로 한 설문조사 결과, ‘A’ 후보자에 대한 지지율은 49%, ‘B’ 후보자의 지지율은 51%로 나타났다. 이번 조사는 신뢰수준 95%, 오차는 ±3.0%포인트다.”

 

신뢰수준 95%는 신뢰계수(confidence coefficient)와 같은 의미이고, +, - 오차범위는 신뢰한계(confidence limits)와 같은 의미로 보면 되겠습니다.  따라서, 위 조사결과의 의미는 동일한 조사를 100번 반복하게 되면, 

'A' 후보자의 지지율(모수 추정치)이 46%~52%, 'B' 후보자의 지지율(모수 추정치)이 48%~54%의 신뢰구간으로 조사 결과가 나올 횟수가 95번 (확률이 95%)이라는 의미입니다.  바꿔 말하면, 다른 결과가 나올 확률이 5% 라는 뜻입니다. 

 

95% 신뢰수준보다는 99%신뢰수준이 좋기는 합니다만, 신뢰수준을 높이려면 표본의 개수를 높여야 하는데요, 표본 = 돈, 시간 이므로 보통은 신뢰수준을 특정 값으로 고정 시켜놓고 샘플 수를 적정한 타협을 보게 됩니다.

 

많은 도움 되셨기를 바랍니다.

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

Posted by R Friend Rfriend

댓글을 달아 주세요

연속형 확률분포 (Continuous probability distribution)에는

 

 - 정규분포 (normal distribution)

   : norm()

 

 - 균등분포 (uniform distribution)

   : unif()

 

  - 지수분포 (exponential distribution)

   : exp()

 

 - t-분포 (t-distribution)

   : t()

 

 - F-분포 (F-distribution)

   : f()

 

  - 카이제곱분포(chisq-distribution)

   : chisq()

 

등이 있습니다.  

 

표본분포를 나타낼 때 t-분포, F-분포, 카이제곱 분포 등을 사용하는데요, 이번 포스팅에서는 이중에서 카이제곱 분포 (chi-squared distribution)에 대해서 소개하도록 하겠습니다.

 

카이제곱 분포 (chi-squared distribution)은 k개의 서로 독립적이고 표준정규분포를 따르는 확률변수 X 를 제곱한 값들을 합하였을 때의 분포이며, 이때 k 는 자유도 (degrees of freedom) 로서 카이제곱 분포의 parameter 가 됩니다.

 

 

 

 

위에서 카이제곱 분포는 표본분포 중의 하나이며, 정규분포와도 관련이 있다고 했는데요, 카이제곱 분포 (chi-squared distribution)은 모집단의 분산에 대한 추정과 검정에 활용됩니다.  그리고 F-분포와도 관련이 있는데요, F-분포가 카이제곱 분포를 따르는 두 개의 확률변수를 가지고 자유도 2개를 parameter로 사용하는 분포이구요, 카이제곱 분포는 정규분포를 따르는 확률변수를 가지고 제곱해서 합한 값의 분포로서 자유도 1개만을 parameter로 사용하는 분포가 되겠습니다. 

 

 

R을 활용한 카이제곱 분포 함수 및 parameter는 다음과 같습니다. 그래프 모양부터 확인하고 함수를 하나씩 예를 들어보겠습니다.   카이제곱 분포를 활용해서 집단 독립성 여부 가설에 대해 검정하는 예도 마지막에 들어보겠습니다(이것이 카이제곱 분포 공부하는 이유! ^^).

 

함수 구분

함수 / parameter

 chisq()

  밀도함수

  (density function) 

d

  dchisq(df)

  누적분포함수

 (cumulative distribution function)

p

  pchisq(df, lower.tail=TRUE/FALSE)

  분위수 함수

 (quantile function)

q

  qchisq(df, lower.tail=TRUE/FALSE)

  난수 발생

  (random number generation)

r

  rchisq(n, df)

 

 

 

(1) 카이제곱 분포 그래프 (chi-squared distribution plot by degrees of freedom)
    : stat_function(fun=dchisq, args=list(df))

 

위의 카이제곱 분포 정의에서 보다시피 확률변수 Xi를 제곱한 값을 summation 한 값의 분포가 카이제곱 분포이므로 모두 + 값만을 가지고 되어 0 보다 오른쪽에 있게 됩니다.  (표준정규분포는 0을 중심으로 + 와 - 가 양분)  그리고 왼쪽으로 심하게 치우쳐 있고, 오른쪽으로 꼬리가 길게 늘어서 있습니다.  오른쪽으로 갈 수록 값이 큰 경우의 도수가 현저하게 줄어드고, 0 에 가까운 왼쪽으로 갈 수록 도수가 급격하게 증가하는, 심하게 왼쪽으로 쏠린 분포라고 하겠습니다.  자유도 k 에 따라서 그래프가 바뀌는 모양은 아래의 3개 case를 통해서 확인해보기 바랍니다. (자유도 3의 경우 F 분포와 모양이 유사하지요?)

 

 

> library(ggplot2)
> 
> 
> ggplot(data.frame(x=c(0,10)), aes(x=x)) +
+   stat_function(fun=dchisq, args=list(df=1), colour="black", size=1.2) +
+   geom_text(x=0.6, y=1, label="df=1") +
+    
+   stat_function(fun=dchisq, args=list(df=2), colour="blue", size=1.2) +
+   geom_text(x=0, y=0.55, label="df=2") +
+      
+   stat_function(fun=dchisq, args=list(df=3), colour="red", size=1.2) +
+   geom_text(x=0.5, y=0.05, label="df=3") +
+      
+   ggtitle("Chisq-Distribution")

 

 

 

 

 

 

(2) 누적 카이제곱 분포 그래프 (cumulative chi-squared distribution plot)

   : stat_function(fun=pchisq, args=list(df))

 

 

> # (2) 누적 카이제곱 분포 그래프 (cumulative chi-squared distribution plot) > > ggplot(data.frame(x=c(0,10)), aes(x=x)) + + stat_function(fun=pchisq, args=list(df=1), colour="black", size=1.2) + + geom_text(x=2.5, y=0.93, label="df=1") + + + stat_function(fun=pchisq, args=list(df=2), colour="blue", size=1.2) + + geom_text(x=2.5, y=0.77, label="df=2") + + + stat_function(fun=pchisq, args=list(df=3), colour="red", size=1.2) + + geom_text(x=2.5, y=0.45, label="df=3") + + + ggtitle("Cumulative Chisq-Distribution")

 

 

 

 

 

(3) 누적 카이제곱 분포 확률 값 계산 (chisq-distribution probability)

    : pchisq(q, df, lower.tail = TRUE)

 

 

> # (3) 누적 카이제곱 분포 확률 값 계산 (cumulative probability at q quantile of chisq-distribution) > pchisq(q=2.5, df=1, lower.tail = TRUE) [1] 0.8861537 >

> pchisq(q=2.5, df=2, lower.tail = TRUE)
[1] 0.7134952
> 
> pchisq(q=2.5, df=3, lower.tail = TRUE)
[1] 0.5247089
 

 

 

 

(4) 카이제곱 분포 분위수 계산 (quantile value for specific probability of chisq-distribution)

 : qchisq(p, df, lower.tail=TRUE) 

 

위 (3)번의 pchisq()와 역의 관계에 있는 아래 (4)번의 qchisq() 함수 예제는 아래와 같습니다.

 

 
> # (4) 카이제곱 분포 분위수 계산 (quantile value for specific probability of chisq-distribution)
> 
> qchisq(p=0.8861537, df=1, lower.tail = TRUE)
[1] 2.5
> 
> qchisq(p=0.7134952, df=2, lower.tail = TRUE)
[1] 2.5
> 
> qchisq(p=0.5247089, df=3, lower.tail = TRUE)
[1] 2.5

 

 

 

 

(5) 카이제곱 분포 난수 발생

 

> # (5) 카이제곱 분포 난수 발생 (chisq-distribution random number generation) : rchisq(n, df)
> 
> rchisq <- rchisq(n=100, df=2)
> rchisq
  [1]  5.55032676  0.11525963  2.07644124  0.51528086  0.11859888  0.42973604  0.70499599  0.14487303
  [9]  5.97267125  1.56690348  0.58191202  0.70909369  3.35827777  0.21920825  0.61245462  0.38701193
 [17]  0.24927629  1.57808212  0.88996040  0.64965766  0.25870279  1.49434540  3.54133658  0.42318565
 [25] 11.70768889  2.13877942  1.06981339  5.29613195  1.51610509  2.70021112  1.01412157  0.64321783
 [33]  0.12791346  0.34348613  3.98654362  0.66337839  5.49603471  2.08228802  0.03906119  1.58919877
 [41]  0.38777794  2.31990975  3.21569380  0.19971089  1.50807461  4.54283466  2.97805146  3.27722149
 [49]  3.22608707  3.59353437  6.00263483  2.61364942  0.40056050  2.41963582  0.33523347  2.36712127
 [57]  8.51839523  0.13128607  1.41192196  5.25834986  0.83043323  5.23781460  1.70021291  3.83561464
 [65]  4.29016862  2.77080963  0.73163409  9.30564564  0.68487159  1.74665897  3.43854539  1.57618785
 [73]  2.15393289  0.56409968  0.03347495  0.21295686  2.33164034  0.12940326  2.08870172  2.88493676
 [81]  0.28369908  6.03088340  0.15760641  0.67655815  1.46900015  0.41068896  3.06536991  0.48259672
 [89]  0.13034650  1.31484537  0.38220237  1.68373899  2.75310928  0.94021333  0.50195281  0.55858578
 [97]  2.15590237 11.30222464  1.45066890  0.09005381
> 
> hist(rchisq, breaks=20)

 

 

 

 

 

 

(6) 범주형 데이터에 대한 분할표와 카이제곱 검정 (contingency table and chisq-test for categorical data)

  : CrossTable(x, y, chisq = TRUE)

 

vcd 패키지에 들어있는, 류머티스 관점염 환자에 대한 신약 테스트 결과를 모아놓은 Arthritis 데이터 프레임을 가지고 예를 들겠습니다.  gmodels 패키지를 사용해서 분할표를 집계하고, chisq=TRUE 옵션을 추가해서 카이제곱 검정 (신약이 류머티즘성 관점염에 효과가 있는지 없는지, 독립적인지 아닌지) 을 해보았습니다.  p-value = 0.001462643  인 것으로  봐서는 유의수준 5% 에서 통계적으로 효과가 있다고 봐야겠네요. (효과가 없다는 귀무가설을 기각하고, 효과가 있다(즉 독립이 아니다)는 대립가설 채택)

 

> # (6) chisq-test
> library(gmodels)
> library(vcd)
> 
> str(Arthritis)
'data.frame':	84 obs. of  5 variables:
 $ ID       : int  57 46 77 17 36 23 75 39 33 55 ...
 $ Treatment: Factor w/ 2 levels "Placebo","Treated": 2 2 2 2 2 2 2 2 2 2 ...
 $ Sex      : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 2 2 ...
 $ Age      : int  27 29 30 32 46 58 59 59 63 63 ...
 $ Improved : Ord.factor w/ 3 levels "None"<"Some"<..: 2 1 1 3 3 3 1 3 1 1 ...
> 
> attach(Arthritis)
> 
> CrossTable(Treatment, Improved, 
+            expected = TRUE, # expected frequency
+            chisq = TRUE) # chisq-test 

 
   Cell Contents
|-------------------------|
|                       N |
|              Expected N |
| Chi-square contribution |
|           N / Row Total |
|           N / Col Total |
|         N / Table Total |
|-------------------------|

 
Total Observations in Table:  84 

 
             | Improved 
   Treatment |      None |      Some |    Marked | Row Total | 
-------------|-----------|-----------|-----------|-----------|
     Placebo |        29 |         7 |         7 |        43 | 
             |    21.500 |     7.167 |    14.333 |           | 
             |     2.616 |     0.004 |     3.752 |           | 
             |     0.674 |     0.163 |     0.163 |     0.512 | 
             |     0.690 |     0.500 |     0.250 |           | 
             |     0.345 |     0.083 |     0.083 |           | 
-------------|-----------|-----------|-----------|-----------|
     Treated |        13 |         7 |        21 |        41 | 
             |    20.500 |     6.833 |    13.667 |           | 
             |     2.744 |     0.004 |     3.935 |           | 
             |     0.317 |     0.171 |     0.512 |     0.488 | 
             |     0.310 |     0.500 |     0.750 |           | 
             |     0.155 |     0.083 |     0.250 |           | 
-------------|-----------|-----------|-----------|-----------|
Column Total |        42 |        14 |        28 |        84 | 
             |     0.500 |     0.167 |     0.333 |           | 
-------------|-----------|-----------|-----------|-----------|

 
Statistics for All Table Factors


Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  13.05502     d.f. =  2     p =  0.001462643 


 
> 
> detach(Arthritis)

 

 

많은 도움 되었기를 바랍니다.

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

 

Posted by R Friend Rfriend

댓글을 달아 주세요

  1. 2020.11.09 16:30  댓글주소  수정/삭제  댓글쓰기

    비밀댓글입니다

    • R Friend Rfriend 2020.11.11 19:48 신고  댓글주소  수정/삭제

      안녕하세요.

      빈도 데이터에 대한 포아송 회귀모형 적합할 때 과대산포 시 3가지 대안에 대한 개념적인 소개글은 아래 포스팅을 참고하세요.

      https://rfriend.tistory.com/490

      로지스틱 회귀모형 시 과대산포 문제는 잘 모르겠습니다.

연속형 확률분포 (Continuous probability distribution)에는

 

 - 정규분포 (normal distribution)

   : norm()

 

 - 균등분포 (uniform distribution)

   : unif()

 

  - 지수분포 (exponential distribution)

   : exp()

 

 - t-분포 (t-distribution)

   : t()

 

 - F-분포 (F-distribution)

   : f()

 

  - 카이제곱분포(chisq-distribution)

   : chisq()

 

등이 있습니다.  

 

표본분포를 나타낼 때 t-분포, F-분포, 카이제곱 분포 등을 사용하는데요, 이번 포스팅에서는 이중에서 F-분포 (F-distribution)에 대해서 소개하도록 하겠습니다. 

 

F-분포는 F-검정 (F-test)과 두 집단 이상의 분산이 같은지 여부를 비교하는 분산분석(ANOVA, Analysis of Variance)에 사용되며, (카이제곱 분포처럼) 분산을 제곱한 값만을 사용하므로 양(+)의 값만을 가지고 되고, 왼쪽으로 치우치고 오른쪽으로 꼬리가 긴 비대칭형 형태를 띠고 있습니다.  (카이제곱 분포와 모양이 유사함)

 

 

 

 

 

 

R에서의 F 분포 함수 및 parameter는 아래와 같습니다.

 

함수 구분

함수 / parameter

 f()

  밀도함수

  (density function) 

d

  df(df1, df2)

  누적분포함수

 (cumulative distribution function)

p

  pf(df1, df2, lower.tail=TRUE/FALSE)

  분위수 함수

 (quantile function)

q

  qf(df1, df2, lower.tail=TRUE/FALSE)

  난수 발생

  (random number generation)

r

  rf(n, df1, df2)

 

 

 

(1) F 분포 그래프 (F-distribution plot) : stat_function(fun=df, args=list(df1=xx, df2=xx)

 

 

> library(ggplot2)
> ggplot(data.frame(x=c(0,5)), aes(x=x)) +
+   stat_function(fun=df, args=list(df1=5, df2=10), colour="blue", size=1) +
+   stat_function(fun=df, args=list(df1=10, df2=30), colour="red", size=2) +
+   stat_function(fun=df, args=list(df1=50, df2=100), colour="yellow", size=3) +
+   annotate("segment", x=3, xend=3.5, y=1.4, yend=1.4, colour="blue", size=1) +
+   annotate("segment", x=3, xend=3.5, y=1.2, yend=1.2, colour="red", size=2) + 
+   annotate("segment", x=3, xend=3.5, y=1.0, yend=1.0, colour="yellow", size=3) + 
+   annotate("text", x=4.3, y=1.4, label="F(df1=5, df2=10)") +
+   annotate("text", x=4.3, y=1.2, label="F(df1=10, df2=30)") + 
+   annotate("text", x=4.3, y=1.0, label="F(df1=50, df2=100)") +
+   ggtitle("F Distribution")

 

 

 

 

 

(2) 누적 F확률분포 그래프 (Cumulative F-distribution plot)
  : stat_function(fun=pf, args=list(df1=xx, df2=xx)

 

 

> # (2) 누적 F분포 그래프 (Cumulative F-distribution plot) : fun=pf > ggplot(data.frame(x=c(0,5)), aes(x=x)) + + stat_function(fun=pf, args=list(df1=5, df2=10), colour="blue", size=1) + + ggtitle("Cumulative F-distribution : F(df1=5, df2=10)")

 

 

 

 

 

 

 

(3) 누적 F분포 확률 값 계산(F-distribution cumulative probability) 
     : pf(q, df1, df2, lower.tail=TRUE)

 

 

> # (3) 누적 F분포 확률 값 계산 : pf(q, df1, df2, lower.tail=TRUE/FALSE)
> pf(q=2, df1=5, df2=10, lower.tail = TRUE)
[1] 0.835805

 

 

 

 

 

(4) F분포 분위수 함수 값 계산 (F-distribution quantile value) : qf(p, df1, df2, lower.tail=TRUE)

 

 

> # (4) F분포 분위수 함수 값 계산 : qf(p, df1, df2, lower.tail=TRUE/FALSE)
> qf(p=0.835805, df1=5, df2=10, lower.tail = TRUE)
[1] 2

 

 

F분포 분위수 함수 값은 (3)번의 F분포 누적 확률분포 값과는 역의 관계에 있습니다.

 

 

 

(5) F분포 난수 발생 (random number generation) : rf(n, df1, df2) 

 

 

 

> rf <- rf(100, df1=5, df2=10)

>
rf [1] 0.62010969 3.88443539 0.49659061 0.80159363 0.17508238 0.27307005 0.26126223 0.73465551 [9] 0.42973215 0.95448854 0.33259309 0.80825296 1.18304459 0.28760500 0.81315967 1.44000014 [17] 0.83041456 0.76020586 2.40585264 1.66243746 0.71626772 4.18541745 1.81731056 0.32127326 [25] 0.59385259 1.48154031 1.28662657 0.51072359 2.95967809 1.48123950 0.85420041 0.29173710 [33] 1.55212508 0.82884045 0.38375658 1.75641755 2.39041073 0.86752014 0.87916583 0.54687123 [41] 0.95937758 1.08620415 2.56198072 3.33297651 0.35509530 1.33342398 0.81011460 0.46871978 [49] 0.67320976 1.82654953 1.57868396 0.51026464 1.81435701 1.59062772 1.79081953 0.26546396 [57] 0.70741755 0.56732326 2.03751027 0.81498081 1.27629405 0.65799392 1.82927296 0.45875179 [65] 3.62942043 0.57416234 0.21089950 1.30947117 1.08285304 0.90621979 1.14831429 0.31727995 [73] 0.11544574 1.53316697 2.16508104 1.15088174 2.26276896 1.50965228 6.53052329 1.96294253 [81] 0.78687566 0.41327969 0.52918777 0.96117638 0.44835433 2.34021493 3.89895973 0.84022777 [89] 1.49158411 1.75840630 0.05668181 0.82050282 1.14190710 0.46880962 0.70338470 1.09613826 [97] 3.95759814 2.81014843 2.08048954 0.30376581 >
> hist(rf, breaks=20)
 

 

 

 

난수 발생은 매번 할 때마다 바뀌므로 위의 결과와는 다르게 보일 것입니다만, 왼쪽으로 치우쳐있고 오른쪽으로 공룡꼬리처럼 가늘고 길게 늘어선 비대칭 모양이 나타날 것입니다.  

 

많은 도움이 되었기를 바랍니다.

아래의 '공감' 꾸~욱 눌러주시면 글 쓰는데 큰 힘이 됩니다. ^^  (로그인 필요 없어요)

 

 

 

Posted by R Friend Rfriend

댓글을 달아 주세요

  1. 통계병아리 2016.10.24 19:26  댓글주소  수정/삭제  댓글쓰기

    그렇다면 다중비교에서 f분포의 p-value를 구할때는 df를 사용해서, 1-df(모집단의 수, 자유도 1, 자유도 2) 의 식으로 계산하면 되나요?

    • R Friend Rfriend 2016.10.25 17:14 신고  댓글주소  수정/삭제

      위 포스팅의 3)번 누적 F분포 확률 값 계산(F-distribution cumulative probability) 내용이 p-value 를 구하는 방법입니다.

      아래 R script 처럼 lower.tail = FALSE로 지정해주시면 p-value를 구할 수 있습니다.

      pf(q=2, df1=5, df2=10, lower.tail = FALSE)

    • 통계병아리 2016.10.25 21:44  댓글주소  수정/삭제

      아, 굳이 1-를 할 필요없이 pf에서 lower.tail을 FALSE로 지정해주면 되는군요! 감사합니다!!

  2. Kang Han 2020.04.08 00:21 신고  댓글주소  수정/삭제  댓글쓰기

    감사합니다. 공부하는데 많은 도움이 되네요. 올려 놓으신 좋은 자료들 잘 보고 있습니다.

 

연속형 확률분포 (Continuous probability distribution)에는

 

 - 정규분포 (normal distribution)

   : norm()

 

 - 균등분포 (uniform distribution)

   : unif()

 

  - 지수분포 (exponential distribution)

   : exp()

 

 - t-분포 (t-distribution)

   : t()

 

 - F-분포 (F-distribution)

   : f()

 

  - 카이제곱분포(chisq-distribution)

   : chisq()

 

등이 있습니다.  

 

표본분포를 나타낼 때 t-분포, F-분포, 카이제곱 분포 등을 사용하는데요, 이번 포스팅에서는 이중에서 t-분포 (Student's t-distribution)에 대해서 소개하도록 하겠습니다. 

 

스튜던트 t-분포 (Student's t-distribution)을 줄여서 t-분포 (t-distribution)라고 보통 부르곤 합니다.  스튜던트 t-분포 (Student's t-distribution)의 유래를 살펴보면 아래와 같습니다.

 

 

[ 스튜던트 t-분포 (Student's t-distribution)의 유래 ]

 

프리드리히 로베르트 헬메르트(독일어: Friedrich Robert Helmert)가 1875년에 도입하였다. 이듬해 야코프 뤼로트(독일어: Jacob Lüroth)도 같은 분포를 재발견하였다. 그러나 헬메르트와 뤼로트의 논문은 영문 학계에 널리 알려지지 않았다.
1908년에 윌리엄 고셋이 "스튜던트"(영어: Student, ‘학생’)라는 필명으로 1908년 재발견하였다. 고셋은 기네스 양조 공장에서 일했고, 맥주에 사용되는 보리의 질을 시험하기 위해 이 분포를 도입하였고, 경쟁사들에게 기네스의 획기적인 통계 기법을 숨기기 위해 필명을 사용하였다고 한다. 이후 저명한 통계학자인 로널드 피셔는 이 분포를 "스튜던트 분포"라고 불렀고, t라는 기호를 사용하였다. 피셔 이후 이 분포는 고셋의 필명을 따 "스튜던트 t 분포"로 알려지게 되었다.

 

   - source : wikipedia

 

 

 

 

정규분포에서는 모분산(σ2)을 알고 있다고 가정하는데요, 실전의 현실세계에서는 모분산(σ2)을 모르는 경우가 대부분이다보니 표본을 추출해서 표본분산(s2)을 계산하여 사용하는 경우가 다반사입니다.  이러다 보니 표준정규분포 통계량 Z 값을 사용할 수 없고 표본확률분포를 사용해야 하는데 그중의 하나가 T통계량을 사용하는 t-분포가 되겠습니다.

 

 

 

 

 

 

t-분포는 평균이 0 이고, 평균을 중심으로 좌우 대칭형태로 되어있으며, 정규분포보다 가운데의 높이가 조금 낮고 좌우의 옆 부분은 정규분포보다 조금 더 높은 형태를 취하고 있습니다.  t-분포는 자유가를 모수로 가지고 있으며, 자유도가 높을 수록, 즉 표본의 갯수가 증가할 수록 중심극한의 정리(central limit theorem)에 의해 정규분포에 근사하게 됩니다.  (보통은 표본의 갯수 30개를 기준으로 이보다 많으면 정규분포, 적으면 t-분포를 사용)

 

표본의 수가 많으면 모집단을 대표할 신뢰도가 높아지게 되어 좋기는 합니다만, 많은 경우는 비용과 시간의 한계로 인해서 한정된 수의 표본만을 추출해서 분석해야 하는 경우가 생기게 됩니다.  표본의 수가 적은 경우에 표본이 모집단의 어느 한쪽으로 쏠려서 추출될 경우 모집단을 잘 대표할 수 없는 신뢰도 이슈를 보완하기 위해서, 정규분포일 때보다 평균에서 양쪽으로 멀어지는 바깥쪽 부분의 확률의 수준을 더 높인 분포가 t-분포(t-distribution)이 되겠습니다.

 

아래에 정규분포와 t-분포를 겹쳐서 그려보았는데요, 위의 설명을 이해하는데 도움이 될 것입니다.

 

 

> install.packages("ggplot2") # install

> library(ggplot2)
> ggplot(data.frame(x=c(-3,3)), aes(x=x)) +
+   stat_function(fun=dnorm, colour="blue", size=1) +
+   stat_function(fun=dt, args=list(df=3), colour="red", size=2) +
+   stat_function(fun=dt, args=list(df=1), colour="yellow", size=3) +
+   annotate("segment", x=1.5, xend=2, y=0.4, yend=0.4, colour="blue", size=1) +
+   annotate("segment", x=1.5, xend=2, y=0.37, yend=0.37, colour="red", size=2) + 
+   annotate("segment", x=1.5, xend=2, y=0.34, yend=0.34, colour="yellow", size=3) + 
+   annotate("text", x=2.4, y=0.4, label="N(0,1)") +
+   annotate("text", x=2.4, y=0.37, label="t(3)") + 
+   annotate("text", x=2.4, y=0.34, label="t(1)") + 
+   ggtitle("Normal Distribution, t-distribution")
 

 

 

 

 

R에서 t-분포 (t-distribution)을 위해 사용하는 함수 및 Parameter는 아래와 같습니다.  정규분포(normal distribution)는 평균(mean)과 표준편차(sd)를 parameter로 사용하고, 균등분포(uniform distribution)는 구간의 최소값(min)과 최대값(max)를 parameter로 사용하며, 지수분포(exponential distribution)는 Lamda(λ, rate)를 parameter로 사용하는데 반해, t-분포(t-distribution)은 자유도(df, degrees of freedom) 를 parameter로 사용합니다.

 

참고로 자유도 (df, degrees of freedom)은 통계량을 구성하는 확률변수들 중에서 자유롭게 선택가능한 확률변수의 개수를 의미합니다. 

 

 

 함수 구분

함수/ Parameter 

t() 

  밀도함수

  (density function) 

  dt(x, df)

  누적분포함수

 (cumulative distribution function) 

p

  pt(q, df, lower.tail = TRUE/FALSE)

  분위수 함수

 (quantile function)

q

  qt(p, df, lower.tail = TRUE/FALSE)

  난수 발생

  (random number generation)

r

  rt(n, df)

 

 

R 함수를 차례대로 하나씩 예를 들어보겠습니다.

 

위에서 t분포 밀도함수 : dt(x, df) 의 그래프를 정규분포와 비교해서 그려보았으며, 아래에는 t-분포의 누적분포함수를 그래프로 그려보고 누적분포확률도 계산해보겠습니다.

 

(1) t분포 누적분포함수 그래프 : stat_function(fun = pt, args=list(df=1))

 

 

> # 누적 t분포 그래프 (Cumulative t-distribution plot) : fun=pt > ggplot(data.frame(x=c(-3,3)), aes(x=x)) + + stat_function(fun=pt, args=list(df=1), colour="brown", size=1.5) + + ggtitle("Cumulative t-Distribution : t(1)")

 

 

 

 

 

(2) t분포 누적분포함수 : pt(q, df, lower.tail = TRUE/FALSE)

 

위의 t(df=1) 그래프에서 x(-inf, 1) 범위의 누적분포확률은 0.75 임을 알 수 있습니다.

 

 

> # 누적 t분포 확률 값 계산 : pt(q, df, lower.tail=TRUE/FALSE) > pt(q=1, df=1, lower.tail = TRUE) [1] 0.75

 

 

 

 

(3) t분포 분위수 함수 : qt(p, df, lower.tail = TRUE/FALSE)

 

누적확률분포와 분위수 함수는 역의 관계에 있는데요, 위의 pt(q=1, df=1, lower.tail=T) = 0.75 였으므로, 누적확률분포값이 0.75가 되는 분위수 q는 '1'이 되는지 아래 R script로 확인해보겠습니다.

 

 

> # t분포 분위수 함수 값 계산 : qt(p, df, lower.tail=TRUE/FALSE) > qt(p=0.75, df=1, lower.tail = TRUE) [1] 1

 

 

 

 

(4) t분포 난수 발생 : rt(n, df)

 

마지막으로 t분포에서 rt(n, df) 함수를 이용해서 난수를 발생시켜보겠습니다.  난수는 매번 생성할 때마다 바뀌므로 매 시행마다 아래의 예제와는 조금씩 달라지게 됩니다만, 평균 0을 중심으로 좌우 대칭형태 (정규분포와 유사)를 띨 것입니다.

 

 

> rt <- rt(50, df=1)
> rt
 [1] -1.89734826 -1.40787172 -2.38022561 -0.21218535  0.71259957 -0.37759519 -0.46671360 -0.20455553
 [9] -0.54603347  2.88325761  0.13593286  0.09719242  1.35843796 -1.86861627  5.69879846  1.23054665
[17] -1.24953468 -0.76327864  1.87092396 -0.39719483 -0.42141574  0.10862682  0.54106231  0.30827837
[25] -1.25508149  0.54352324 -1.92825750  0.22491497 -0.63797793 -0.37089263  7.31876302  2.25023970
[33] -1.60455685 -1.64779189 -5.54583982 -8.82959795  0.53445244 -0.47451960 -2.52582931  1.57372391
[41]  2.30557669 -0.04118914 -0.71146732 -0.27621122 -7.29220086 -0.52472800 -0.78465027 -2.07649916
[49]  0.38322764 -1.71782797
>
> hist(rt, breaks=20)
 

 

 

 

많은 도움이 되었기를 바랍니다. 

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

 

Posted by R Friend Rfriend

댓글을 달아 주세요

R 프로그래밍 중에 "target of assignment expands to non-language object" error 메시지가 발생하면 assign() 함수를 사용하여 문제를 해결할 수 있습니다.

 

"target of assignment expands to non-language object" error 가 발생하는 이유는 할당하려는 목표 객체(Target object) 를 R이 인식하지 못해 NULL 값으로 처리되어 있는 상태이기 때문입니다. 

 

아래에 이런 에러가 발생하는 간단한 예와 해결 방법을 소개해보겠습니다.

 

 

x변수는 1부터 30까지의 정수, y변수는 1부터 5까지의 정수이고, 이 두변수를 묶어서 xy라는 데이터 프레임을 만들어보겠습니다.

 

 

> x <- c(1:30) >

> y <- c(1:5)
>
> xy <- data.frame(x, y)
>
> str(xy)
'data.frame':	30 obs. of  2 variables:
 $ x: int  1 2 3 4 5 6 7 8 9 10 ...
 $ y: int  1 2 3 4 5 1 2 3 4 5 ...
>
> xy
    x y
1   1 1
2   2 2
3   3 3
4   4 4
5   5 5
6   6 1
7   7 2
8   8 3
9   9 4
10 10 5
11 11 1
12 12 2
13 13 3
14 14 4
15 15 5
16 16 1
17 17 2
18 18 3
19 19 4
20 20 5
21 21 1
22 22 2
23 23 3
24 24 4
25 25 5
26 26 1
27 27 2
28 28 3
29 29 4
30 30 5

 

 

위에서 생성한 xy 데이터프레임을 가지고, y변수 1, 2, 3, 4, 5 를 포함한 행(row) 별로 각 각 개별 데이터 프레임을 Loop 를 사용해서 만들어보겠습니다.  아래 Loop program에서 목표 객체(Target object)에 paste() 함수를 사용해서 공통 접두사(predix)로 'xy_'를 사용하고 뒤에 'i'로 loop 를 돌면서 1, 2, 3, 4, 5 를 붙여주려고 프로그램을 짰습니다만, "target of assignment expands to non-language object" 라는 오류 메시지가 떴습니다.  아래 프로그램에 대해 R은 [ paste("xy_", i, sep="") ]부분을 NULL 값으로 인식하기 때문에 이런 오류가 발생하게 됩니다.

 

 

> for (i in 1:5) {
+   paste("xy_", i, sep="") <- subset(xy, subset = (y == i))
+ }
Error in paste("xy_", i, sep = "") <- subset(xy, subset = (y == i)) : 
  target of assignment expands to non-language object

 

 

 

 

위 프로그램에서 하고자 했던 의도대로 R이 이해하고 실행을 하게끔 하려면 assign() 함수를 사용해야만 합니다.  assign() 함수를 사용할 때는 아래 처럼 '<-' 대신에 ',' 가 사용되었음에 주의하시기 바랍니다.

 

 
> for (i in 1:5) {
+   assign(paste("xy_", i, sep=""), subset(xy, subset = (y == i)))
+ }
>
> xy_1
    x y
1   1 1
6   6 1
11 11 1
16 16 1
21 21 1
26 26 1
>
> xy_2
    x y
2   2 2
7   7 2
12 12 2
17 17 2
22 22 2
27 27 2
>
> xy_3
    x y
3   3 3
8   8 3
13 13 3
18 18 3
23 23 3
28 28 3
>
> xy_4
    x y
4   4 4
9   9 4
14 14 4
19 19 4
24 24 4
29 29 4
>
> xy_5
    x y
5   5 5
10 10 5
15 15 5
20 20 5
25 25 5
30 30 5

 

 

많은 도움이 되었기를 바랍니다.

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡' 단추를 꾸욱 눌러주세요.^^

 

 

Posted by R Friend Rfriend

댓글을 달아 주세요

  1. 산낙지 2016.10.07 15:42  댓글주소  수정/삭제  댓글쓰기

    안녕하세요! 오늘도 포스팅 잘 보았습니다. 제가 본문을 참고해서 명령어를 실행해보았습니다. 그런데 문제가 발생해서 하나를 여쭙고 싶습니다.
    우선 저는 다음과 같은 명령어를 실행하고 싶은데요.

    pdt$birth1 <- c(1, 0, 0, 0, 0)[match(pdt$r_birth1, c(1, 2, 3, 4, 5))]
    pdt$birth2 <- c(1, 0, 0, 0, 0)[match(pdt$r_birth2, c(1, 2, 3, 4, 5))]
    pdt$birth3 <- c(1, 0, 0, 0, 0)[match(pdt$r_birth3, c(1, 2, 3, 4, 5))]

    그래서 assign을 이용해서 다음과 같은 명령어를 만들어 보았습니다.

    for (i in 1:3) {
    assign(paste("pdt$birth", i, sep=""), c(1, 0, 0, 0, 0)[match(paste("pdt$r_birth", i, sep=""), c(1, 2, 3, 4, 5))])
    }

    그런데 명령어가 오류 없이 실행은 되는데, pdt라는 객체에 birth1~3이 아예 안 생기더라고요.
    혹시 제 명령어에는 어떤 문제가 있는 건지 조언을 구할 수 있을까요? ㅠㅠ

    • R Friend Rfriend 2016.10.08 15:29 신고  댓글주소  수정/삭제

      산낙지님, 반갑습니다.

      assign() 함수 연습하는 용도라면요, birth1, birth2, birth3으로 개별 객체를 loop 연산으로 assign()으로 생성 후에 data.frame(cbind()) 하면 되겠네요.

      근데, 문의주신 형태의 output을 얻는데는 assign() 함수를 안쓰고 그냥 transform() 함수에 loop 돌리면 될거 같습니다.

연속형 확률분포 (Continuous probability distribution)에는

 

 - 정규분포 (normal distribution)

   : norm()

 

 - 균등분포 (uniform distribution)

   : unif()

 

 - 지수분포 (exponential distribution)

   : exp()

 

 - t-분포 (t-distribution)

   : t()

 

 - F-분포 (F-distribution)

   : f()

 

 - 카이제곱분포(chisq-distribution)

   : chisq()

 

등이 있습니다.  

 

 

이번 포스팅에서는 지수분포 (exponential distribution)에 대해서 소개하도록 하겠습니다.  지수분포 (exponential distribution)는 어떤 특정 사건이 발생하기 전까지 걸리는 시간을 나타내기 위해 많이 사용하는 확률분포 입니다. 

 

지수분포(exponential distribution)의 예로는 전자레인지의 수명시간, 콜센터에 전화가 걸려 올 때까지 걸리는 시간, 경부고속도로 안성나들목에서 다음번 교통사고가 발생할 때까지 걸리는 시간, 은행 지점에 고객이 내방하는데 걸리는 시간 등이 있겠습니다.

 

참고로, 이산형 확률분포 중에서 포아송 분포 (Poisson distribution)이 단위 시간 혹은 단위 공간에서 특정 사건이 발생하는 횟수에 대한 분포를 나타낼 때 주로 사용한다고 했는데요, 헷갈리지 않도록 주의를 해야겠습니다.

 

 

 

확률변수 X의 확률밀도함수가 위와 같을 때, 확률변수 X는 모수 λ인 지수분포를 따른다고 말합니다.

 

R에서 지수분포(exponential distribution)와 관련된 함수 및  파라미터는 다음과 같습니다.

 

 함수 구분

지수분포 함수/파라미터 

 exp()

 밀도 함수
 (Density function)

 d

  dexp(x, rate)

 누적분포 함수

 (Cumulative distribution function)

 p

  pexp(q, rate, lower.tail=TRUE/FALSE)

 분위수 함수

 (Quantile function)

 q

  qexp(p, rate, lower.tail=TRUE/FALSE)

 난수 생성

 (Random number generation)

 r

  rexp(n, rate)

 

 

(1) 지수분포 그래프 (Exponential distribution plot) : fun=dexp

 

 
> library(ggplot2)

> # (1) 지수분포 그래프 (Exponential distribution plot) : fun=dexp > ggplot(data.frame(x=c(0,10)), aes(x=x)) + + stat_function(fun=dexp, args=list(rate=1), colour="brown", size=1.5) + + ggtitle("Exponential Distribution") > + ggtitle("Exponential Distribution")

 

 

 

 

 

(2) 누적 지수분포 그래프 (Cumulative exponential distribution plot) : fun = pexp

 

 

> # (2) (Cumulative exponential distribution plot) : fun=pexp > ggplot(data.frame(x=c(0,10)), aes(x=x)) + + stat_function(fun=pexp, args=list(rate=1), colour="brown", size=1.5) + + ggtitle("Cumulative Exponential Distribution")

 

 

 

 

(3) 누적 지수분포 확률 값 계산 : pexp(q, rate, lower.tail=TRUE/FALSE)

 

 

 
> pexp(q=2, rate=1, lower.tail = TRUE)
[1] 0.8646647
 

 

λ (rate) = 1 인 지수분포에서 0~1까지의 누적 확률 값은 0.8646647 임을 알 수 있습니다.  

 

 

(4) 지수분포 분위수 함수 값 계산 : qexp(p, rate, lower.tail=TRUE/FALSE)

 

 

> qexp(p=0.8646647, rate=1, lower.tail = TRUE)
[1] 2
 

 

 

qexp()는 pexp()와 역의 관계에 있다고 보면 되겠습니다.

 

 

(5) 지수분포 난수 발생 (exponential distribution random number generation) : rexp(n, rate)

 

 

> rexp(100, rate=1)
  [1] 0.805385854 1.077017598 0.941678341 2.059229603 0.517943248 0.955476408 0.575837716
  [8] 0.851462637 0.086982322 1.243358626 1.077268675 2.604957888 0.007571515 1.793674221
 [15] 0.118729103 0.096055712 0.015758928 0.201158101 0.914114063 0.130984491 2.752139235
 [22] 0.829986667 0.651976457 1.265562156 2.635988993 1.190808342 0.444055191 1.480476206
 [29] 1.741018226 2.692880185 0.804053361 3.127147071 0.902618388 1.432761851 0.369694262
 [36] 0.290926187 0.576759913 0.827636680 1.634353038 2.113214617 0.570110160 0.609782309
 [43] 1.985241502 1.067016441 0.098556668 3.326005637 2.261946740 6.395236475 1.906314444
 [50] 0.503994692 1.578938061 0.144050682 0.361734510 1.495605791 1.167056286 1.397221429
 [57] 1.598533234 0.370363955 0.153343928 0.351399011 0.957647500 1.053767695 0.272256882
 [64] 1.176451771 0.222995248 1.125289913 0.076051627 3.489328747 0.199440748 3.143880822
 [71] 1.640546855 4.492400575 1.102261695 0.189814932 0.222941682 2.305212835 0.069710370
 [78] 1.972949872 0.201703040 1.783014953 1.297271054 2.173743544 2.350792197 2.018307233
 [85] 0.417343667 2.533255698 0.522270208 2.068958899 1.062778262 0.210765630 0.804149691
 [92] 1.261281259 0.006859250 0.620238345 3.478939042 2.692230696 0.557543887 1.830330845
 [99] 0.478452368 0.904496278
> hist(rexp(100, rate=1), breaks=10)
 

 

 

 

 

(6) dexp(x, rate, log=TRUE)

 

log=TRUE 옵션을 설정하면 지수분포의 확률밀도함수값에 밑이 e (약 2.17)인 자연로그 ln 을 취한 값을 계산합니다. (참고 : 자연로그인 ln 의 역함수는 지수함수인 exp() )

 

아래에는 X가 모수 λ가  1인 지수분포를 따른다고 했을 때, X: 1, 2, .., 10 의 확률밀도함수값을 계산한 것입니다.  (지수분포는 연속형 확률분포이지만, log 취한거와 비교를 이해하기 쉽도록 1, 2, ..., 10 의 값을 가지고 계산했음)  'log=TRUE' 라는 옵션을 취한 값과, log(dexp) 라는 수식을 직접 입력해서 계산한 값이 서로 정확히 일치함을 알 수 있습니다.

 

> # dexp(x, rate, log=TRUE)
> dexp <- dexp(c(1:10), rate=1)
> dexp_log <- dexp(c(1:10), rate=1, log=TRUE)
> 
> exp_df <- data.frame(cbind(c(1:10), dexp, dexp_log))
> exp_df
   V1         dexp dexp_log
1   1 3.678794e-01       -1
2   2 1.353353e-01       -2
3   3 4.978707e-02       -3
4   4 1.831564e-02       -4
5   5 6.737947e-03       -5
6   6 2.478752e-03       -6
7   7 9.118820e-04       -7
8   8 3.354626e-04       -8
9   9 1.234098e-04       -9
10 10 4.539993e-05      -10
> 
> exp_df <- transform(exp_df, dexp_logarithm = log(dexp))
> exp_df
   V1         dexp dexp_log dexp_logarithm
1   1 3.678794e-01       -1             -1
2   2 1.353353e-01       -2             -2
3   3 4.978707e-02       -3             -3
4   4 1.831564e-02       -4             -4
5   5 6.737947e-03       -5             -5
6   6 2.478752e-03       -6             -6
7   7 9.118820e-04       -7             -7
8   8 3.354626e-04       -8             -8
9   9 1.234098e-04       -9             -9
10 10 4.539993e-05      -10            -10

 

 

 

 

디폴트인 dexp(x, log=FALSE) 와 dexp(x, log=TRUE) 옵션 설정해서 나온 값을 가지고 그래프로 그려서 비교해보면 아래와 같습니다.  log=TRUE 설정을 해서 자연로그를 취했더니 원래 밑으로 축 쳐졌던 지수분포 그래프가 곧은 직선으로 변환되었음을 알 수 있습니다.

 

> my_par = par(no.readonly = TRUE)
> par(oma = c(0, 0, 1, 0))
> par(mfrow = c(1, 2))
> 
> plot(dexp, main = "dexp : log = F")
> plot(dexp_log, main = "dexp : log = T")
> 
> mtext("density function of exponential distributin : log = FALSE vs. log = TRUE", 
outer = TRUE, cex = 1.2)

 

 

 

 많은 도움 되었기를 바랍니다.

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요.

 

Posted by R Friend Rfriend

댓글을 달아 주세요

  1. 2016.03.28 15:17  댓글주소  수정/삭제  댓글쓰기

    비밀댓글입니다

    • R Friend Rfriend 2016.03.28 23:41 신고  댓글주소  수정/삭제

      (1) Histogram 을 그릴 때 freq = FALSE 옵션을 설정해서 y축이 frequency 가 아니라 density가 되게 해야만 합니다.

      (2) lines(density(x)) 함수를 써서 커널 밀도 곡선을 (1)에서 만든 밀도 기반의 히스토그램 위에 덧쓰우면 됩니다.

      아래의 두개의 프로그램을 순차적으로 실행하시면 됩니다.

      # (1) exponential distribution, Histogram + Kernel Density Curve
      set.seed(12345)
      x <- rexp(1000, rate=5)
      hist(x,
      freq = FALSE, # density (not frequency)
      breaks = 50,
      col = "blue",
      main = "Histogram + Kernel Density Curve" )

      # (2) adding kernel density curve
      lines(density(x),
      col = "red",
      lwd = 3, # line width =2
      lty = 3) # dotted line

  2. 2016.05.29 17:38  댓글주소  수정/삭제  댓글쓰기

    비밀댓글입니다

    • R Friend Rfriend 2016.05.31 00:35 신고  댓글주소  수정/삭제

      dexp(x, rate, log=TRUE) 에서 log=TRUE 는 지수분포를 따르는 확률밀도함수값에 자연로그 ln을 취하는 옵션입니다.

      블로그의 마지막 부분 (6)번에 log=TRUE 옵션 설정에 대해서 추가로 내용을 보충해서 포스팅 업데이트 하였습니다. 참고하시기 바랍니다.

  3. 2020.05.06 21:21  댓글주소  수정/삭제  댓글쓰기

    비밀댓글입니다

    • R Friend Rfriend 2020.05.06 23:11 신고  댓글주소  수정/삭제

      안녕하세요.
      이중지수분포 DE(x, mu, lambda) 사용자 정의 함수는 아래 코드를 참고하세요.

      # create DE UDF
      DE <- function (x, mu = 0, lambda = 1) {
      exp(-abs(x - mu)/lambda)/(2 * lambda)
      }

      # generate x
      x <- seq(-4,4, length.out = 1000)
      print(DE(x, mu=0, lambda = 1))

      # plot DE
      curve(DE, -4, 4)