통계적 검정 (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

 

 

 

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

 

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

 


728x90
반응형
Posted by Rfriend
,

이전 포스팅들에서 이산형 확률분포, 연속형 확률분포에 대해서 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%신뢰수준이 좋기는 합니다만, 신뢰수준을 높이려면 표본의 개수를 높여야 하는데요, 표본 = 돈, 시간 이므로 보통은 신뢰수준을 특정 값으로 고정 시켜놓고 샘플 수를 적정한 타협을 보게 됩니다.

 

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

 

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

 

728x90
반응형
Posted by 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)

 

 

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

 

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

 

 

728x90
반응형
Posted by 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-분포, 카이제곱 분포 등을 사용하는데요, 이번 포스팅에서는 이중에서 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)
 

 

 

 

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

 

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

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

 

 

 

728x90
반응형
Posted by 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-분포, 카이제곱 분포 등을 사용하는데요, 이번 포스팅에서는 이중에서 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)
 

 

 

 

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

 

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

 

 

728x90
반응형
Posted by 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

 

 

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

 

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

 

 

728x90
반응형
Posted by Rfriend
,

연속형 확률분포 (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)

 

 

 

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

 

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

 

728x90
반응형
Posted by Rfriend
,

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

 

 - 정규분포 (normal distribution)

   : norm()

 

 - 균등분포 (uniform distribution)

   : unif()

 

  - 지수분포 (exponential distribution)

   : exp()

 

 - t-분포 (t-distribution)

   : t()

 

 - F-분포 (F-distribution)

   : f()

 

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

   : chisq()

 

등이 있습니다.  

 

 

이번 포스팅에서는 균등분포(uniform distribution)에 대해서 알아보겠습니다.  균등분포(uniform distribution)은 연속형 확률 분포 중에서 가장 간단한 형태로서, 구간 [mi=a, max=b]에서 값이 균등하게 퍼져 있는 집단, 일어날 확률이 균등한 분포를 말합니다. 

 

가령, 김포공항에서 제주도 공항까지 비행기로 이륙에서 착륙까지 걸리는 총 비행시간이 1시간~1시간5분 사이라고 하면, 0시~59분59초까지는 비행기가 도착할 확률이 0, 1시간~1시간5분 사이에 도착할 확률은 1, 1시간 5분 이후는 다시 확률이 0이 되는 균등분포를 따른다고 할 수 있겠습니다.

 

 

 

 

 

 

 

R에서 사용하는 균등분포 함수(uniform distribution function) 및 파라미터(parameter)들은 아래와 같으며, 필요한 함수, 파라미터를 가져다 사용하면 되겠습니다.

 

 

 함수 구분

 균등분포 함수/파라미더

 unif()

  밀도함수

  (density function)

 d

  dunif(x, min, max)

  누적분포함수

 (cumulative distribution function)

 p

  punif(q, min, max, lower.tail=TRUE/FALSE)

  분위수 함수

 (quantile function)

 q

  qunif(p, min, max, lower.tail=TRUE/FALSE)

  난수 발생

 (random number generation)

 r

  runif(n, min, max)

 

 

 

 

(1) 균등분포 그래프(uniform distribution plot) : fun = dunif

 

ggplot2의 fun= dunif() 함수를 사용해서 균등분포를 그래프로 그려보면 아래와 같이 특정 구간 [a, b]에서 확률이 균등함을 알 수 있습니다.

 

 

> library(ggplot2)

> # uniform distribution plot (min=0, max=10) > # 균등분포 : fun = dunif > ggplot(data.frame(x=c(-2,20)), aes(x=x)) + + stat_function(fun=dunif, args=list(min = 0, max = 10), colour="black", size=1) + + ggtitle("Uniform Distribution of (min=1, max=10)")

 

 

 

 

 

 

 

(2) 누적 균등분포 그래프(cumulative uniform distribution plot) : fun = punif

