전산통계학 7주차 강의 - Attributes, class, do.call, RCBD, RCD, scheffe, bonferroni
코딩 공부/R-전산 통계학 2019. 11. 17. 18:08Inf or NaN
1/0
1/Inf
0/0
R에서도 infinite , not a number 이 존재함을 볼 수 있다
is
is.na(0/0)
"na == Not available ; 모른다 , 사용할 수 없다"
is.nan(0/0)
"na > nan ; na가 포함하는 + 더 큰 개념"
x <- c(1, 2, NA, 10, 3)
is.na(x)
is.nan(x)
"두 결과가 다른 이유는? , NA와 NaN의 차이로 생각하면 됨"
<첨언>
명령어를 입력해보고 어떤 결과가 나오는지 확인해보자
x <- c(1, 2, NaN, NA, 4)
is.na(x)
is.nan(x)
<문제>
! 가 붙고, 붙지않고 의 의미를 파악해보자
x <- c(1, 2, NA, 4, NA, 5)
bad <- is.na(x)
bad
x[bad]
x[!bad]
data type (numeric, character, logical)
x <- vector("numeric",length = 10)
x
x <- vector("character",length = 10)
x
x <- vector("logical",length = 10)
x
세 타입의 상하관계 ?
y <- c(1, "a")
y ; "벡터에 숫자(뉴메릭) 와 문자(캐릭터)가 들어가면 문자로 다 인식해버림"
class(y)
y <- c(FALSE, 2)
y ;"로지칼과 뉴메릭이 같이 들어가면 뉴메릭으로 다 인식해버림 , TRUE=1, FALSE=0"
class(y)
y <- c("a", TRUE)
y ; "로지칼 < 캐릭터"
class(y)
as
x <- 0:6
class(x) x의 타입을 확인할 수 있다
as.numeric(x)
as.logical(x)
as.character(x)
as.complex(x)
class(as.complex(x))
<첨언>
다음 명령어를 입력해보고, 실행 여부를 확인 후에 의미를 파악해보자
x <- c("a", "b", "c")
as.numeric(x)
as.logical(x)
attributes
x <- cbind(a = 1:3, pi = pi)
x
attributes(x)
dim(m) "attributes 에 비해서 몇행,몇열인지만 가르쳐줌"
attributes(m) "m 에 대해서 더 자세한 내용"
m <- matrix(nrow = 2, ncol = 3)
m
"missing value ; 빈 공간이다"
<첨언>
캐릭터 사이에도 비교(상하관계) 를 할 수 있을까?
x <- c("a", "b", "c", "c", "d", "a")
x[x > "a"]
u <- x > "a"
u ; "캐릭터 사이도 비교를 해주었다"
<문제>
다음 명령어를 입력하면, 값은 동일하게 나옴을 확인할 수 있다. 도대체 무슨 차이가 있는걸까?
(사진을 보면 알 수 있듯이, 스칼라와 매트릭스 꼴로 나와있다, 어떤 명령어가 그렇게 만들었을까?)
do.call(funname, args) (다른 함수에서 이미 썻던 함수를 불러오는 것)
f<- function(x) print (x^2)
f(5)
do.call(f,list(5))
"만약 f 함수의 변수가 2개였다면 list 안에 값도 2개를 넣어여함"
g<-function(x,y) print(x^2+y)
g(1,5)
do.call(g,list(1,5))
R.C.D ANOVA (randomized complete design)
DATA는 교재의 '514p. ANOVA - RCD ; 석탄문제' 로 할 것이다
#DATA
y=list()
y[[1]]<-c(1.51, 1.92, 1.08, 2.04, 2.14, 1.76, 1.17)
y[[2]]<-c(1.69, 0.64, 0.90, 1.41, 1.01, 0.84, 1.28, 1.59)
y[[3]]<-c(1.56, 1.22, 1.32, 1.39, 1.33, 1.54, 1.04, 2.25, 1.49)
y[[4]]<-c(1.30, 0.75, 1.26, 0.69, 0.62, 0.90, 1.20, 0.32)
y[[5]]<-c(0.73, 0.80, 0.90, 1.24, 0.82, 0.72, 0.57, 1.18, 0.54, 1.30)
y
smanova<-function(y) {
k<-length(y)
n<-c()
t<-c()
for (i in 1:k) {
n[i] <- length (y[[i]])
t[i] <- sum(y[[i]])
}
N<-sum(n)
t..<-sum(t)
tn<-c()
for (i in 1:k) {
tn[i]<- sum(t[i]^2/n[i])
}
sstrt<-sum(tn)-(t..^2)/N
yij<-0
for (i in 1:k) {
for (j in 1:n[i]) {
yij<-yij+(y[[i]][j])^2
}
yij
}
sstot <- yij-(t..^2)/N
sse <- sstot-sstrt
mstrt <- sstrt/(k-1)
mse <- sse/(N-k)
fstar <- (mstrt/mse)
pvalue <- pf(fstar, (k-1), (N-k), lower.tail = F)
cat ("fstar=",fstar," pvalue=",pvalue)
}
smanova(y)
<첨언>
위에서는 데이터를 list를 채워서 넣는 방법이였더라면 (2019년 방식)
, 데이터 안에서 나누는 방법도 존재합니다 (2018년 방식)
x <-c(rep(1,7),rep(2,8),rep(3,9),rep(4,8),rep(5,10))
y<-c(1.51, 1.92, 1.08, 2.04, 2.14, 1.76, 1.17,
1.69, 0.64, 0.9, 1.41, 1.01, 0.84, 1.28, 1.59,
1.56, 1.22, 1.32, 1.39, 1.33,1.54, 1.04, 2.25, 1.49,
1.30, 0.75, 1.26, 0.69, 0.62, 0.90, 1.20, 0.32,
0.73, 0.80, 0.90, 1.24, 0.82, 0.72, 0.57, 1.18, 0.54, 1.30)
z <- data.frame(y,factor(x))
x를 factor 로써 y안의 데이터를 1~5까지 나누어서 분류하고 있음을 알 수 있다
R.C.B.D ANOVA (randomized complete block design)
DATA는 교재의 '535p. 예제 13.5.1 ; ANOVA - RCBD ; 버스+타이어사 문제' 로 할 것이다
#DATA
x<-matrix(c(9.1,17.1,20.8,11.8,
13.4,20.3,28.3,16.0,
15.6,24.6,23.7,16.2,
11.0,18.2,21.4,14.1,
12.7,19.8,25.1,15.8) , nrow = 5, ncol = 4,byrow = T)
smanova<- function (x,alpha1, alpha2) {
b<-nrow(x) ; k<-ncol(x)
t..<-sum(x)
ti.<-apply(x,2,sum)
sstrt<-sum(ti.^2/b) - t..^2/(k*b)
t.i<-apply(x,1,sum)
ssblock<-sum(t.i^2/k) - t..^2/(k*b)
sstotal <- sum(x^2) -t..^2/(k*b)
sserror <- sstotal-(sstrt+ssblock)
"Mean Square"
mstrt <- sstrt/(k-1)
msblock <- ssblock/(b-1)
mserror <- sserror/((k-1)*(b-1))
"F star"
f1star <- mstrt/mserror
f2star <- msblock/mserror
pvalue1 <- pf(f1star,(k-1),(k-1)*(b-1),lower.tail = F)
pvalue2 <- pf(f2star,(b-1),(k-1)*(b-1),lower.tail = F)
ifelse (pvalue1<alpha1, A<-"test1 is reject H0" , A<-"test1 is accept H0")
ifelse (pvalue2<alpha2, B<-"test2 is reject H0" , B<-"test2 is accept H0")
cat("f1star=",f1star, "#1.pvalue=",pvalue1, "result:",A,
"\nf2star=",f2star, "#2.pvalue=",pvalue2, "result:",B)
}
smanova(x,0.1,0.05)
아마 필자의 코딩을 보고 좀 어색하거나, 찝찝~한 명령문들이 몇 개 보였다면 알아차린 분들은 정말 수리통계학을 열심히 배우고 공부하는 사람들이라고 말씀드리고 싶군요... 사실 R상에서는 큰 문제 없이 돌아가지만, 그래도 독자들이 직접 만들때는 올바른 문자를 쓰시길....
(못 알아차렸다면, R.C.B.D 의 데이터 구성은 어떻게 되있고, treatment 와 block 의 범위를 생각하고 T 의 문자를 다시한번 봅시다)
필자가 했던 밑의 <문제>도 같은 문자를 썼으니, 독자들은 꼭 올바르게 고쳐서 코딩을 해봅시다
<문제>
수리 통계학 시간 때, 아마도 reject H0 이후에 각 그룹별로 어떻게 차이나는지 검증하는 '1. Bonferroni' 와 '2. Scheffe' 를 배웠을 것이다
이 것도 프로그램화 시킬 수 있을 것이다
2019년에는 '1. Bonferroni' 를 했었으니, 예시로 보여주겠다 (그러니까, 2. Scheffe 를 해보도록 합시다~ ^^)
smanova<- function (x,overallalpha,c) {
b<-nrow(x) ; k<-ncol(x)
t..<-sum(x)
ti.<-apply(x,2,sum)
sstrt<-sum(ti.^2/b) - t..^2/(k*b)
t.i<-apply(x,1,sum)
ssblock<-sum(t.i^2/k) - t..^2/(k*b)
sstotal <- sum(x^2) -t..^2/(k*b)
sserror <- sstotal-(sstrt+ssblock)
"Mean Square"
mstrt <- sstrt/(k-1)
msblock <- ssblock/(b-1)
mserror <- sserror/((k-1)*(b-1))
"F star"
f1star <- mstrt/mserror
f2star <- msblock/mserror
pvalue1 <- pf(f1star,(k-1),(k-1)*(b-1),lower.tail = F)
pvalue2 <- pf(f2star,(b-1),(k-1)*(b-1),lower.tail = F)
ifelse (pvalue1<(overallalpha/c), A<-"test1 is reject H0" , A<-"test1 is accept H0")
ifelse (pvalue2<(overallalpha/c), B<-"test2 is reject H0" , B<-"test2 is accept H0")
cat("f1star=",f1star, "#1.pvalue=",pvalue1, "result:",A,
"\nf2star=",f2star, "#2.pvalue=",pvalue2, "result:",B)
}
smanova(x,0.1,2)
TAG
