전산통계학 3주차 강의 - dnorm, pnorm, qnorm, rnorm, F test, T, var test
코딩 공부/R-전산 통계학 2019. 11. 14. 15:37Statistical Distributions
dnorm
dnorm(x,mean,sd) sd는 standard devaition 이므로 sigma^2 의 분포를 따르면 sigma 값만 적어줘야만 한다
dnorm(8,5,10) 교수님이 적어주신 f(x=n)=~~ 식만 보고도 구할 줄 알아야 한다"
pnorm
"X~N(5,10^2)을 따를 때, 분포표 상에서 우측에서부터 0.025의 영역 너비를 가지게끔 하는(=P(X,x0)=0.975)) 좌표 x0를 구할수 있는가?"
위의 0.975 식을 가지게끔 하는 식을 다시 변경하면, P(z < (x0-5)/10 )=0.975
그런데 (x0-5)/10 = 1.96 <- standard normal 에서 찾는 거고, 이 값을 식에 대입해서 구하면 된다.
이거를 norm 식으로 구할수 있을까???
pnorm = "스탠다드 노말 분포표에서 왼쪽 끝에서부터 1.96까지의 좌표 영역너비를 뜻한다"
pnorm(1.96,0,1)
pnorm(1.96,0,1, lower.tail = T)
pnorm(1.96,0,1, lower.tail = F) "lower.tail = 오른쪽 끝에서부터 영역 계산 ; 0.025값이 나온다."
<첨언>
응용문제로 -1.28~2.3 까지의 영역너비도 구할 수 있을까?
pnorm(2.3,0,1)-pnorm(-1.28,0,1) "0.9893-0.1003 을 직접적으로 계산해서 구함"
qnorm (pnorm 과 짝꿍이다)
qnorm = "식에 너비값을 대입하면 그만큼 차지하게끔 하는 좌표가 나온다"
qnorm(0.88)
<문제>
응용문제로 가운데에서 부터 왼쪽 너비는 0.4 ; 오른쪽 너비는 0.32가 되게끔 하는 좌표를 찾고 , 두 점 사이의 거리도 같이 구하여라
x1<-qnorm(0.82)
x2<- -qnorm(0.9) (또는 qnorm(0.1) 에 대해서 생각해도 상관없다))
x1-x2
<문제>
좀 더 심화시켜서 응용문제로 이제는 standard normal 이 아닌 N(5,10^2) 을 따를 때,
가운데 영역이 0.94가 되게끔 되는 양쪽 점 좌표 및 거리도 구하여라
qnorm(0.94/2+0.5)
qnorm(0.94/2+0.5,lower.tail=F)
rnorm
rnorm(n,mean,sd) <- "X1...Xn ~ N(u,sigma^2) 을 의미함
rnorm(10,5,10)
t (t distribution)
qt ="t분포에서 너비를 입력하면 그만큼의 좌표가 나옴"
qt(0.023,21)
qt(0.99,21)
qt(0.5+0.31,15) ; qt(0.5-0.363,15) 이 두 명령어는 그럼 어떤 것을 구한 것인지 쉽게 알 수 있다
<첨언>
pt 는 어떤 의미일까? (norm 에서 이미 다뤘으므로, 추측 가능)
f distribution
"f분포는 degree of freedom 이 2개임을 수리통계학 시간에 배웠다"
qf(0.85,10,15)
qf(0.035,10,15) <-"qt(1-0.035,10,15, lower.tail=F) 와 같이 써도 무관"
<문제>
응용문제로 시작점부터 0.1만큼의 면적을 가지게끔 하는 점에서 1.5의 거리뒤에 있는 좌표 x1을 구하고, x1의 오른쪽 영역 너비를 구하여라
x1<-qf(0.1,15,20)+1.5
x1
1-pf(x1,15,20)
alpha & d.f ?
alpha = 0.1 이라고 가정해보자
"수업시간에 반대영역은 역수하고 degree of freedom 을 뒤엎은 것임을 알고있다."
a<-qf(0.9,10,15)
b<-qf(0.1,15,10)
all.equal(a,1/b) "두 값이 같은지를 확인"
<문제>
"X~B(n=100, P=0.1" 을 따를 때," P(8<=X<=12) = f(X=8)+ ~ + f(x=12) = 바이노미얼 식" 을 함수로 만들 수 있겠는가?
choose(3,2)
a=dbinom(8,100,0.1)
b=choose(100,8)*(0.1)^8*(0.9)^92
all.equal(a,b)
pbinom(12,100,0.1) - pbinom(7,100,0.1)
함수식 만들기
sa<-function(n,p,x1,x2){
yo=0
for(i in x1:x2) {
yo=yo+choose(n,i)*(p^i)*(1-p)^(n-i)
}
cat(yo)
}
sa(100,0.1,8,12)
( Q ) cat은 생성한 함수의 결과값을 어떤 것으로 정리해서 나오게 할지 정하는 것이다, 따라서 cat을 안 쓴다면 함수의 마지막에 'yo'를 적자)
making function (program)
우선 x와 y의 데이터를 다음과 같이 주도록 하자 (교재 361p. 연습문제 13번)
#DATA
x<-c(1.79,1.87,1.62,1.96,1.75,1.74,2.06,1.69,1.67,1.94,1.33,1.70,1.65)
y<-c(2.39,2.56,2.36,2.62,2.51,2.29,2.58,2.41,2.86,2.49,2.33,1.94,2.14)
<문제>
"x가 1.5보다 작으면 0 ; 1.5이상~1.8미만이면 1 ; 1.8이상이면 2 (3개 조건부)" 가 되게끔 하는 프로그램을 짜보아라
x<-c(1.79,1.87,1.62,1.96,1.75,1.74,2.06,1.69,1.67,1.94,1.33,1.70,1.65)
for(i in 1:length(x)) {
if(x[i]<1.5) {x[i]<-0}
else if (x[i]>=1.5 & x[i]<1.8) {x[i]<-1}
else {
x[i]<-2
}
}
x
이번에는 "mm(x) 만 넣으면 mean값이 나오게끔 되는 함수 만들기"
mm=function(x){
pro=0
n=length(x)
for (i in 1:n) {
pro=pro+x[i]
}
cat(pro/n)
}
mm(x)
"이번에는 s^2 이 나오게끔 하는 함수 식 만들기"
omg=function(x){
hot=0
for(i in 1: length(x)) {
hot=hot+(x[i]-mean(x))^2
}
cat(hot/(n-1))
}
omg(x) ; var(x)
Hyper Testing
"우선 두 그룹의 분산이 같은지 sigma 1^2 == sigma 2^2 임을 확인해보자
,그리고 f* 를 구해서 C.R 에 넣어 H_o 에 부정하는지 수긍하는지 에 대한 프로그램을 짤것이다"
1. "f*를 구하는 알고리즘"
fstar<-(var(x)/var(y)) x,y 데이터는 위에 적혀있으니, 다시 입력해주어도 좋다
fstar
2. "Let a=0.05, f분포표에서 양쪽의 값을 구해보자"
alpha=0.05
q1=qf((alpha/2), (length(x)-1),length(y)-1) ; qf(0.025,12,12)
q2=qf((alpha/2), (length(x)-1),length(y)-1, lower.tail = F) ; qf(0.975,12,12)
q1
q2
3. "0.025~0.975 구간 안에 들어가면 accept H_o , 아니면 reject H_o 이게끔 하자"
ifelse (fstar>=q1 & fstar<=q2, "accept H0", "reject H0")
<문제>
위의 1~3 단계를 function 으로 만들어보자
sm <- function(x1,x2,alpha) {
fstar<-(var(x1)/var(x2))
q1=qf((alpha/2), (length(x1)-1),length(x2)-1)
q2=qf((alpha/2), (length(x1)-1),length(x2)-1, lower.tail = F)
ifelse (fstar>=q1 & fstar<=q2, "accept H0", "reject H0")
}
sm(x1,x2,0.05)
var.test
"위에서는 (sigma 1^2 / sigma 2^2) = 1 (two-sided) 을 확인한 것이라면, 통상적으로 배울때 쓰는 l('엘')로 두고 한번 해보자"
?var.test ?+명령어(문) = R studio 내부의 가이드북을 볼 수 있도록 하니, 많이많이 쓰도록 하자
var.test(x, y, ratio = 1, alternative = c("two.sided", "less", "greater") 명령어에서 각각 테스트의 1,2,3 중에 골라서 쓰는 것이다.
conf.level = 0.95, ...) 은 유의 수준을 말한다.
var.test(x1,x2,ratio=1,alternative=c("two.sided"))
명령어를 입력하면 fstar, degree of freedom, p-value, 1-a 에 해당하는 점 값 들이 나온다
뒤에 $를 붙이고 추가명령어를 쓰면 그 값만 나온다.
var.test(x1,x2,ratio=1,alternative=c("two.sided"))$'statistic'
var.test(x1,x2,ratio=1,alternative=c("two.sided"))$'p.value'
<문제>
"f분포에서 0.5영역에 해당하는값 = b 로 두고 two-sided 에서 작은지 큰지를 구분 후, p.value 를 *2 해서 찾는 함수를 생성해라"
psm <- function(x1,x2,alpha) {
fstar<-(var(x1)/var(x2))
b <- qf(0.5,(length(x1)-1),length(x2)-1)
if(fstar < b) {oh<- pf(fstar,(length(x1)-1),length(x2)-1)}
else if (fstar > b) {oh<-(1- pf(fstar,(length(x1)-1),length(x2)-1))
}
oh*2
}
psm(x1,x2,0.05)
all.equal (var.test(x1,x2,ratio=1,alternative=c("two.sided"))$'p.value', psm(x1,x2,0.05) )
그렇다면 이번에는 confidene interval (= 신뢰구간) 양쪽을 구하는 함수도 만들어 볼 수 있겠죠?
TAG