누적 균등분포 그래프를 그려보면 아래와 같습니다.

 

 

> # (2) 누적균등분포 함수 그래프 (Cumulative Uniform distribution plot) : fun = punif > ggplot(data.frame(x=c(-2,20)), aes(x=x)) + + stat_function(fun=punif, args=list(min = 0, max = 10), colour="black", size=1) + + ggtitle("Cumulative Uniform Distribution of (min=0, max=10)")

 

 

 

 

 

 

(3) 누적 균등분포 함수의 확률 값 계산 : punif()

 

 

> # (3) 누적 균등분포함수(cumulative uniform distribution function) 확률 값 계산 : punif() > # : punif(q, min, max, lower.tail = TRUE/FALSE) > punif(3, min=0, max=10, lower.tail=TRUE) [1] 0.3 > >

> # Uniform Distribution of (min=1, max=10), x from 0 to 3"

> ggplot(data.frame(x=c(-2,20)), aes(x=x)) + + stat_function(fun=dunif, args=list(min = 0, max = 10), colour="black", size=1) + + annotate("rect", xmin=0, xmax=3, ymin=0, ymax=0.1, alpha=0.2, fill="yellow") + + ggtitle("Uniform Distribution of (min=1, max=10), x from 0 to 3")

 

 

 

 

 

(4) 균등분포 분위수 함수 값 계산 : qunif(p, min, max, lower.tail=TRUE/FALSE)

 

이전 포스팅의 정규분포와는 함수는 qunif()로 동일하지만 괄호 안의 parameter 들은 다릅니다.

(참고: 정규분포에서는 qunif(p, mean, sd, lower.tail=T/F)

 

 
> # (4) 분위수 함수 : qunif(p, min, max, lower.tail=TRUE/FALSE)
> qunif(0.3, min=0, max=10, lower.tail = TRUE)
[1] 3

 

 

 

 

(5) 난수 발생 : runif(n, min, max)

 

난수는 매번 실행할 때마다 바뀌므로 제가 아래에 제시한 것과는 다른 숫자, 다른 그래프가 그려질 것입니다만, 형태는 균등분포를 띠는 유사한 모양이 될 것입니다.

 

 

> ru_100 <- runif(n=100, min=0, max = 10)
> ru_100

 

 [1] 7.33957568 2.78596723 6.30797744 5.01438337 6.57949706 5.90883342 3.51446293 9.28736811
  [9] 9.55213668 5.59377524 4.71003185 3.29525512 0.25759555 9.40326151 6.56466466 2.44973803
 [17] 4.88714900 3.10710648 3.84375758 8.55017741 3.09487276 0.13411621 0.44285713 8.90632265
 [25] 0.07968823 5.03465390 4.64601169 1.23565062 4.81310463 1.59225023 7.03799510 0.68870704
 [33] 4.03014086 9.97756283 5.55815726 2.01819345 7.00497545 8.50399118 2.29608430 2.92359120
 [41] 0.85656712 6.52544881 6.37193951 6.15247601 5.29502105 7.68988134 6.37691223 0.37387705
 [49] 6.89023959 1.65049129 3.75195268 7.97220092 6.50160025 9.52491436 1.70569894 9.80475205
 [57] 0.24770673 8.47412000 4.66718922 2.52269224 2.81985175 8.79845402 6.03852213 8.10848875
 [65] 1.10510449 9.35548906 1.83535387 0.47889795 6.54578585 1.61742080 4.51840400 3.99912651
 [73] 4.82545376 4.04589108 0.71750065 7.56085867 1.22887762 2.97822070 5.14541682 3.59126885
 [81] 5.00911758 1.02152702 7.78324707 4.69437196 1.13090493 3.70933500 0.03173870 5.74159309
 [89] 2.68879279 3.36398725 9.34593590 6.18818473 9.43490689 5.82578697 4.49576854 2.90029081
 [97] 3.34726356 7.19013351 9.97276521 9.39421932

 

 

> # density plot of runif(n=100, min=0, max = 10) & adding line of 0.1 uniform probability

> hist(ru_100, freq=FALSE, breaks=10, col="yellow")

> abline(h=0.1, lty=3, lwd=3, col="red")

 

 

 

 

 

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

 

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

 

 

728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 지혜 (경험 + 지식)와 자신감의 관계 (Relationship between Wisdom (Knowledge + Experience) and Confidence) 에 대해 탁월한 통찰력이 녹아있는 Dunning-Kruger Effect 에 대해서 소개해 보겠습니다.

 

아래에 소개하는 주요 내용은 WIKIPEDIA 의 내용을 많이 참고/번역하였으며, 개인적인 경험과 소견을 좀더 붙여 보았습니다.

 

코넬 대학교의 교수였던 David Dunning과 Juskin Kruger 는 1999년에 Journal of Personality and Social Psychology에 『Kruger, Justin; Dunning, David (1999). "Unskilled and Unaware of It: How Difficulties in Recognizing One's Own Incompetence Lead to Inflated Self-Assessments". 』라는 논문을 발표합니다. 

 

핵심 내용만 짧게 소개하자면, "상대적으로 지식과 경험이 부족한 사람은 자신의 역량을 실제보다 훨씬 높다고 평가하는 인지 오류 경향이 있으며, 반대로 지식과 경험이 많은 사람일 수록 자신에게 쉬운 일이 타인에게도 쉬울 것이라고 잘못 평가하는 인지오류에 빠져서 자신의 역량을 과소평가하는 경향이 있다"는 것입니다.

 

결국 자신의 역량을 제대로, 있는 그대로 평가할 수 있는 사람은 없다는 얘기가 되겠습니다. 

 

 

아래의 Dunning-Kruger Effect 그래프를 한번 보시지요.  Know-nothing에서 살짝 오른쪽으로 간 Know-a-little 단계에서 Confidence level이 Peak of "Mt. Stupid" 로 최고 수준을 보인다는 점, 그러다가 조금씩 지식과 경험이 쌓여갈 수록 자신의 역량 수준이 사실은 과대평가된 것이었음을 깨닫고 겸손해진다는점, 자신감의 수준이 지식과 경험이 쌓인다고 선형적으로 급상승하는 것이 아니라 완만한 언덕 형태를 띠면서 오랜 시간을 두고 조금씩 증가한다는 점을 알 수가 있습니다.

 

 

 

 

 

 

제가 분석 업무를 하는데 있어서도 위의 Dunning-Kruger Effect 가 잘 들어맞는 것 같다는 생각을 해보았습니다.  통계책 한권 띠면 분석 문제는 그게 무엇이든지 다 풀 수 있을 것 같은 자신감이 분기탱천하였던 저의 초짜 때의 모습이 떠오릅니다. 

 

그러다가 다변량통계분석으로 넘어가다보니 선형대수, 미적분을 모르고서는 제대로 이해할 수 없다는 것을 알고 공부할 거리가 차고도 넘친다는 것, 저의 수준이 정말 왕초짜라는 것을 알고 '이 길이 내가 가야할 길이 맞나?'라고 겁도 먹었던 시기도 있었습니다.  이론만 가지고는 성과, 결과를 낼 수 없고 R이든 SAS든 SQL이든 HIVE 든 분석 툴도 자유자재로 능수능란하게 사용할 수 있어야 한다는 것을 알고 한번 더 한숨을 내쉬며 이걸 언제 다 배우냐며 부담을 확 느끼기도 했구요.  이랬던 단계가 "Valley of Despair" 에 해당될 거 같아요.

 

퇴근 후 밤마다 12시, 1시까지 공부하고, 주말마다 도서관 가서 책보고 공부하고, R로 데이터 분석 연습도 해보면서 그동안 소련말처럼 보였던 통계, 수학 공식들이 하나, 둘씩 눈에 들어오기 시작하고 이해되기 시작했을 때의 짜릿한 지적 희열, 분석 역량 두뇌 근육이 조금씩 강해지고 있다고 느끼면서 조금씩 오르는 자신감, 더 공부해보고 싶고 더 분석해보고 싶은 욕심과 욕구... 이런 단계가 "Slope of Enlightment" 이지 않을까 싶습니다.

 

Biz. 문제를 접하면 이건 어떤 문제에 속하고, 무슨 분석 기법을 사용해서 이런 저런 데이터를 같이 활용해 어떻게 approach하면 될것 이라는 것을 통찰력을 가지고 제시할 수 있는 단계, 그래서 데이터 기반의 Biz. 문제 해결 및 솔루션까지 제시해줌으로써 성과를 이끌어 낼 수 있는 단계가 "Plateau of Sustainability"가 아닐까 싶은데요, 저도 계속 공부하고 노력하면 어느날은 이 단계에 다다를 날이 오리라 믿고 오늘도 공부하렵니다.  겸손한 맘과 지적 호기심을 채울 때의 기쁜 맘을 같이 가지고서요.

 

 

 

EBS 다큐 중에서 "학교란 무엇인가?" 방송 내용 중 "최상위 0.1%의 공부잘하는 학생들은 무엇이 다른가?"라는 부분에서 실험을 통해 밝히기를, "최상위권 아이들은 자신이 무엇을 알고, 무엇을 모르는지 정확히 파악하고 있었다" 라고 말하고 있습니다.  Dunning-Kruger Effect 의 인지 오류 함정에 빠지지 않는 학생이 최상위 0.1%의 학생이 될 수 있다는 것입니다. 

 

 

● 0.1% 학생들은 뭐가 다르던가.

 “IQ가 평균 134로 일반 학생(125)보다 약간 높긴 했다. 하지만 이 정도 IQ 차이로는 도저히 설명이 불가능했다. 실제로 IQ 상위 0.1% 영재들의 학교 성적은 상위 0.1%에 크게 못 미치는 경우가 대부분이다. 하루 일과도 0.1% 아이들이나 보통 학생들이나 크게 다르지 않았다. 기억력이 딱히 더 좋다고 하기도 어려웠다. ‘0.1%의 비밀’은 이들이 자신이 뭘 알고, 뭘 모르는지 정확히 파악하고 있었다는 점이다. 자신의 생각에 대해 비판적 사고를 하고, 스스로를 객관적으로 바라볼 줄 아는 능력을 ‘메타 인지’라고 한다. 최상위권 아이들은 메타 인지 능력이 아주 뛰어났다.”

 

 

● 메타 인지는 어떻게 길러지나.

 “서울대에 들어간 학생들의 공부법을 보면 스스로 공부하는 시간이 많았다. 0.1% 학생들은 혼자 차분히 복습하는 시간이 월등히 길었다.”

 

 

● 다른 특징은 없었나.

 “학습법뿐 아니라 부모와의 관계에도 차이가 있었다. 0.1% 아이들 중에는 부모와 대화하면 ‘편안함을 느낀다’ ‘즐겁고 유쾌하다’는 등의 긍정적 반응을 보인 비율이 74%나 됐다. 또 하나 차이는 아침식사를 하는 비율이 92%로 압도적이었다는 점이다. 아침을 먹어서 두뇌 활동이 활발해졌기 때문만은 아닐 것이다. 이보다는 아침 식탁에서 가족과 대화하면서 얻게 되는 것이 훨씬 많아 보였다.”

 

 

* 출처 : EBS, 학교란 무엇인가?  PD 인터뷰 中에서 (최낙언의 자료보관소)

 

 

서울대 들어간 학생들 중 차분히 스스로 공부하고 복습하는 시간을 가진 경우가 많았다고 하는데요, 부산한 가운데서도 혼자만의 시간을 내서 복습, 복기, 자기성찰을 해보는 시간을 꾸준이 가져야겠다는 생각을 해봅니다. 

 

 

그리스의 선인들이 했다던 말 "너 자신을 알라 (γνῶθι σεαυτόν 그노티 세아우톤, Know Yourself)"가 얼마나 어려운지, 얼마나 중요한지를 새삼 깨닫게 됩니다.

728x90
반응형
Posted by Rfriend
,

한두개 정도 일회성으로 그래프 그리고 말거면 그냥 화면 캡쳐하는 프로그램 사용하거나 아니면 RStudio의 파일 내보내기를 사용하면 됩니다. 

 

한두번 분석하고 말거면 그냥 마우스로 Console 창 분석결과에 블럭 설정하고 Copy & Paste 하면 됩니다.

 

하지만, 수백개, 수천개의 그래프를 그리고 이를 파일로 저장해야 하고 자동화(사용자 정의 함수, 루프) 해야 한다거나, 분석이나 모형개발을 수백개, 수천개 해야 하고 이의 결과를 따로 저장해야 한다면 이걸 수작업으로 매번 할 수는 없는 노릇입니다.  시간도 많이 걸리고, 아무래도 사람 손이 자꾸 타다 보면 실수도 하기 마련이기 때문입니다. 

 

이에 이번 포스팅에서는

 

(1) ggplot2로 그린 그래프를 jpg 나 pdf 파일로 저장하는 방법

     : ggsave() 

 

(2) Console 창에 나타나는 분석 결과, 모형 개발 결과를 text 파일로 저장하는 방법

     : capture.output()

 

에 대해서 소개하겠습니다.  R script로 위 작업을 수행할 수 있다면 프로그래밍을 통해 자동화도 할 수 있겠지요.

 

 

예제로 사용할 데이터는 MASS 패키지 내 Cars93 데이터프레임의 고속도로연비(MPG.highway), 무게(Weight), 엔진크기(EngineSize), 마련(Horsepower), 길이(Length), 폭(Width) 등의 변수를 사용해서 선형 회귀모형을 만들고 이의 적합 결과를 text 파일로 내보내기를 해보겠습니다.

 

 
> # dataset : Cars93 dataframe,  Weight, MPG.highway variable
> library(MASS)
> str(Cars93)
'data.frame':	93 obs. of  27 variables:
 $ Manufacturer      : Factor w/ 32 levels "Acura","Audi",..: 1 1 2 2 3 4 4 4 4 5 ...
 $ Model             : Factor w/ 93 levels "100","190E","240",..: 49 56 9 1 6 24 54 74 73 35 ...
 $ Type              : Factor w/ 6 levels "Compact","Large",..: 4 3 1 3 3 3 2 2 3 2 ...
 $ Min.Price         : num  12.9 29.2 25.9 30.8 23.7 14.2 19.9 22.6 26.3 33 ...
 $ Price             : num  15.9 33.9 29.1 37.7 30 15.7 20.8 23.7 26.3 34.7 ...
 $ Max.Price         : num  18.8 38.7 32.3 44.6 36.2 17.3 21.7 24.9 26.3 36.3 ...
 $ MPG.city          : int  25 18 20 19 22 22 19 16 19 16 ...
 $ MPG.highway       : int  31 25 26 26 30 31 28 25 27 25 ...
 $ AirBags           : Factor w/ 3 levels "Driver & Passenger",..: 3 1 2 1 2 2 2 2 2 2 ...
 $ DriveTrain        : Factor w/ 3 levels "4WD","Front",..: 2 2 2 2 3 2 2 3 2 2 ...
 $ Cylinders         : Factor w/ 6 levels "3","4","5","6",..: 2 4 4 4 2 2 4 4 4 5 ...
 $ EngineSize        : num  1.8 3.2 2.8 2.8 3.5 2.2 3.8 5.7 3.8 4.9 ...
 $ Horsepower        : int  140 200 172 172 208 110 170 180 170 200 ...
 $ RPM               : int  6300 5500 5500 5500 5700 5200 4800 4000 4800 4100 ...
 $ Rev.per.mile      : int  2890 2335 2280 2535 2545 2565 1570 1320 1690 1510 ...
 $ Man.trans.avail   : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 1 1 1 1 ...
 $ Fuel.tank.capacity: num  13.2 18 16.9 21.1 21.1 16.4 18 23 18.8 18 ...
 $ Passengers        : int  5 5 5 6 4 6 6 6 5 6 ...
 $ Length            : int  177 195 180 193 186 189 200 216 198 206 ...
 $ Wheelbase         : int  102 115 102 106 109 105 111 116 108 114 ...
 $ Width             : int  68 71 67 70 69 69 74 78 73 73 ...
 $ Turn.circle       : int  37 38 37 37 39 41 42 45 41 43 ...
 $ Rear.seat.room    : num  26.5 30 28 31 27 28 30.5 30.5 26.5 35 ...
 $ Luggage.room      : int  11 15 14 17 13 16 17 21 14 18 ...
 $ Weight            : int  2705 3560 3375 3405 3640 2880 3470 4105 3495 3620 ...
 $ Origin            : Factor w/ 2 levels "USA","non-USA": 2 2 2 2 2 1 1 1 1 1 ...
 $ Make              : Factor w/ 93 levels "Acura Integra",..: 1 2 4 3 5 6 7 9 8 10 ...

 

 

 

 

고속도로연비(PMG.highway)와 차 무게(Weight) 간의 산포도를 ggplot으로 그리면 아래와 같이 RStudio 우측 하단의 Plot 창에 그래프가 생성됩니다.

 

 

> # Scatter Plot of Weight & MPG.highway
> library(ggplot2)
> 
> ggplot(Cars93, aes(x=Weight, y=MPG.highway)) +
+   geom_point(shape=19, size=3, colour="blue") +
+   ggtitle("Scatter Plot of Weight & MPG.highway")
 

 

 

 

 

이를 ggsave() 함수를 사용해서 jpg 파일로 저장해서 내보내기를 해보겠습니다.  pdf 파일로 저장하려면 jpg 대신에 pdf 를 사용하면 됩니다.

 

 

> # saving ggplot with jpg format file : ggsave() > ggsave(file="C:/Users/user/Documents/R/scatter_plot.jpg", # directory, filename + width=20, height=15, units=c("cm")) # width, height, units

 

 

 

 

 

단순 선형회귀모형과 다변량 선형회귀모형을 각각 적합시켜보면 아래와 같습니다.

 

 
> # linear regression modeling
> fit_1 <- lm(MPG.highway ~ Weight, Cars93)
> summary(fit_1)

Call:
lm(formula = MPG.highway ~ Weight, data = Cars93)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.6501 -1.8359 -0.0774  1.8235 11.6172 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 51.6013654  1.7355498   29.73   <2e-16 ***
Weight      -0.0073271  0.0005548  -13.21   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.139 on 91 degrees of freedom
Multiple R-squared:  0.6572,	Adjusted R-squared:  0.6534 
F-statistic: 174.4 on 1 and 91 DF,  p-value: < 2.2e-16

> 
> 
> fit_2 <- lm(MPG.highway ~ Weight + EngineSize + Horsepower + Length + Width, Cars93)
> fit_3=stepAIC(fit_2, direction="both")
Start:  AIC=210.77
MPG.highway ~ Weight + EngineSize + Horsepower + Length + Width

             Df Sum of Sq     RSS    AIC
- Horsepower  1      1.38  789.67 208.93
- EngineSize  1      2.62  790.91 209.07
- Width       1      7.03  795.33 209.59
<none>                     788.30 210.77
- Length      1     44.23  832.53 213.84
- Weight      1    562.35 1350.65 258.84

Step:  AIC=208.93
MPG.highway ~ Weight + EngineSize + Length + Width

             Df Sum of Sq     RSS    AIC
- EngineSize  1      1.62  791.29 207.12
- Width       1      7.95  797.62 207.86
<none>                     789.67 208.93
+ Horsepower  1      1.38  788.30 210.77
- Length      1     48.74  838.41 212.50
- Weight      1    699.19 1488.87 265.90

Step:  AIC=207.12
MPG.highway ~ Weight + Length + Width

             Df Sum of Sq     RSS    AIC
- Width       1     13.71  805.00 206.72
<none>                     791.29 207.12
+ EngineSize  1      1.62  789.67 208.93
+ Horsepower  1      0.38  790.91 209.07
- Length      1     52.31  843.61 211.07
- Weight      1    749.41 1540.70 267.09

Step:  AIC=206.72
MPG.highway ~ Weight + Length

             Df Sum of Sq     RSS    AIC
<none>                     805.00 206.72
+ Width       1     13.71  791.29 207.12
+ EngineSize  1      7.38  797.62 207.86
+ Horsepower  1      0.21  804.79 208.69
- Length      1     91.62  896.62 214.74
- Weight      1   1039.48 1844.48 281.82
> summary(fit_3)

Call:
lm(formula = MPG.highway ~ Weight + Length, data = Cars93)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.0988 -1.8630 -0.2093  1.4199 11.3613 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 37.5217809  4.6998055   7.984 4.41e-12 ***
Weight      -0.0096328  0.0008936 -10.780  < 2e-16 ***
Length       0.1155263  0.0360972   3.200   0.0019 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.991 on 90 degrees of freedom
Multiple R-squared:  0.6922,	Adjusted R-squared:  0.6854 
F-statistic: 101.2 on 2 and 90 DF,  p-value: < 2.2e-16

 

 

 

 

이를 capture.output() 함수를 사용하여 text 파일로 내보내서 차곡 차곡 쌓아가면서 저장하는 방법은 아래와 같습니다.  결과 text 파일을 화면캡채해서 같이 올립니다. 

 

> 
> # capture.output()
> # (1) Simple Linear Regression Model (y=MPG.highway, x=Weight)
> cat("\n", 
+     "\n",
+     "==============================================================", "\n", 
+     " [ Simple Linear Regression Model (y=MPG.highway, x=Weight)]  ", "\n", 
+     "==============================================================", "\n", 
+     file="C:/Users/user/Documents/R/lm_MPG_highway.txt", append = TRUE)
> 
> capture.output(summary(fit_1), 
+                file="C:/Users/user/Documents/R/lm_MPG_highway.txt", append = TRUE)
> 
> 
> # (2) Multivariate Linear Regression Model (y=MPG.highway, x1~x5)
> cat("\n", 
+     "\n",
+     "===============================================================", "\n", 
+     " [ Multivariate Linear Regression Model (y=MPG.highway, x1~x5)]  ", "\n", 
+     "===============================================================", "\n", 
+     file="C:/Users/user/Documents/R/lm_MPG_highway.txt", append = TRUE)
> 
> capture.output(summary(fit_3), 
+                file="C:/Users/user/Documents/R/lm_MPG_highway.txt", append = TRUE)

 

 

 

 

 

 

 

 

노가다 하기 싫다면 프로그래밍하고 자동화하는 것이 정답이지요. 

위의 작업을 반복해야 한다면 사용자 정의 함수를 덧붙여서 파일 경로 끝 부분에 파일 이름 부분을 paste() 함수를 써서 사용자 정의 함수에 입력한 값으로 매 루프 돌때 마다 바꿔치기 될 수 있도록 프로그래밍을 해주면 되겠습니다.

 

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

 

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

 

 

728x90
반응형
Posted by Rfriend
,