$$config title: R keywords: - IT - ... $$ # 3/20 ``` Sys.getenv("jAVA_HOME") Sys.getlocale() ls() # 대입(<-), equal(=) x <- 5 x typeof(x) # double mode(x) # 데이터타입 class(x) # 자료구조 # 정수형으로 변환 x = as.integer(x) # R에서는 .이 이름으로 사용(.앞이 클래스이름 아님 뒤에 함수이름 아님님) y <- x; typeof(y) x <- 10 # 값에 의해 전달을 하고 있음(주소X) y #5 x #10 z="대한민국" print(z) z x+y x-y x*y y/x y%/%x # 몫연산자 y%%x # 나머지 연산자 y^x x+z # 자바라면 암묵적 캐스팅으로 문자로 출력, R은 연산이 우선선 # 비교(관계)연산자 (TRUE/FALSE) x <- 5 y <- 16 xy x<=5 y>=20 y==16 x!=5 x <- c(TRUE,FALSE,0,6) # combine -> vector사용 # 숫자형 #1 1 0 6 x (x <- c(TRUE,FALSE,0,6)) # 가장 바깥쪽에 괄호 붙이면 바로 출력 (y <- c(FALSE,TRUE,FALSE,TRUE)) # boolean형 #FALSE TRUE FALSE TRUE typeof(x) #double mode(x) class(x) typeof(y) #logical mode(y) class(y) x+1 #2 1 1 7 # for문을 사용하지 않음, vector화 연산 !x #FALSE TRUE TRUE FALSE # 데이터를 TRUE/FALSE로 자동 변환하고 부정 x&y #FALSE FALSE FALSE TRUE # 개별 연산 x&&y # 통합연산 -> 모두가 참이면 참인데 하나라도 FALSE면 FALSE x y # 숫자를 벡터화할 때 v <- 2:8 print(v) #2 3 4 5 6 7 8 v1 <- 8 v2 <- 12 t <- 1:10 v1 %in% t # 매칭연산자 #TRUE print(v1 %in% t) #TRUE print(v2 %in% t) #FALSE # 행렬(vector로부터 생성) # 2x3 행렬 (M = matrix(c(2,6,5,1,10,4), nrow = 2, ncol = 3, byrow = TRUE)) # 행방향으로 t(M) # transpose 전치행렬-행과 열의 순서를 변경(전치행렬을 사용하는 이유-행렬곱 때문에) # 행렬곱( %*% ) -> 앞의 행렬의 열 수와 뒤의 행렬의 행 수가 일치해야 함 # 결과는 앞의 행렬의 행수 * 뒤의 행렬의 열수 # 앞의 행렬(2x3) ->2 # 뒤의 행렬(3x2) ->2 t = M %*% t(M) print(t) # 결과행렬은 2x2 # 행렬의 거듭제곱을 하려면 행렬의 전치행렬을 구해서 행렬곱을 해야한다 # M%*%M # 2x3 ->2 2x3 ->3 => 곱할 수 없음 num <- 6 num ^ 2 num num <- num^2*5+10/3 num # scalar량 c(1,3,5,7,9) * 2 #2 6 10 14 18 c(1,3,5,7,9) * c(2,4) # recycling 배수가 되어야 함 (짝이 맞아야 함) c(1,3,5,7,9,10) * c(2,4) #2 12 10 28 18 40 c(2,4) * c(1,3,5,7,9,10) # 리싸이클링을 통해서 개수를 맞춤 factorial(1:5) #1 2 6 24 120 # 함수의 매개변수는 기본이 벡터 exp(2) # 밑을 자연대수로 하는 지수 exp(2:10) sqrt(c(1,4,9,16)) cos(c(0, pi/4)) ?exp 1/0 #Inf 무한대 0/0 #NaN 수가 정의되지 않음, NA는 결측치 sum(c(1,2,NA,3)) # 연산에서 결측치가 하나라도 있으면 결과는 NA # 반드시 결측치 처리를 해야함 # 기본이 수치는 double vec <- c(0, Inf, NaN, NA) typeof(vec) #double mode(vec) #numeric class(vec) #numeric # ↓확인함수 is.finite(vec) #TRUE FALSE FALSE FASLE is.nan(vec) #FALSE FALSE TRUE FALSE is.na(vec) #FALSE FALSE TRUE TRUE is.infinite(vec) #FALSE TRUE FALSE FALSE sum(vec) #NaN # 연산자함수 작성 '%divisible%' <- function(x,y) { if(x%%y == 0) return(TRUE) else return(FALSE) } 10 %divisible% 3 10 %divisible% 2 '%divisible%'(10,5) x <- c(5,6,7) # vector y <- c(1,2,3) x%%y #0 0 1 x%*%y #38 (5+12+21) # 행렬곱,,벡터곱을 내적이라고 함 # 내적의 의미: 크기값을 고려한 사이각 # 같은 위치의 요소를 곱한 다음 다 더해주면 내적이다. # 내적이 0이면 직교 x*y #5 12 21 sum(x*y) #38 # R은 대소문자를 구분 LETTERS letters month.abb #Jan Feb Mar ... # 축약어 month.name pi ls() # 메모리에 있는 변수 확인 int <- 50 # int를 변수명으로 사용해도 상관x (데이터타입이 자료형으로 사용가능) # 통계함수 score <- c(85,95,75,65) score mean(score) # 평균 (avg <- sum(score)/4) sum(score) # 합 var(score) # 분산 = sum((score - avg)^2)/3 => 분산이 큰 변수가 중요한 변수(주성분) sd(score) # 표준편차 = sqrt(sum((score - avg)^2)/3) # recycling score -avg sum(score - avg) #0 (score - avg)^2 # 평균값의 차의 제곱이 분산(부호를 떼기위해서 제곱함) var_value = sum((score - avg)^2)/3 #166.6667 sqrt(var_value) #12.90994 표준편차 -> 표준오차 => 신뢰구간 score <- c(85,95,NA,75,65) score sum(score) #NA mean(score,na.rm = TRUE) #80 sum(score, na.rm = T) #320 # TRUE를 T로 쓸 수도 있음 x <- c(10,20,30,'40') mode(x) #character typeof(x) (xx <- as.numeric(x)) #10 20 30 40 mode(xx) #numeric # 정성적 데이터(질적), 정량적 데이터(양적) : 이산적, 연속적 # factor형 데이터, 범주형 데이터 gender <- c('M','F','M','F','M') # male, female gender mode(gender); class(gender) plot(gender) fgender <- as.factor(gender) fgender # Levels: F M fgender[2] # 범주형은 그래프 출력 가능능 mode(fgender); class(fgender) levels(fgender) plot(fgender) # 숫자로 매핑됨 Sys.Date() Sys.time() today <- '2017-03-11 10:15:24' # 입출력은 문자열 계산은 특별숫자계산 (today2 <- as.Date(today)) ctoday <- strptime(today, '%Y-%m-%d %H:%M:%S') # 문자열을 날짜로 분해해 ctoday args(sum) example(sum) Nile # 시계열 데이터 length(Nile) plot(Nile) ``` ``` # vector : 배열개념으로 저장(ArrayList) # 왜 vector인가 - 데이터를 선형대수로 처리 # Vector : 동질적 데이터를 담음 # vector의 기본연산(+,-,*,/:요소별로 연산) # 내적과 외적: 분석에서는 외적은 사용하지 않음(게임에서:면에 직교하는 두 벡터가 이루는 면에 수직벡터를 이용해서 음영을 줄 때 사용) # 내적은 요소끼리 곱한다음 다 더하면 크기가 고려된 사이각으로 표현할 수 있다 # A,B = |A| |B| cos theta # norm : 벡터의 크기 # normalize : 정규화, 방향값으로 표현되고 사이즈는 1이다. # 두 개의 벡터인 경우 방향은 내적, 크기는 피타고라스 정리 # 벡터:크기와 방향을 표현 # 가격벡터 수량벡터의 내적 매출이 됨 # 유사도(추천)와 거리(군집분석) # 사이각이라면 각도는 사이즈가 1이다. -> 곱셈의 항등원 # 방향을 표시할 때는 사이즈 1인 것으로 표현 # 덧셈의 항등원 0 x <- c(1,5,4,9,0) # 벡터는 위아래로 저장 class(x) typeof(x) length(x) NROW(x) #5 NCOL(x) #1 # vector -> matrix(2차원) -> array:동질적인 데이터 # vector -> 이질적인 요소로 구성된 dataframe # 열로 열내 동질적, 열간 이질적 # matrix, array, dataframe은 사각형으로 이빨이 빠지면 안됨(사각형 모양 유지해야됨) # list는 이빨이 빠져도 됨 (x=vector("list",10)) (x=vector("numeric",10)) rm(x) (y <- 2:-2) seq(1, 3.2, by=0.2) # 간격을 지정 seq(1, 5, length.out=4) # 개수를 지정정 seq(1, 6, length.out=4) (b <- rep(1:4,2)) #1 2 3 4 1 2 3 4 # 반복 (d <- rep(1:3, each=3)) #1 1 1 2 2 2 3 3 3 (d <- rep(1:3, 2, each=3)) #1 1 1 2 2 2 3 3 3 1 1 1 2 2 2 3 3 3 # each가 먼저 작동 (x <- c(1:9)) x[3] #3 x[c(2,4)] #2 4 # 벡터를 이용해서 인덱싱 x[-1] #2 3 4 5 6 7 8 9 # 첫 번째를 제외하고 x[c(2, -4)] # 복합 안됨 x[c(2.4, 3.54)] # 소수점은 절사 x[c(TRUE, FALSE, FALSE, TRUE)] #1 4 5 8 9 # Boolean indexing # 숫자가 안맞으면 recycling x[x <= 3] #1 2 3 x[x > 3] #4 5 6 7 8 9 # 관계연산자를 이용한 인덱싱 (x <- c("first"=3, "second"=0, "third"=9)) x x["second"] # 키데이터 접근방법 x[c("first","third")] ``` ``` # 비교 c(1,2,3) == c(1,2,100) x <- 1.5 all(x>2) # 모두 조건을 만족하는가 any(x>2) # 하나라도 조건을 만족하는가 # 정렬 : 원본에 영향을 미치지 않음 # 데이터를 정렬 x <- c(3,2,6,4,5,8,1) x sort(x) sort(x, decreasing=TRUE) # 내림차순 x # 인덱스 정렬(하나의 열을 기준으로 다른 데이터를 정렬하고 싶은 경우) # 두 개의 열에 대하여 하나는 기준 하나는 정렬된 것을 참조 # 국어 점수를 기준으로 이름을 참조 order(x) order(x,decreasing=TRUE) # 결과가 인덱스 값으로 나옴옴 x[order(x)] # 값으로 나오고고 # NA결측치 sum(2,7,5) x = c(2,NA,3,1,4) sum(x) sum(x, na.rm=TRUE) mean(x, na.rm=TRUE) median(x, na.rm=TRUE) # 중위수 - 통걔학(이상치를 고려하기 때문에) prod(x, na.rm=TRUE) # 벡터에 이름 정하기 vectorA <- c(1,2,3,4) names(vectorA) <- c("국어","영어","수학","과학") vectorA vectorA["국어"] vectorA[2] vectorA[-1] vectorA[c(1,3)] # 벡터 인덱싱(subset) vectorA[vectorA > 5] vectorA[c(FALSE,FALSE,TRUE)] append(vectorA, c(3,4,5)) vectorA # 원래 데이터에 영향을 미치지 않습니다. : 값에 의한 전달 (vectorB = append(vectorA, c(3,4,5))) # 집합연산(set:중복이 없음음) union(vectorA, vectorB) # 합집합 intersect(vectorA, vectorB) # 교집합 setdiff(vectorA, vectorB) #numeric(0) # 차집합 - 순서가 있음 setdiff(vectorB, vectorA) #5 setequal(vectorA, vectorB) # 체크 함수 # 인덱스를 리턴하는 것들은 하나를 기준으로 다른 데이터를 참조할 때 (x <- c(3,2,6,4,5,8,1)) subset(x, x>3) #6 4 5 8 which(x*x > 8) #1 3 4 5 6 # 인덱스로 찾을 경우 which x[which(x*x > 8)] #3 6 4 5 8 # 값으로 찾고 싶을 때 # 확률(표본을 뽑는 행위) # 우리나라의 키의 평균 # 전국민을 대상으로 키를 측정 => 비용과 시간이 많이 들어 # 그래서 표본을 추출(확률) # 복원추출 (x <- c(sort(sample(1:3,3,replace=T,prob=c(0.1,0.3,0.2))))) # 1부터3중에 3개를 뽑아라 nums <- c(5,8,10,NA,3,11) nums which.min(nums) # 결과가 인덱스 which.max(nums) nums[which.min(nums)] # 결과가 값 nums[which.max(nums)] # 문제 (for문을 사용해서 프로그래밍 => 인덱스):if문 => [], for문은 자동동 # vector1 벡터 변수를 만들고,"R"문자가 5회 반복되도록 하시오 (vector1 <- rep("R",5)) # vector2 벡터 변수에 1~20까지 3간격으로 연속된 정수를 만드시오. (vector2 <- seq(1,20,by=3)) # (vector2 <- seq(1,20,3)) # vector3에는 1~10까지 3간격으로 연속된 정수가 3회 반복되도록 만드시오. (vector3 <- rep(seq(1,10,by=3),each=3)) # (vector3 <- rep(seq(1,10,3),3)) # vector4에는 vector2~vector3가 모두 포함되는 벡터를 만드시오 (vector4 <- append(vector2,vector3)) # 25~ -15까지 5간격으로 벡터 생성 - seq()함수 이용 (seq(-15, 25, by=5)) # vector4에서 홀수번째 값들만 선택하여 vector5에 할당하시오(첨자[] 이용) (vector5 <- vector4[seq(1,length(vector4),2)]) # vector5에서 짝수 번째 값들만 선택하여 vector6에 할당하시오(첨자 이용) (vector6 <- vector5[seq(2,length(vector5),2)]) # vector5의 데이터를 문자로 변경하여 vector7에 저장하시오(as.character()) (vector7 <- as.character(vector5)) # vector5와 vector6을 벡터 연산하여 더하시오 vector5 + vector6 # vector5와 vector7을 벡터 연산하여 더하시오 (as.numeric()) vector5 + as.numeric(vector7) # package : base packages(R studio가 로딩할 때 자동으로 로딩), # recommended packages(설치는 되어 있는데 로딩이 안된 패키지), # others(설치도 안되어 있고 로딩이 안된 것들) #install.packages("NISTunits") # 설치 library(NISTunits) # 로딩 NISTdegTOradian(180) # 각도를 라디안 바꿔줘야함 # 각도에서 sin(x), cos(x)의 값의 범위 -1 ~ 1 사이의 값은 가짐 # 반지름이 1인 원에서의 계산 NISTradianTOdeg(pi) #180 # 문제 # 45도를 라디안으로 변경해보시오 ang <- 45 rad1 <- NISTdegTOradian(ang) # 밑변의 길이와 높이를 구하시오(반지름이 1) : 곱할 때의 항등원원 cosbi <- cos(rad1) # 밑변 (대각선과 밑변의 비) sin(rad1) # 높이 # 대각선과 밑변의 비를 각도로 변경하시오 cosrad <- acos(cosbi) # 라디안으로 변경 NISTradianTOdeg(cosrad) # 삼각함수, 지수함수, 로그함수 -> 초월함수 ( +,-,*,/ 로 값을 구할 수 없음:즉 대수적으로 정확한 값이 없음) # 극한 # 회전, 주기적으로 값이 변화하는 경우 # 독립변수의 영향을 균등하게 하기 위해서 정규화 함 # 벡터 -> 정규화(nomalize) # min-max정규화, normalize, z점수 정규화 # 모든 값을 0 ~ 1사이의 값으로 변환 # 벡터의 크기값을 구해서 벡터의 요소값을 나누어 주면 사이즈 1인 벡터 # -1, 1값을 가짐 # -oo ~ oo 무한대 ``` # 3/21 ``` # matrix # 동질적인 vector가 열로 여러개 모여있는 행렬 # 데이터 변환, 데이터 처리 # 공간매핑 : MRS # 데이터 처리: 차원감소(시각화:3차원(컬리값, 선타입)), 특징추출 # 시각화 없이 데이터를 이해하는 것은 불가능 ( 2차원 3차원 ) # 특징이 많으면 복잡 -> 특징 추출 # 클래스를 함수로 만들어 함수베이스 프로그래밍이다. # 사각형을 유지해야 함함 (a = matrix(1:9, nrow = 3, ncol = 3, byrow = TRUE)) matrix(1:9, nrow = 3) // ncol은 자동으로 계산이 되어짐 class(a) (b = matrix(1:9, nrow = 3, ncol = 3, byrow = FALSE)) attributes(a) # 차수(dimension) x <- (matrix(1:9, nrow=3)) colnames(x) colnames(x) <- c("C1","C2","C3") rownames(x) <- c("R1","R2","R3") x #x["C1"] #x["R1"] # 왜 dimnmaes는 list로 전달이 될까 # list는 열별 행별 개수가 달라도 됨 (x <- matrix(1:9, nrow=3, dimnames = list(c("X","Y","Z"),c("A","B","C")))) # 선형대수를 동질적 데이터에만 가능 - 이질적인 불가능 # dataframe에는 행렬 연산이 적용되지 않음 # 빈 매트릭스를 만들고 값을 지정 y <- matrix(nrow=2, ncol=2) y[1,2] y[1,1] <- 1 y[2,1] <- 2 y[1,2] <- 3 y[2,2] <- 4 y mode(y) class(y) # 두 개의 매트릭스를 결합할 때 (cbind(c(1,2,3),c(4,5,6))) # column (rbind(c(1,2,3),c(4,5,6))) # row bind # matrix의 attribute는 dim이 있음 (x <- c(1,2,3,4,5,6)) # vector dim(x) # dim정보가 없음 (dim(x) <- c(2,3)) x (x <- matrix(1:9, nrow=3, byrow=TRUE)) (y <- matrix(11:19, nrow=3, byrow=TRUE)) (c <- rbind(x,y)) (d <- cbind(x,y)) (x <- matrix(1:9, nrow=3, byrow=TRUE)) x[c(1,2),c(2,3)] x[c(2,3),] # 없으면 모두 다 출력력 x[,] x[-1,] # 생략 a = x[1,] class(a) typeof(a) x[c(TRUE,FALSE,TRUE),c(TRUE,TRUE,FALSE)] x[c(TRUE,FALSE),c(2,3)] x[c(TRUE,FALSE)] # 하나의 벡터로 보고 데이터 추출 x[x>5] # 데이터를 1차원으로 보고 출력 x[x%%2 == 0] (mdat <- matrix(seq(20,4,-2), nrow=3, ncol=3, byrow=TRUE, dimnames=list(c("a","b","c"),c("x","y","z")))) t(mdat) nrow(mdat) ncol(mdat) dim(mdat) # 차수 rowMeans(mdat) # 행방향으로 평균 rowSums(mdat) # 열개수 만큼 (행이 없어지는 것) colMeans(mdat) # 행개수 만큼 colSums(mdat) diag(mdat) # 행렬의 대각선에 있는 값의 의미:크기를 결정 diag(20,12,4) # diag(값, 행, 열) => 대각선으로 값이 써짐 # 벡터는 크기와 방향값을 표현 eigen(mdat) # 고유값 분해(크기와 방향을 분해) # 고유치 + 고유벡터 # 고유치는 축 방향으로의 크기 # 고유벡터는 축의 방향을 의미 # 정방행렬이면서 대칭행렬(상관계수행렬,거리값행렬)에 대하여 고유값 분해를 하면 고유치와 고유벡터 출력 # 고유벡터는 othogonal한다.( 정직교한다.) -> 고유값 분해라고 함 # 고유값 분해를 하여 생성되는 고유벡터에 데이터를 재표현하면 noise가 사라진다. # 주성분 분석(Principal Component Analysis) svd(mdat) # singular value decomposition(특이행렬분해:희소행렬에서 특징추출) # 비정방행렬에 사용 # 비정방행렬에는 정방행렬도 포함 # 희소행렬(sparce matrix : 요소가 0이 많은 경우) # 텍스트 분석에서 사용하는 단어 벡터 # y = ax # y = 10x 일 때 y값이 20이라면 x값은 solve(10,20) # solve(a,y) # 2x + 3y = 5 # 3x + 5y = 6 (mdat <- matrix(c(2,3,3,5), nrow=2, ncol=2, byrow=TRUE)) (c=c(5,6)) solve(mdat,c) # x = 7, y = -3 2*7 + 3*-3 3*7 + 5*-3 # 문제 # 연립방정식의 해를 구하시오. 계산과정에서 발생하는 매트릭스의 # 행의 합, 행의 평균, 열의 합, 열의 평균을 구하시오. # x+3y+5z=10 # 2x+5y+z=8 # 2x+3y+8z=3 #c(1,3,5,2,5,1,2,3,8) [x,y,z] = [10,8,3] #[x,y,z] = [10,8,3]/c(1,3,5,2,5,1,2,3,8) #[x,y,z] = [10,8,3] %*% inverse matrix # 행렬은 나눗셈이 없다 -> 역행렬을 이용해서 곱하면 됨. x=c(1,3,5,2,5,1,2,3,8) b = c(10,8,3) (mat <- matrix(x, nrow=3, byrow=T)) solve(mat,b) # -9.28 5.16 0.76 첫매개변수는 역행렬로 구하고 다음 매개변수와 행렬곱 1*-9.28 + 3*5.16 + 5*0.76 (inv_m <- solve(mat)) # 역행렬 inv_m %*% b (identity_m <- inv_m %*% mat) # 단위행렬 1과 마찬가지 identity_m %*% inv_m # 행렬의 항등원은 단위행렬 # 덧셈의 항등원은 0이고, 곱셈의 항등원은 1이고, 행렬의 항등원은 단위행렬이다. rowMeans(mat) rowSums(mat) colMeans(mat) colSums(mat) # Matric Operations 문제 # 다음 두 매트릭스를 rbind, cbind를 이용하여 묶으시오 # 다음 두 매트릭스의 사칙연산을 구하시오 # 다음 두 매트릭스의 행렬내적(곱: %*%)를 구하시오 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 (a <- matrix(1:15, nrow=3, byrow=T)) (b <- matrix(16:30, nrow=3, byrow=T)) rbind(a,b) # 6x5 cbind(a,b) # 3x10 a+b # 행렬 사칙연산(요소끼리) a-b a*b a/b # a%*%b -> 앞에 있는 행렬의 열 수와 뒤에 있는 행렬의 행 수가 일치하지 않아서 에러남 # 내적을 내는 의미는 투영의 의미가 있음 (a %*% t(b)) # 앞에 있는 행렬의 행벡터에 대하여 뒤에 있는 행렬의 열개수만큼 별도의 특징 추출을 해서 데이터를 변환 # 2000 x 3 => character # 앞에 있는 행렬이 캐릭터, 뒤에 있는 행렬이 회전값(변환값) # 3 x 3 => 결과: 2000 x 3 # cos theta 0 sin theta # 새로운 축에 투영된 값을 다시 구하면 # 0 0 0 # -sin theta 0 cos theta ``` ``` # Array : 기준은 행렬 vector1 <- c(5,9,3) vector2 <- c(10,11,12,13,14) # combine해서 두개를 묶었엉 # 24개 데이터 필요함 -> 근데 실제로는 8개만 들어왔으니 recycling 해야함 result <- array(c(vector1,vector2),dim = c(4,3,2)) # 면행렬이 행렬면으로 표현 # dim = c(행,열, 개수) print(result) # 3차원 데이터 column.names <- c("COL1","COL2","COL3") row.names <- c("ROW1","ROW2","ROW3") matrix.names <- c("Matrix1","Matrix2","Matrix3","Matrix4","Matrix5") result <- array(c(11:19, 21:29, 31:39, 41:49, 51:59),dim = c(3,3,5), dimnames=list(row.names, column.names,matrix.names)) print(result) print(result[3,,2]) # 2번째 매트릭스에 3번째 행 print(result[,"COL1","Matrix1"]) # 1번째 매트릭스에 1번째 열열 # 2면에 있는 데이터를 모두 출력하시오 print(result[,,"Matrix2"]) # 2번 매트릭스의 3행2열에 있는 데이터를 추상해보고 확인하시오 print(result[3,2,2]) # 1번 매트릭스의 2행3열에 있는 데이터를 추상해보고 확인하시오 print(result[2,3,1]) # 2번 매트릭스의 1행3열에 있는 데이터를 추상해보고 확인하시오 print(result[1,3,2]) # 3번 매트릭스의 2행2열에 있는 데이터를 추상해보고 확인하시오 print(result[2,2,3]) # 3번 매트릭스의 2행3열에 있는 데이터를 추상해보고 확인하시오 print(result[2,3,3]) # 4번 매트릭스의 3행2열에 있는 데이터를 추상해보고 확인하시오 print(result[3,2,4]) # 5번 매트릭스의 3행1열에 있는 데이터를 추상해보고 확인하시오 print(result[3,1,5]) # data.frame # 사각형을 유지해야 하고 열내 동질적이고 열간 이질적인 데이터 # R은 데이터를 표현할 때 열중심으로 데이터를 저장하고 통계처리에 유리(10배정도 빠르다) # 데이터를 실시간으로 추가할 수 있음(javascript) x <- c(10,20,30,40) y <- c(6,7,8,9) # 변수의 개수가 특성의 개수를 의미 (data <- data.frame(width=x,height=y)) # 처음부터 열이름 붙여서 사용 str(data) data$width # 열이름으로 data[,1] # matrix방식으로 접근가능 head(data) # 상위 6개만 출력 data$delinetr = c(30,40,50,60) # 실시간으로 데이터 추가 data data <- edit(data) # 데이터 편집기 나옴 data L3 = LETTERS[1:3] # A,B,C d <- data.frame(cbind(x=1,y=1:10, fac=sample(L3,10,replace = T))) d d$fac (d$yhat <- d$fac) # 파생변수 str(d) head(d) # 범주형 변수 d$fac = factor(d$fac) rownames(d) = c("일","이","삼","사","오","육","칠","팔","구","십") d # factor형으로 자동으로 변경되는 현상이 버전업되면서 해결됨 x <- data.frame("SN"=1:2, "Age"=c(21,15),"Name"=c("에어컨","삼성SDS")) str(x) x$Name # 문자열로 참조 x["Name"] # 데이터 프레임으로 참조 x$Name = as.character(x$Name) ncol(x) nrow(x) length(x) # 2차원에서는 열의 개수의 의미 x["Age",1] <- 20; # Age행이 생성 x x[1,"Age"] <- 21; x x$SN <- NULL data() str(trees) head(trees) # 위에서부터 6개 tail(trees) # 뒤에서부터 6개 trees[2:3,] trees[trees$Height > 82,] # 조건식으로 boolean indexing trees[10:12,2] # numeric trees[10:12,2,drop=FALSE] # data.frame summary(trees) # 데이터를 열로 봄 # Min :최소값, 1st quantile(1사분위수), Median(중위수), Mean 산술평균 # 3rd quantile (3사분위수), Max 최대값 boxplot(trees) # 이상치 # 3사분위수 - 1사분위수 : inter Quantile range(IQR) # 3사분위수 + IQR*1.5 # 1사분위수 - IQR*1.5 # Data Cleaning(데이터 전처리) # -결측치 처리 NA : 제거, 평균으로 대치, 비슷한 데이터로 대치 # -이상치 처리 boxplot(trees) # 신뢰구간을 구한다음에 2사분위수를 넘어가는 것을 제외 # -정규화 min-max(0 ~ 1), normalize(-1 ~ 1), z점수 정규화 # -사이즈를 1로 # -평균을 0으로 하고 +,-값이 존재할 때 모델학습이 잘 된다. # -범주화: 정성적 데이터 ``` ``` pairs(trees) # 변수간의 상관도를 출력 :대칭행렬 # 상관계수 행렬은 대칭행렬이다. # Excel파일 열기 "C:\\Users\\401-26\\Documents" data <- read.csv("C:/Users/401-26/Documents/input.csv",fileEncoding="CP949",header=TRUE,) data colnames(data)[1] <- "id" print(data) print(is.data.frame(data)) #TRUE ncol(data) nrow(data) (sal <- max(data$salary)) (retval <- subset(data, salary == max(salary))) install.packages("stringr") library(stringr) (retval <- subset(data, str_trim(dept)=="전산")) (info <- subset(data, salary > 620 & str_trim(dept)=="전산")) (retval <- subset(data, as.Date(start_date)>as.Date("2014-01-01"))) write.csv(retval,"output.csv") # 저장할 때 (newdata <- read.csv("output.csv")) # 문제 # 다음을 데이터 프레임에 입력하시오 # 영어 등급 #퀴즈 67 "C" #중간 92 "A" #기말 89 "B" eng <- c(67,92,89) grade <- c("C","A","B") name <- c("퀴즈","중간","기말") data <- data.frame("영어"=eng,"등급"=grade, row.names=name) data # 수학점수를 50,100,80으로 입력하시오(열로 추가) data$"수학" <- c(50,100,80) # 보충이라는 이름으로 dataframe을 만들어 rbind시키시오(행으로 추가) da <- data.frame("영어"=c(90),"등급"=c("A"),"수학"=c(100), row.names="보충") # 영어점수 90, 등급 "A", 수학 100 data2 <- rbind(data,da) data2 # 열별 합계를 내서 데이터에 추가하시오(colSums()사용 못함,이질적인 데이터가 들어와있어서) hab <- colSums(as.matrix(cbind(data2$"영어",data2$"수학"))) da <- data.frame("영어"=hab[1],"등급"=0,"수학"=hab[2], row.names="합계") data2 <- rbind(data2, da) data2 # 국어점수를 추가하시오(90,80,100,77) kor <- c(90,80,100,77) (kor <- c(kor,sum(kor))) data2$"국어" <- kor write.table(data2,file="sungjuk.csv", sep=",",col.names = NA) (data3 <- read.table("sungjuk.csv", header=TRUE, sep=",")) ``` ``` # 정규분포:평균이 0이고 전체 확률이 1인 분포표 (x <- rnorm(400,mean=50,sd=10)) # random # 표준정규분포표 (확률을 확인하는 것이 가능) - 확률표를 단일작성 # z점수 = (데이터-평균)/표준편차 hist(x) x <- runif(400,min=0,max=100) # uniform 균등분포 hist(x) x <- floor(runif(400,min=0,max=100)) plot(x) # type= "b" => both(line+dot) plot(c(data3$"수학"[1:4]),main="그래프",col="red", ylab="수학점수",ylim=c(0,100),type="b") # 문제 # 수학의 평균, 중위수, 분산, 표준편차를 구하시오 mean(data3$"수학"[1:4]) median(data3$"수학"[1:4]) var(data3$"수학"[1:4]) sd(data3$"수학"[1:4]) # 다음과 같은 벡터를 컬럼으로 갖는 데이터프레임을 생성하시오 name <- c("대한이","민국이","만세야","희망이","다함께") age <- c(55,45,45,53,15) gender <- c(1,1,1,2,1) # 1:남자 2:여자 job <- c("연예인","주부","군인","직장인","학생") sat <- c(3,4,2,5,5) # 만족도 grade <- c("C","C","A","D","A") total <- c(44.4,28.5,43.5,NA,27.1) # 성별 변수를 이용하여 히스토그램 그리기(hist) # 만족도 변수를 이용하여 산점도를 그리시오(plot) # 총 구매금액에 대하여 평균과 분산 그리고 표준편차를 구하시오 user <- data.frame(name=name, age=age, gender=gender, job=job, sat=sat, grade=grade, total=total) user str(user) user$gender <- as.factor(user$gender) str(user) # histogram : 연속된 수치에 대한 카운트 hist(as.integer(user$gender),main="성별분포",xlab="성별(남자,여자)",xlim=c(1,2),ylab="인원",col="lightblue") plot(user$sat, main="만족도 산점도",col="red", ylab="만족도",ylim=c(0,5),type="b") mean(user$total, na.rm=T) var(user$total, na.rm=T) sd(user$total, na.rm=T) # 문제 # 다음 데이터를 data.frame으로 만드시오 name <- c("유관순","홍길동","이순신","신사임당") gender <- c("F","M","M","F") price <- c(50,65,45,75) client <- data.frame(name, gender, price) # 조건1 price가 65 이상인 고객은 "Best", 미만이면 "Normal" 문자열을 result변수에 넣고 client에 파생변수로 추가하시오 # 3항연산자:ifelse(조건, 참, 거짓) client$result = ifelse(client$price >= 65,"Best","Normal") client # 조건2 result의 빈도수를 구하시오 # 도수분포표(table) table(client$result) # 조건3 gender가 'M'이면 Male,'F'면 Female형식으로 client의 파생변수로 추가하시오 client$gender_origin <- ifelse(client$gender=="M",'Male','Female') client ``` ``` # list # 이질적인 데이터를 가질 수 있고 사각형을 유지하지 않아도 됨 x <- c(82,95,78) y <- c(84,92) z <- c(96,86,82,77) (core <- list("수학"=x, "영어"=y,"국어"=z)) # $로 접근하는 것은 동일 core typeof(core) length(core) str(core) # data.frame과 형식이 비슷함:사이즈가 다름름 # 리턴데이터로 참조해야되는 경우가 많음 (x <- list("name"="에어컨","age"=19,"speaks"=c("한국어","영어"))) x$name x[c(1:2)] x[-2] x[c(T,F,F)] x[c("age","speaks")] x["age"] typeof(x["age"]) #list x["age"] # 요소로 데이터를 출력 typeof(x[["age"]]) #double a <- list(c(1:5)) b <- list(6:10) lapply(c(a,b),max) # 함수 호출: 함수의 매개변수로 함수를 전달할 때는 ()안붙임 # factor형(범주형데이터: 무순서, 순서범주데이터) : 정성적 데이터 -> 범주형 gender <- c('M','F','M','F','M') gender plot(gender) # 문자는 그래프 불가 fgender <- as.factor(gender) fgender fgender[2] mode(fgender); class(fgender) plot(fgender) (x <- factor(c("single","married","married","single"),levels=c("single","married","divorced"))) str(x) levels(x) x[3] x[c(2,4)] x[-1] x[2] <- "divorced" #x[3] <- "widowed" # levels에 없기 때문에 입력불가 # 연속형 숫자의 범주화 (x = c(1:100)) cut(x,breaks=4) # 구간범주화 : 분석은 범주형 : 구간개수수 # 함수 충돌을 방지하는 패키지:: z <- stats::rnorm(10000) # rnorm을 패키지마다 다르게 정의할 수 있어서 중복되지 말라고 앞에 패키지이름(stats::) 붙임 z table(cut(z,breaks=-6:6)) # 일정한 구간으로 graphics::hist(z,breaks=-6:6) IQR(x) # Z점수를 구하는 방법: 평균을 빼주고 표준편차로 나누어주면 # 표준정규분포표에 접근할 수 있음 1sigma: 0.675 2sigma: 95% 3sigma: 99.7% # 1.5배 하는 이유 # 표준 정규분포에서 3sigma : 99.7% 일반적으로 정상데이터라고 본다 # IQR 1배수로 하면 너무 작은 데이터를 정상적으로 보기 때문에 이상치 많고 # IQR 2배수로 하면 너무 큰 범위를 정상으로 보기 때문에 이상치 없고 # 그래서 1.5배를 사용한다. ``` # 3/22 ``` (x <- 1:10) (x= rnorm(100)) # random normal 정규분포로부터 랜덤하게 선택 # vector -> matrix, dataframe, array, list var(x) # 기본이 벡터 sd(x) range(x) # 상한선 하한선 diff(range(x)) quantile(x,0.75) - quantile(x, 0.25) # 사분위를 구하는 함수 IQR(x) IQR(rnrom(10000)) summary(x) boxplot(x) summary(x) vec<-1:10 (vec<-sample(vec)) # 랜덤선택 summary(vec) table(vec) # 도수분포표 directory<-getwd() # "C:/Users/401-26/Documents" (현재 작업 경로 위치) setwd(directory) getwd() x <- c(1.5,2.5,-1.3,2.5) x round(mean(x)) # 반올림 ceiling(mean(x)) # 올림 floor(mean(x)) # 내림 #install.packages("mlbench") library(mlbench) data("BostonHousing", package="mlbench") # BostonHousing : 보스턴 주택가격 데이터셋 original <- BostonHousing # 예측데이터(선형회귀:종속변수+독립변수) ?BostonHousing # 506x14 (14개중 1개는 종속변수:medv, 13개는 독립변수) head(BostonHousing) set.seed(100) # 컴퓨터의 난수는 의사난수 -> 똑같은 결과가 나오기 위해서 # 3(평가):7(학습) # 컴퓨터의 모델은 난수에 의해 작동 # 참고로 선형회귀는 다 수치데이터이어야 함 # 범주형데이터는 2개 이상 사용못함 str(BostonHousing) sum(is.na(BostonHousing)) #0 (na없음) # 값에 의해서 전달 BostonHousing <- na.omit(BostonHousing) # 506개에서 40개 선택해라 BostonHousing[sample(1:nrow(BostonHousing),40),"rad"] <-NA # rad:BostonHousing 변수 중 하나 BostonHousing[sample(1:nrow(BostonHousing),40),"ptratio"] <-NA # patratio:BostonHousing 변수 중 하나(학생-교사 비율) BostonHousing$"rad" sum(is.na(BostonHousing)) #80 # 어떤 변수에 na가 있는지 모름 # apply 데이터에 함수를 적용 # apply, lapply(리스트인 경우), sapply(벡터나 행렬), tapply(그룹핑한 다음 함수 적용) # 함수적 프로그래밍(반복문x - 자동으로 모든 데이터에 대하여 적용) # 열별로 적용 # 데이터프레임은 열우선 sapply(BostonHousing, function(x) sum(is.na(x))) # x는 열벡터로 입력 # 종속변수(medv) ~ 독립변수(ptratio, rad) # 506 - 80 # 426 lm(medv ~ ptratio + rad, data = BostonHousing, na.action = na.omit) #install.packages("Hmisc") # 함수제공하는 패키지 다운로드 library(Hmisc) # Hmisc 패키지 불러오는 코드드 impute(BostonHousing$ptratio, mean) # NA가 있으면 평균으로 대체해라 impute(BostonHousing$ptratio, median) # NA가 있으면 중위수로 대체해라 impute(BostonHousing$ptratio, 20) # NA가 있으면 20으로 대체해라 # if문 BostonHousing$ptratio[is.na(BostonHousing$ptratio)] <- mean(BostonHousing$ptratio, na.rm=T) # na.rm은 함수에서 계산이나 요약할 때 누락된 값 (NA)을 처리하는 방법. TRUE로 설정하면, 함수는 NA 값을 무시하고 계산 또는 요약을 수행 ``` ``` # byrow 디폴트 값은 F (m1 <- matrix(c <-(1:10),nrow=5,ncol=6)) # 5x6 > 30개 a_m1 <- apply(m1,2,sum) # 1이면 행으로 계산, 2이면 열방향으로 계산 a_m1 vec <- c(1:10) vec apply(vec,2,sum) # 전제조건: 최소 dim 2차원 sapply(vec,sum) #1 2 3 4 5 6 7 8 9 10 # 행렬상관없이 열별로 계산산 avg <- function(x){ (min(x)+max(x))/2 } cars head(cars) # 속도에 따른 제동거리 dt <- cars sapply(dt,avg)# 열별로 적용 data(iris) iris head(iris) dim(iris) # 150x5 (데이터 150개, 변수 5개(종속1,독립4)) => Species가 종속변수 # iris 종속변수 Species는 setosa versivolor virginica 3가지 그룹으로 되어있다. tapply(iris$Sepal.Width, iris$Species, median) # 그룹별 Width의 중위수 tapply(iris$Sepal.Width, iris$Species, mean) # 그룹별 width의 평균 tapply(iris$Sepal.Width, iris$Species, sum) # 그룹별 width의 합 ?mtcars # 32x11 연비테스트:예측 head(mtcars) tail(mtcars) summary(mtcars) str(mtcars) (exdata <- mtcars[1:10,1:3]) class(exdata) apply(exdata,2,mean) # 열방향으로 변수마다 평균 (x=apply(exdata,1,mean)) boxplot(mtcars) # 열별 ``` ``` list.o.nums <- list(runif(100), rnorm(100), rpois(100, lambda=1)) lapply(list.o.nums, sd) # 결과 list sapply(list.o.nums, sd) # 결과 벡터 apply(list.o.nums, 2, sd) # 이빨빠진 2차원은 2차원으로 안본다. # 1부터 10의 데이터에 대하여 2승 값을 구해서 출력하시오(apply) num <- 1:10 lapply(num, function(x){x^2}) sapply(num, function(x){x^2}) # 문제 # iris의 4개의 독립변수에 대하여(열별) 평균을 계산하시오 sapply(iris[,1:4],mean) # 정규분포에서 랜덤하게 데이터를 20개를 선택한 다음 5x4행렬을 만들고 # 행합계와 열합계를 구하시오 data <- rnorm(20) (mat <- matrix(data, nrow=5, ncol=4)) apply(mat,1, sum) apply(mat,2, sum) colSums(mat) colMeans(mat) # 1사분위수와 2사분위수 3사분위수를 수하시오 quantile(mat, probs=0.25) # 1사분위수 quantile(mat, probs=0.5) # 2사분위수 quantile(mat, probs=0.75) # 3사분위수 apply(mat,2,quantile, probs=c(0.25,0.5,0.75)) (data <- c(1:10, rnorm(10,2), runif(10))) groups <-gl(3,10) tapply(data, groups, mean) ``` ``` # factor # 구간범주화 #install.packages("ggplot2") library(ggplot2) library(MASS) ?Cars93 str(Cars93) hist(Cars93$MPG.highway) disc_1 <- Cars93[,c("Model","MPG.highway")] # highway:고속도로연비(수치형이고 범주형데이터 아님):6가지 구간화해야함 head(disc_1) range(disc_1["MPG.highway"]) dim(Cars93) # 93x27 sum(disc_1$MPG.highway >= 25 & disc_1$MPG.highway < 30) disc_1 <- within(disc_1,{ # 연비 구간화 시킴 MPG.highway_cd = character(0) # 변수추가하고 문자타입으로 초기화 MPG.highway_cd[ MPG.highway >= 20 & MPG.highway < 25 ] = "20~25" # 구간이름 = 그룹이름 MPG.highway_cd[ MPG.highway >= 25 & MPG.highway < 30 ] = "25~30" #41 MPG.highway_cd[ MPG.highway >= 30 & MPG.highway < 35 ] = "30~35" #27 MPG.highway_cd[ MPG.highway >= 35 & MPG.highway < 40 ] = "35~40" MPG.highway_cd[ MPG.highway >= 40 & MPG.highway < 45 ] = "40~45" MPG.highway_cd[ MPG.highway >= 45 & MPG.highway < 50 ] = "45~50" MPG.highway_cd = factor(MPG.highway_cd, level = c("20~25","25~30","30~35","35~40","40~45","45~50")) }) attributes(disc_1$MPG.highway_cd) # 문제 # 범주화된 구간값에 따라 MPG.highway의 합계, 평균, 표준편차를 구하시오. tapply(disc_1$MPG.highway, disc_1$MPG.highway_cd, sum) # 합계 tapply(disc_1$MPG.highway, disc_1$MPG.highway_cd, mean) # 평균 tapply(disc_1$MPG.highway, disc_1$MPG.highway_cd, sd) # 표준편차 # ggplot은 layer로 적용 ggplot(disc_1, aes(x=MPG.highway_cd, fill=MPG.highway_cd))+geom_dotplot(binwidth=0.3) ggplot(disc_1, aes(x=MPG.highway_cd, fill=MPG.highway_cd))+geom_bar() # 막대그래프 # 사분위수 구간 quantile(disc_1$MPG.highway, 0.25) quantile(disc_1$MPG.highway, 0.5) quantile(disc_1$MPG.highway, 0.75) # dummy 변수화 cust_id <- c("c01","c02","c03","c04","c05","c06","c07") age <- c(25,45,31,30,49,53,27) cust_profile <- data.frame(cust_id, age, stringsAsFactors = F) cust_profile # 범주 수만큼 변수가 증가하고 데이터는 0과 1로 표현된다. cust_profile <- transform(cust_profile, age_20 = ifelse(age<=20 & age<30,1,0), age_30 = ifelse(age<=30 & age<40,1,0), age_40 = ifelse(age<=40 & age<50,1,0), age_50 = ifelse(age<=50 & age<60,1,0)) cust_profile ``` ``` # 정규화(min-max) head(iris) min_max_norm <- function(x){ # 열별 (x - min(x)) / (max(x) - min(x)) } lapply(iris[1:4], min_max_norm) # 리스트로 출력 iris_norm <- as.data.frame(lapply(iris[1:4], min_max_norm)) # 열별적용 head(iris_norm) iris_norm$Species <- iris$Species head(iris_norm) # 종속변수에 동일한 영향을 미치기 위해서 정규화 함 z_normal <- function(x){ (x-mean(x))/sd(x) # 표준편차:sqrt(var) # 분산은 평균과의 차의 제곱의 합 } iris_z_norm <- as.data.frame(lapply(iris[1:4],z_normal)) iris_z_norm class(iris) (iris_standardize <- as.data.frame(scale(iris[1:4]))) # 열별로 head(iris) re = boxplot(iris) re # 리스트 fivenum(iris[,1], na.rm = TRUE) # minimum, lower, median, upper, maiznum #install.packages("dplyr") library(dplyr) outdata <- iris[1:4] out <- outdata for(i in 1:(ncol(outdata)-1)){ # 사분위수:열수 uppercut = fivenum(outdata[,i])[4]+1.5*IQR(outdata[,i]) lowercut = fivenum(outdata[,i])[2]-1.5*IQR(outdata[,i]) # ,로 분리된 조건 out <- filter(out, out[,i] <= uppercut, out[,i] >= lowercut) } str(out) for(i in 1:(ncol(outdata)-1)){ uppercut = quantile(outdata[,i],0.75)+1.5*IQR(outdata[,i]) lowercut = quantile(outdata[,i],0.25)-1.5*IQR(outdata[,i]) out <- filter(out, out[,i] <= uppercut, out[,i] >= lowercut) } str(out) ``` ``` # 위치 데이터 state_table <- data.frame(key=c("SE","DJ","DG","SH","QD"), name=c("서울","대전","대구","상해","칭따오"), country=c("한국","한국","한국","중국","중국")) state_table #년도 month_table <- data.frame(key=1:12, desc=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Agu","Sep","Oct","Nov","Dec"), quarter=c("Q1","Q1","Q1","Q2","Q2","Q2","Q3","Q3","Q3","Q4","Q4","Q4")) month_table # 상품 테이블 prod_table <- data.frame(key=c("Printer","Tablet","Laptop"), price=c(225,570,1120)) prod_table names(prod_table) row.names(prod_table) (prod_table1 <- data.frame(Printer = 255, Tablet = 570, Laptop=1120)) (prod_table1 <- t(prod_table1)) gen_sales <- function(no_of_recs){ # 가짜 데이터 생성 loc <- sample(state_table$key, no_of_recs, replace=T, prob = c(2,2,1,1,1)) time_month <- sample(month_table$key, no_of_recs, replace = T) time_year <- sample(c(2012,2013),no_of_recs,replace=T) prod <- sample(prod_table$key, no_of_recs, replace=T, prob=c(1,3,2)) unit <- sample(c(1,2),no_of_recs, replace=T, prob=c(10,3)) price <- unit*prod_table1[prod,1] sales <- data.frame(month=time_month, year=time_year, loc=loc, prod=prod, unit=unit, price=price) sales <- sales[order(sales$year,sales$month),] row.names(sales) <- NULL return(sales) } sales_fact <- gen_sales(500) dim(sales_fact) str(sales_fact) # product별로 월,년으로 location별로 합계를 낼 때 # 범주형 3개, 월은 범주형 3x12x2x5 (revenue_cube <- tapply(sales_fact$price, sales_fact[,c("prod","month","year","loc")], # 데이터를 4차원으로 확인 FUN=function(x){return(sum(x))})) # 문제 # 제품별 판매액을 계산하시오 tapply(sales_fact$price,sales_fact[,c("prod")], FUN = function(x){return(sum(x))}) # 지역별, 제품별 판매액을 계산하시오 tapply(sales_fact$price, sales_fact[,c("loc","prod")], FUN = function(x){return(sum(x))}) # 2012년도 1월 태블릿에 대한 판매현황을 출력하시오 (revenue_cube["Tablet","1","2012",]) revenue_cube apply(revenue_cube, c("year","prod"),FUN=function(x){return(sum(x, na.rm=TRUE))}) dim(revenue_cube) #3 12 2 5 # 문제 # 제품별로 년도별 월별 매출액을 확인하시오 apply(revenue_cube, c("prod","month","year"), FUN=function(x){return(sum(x, na.rm=T))}) # 연도별 월별로 작업하시오 apply(revenue_cube, c("year","month"), FUN=function(x){return(sum(x, na.rm=T))}) # 매장별 제품별로 월별 매출액을 확인하시오 apply(revenue_cube, c("loc","month","prod"), FUN=function(x){return(sum(x, na.rm=T))}) # 3 * 12 * 2 * 5 경우의 수가 360가지를 입체적으로 표현(차원확대(그룹의 개수),차원축소(그룹의 개수를 줄임)) ``` ``` # 로딩된 패키지 확인 search() mean(iris$Sepal.Length) # with vs within(대입가능) with(iris,{ mean(Sepal.Length) }) head(iris) iris[1,1]=NA head(iris) # vector나 metrix # 종별로 중위수를 구함 median.per.species <- sapply(split(iris$Sepal.Length, iris$Species), median, na.rm=T) median.per.species iris <- within(iris,{ # 중위수로 대체하거나 아니면 자기자신의 값으로 Sepal.Length <- ifelse(is.na(Sepal.Length), median.per.species[Species],Sepal.Length) # NA를 중위수로 대체 }) head(iris) mtcars str(mtcars) # 실린더가 4 또는 5개인 자동차들 중 연비, 실리더수, 변속기 종류를 출력하시오 mtcars[which(mtcars$cyl == c(4,6)), c("mpg","cyl","am")] # 행은 조건이 오고, # 연비가 20보다 큰 데이터 중 연비, 실린더수, 변속기 종류를 출력하시오 mtcars[which(mtcars$mpg >= 20), c("mpg","cyl","am")] # 평균 마력 이상인 자동차중 연비, 실린더, 변속기를 출력하시오 mtcars[mtcars$hp >= mean(mtcars$hp), c("mpg","cyl","am")] # mtcars의 연비에 대하여 min/max 정규화해 보시오 ( mtcars$mpg - min(mtcars$mpg) )/( max(mtcars$mpg)-min(mtcars$mpg) ) # 위를 with문을 이용해서 수정해보시오 with(mtcars, (mpg-min(mpg)) /(max(mpg)-min(mpg))) # merge # merge가 rbind,cbind와 다른점: 양쪽에 키가 있어야 함 x <- data.frame(name=c("a","b","c"),math=c(1,2,3)) y <- data.frame(name=c("a","b","d"), english=c(4,5,6)) merge(x,y) # 키가 일치하는 경우에 출력 merge(x,y,all=TRUE) # 모든 경우, 없으면 NA로 출력 merge(x,y,all.x=T) # left join merge(x,y,all.y=T) # right join merge(x,y,all=T) # full join # search search() attach(airquality) # 데이터를 패키지처럼 로딩 transform(Ozone, logOzone=log(Ozone)) search() detach(airquality) search() ls() # subset subset(iris,subset=(Species=="setosa")) # 조건 subset(iris,subset=(Species=="setosa" & Sepal.Length>=5.0)) subset(iris,select=c(Sepal.Length, Species)) # 열선택 subset(iris,select=c(Sepal.Length, Species)) # 제외하고 subset(iris,subset=(Species=="setosa"), select = c(Sepal.Length, Species)) # 실린더가 4~6개인 자동차들 중 연비, 실린더수, 변속기 종류를 출력하시오 subset(mtcars,select=c(mpg,cyl,am),subset= cyl>=4 & cyl<=6) subset(mtcars,select=c(mpg,cyl,am),subset=(cyl %in% c(4,5,6))) # 실린더가 4또는6개인 자동차들 중(연비mpg, 실린더cyl, 변속기am, 마력hp) 종류를 출력하시오 subset(mtcars, select=c(mpg,cyl,am), subset=(cyl%in%c(4,6))) # 문제 # 연비가 20보다 큰 자동차들 중 연비, 실린더수, 변속기 종류를 출력 subset(mtcars, select=c(mpg,cyl,am), subset= mpg>20) # 변속기가 자동이고 실린더가 4,6개인 자동차들의 연비 평균은? colMeans(subset(mtcars, select=c(mpg), subset=(cyl %in% c(4,6))&(am %in% 0))) ``` # 3/23 ``` #install.packages("tidyverse") library(tidyverse) exp(1) 1 %>% exp head(mtcars,4) mtcars %>% head(4) iris # 150x5 iris %>% subset(Sepal.Length > 5,c(Species, Sepal.Length, Sepal.Width)) %>% aggregate(. ~ Species, ., mean) # . 이 데이터가 들어갈 위치 str(mtcars) # 32x11 # 절차적 프로그래밍이 가능 mpg_inds = order(mtcars$mpg, decreasing=T) # 정렬된 데이터의 인덱스 data <- mtcars[mpg_inds,] head(data, 4) mtcars %>% arrange(desc(mpg)) %>% # 정렬 head(10) # 문제 # mtcars(gear, hp) 순으로 정렬해보시오 mtcars %>% arrange(desc(gear),desc(hp)) %>% # 정렬 head(10) mtcars %>% select(cyl, disp, hp) %>% # 변수선택(속성선택, 열선택) head(2) mtcars %>% select(starts_with("d")) %>% # starts_with 문자열(선두 문자열을 체크), ends_with head(2) # 문제 # iris데이터 중에 Sepal로 시작하는 데이터만 선택하시오 iris %>% select(starts_with(("Sepal"))) # iris데이터 중에 Petal이나 Sepal로 시작하는 데이터를 선택하시오 iris %>% select(starts_with(c("Petal","Sepal"))) iris %>% select(matches(("[pt]al"))) # 문자열 정규표현식 함수 # pal이나 tal을 의미 mtcars %>% select(contains('ar')) %>% head(2) # 파생변수 mtcars = mtcars %>% mutate(hp_wt = hp/wt, mpg_wt = mpg/wt) # 파생변수로 추가 mtcars # 여러 열에 대하여 함수를 적용 mtcars = mtcars %>% mutate_at(c("hp_wt","mpg_wt"),log) #install.packages("EDAWR") -- version을 타서 안됨 #library(EDAWR) #EDAWR::cases ?diamonds head(diamonds) (df.diamonds_ideal <- filter(diamonds, cut=="Ideal")) head(df.diamonds_idea) (df.diamonds_ideal <- dplyr::select(df.diamonds_ideal, carat, cut, color, price, clarity)) head(df.diamonds_ideal) (df.diamonds_ideal <- mutate(df.diamonds_ideal, price_per_carat = price/carat)) head(df.diamonds_ideal) # 위의 명령을 pipeline 명령으로 변경해보시오 df.diamonds <- diamonds %>% filter(cut=="Ideal") %>% dplyr::select(carat, cut, color, price, clarity) %>% mutate(price_per_carat = price/carat) diamonds %>% filter(cut=="Ideal") %>% ggplot(aes(x=color,y=price))+ geom_boxplot() diamonds %>% filter(cut=="Ideal") %>% ggplot(aes(price))+ geom_histogram() + facet_wrap(~color) # 컬러로 범주화한 다음 시각화 # tapply : stat="identity" y변수에 대하여 합계를 계산하라 iris %>% group_by(Species) %>% summarize(meanSepLength=mean(Sepal.Length)) %>% ggplot(aes(Species, meanSepLength)) + geom_bar(stat="identity") # 문제 # iris 데이터에 대하여 Sepal.Length에 10000을 곱해서 x변수로 추가하고 Species와 x변수를 선택해서 result변수에 대입하시오 (result <- iris %>% mutate(x = Sepal.Length*10000) %>% dplyr::select(Species, x) %>% map_dfr(., ~format(.x, big.mark=','))) head(mtcars,4) (by_cyl <- mtcars %>% group_by(cyl)) by_cyl %>% summarise( # 그룹표시가 있는 데이터에 대하여 그룹별로 요약 disp = mean(disp), hp=mean(hp) ) # 문제 # mtcars데이터에 대하여 vs, am으로 그룹화하고 disp의 평균을 구해보시오 by_vs_am <- mtcars %>% group_by(vs,am) %>% dplyr::summarise(disp=mean(disp)) # 명령어 이름이 동일해서 충돌이 일어나니까 dplyr 패키지 써주는거 result <- mtcars %>% group_by(vs,am) %>% # ungroup() %>% -> 그룹 푸는거 dplyr::summarise(n=n()) # 개수를 세는 경우 result %>% dplyr::summarise(n=sum(n)) (by_cyl %>% group_by(vs, am) %>% group_vars()) #install.packages("hflights") library(hflights) ?hflights # 227496x21 # SQLife 경량 데이터베이스(웹) head(hflights) str(hflights) dim(hflights) library(data.table) # 고속으로 데이터처리를 위해서 hflights_df <- as_tibble(hflights) # "tbl" # df, dt : data.table class(hflights) class(hflights_df) # 정렬 # %>% 파이프 없이 단독으로 사용이 될 때는 데이터를 입력 arrange(hflights_df, Month, DayofMonth, desc(AirTime)) summarise(hflights_df, delay=mean(ArrDelay, na.rm=T)) # 목적지(Dest) 별로 그룹핑한 다음 비행기(TailNum) 별로 요약한 것을 planes라는 변수에, 횟수(n())를 flights라는 변수로 요약해보시오 hflights_df %>% group_by(Dest) %>% summarise(planes = n_distinct(TailNum), flights = n()) # 년,월,일로 그룹핑해서 년,월,일 별 비행횟수를 확인해보시오 hflights_df %>% group_by(Year, Month, DayofMonth) %>% dplyr::summarise(number_flights=n()) %>% dplyr::summarise(number_flights=n()) # 년,월,일로 그룹핑한 다음 년,월,일과 ArrDealy, DepDelay만 선택하고(select) ArrDelay,DepDelay 평균을 구하시오(summarise) # ArrDelay, DepDelay의 합계가 30이 넘는 경우만 출력하시오(filter) hflights_df %>% group_by(Year, Month, DayofMonth) %>% dplyr::select(Year:DayofMonth, ArrDelay, DepDelay) %>% # 집계함수 결과 summarise(arr = mean(ArrDelay, na.rm=T), dep = mean(DepDelay, na.rm=T)) %>% filter(arr + dep>30) # dplyr명령은 column해주는 명령이 있고, row해주주는 명령, dataframe하는 명령, grouping하는 명령이 있다. # 데이터는 몇년의 데이터인가 (년도데이터 중복을 제거하고 카운트) distinct, n_distinct n_distinct(hflights_df$Year) # 매개변수가 vector # 1월1일의 데이터를 추출하시오(filter) filter(hflights_df, Month=="1",DayofMonth=="1") filter(hflights_df, Month==1 & DayofMonth==1) # 1월 혹은 2월 데이터를 추출하시오(filter) filter(hflights_df, Month==1 | Month==2) # arrange를 이용해서 년,월,도착시간으로 정렬하시오(arrange) arrange(hflights_df, Year, Month, ArrTime) # Month를 기준으로 내림차순하시오(arrange) arrange(hflights_df, desc(Month)) # select를 이용해서 년,월,일을 제외하고 선택하시오 select(hflights_df, -(Year:DayofMonth)) # 비행편수를 구하기 위해서 항공기별(TailNum)로 그룹화하시요 # 그룹화 한 결과를 항공기별 횟수, Distance의 평균, Arrdelay의 도착지연 시간 평균을 요약하시요 # 위의 결과에서 횟수가 20 초과 거리의 평균이 2000 이하인 데이터만 출력하시요 # 위의 요약결과의 결측치 처리 : 위의 결과값이 dist(distance의 평균)와 delay(ArrDelay의 평균)이 # 모두 1 이상인 경우만 필터링하시요 hflights_df %>% group_by(TailNum) %>% summarise(count=n(),dist=mean(Distance, na.rm=T), delay=mean(ArrDelay, na.rm=T)) %>% filter(count>20, dist<=2000) %>% filter(dist>1 & delay>1) # 다음 데이터를 데이터프레임으로 생성하시오 # empno pay bonus # 101 250 0.10 # 102 180 0.10 # 103 200 0.12 # 104 300 0.15 # 105 1000 0.00 empno <- 101:105 pay <- c(250,180,200,300,1000) bonus <- c(0.10,0.10,0.12,0.15,0.00) (data <- data.frame(empno=empno,pay=pay,bonus=bonus)) # total이라는 이름으로 보너스가 고려된 급여를 계산해서 추가하시오 (data <- data %>% mutate(total=(pay*(1+bonus)))) # total의 급여가 300이상인 사원번호와 총급여를 출력하시오 data %>% filter(total >= 300) %>% select(empno, total) # 다음 deptno변수의 값을 추가하시오 deptno <- c(1,2,1,2,2) data$deptno = deptno # 부서별 평균급여를 출력하시오 data %>% dplyr::select(deptno, total) %>% group_by(deptno) %>% summarise(mean_total=mean(total, na.rm=T), sum_total= sum(total, na.rm=T), count=n()) ``` ``` library(dplyr) library(dbplyr) library(DBI) library(rstudioapi) install.packages('RMySQL') library(RMySQL) conn=dbConnect(MySQL(),user='root',password='acorn1234',dbname='Rtest',host='192.168.41.183') dbListTables(conn) mtcars dbWriteTable(conn,"mtcars",mtcars) dbListTables(conn) dbListFields(conn,"mtcars") db_mtcars <- dbReadTable(conn,"mtcars") dim(db_mtcars) res <- dbSendQuery(conn,"SELECT * FROM mtcars WHERE cyl=4") res # MySQLResult (java에서는 ResultSet -> Class로 변환) query_result <- dbFetch(res) dim(query_result) query_result dbClearResult(res) dbDisconnect(conn) # 문제 # daejeon 데이터베이스에 접속한 다음 테이블 리스트를 확인해보시고 emp의 데이터를 "select * from emp"로 쿼리를 실행한 다음 결과를 출력하시오 conn=dbConnect(MySQL(),user='root',password='acorn1234',dbname='daejeon',host='192.168.41.183') dbListTables(conn) res <- dbSendQuery(conn,"SELECT * FROM emp") query_result <- dbFetch(res) dim(query_result) query_result sum(is.na(query_result)) sapply(query_result, function(x) sum(is.na(x))) #drop_na(query_result) fill(query_result, COMM, .direction='up') replace_na(query_result,list(COMM=0)) dbClearResult(res) dbDisconnect(conn) ``` ``` # install.packages("nycflights13") library(nycflights13) nycflights13::flights names(flights) # conn = dbConnect(RMySQL::MySQL(), user='root', password='acorn1234', dbname='Rtest', # host = '192.168.41.183') conn <- DBI::dbConnect(RMySQL::MySQL(), dbname = "rtest", host = "192.168.41.183", user = "root", password = rstudioapi::askForPassword("Database password")) dbListTables(conn) dbWriteTable(conn, "flights", nycflights13::flights) dbWriteTable(conn, "airlines", nycflights13::airlines) dbWriteTable(conn, "planes", nycflights13::planes) dbWriteTable(conn, "airports", nycflights13::airports) dbWriteTable(conn, "weather", nycflights13::weather) # copy_to(conn, nycflights13::flights, "flights", # temporary = FALSE, # indexes = list( # "year", "month", "day", # "carrier", # "tailnum", # "dest" # ) # ) # # nycflights13::airlines # copy_to(conn, nycflights13::airlines, "airlines", # temporary = FALSE, # indexes = list("carrier"), overwrite=T # ) # # copy_to(conn, iris, overwrite=TRUE) (mtcars<-tbl(conn, "mtcars")) (flights_db <- tbl(conn, "flights")) dim(flights_db) flights_db %>% select(year:day, dep_delay, arr_delay) db_test <- flights_db %>% filter(dep_delay > 240) db_test %>% show_query() gb_test <- flights_db %>% group_by(dest) %>% summarise(delay = AVG(dep_time)) gb_test %>% show_query() dbDisconnect(conn) ``` 과제 ``` library(dplyr) library(dbplyr) library(DBI) library(rstudioapi) #install.packages('RMySQL') library(RMySQL) # 종합적으로 연결하는 # 엑셀에서 성적 데이터를 입력(국어, 영어, 수학, 과목과 3인분의 데이터)합니다. (함계,평균 제외) - csv로 저장 # R에서 데이터를 로딩하여 합계 평균 및 학점을 구합니다. # 계산된 결과 data.frame을 sungjukresult데이터 베이스의 student테이블로 저장합니다.(RMySQL사용) # 저장된 데이터를 dbplytr을 이용하여 데이터베이스로부터 로딩합니다. # 국어 영어를 선택합니다. # 국어 점수가 80점 이상인 데이터만 필터링 합니다. # 처리 결과를 sungjukresult 데이터 베이스의 calcsungjuk 테이블에 저장합니다. "C:\\Users\\401-26\\Documents\\R" data <- read.csv("C:/Users/401-26/Documents/R/sungjuk.csv", fileEncoding = "CP949", header=FALSE,) data colnames(data)[c(1:4)] <- c("이름","국어","영어","수학") data$"합계" <- data$"국어"+data$"수학"+data$"영어" data$"평균" <- round(data$"합계"/3,1) data$"학점" <- ifelse(data$"평균" >= 90,"A",ifelse(data$"평균" >= 80,"B",ifelse(data$"평균" >= 70,"C",ifelse(data$"평균" >= 60,"D","F")))) conn <- DBI::dbConnect(RMySQL::MySQL(), dbname="sungjukresult", host="192.168.41.183", user="root", password=rstudioapi::askForPassword("Database password")) dbListTables(conn) dbWriteTable(conn,"student",data) res <- dbSendQuery(conn, "select 국어,영어 from student") query_result <- dbFetch(res) query_result data %>% select(국어,영어) res2 <- dbSendQuery(conn, "select * from student where 국어>=80") query_result2 <- dbFetch(res2) query_result2 data %>% filter(국어>=80) da <- data %>% select(국어,영어) %>% filter(국어>=80) da dbListTables(conn) dbWriteTable(conn,"calcsungjuk",da) ``` # 3/24 ``` # graphic data(cars) str(cars) # 50x2행렬:data.frame head(cars) plot(cars$speed, type="l") # 2차원(x,y) 데카르트 좌표: 90도 직교 plot(cars$dist, type="l") plot(cars,type="l") # line plot(cars,type="p") # point plot(cars,type="b") # both # point shapes is pch plot(cars,type="b", pch=20) # plotting character plot(cars,type="b", pch=22) plot(cars,type="b", pch="+") plot(cars,type="b", pch=2, cex=2) plot(cars,type="l", pch="+", cex=2, col='red', lty="dotdash") # lty(line type)0~6번까지 plot(cars, main="그래픽", sub="옵션확인을 위해서",col='red',xlab="speed",ylab="dist",ylim=c(0,40),xlim=c(1,10),type="l",asp=1) # aspect # 문제 # 라인타입의 종류에 대하여 확인해보시오 # 0:blank, 1:solid(실선), 2:dashed(대시선), 3:dotted(점선), 4:dotdash(점-대시-점), 5:longdash(긴 대시선), 6:twodash(두 줄 대시선) # plotting character: 0~25 :26개 # cex: character expansion plot(cars, cex=0.5) identify(cars$speed, cars$dist) plot(cars, cex=0.5) text(cars$speed, cars$dist, pos=1) # 앞면하고 뒷면이 있는 코인에 대하여 5회 던졌을 때 앞면이 나온 횟수 coin <- c(2,2,0,1,1,1,2,3,1,2,2,3,3,4,1,3,2,3,3,1,2,2,1,2,2) (frequency <- table(coin)) # 도수분포표 (relative <- frequency/length(coin)) # 상대 도수분포표 (cumulative <- cumsum(relative)) # 누적 상대 도수분포표 # par : margins, fonts, colors 지정 opar <- par(mfrow=c(1,4)) # 전의 상태값을 리턴턴 plot(frequency, xlab="값",ylab="도수",type="b", col='red', main="도수", sub="순수도수",frame.plot = F,panel.first=T) plot(1:5, frequency, xlab="값",ylab="도수",type="b", col='red',frame.plot = F) plot(round(relative,2),type="b",pch=23,col='red') plot(round(cumulative,2),type="b",col='red',axes=F) # 축출력 안됨됨 par(opar) # 상태복원 opar <- par(mfrow=c(1,1)) plot(round(cumulative,2),type='b',col='red',axes=F) plot(c(1:10), xlab=expression(xLab ~ x^2 ~ m^~2 ~ mu ~ sqrt(2) ~ pi), ylab= expression(yLab ~ y^2 ~ m^~2),main="Plot 1") max.temp <- c(22,27,26,24,23,26,28) names(max.temp) <- c("Sun","Mon","Tue","Wen","Thu","Fri","Sat") par(mfrow=c(1,2)) barplot(max.temp, main="Barplot") pie(max.temp, main="Piechart", radius=1) par(mfrow=c(1,1)) ?faithful str(faithful) with(faithful,{ # bw: bandwidth plot(density(eruptions, bw=0.15)) # 밀도 -> 보간법을 이용해서 연속적으로 표현해라 rug(eruptions) # jitter :겹치는 데이터를 더 잘 구분할 수 있도록 무작위로 데이터를 조정하는 방법 rug(jitter(eruptions, amount = 0.01), side=3, col="light blue") }) library(MASS) Boston$tax plot(density(Boston$tax, bw=5)) rug(Boston$tax + rnorm(length(Boston$tax), sd=5), col=2, lwd=3.5) (x=(0:100)*pi/10) plot(x, sin(x), main="cos,sin 출력", ylab="", type="l", col='purple') plot(x, sin(x)+exp(x/10), main="cos,sin 출력", ylab="", type="l", col='blue') amp.1 <- 2 # 진폭을 결정 amp.2 <- 2 amp.3 <- 5 amp.4 <- 5 wav.1 <- 1 # 주기를 결정정 wav.2 <- 2 wav.3 <- 3 wav.4 <- 7 signal.1 <- amp.1*sin(wav.1*x) # 라디오파는 장파, 저주파 signal.2 <- amp.2*sin(wav.2*x) signal.3 <- amp.3*sin(wav.3*x) signal.4 <- amp.4*sin(wav.4*x) # 단파, 고주파 par(mfrow = c(1,4)) plot(x, signal.1, type="l",ylim=c(-5,5));abline(h=0,v=0,lty=3) plot(x, signal.2, type="l",ylim=c(-5,5));abline(h=0,lty=3) plot(x, signal.3, type="l",ylim=c(-5,5));abline(h=0,lty=3) plot(x, signal.4, type="l",ylim=c(-5,5));abline(h=0,lty=3) par(mfrow = c(1,1)) par(mfrow=c(1,2)) str(cars) plot(cars$dist, type="o", cex=0.5, xlab="speed",ylab="dist") tapply(cars$dist, cars$speed, mean) plot(tapply(cars$dist, cars$speed, mean),type="o",cex=0.5, xlab="speed",ylab="dist") par(mfrow=c(1,1)) ``` ``` # 통계학: 추정과 검정 # 우리나라 전국민의 키 평균 구하려면 시간, 비용 걸리니까 표본을 구해서 # 키의 평균 -> 모집단 키의 평균으로 추정 # 추정:점추정,구간추정 (구간추정이 유리) # 구간은 표준오차로 구한다. -> 표준편차를 sqrt(표본의 개수)해서 표준오차를 구함 # 표본 많아지면 분산도 커짐 # 1.96 : 표준편차의 2배수 -> 95%신뢰구간 # 평균 - 2*표준오차 ~ 평균 + 2*표준오차 : 신뢰구간 # 선형회귀 -> 점추정과 구간추정 -> 신뢰구간 # 데이터를 기울기와 절편을 구해서 추정 # lm:linear model (m <- lm(dist ~ speed, data=cars)) # 기울기, 절편 plot(cars) abline(m) # 직선을 작도할 때 (p <- predict(m, interval="confidence")) # 신뢰구간을 가지고 예측 head(p) x <- c(cars$speed, tail(cars$speed,1),rev(cars$speed),cars$speed[1]) y <- c(p[,"lwr"],tail(p[,"upr"],1), rev(p[,"upr"]),p[,"lwr"][1]) # 폐다각형 polygon(x,y,col=rgb(.7,.7,.7,.5)) # 비선형회귀 opar <- par(mfrow=c(1,1)) plot(cars) lines(lowess(cars)) # 딥러닝은 비선형 str(women) # 키와 체중간의 관계 head(women) # 문제 # 위의 데이터에 대하여 선형회귀를 하고 회귀선을 출력해보시오 summary(women) (res=lm(height ~ weight, data=women)) # 독립변수가 x축, 종속변수가 y축 plot(women$weight, women$height) plot(women) abline(res) # 관측값 뿐만 아니라 다양한 파생변수와와 관계고려 str(trees) summary(trees) Height2 <- trees$Height^2 # 파생변수 trees2 <- cbind(trees, Height2) attach(trees2) (test2 <- lm(Girth ~ Height + Height2)) plot(Girth~Height) fitted(test2) lines(sort(Height),fitted(test2)[order(Height)],col='red',type='l') detach(trees2) par(mfrow=c(1,1)) honeymoon <- c(12,14,8,9) names(honeymoon) <- c("하와이","동남아","유럽","기타") pie(honeymoon, col=rainbow(4)) per <- round(100*honeymoon/sum(honeymoon),1) # 비율 lab <- paste(names(honeymoon),"\n",per,"%") # java(cat) pie(honeymoon, labels=lab,col=rainbow(4)) # color palette # pie의 원점은 가운데 지점이다. text(-0.5,0.2,"35%",col="black") text(0.3,0.3,"30%",col="black") text(-0.1,-0.5,"20%",col="black") text(0.5,-0.3,"15%",col="black") # 정규분포와 데이터를 비교한 결과 출력 x <- rnorm(1000, mean=100, sd=1) qqnorm(x) qqline(x, lty=2) x <- runif(1000) # 균등분포 => 정규분포가 아님 qqnorm(x) qqline(x, lty=2) # lattice 패키지(범주형을 고려해서 만든 패키지) library(lattice) str(airquality) xyplot(Ozone ~ Wind, data=airquality) xyplot(Ozone ~ Wind | Month, data=airquality) # int형에 그룹별로 xyplot(Ozone ~ Wind | Month, data=airquality, layout=c(3,2)) convert <- transform(airquality, Month=factor(Month)) xyplot(Ozone ~ Wind | Month, data=convert, layout=c(5,1)) # qplot, ggplot library(ggplot2) str(mpg) qplot(displ, hwy, data=mpg) qplot(displ, hwy, data=mpg, color=drv) # 범례 # drive마다 color값을 줘라 qplot(hwy, data=mpg, fill=drv, binwidth=2) # 범주형 출력 qplot(hwy, data=mpg, fill=drv, facets=.~ drv, binwidth=2) # 중량별 연비 -> carb 종류별로 qplot(wt, mpg, data=mtcars, color=factor(carb), shape=factor(cyl)) # 비선선형회귀 qplot(wt, mpg, data=mtcars, color=factor(cyl), geom=c("point","smooth")) qplot(mpg, wt,data=mtcars, color=factor(cyl),geom="point")+geom_line() # ggplot : layer # Data, Function(요약), Coodinates(축을 지정:aes), # Mapping(color,size,shape,transparency,fill), # Geometries(geom_point:출력 모양), Scales, # Facets(범주별로 출력), Themes(전체적인 배경이나 글꼴) p <- ggplot(mtcars, aes(wt, mpg)) # 기본형태를 지정 p + geom_point() p + geom_point(aes(color=factor(cyl))) p + geom_point(aes(shape=factor(cyl))) p + geom_point(aes(size=qsec)) ggplot(mtcars, aes(wt,mpg)) + geom_point(color='red', size=3) # y축이 지정되지 않음: 자동으로 카운트라는 통계처리를 진행 ggplot(data=mtcars) + geom_bar(mapping = aes(x=cyl, fill=as.factor(am)),position='dodge') ggplot(mpg,aes(x=hwy, fill=drv))+geom_dotplot(binwidth = 2) ggplot(mpg,aes(x=hwy))+geom_histogram(binwidth = 2) # 연속된 수치에 대해 구간별로 카운트 ggplot(mpg, aes(x=hwy)) + geom_density() # 보간법에는 다양한 kernel(선형, 비선형) ggplot(mpg, aes(x=hwy)) + geom_density() + coord_flip() # 문제 # iris 데이터에 대하여 x=Sepal.Length, y=Sepal.width로 해서 geom_point 기하형태로 출력하시오 # point size는 3으로 하시오 # 위의 그래프에서 Species를 구분하도록 모양(shape)을 달리해서 출력하시오 # Species를 범주형으로 3개의 그래프로 출력하시오오 ggplot(iris, aes(Sepal.Length, Sepal.Width, shape=Species))+geom_point(size=3)+facet_wrap( ~ Species, ncol=3) res <- ggplot(iris, aes(Sepal.Length, Sepal.Width, shape=Species)) res <- res+geom_point(size=3) res + facet_wrap( ~ Species, ncol=3) ``` ``` # 도수분포표, 상대도수분포표, 누적도수 분포표를 만든다음 dataframe에 하나의 데이터로 묶고 coin <- (c(2,2,0,1,1,1,2,3,1,2,2,3,3,4,1,3,2,3,3,1,2,2,1,2,2) #ggplot을 이용해서 시각화해보시오.(하나의 화면에 모두 출력되게 하는게 좋습니다.:geom_point,geom_line) # 열 0,1,2,3,4, 열 value, 열 type (coin.freq <- table(coin)) class(coin.freq) # table (coin.rel = round(coin.freq/length(coin),2)) (coin.cum = cumsum(coin.rel)) (coin.num = sort(unique(coin))) # 경우의 수 만들기 # data.frame의 열은 vector로 생성 (coin.freq <- data.frame(coin.num = coin.num, coin.freq=as.vector(coin.freq))) names(coin.freq)[2] <- paste("val") coin.freq$type = rep("freq",length(coin.num)) coin.freq (coin.rel = data.frame(coin.num,coin.rel = as.vector(coin.rel))) (names(coin.rel)[2] <- paste("val")) (coin.rel$type = rep("rel",length(coin.num))) coin.rel (coin.cum <- data.frame(coin.num, coin.cum=as.vector(coin.cum))) names(coin.cum)[2] <- paste("val") (coin.cum$type = rep("cum",length(coin.num))) class(coin.cum) coin.cum # long형으로 묶음 (coin.graph <- rbind(coin.freq,coin.rel,coin.cum)) ggplot(coin.graph, aes(coin.num, val, group=type, col=type))+ geom_point()+ geom_line() ``` ``` #getwd() #dataset <- read.csv("dataset.csv", header=TRUE) #str(dataset) # 다음 코드를 역코딩하시오(주로 범주형 데이터의 값의 변환이고 요약변수를 생성하는 것) # resident2 : 1서울특별시, 2인천광역시, 3대전광역시, 4대구광역시, 5시구군 # job2: 1공무원, 2회사원, 3개인사업 # age2: 청년층, 중년층, 장년층 # position2: 1급,2급,3급,4급,5급 # gender2: 성별 1남자, 2여자 # age3: 청년층, 중년층, 장년층의 level번호로 저장 "C:\\Users\\401-26\\Documents\\R" dataset <- read.csv("C:/Users/401-26/Documents/R/dataset.csv", fileEncoding = "CP949", header=TRUE,) dataset library(tidyverse) library(dplyr) library(dbplyr) library(DBI) library(rstudioapi) library(RMySQL) # 1. 데이터구조 확인 str(dataset) dim(dataset) nrow(dataset) ncol(dataset) length(dataset) summary(dataset) # 숫자변수만 확인인 names(dataset) attributes(dataset) # 2. 다양한 방법으로 데이터 조회 print(dataset) view(dataset) head(dataset) tail(dataset) dataset$age dataset[1,5] dataset[2] dataset[c('age','gender')] # 3. 결측치를 확인하고 - 결측치 갯수를 확인하고 제거하시오 - 결측치의 값을 0, 또는 평균값으로 대체하시오 sum(is.na(dataset)) sapply(dataset, function(x) sum(is.na(x))) nrow(dataset) na.omit(dataset) nrow(dataset) dataset %>% is.na() dataset %>% na.omit() x <- dataset$price sum(is.na(x)) price2 <- na.omit(dataset$price) dataset$price2 <- ifelse(!is.na(x),x,0) # 전체 데이터를 처리하는 방식 dataset$price2 <- ifelse(is.na(x),0,x) dataset$price3 <- ifelse(!is.na(x),x,round(mean(x,na.rm=T),2)) # 4. 이상치를 확인하고 - 변수별 이상치를 확인하시오 - 이상치를 제거하시오 # 5. gender변수를 정제해보시오 (as.factor) # 범주형데이터 -> 숫자매핑(데이터베이스) # boxplot은 정량적 데이터에만 가능 # data.frame:범주형인지 정량적 데이터인지 pie(table(dataset$gender)) # 범주형일 때는 table로 확인해봐야함 dataset<- subset(dataset, gender==1 | gender==2) table(dataset$gender) # 확인 dataset %>% select(gender) %>% filter(gender==1|gender==2) %>% table() # 처리 dataset <- dataset %>% filter(gender==1|gender==2) # 정량적 데이터 boxplot(dataset) boxplot(dataset$price) # 3사분위수+IQR*1.5 dataset2 <- subset(dataset, price >= -100 & price <= 100) # 일정한 값으로 제한 # 6. price변수를 확인하고 범위값으로 제한하시오(2~10) dataset2 <- subset(dataset, price >= 2 & price <= 10 ) # 7. age변수를 20~69사이의 값으로 제한하고 시각화해보시오 dataset2 <- subset(dataset, age >= 20 & age <= 69 ) # 8. resident2변수를 resident변수의 값에 따라 값을 입력하시오(1:서울특별시, 2:인천광역시,3:대전광역시,4:대구광역시,5:시구군군) dataset2$resident2[dataset2$resident == 1] <- '1.서울특별시' dataset2$resident2[dataset2$resident == 2] <- '2.인천광역시' dataset2$resident2[dataset2$resident == 3] <- '3.대전광역시' dataset2$resident2[dataset2$resident == 4] <- '4.대구광역시' dataset2$resident2[dataset2$resident == 5] <- '5.시구군군' # 9. job2변수에 job에 따라서 값을 입력하시오(1:공무원,2:회사원,3:개인사업) dataset2$job2[dataset2$job == 1] <- '공무원' dataset2$job2[dataset2$job == 2] <- '회사원' dataset2$job2[dataset2$job == 3] <- '개인사업' head(dataset2) # 10. 연속형 변수인 age를 age2에 범주화하시오(age<=30:청년층,age<=50:중년층,age>55:장년층) dataset2$age2[dataset2$age <= 30] <- '청년층' dataset2$age2[dataset2$age > 30 & dataset2$age <= 55] <- '중년층' dataset2$age2[dataset2$age > 55] <- '장년층' head(dataset2) # 11. 긍정점수(survey)를 5점 척도를 역순으로 다시 입력하시오(1매우만족, 5매우 불만족) unique(dataset$survey) csurvey <- 6-dataset$survey unique(csurvey) dataset$survey <- csurvey # 12. 거주지역별 성별분포를 시각화하시오 resident_gender <- table(dataset2$resident2, dataset2$gender) barplot(resident_gender, beside=T, horiz=T, col=rainbow(5), legend=row.names(resident_gender),main='성별에 따른 거주지역 분포 현황') # 13. 성별에 따른 거주지역 분포현황을 시각화하시오 resident_gender <- table(dataset2$gender, dataset2$resident2) barplot(resident_gender, beside=T, horiz=T, col=rep(c(2,4),5), legend=c("남자","여자"),main='거주지역별 성별 분포 현황') # 14. 직업 유형에 따른 나이 분포현황을 시각화하시오 densityplot( ~ age, data=dataset2, groups=job2, plot.points=T, auto.key = T) # 15. 성별에 따른 직급별 구매비용을 시각화하시오 densityplot(~price | factor(position), data=dataset2, groups=gender, plot.points=T, auto.key=T) ``` # 3/27 ``` data(iris) head(iris) str(iris) summary(iris) dim(iris) levels(iris$Species) unique(iris$Species) score4 <- c(3,3,6,7,7,10,10,10,11,13,30) (n <- length(x)) # 대표값값 min(score4) max(score4) range(score4) diff(range(score4)) mean(score4) median(score4) sum(score4) sort(score4) order(score4) summary(score4) var(score4) # 분산 sd(score4) # 문제 # var를 수식으로 구해보시오 (total <- sum((score4 - mean(score4))^2)) #542 # 표준편차를 수식으로 구해보시오 (total <- sum((score4 - mean(score4))^2) / length(score4)) #49.2 # 자유도 (표본에서는 전체개수-1 : n-1) (total <- sum((score4 - mean(score4))^2) / (length(score4)-1)) #54.2 sqrt(total) range(score4) quantile(score4) quantile(score4, probs=c(1:100)/100) (Iqr <- quantile(score4)[4] - quantile(score4)[2]) IQR(score4) # 정규분포인 경우 # 1차 원점에 대한 적률 (평균을 빼는가 빼지 않는가에 따라서):누적했다는 말 # 2차 중심적률: 분산 (중심:평균) # 3차 중심적률: 왜도 - 0:좌우 대칭 # 4차 중심적률률: 첨도 - 3:뾰족한 정도 install.packages("moments") library(moments) set.seed(1234) n.sample <- rnorm(n=10000, mean=55,sd=4.5) skewness(n.sample) # 왜도:정규분포이면0 kurtosis(n.sample) # 첨도:정규분포이면3 # 정규분포화하여 분석 => 변환 install.packages("UsingR") library(UsingR) data(cfb) head(cfb) summary(cfb$INCOME) hist(cfb$INCOME, breaks=500, freq=T) # log, sqrt, exp, 제곱 => 변환 cfb <- transform(cfb, INCOME_log = log(INCOME+1)) # 로그를 취하면 0일 때 무한대 hist(cfb$INCOME_log, breaks=500, freq=T) cfb <- transform(cfb, INCOME_sqrt = sqrt(INCOME+1)) hist(cfb$INCOME_sqrt, breaks=500, freq=T) par(mfrow=c(1,3)) qqnorm(cfb$INCOME, main="INCOME") qqline(cfb$INCOME) qqnorm(cfb$INCOME_log, main="INCOME_log") qqline(cfb$INCOME_log) qqnorm(cfb$INCOME_sqrt, main="INCOME_sqrt") qqline(cfb$INCOME_sqrt) par(mfrow=c(1,1)) skewness(cfb$INCOME) kurtosis(cfb$INCOME) skewness(cfb$INCOME_log) kurtosis(cfb$INCOME_log) skewness(cfb$INCOME_sqrt) kurtosis(cfb$INCOME_sqrt) # 키의 집단, 체중 집단:사이즈가 다른 집단을 비교할 때는 표준화 # 대표치:평균, 분산, 표준편차, 표준오차(표본개수고려):구간추정, 변동계수(평균사이즈고려) : 비교 data(trees) dim(trees) head(trees) summary(trees) sd(trees$Height) #6 sd(trees$Girth) #3 sd(trees$Volume) #16 -> 변동이 심하다 # 변동 여부를 변동계수를 만들어서 이질적인 집단에 대한 비교를 해야한다. sd(trees$Volume)/mean(trees$Volume) #0.54 sd(trees$Girth)/mean(trees$Girth) #0.23 sd(trees$Height)/mean(trees$Height) #0.08 -> 변동이 제일 작다 # 분포: 정규분포, 카이제곱분포(범주형데이터), t분포(평균비교), f분포(3집단 이상인 경우) # 신뢰구간 # 표준오차(standard error) x = c(0,1,2,3,4,5,6,7,8,9,10,11) (sem <- sd(x)/sqrt(length(x))) c(mean(x)-2*sem, mean(x)+2*sem) # 95% -> 1.96 : 모집단의 정규분포 # 표본 -> 1.97 : t분포 c(mean(x)-3*sem, mean(x)+3*sem) sqrt(sum((x-mean(x))^2)/(length(x)-1))/sqrt(length(x)) # => 표준오차 c(mean(x) - 2*sem, mean(x) + 2*sem) # 95% # 표준편차의 2배수:95% 신뢰구간 c(mean(x) - 3*sem, mean(x) + 3*sem) # 99% 신뢰구간(국방, 의료) 1% 유의수준을 사용 # 문제 # iris에서 Sepal.Length에 대하여 표준오차를 구하고 95% 신뢰구간 값을 구해보시오 # 구해야할 값: 표준편차, 데이터개수-1 == 자유도, 평균값, 표준오차차 mean_value <- mean(iris$Sepal.Length) n <- length(iris$Sepal.Length) std <- sd(iris$Sepal.Length) st_error <- std/sqrt(n) # 표준오차 degrees_of_freedom = n-1 # 분위수 (0.05유의수준):95% alpha <- 0.05 # 정규분포(rnorm, dnorm(발생빈도), qnorm(확률에 대응), pnorm(누적확률)) # t분포:표본에 대해서는 t분포를 사용 : 자유도에 따라서 분포도가 달라짐 # 확률에 대해서 분위수(확률을 입력하면 -> 분위수) # 자유도에 따라서 분포도가 달라짐(정규분포는 표준정규분포포) (t_score = qt(p=alpha/2, df=degrees_of_freedom, lower.tail = F)) #1.976013, 정규분포에서는 1.96/t분포에서는 1.976 margin_error <- t_score * st_error lower_bound <- mean_value-margin_error upper_bound <- mean_value+margin_error # 신뢰구간 print(c(lower_bound, mean_value, upper_bound)) ``` ``` # 모집단 분포 -> 검정이 가능(표본이 어느지점에서 뽑혔는지를 알 수 있다):유의수준에 따라서(5%(사회적문제),1%(국방,의료)) # 모수 분석(정규분포인 경우), 미모수분석(정규분포가 아닌 경우) => Machine Learning # 선형 회귀 -> tree분석 -> ML rnorm(5,64.5,2.5) # random으로 (5개, 평균, 표준편차) -> 정규분표에서 확률이 가장 높은 분위수는 (평균) # z점수 = (x-xbar)/sd(x) (x <- rnorm(1000,0,1)) (x <- sort(x)) (d <- dnorm(x,0,1)) plot(x,d,main="표준정규분포") curve(dnorm(x),-3,3) plot(density(rnorm(10000,0,1)), main="바보") # interpolation(보간법을 적용) # 확률은 0~1사이의 값 pnorm(0) #0.5 pnorm(1) #0.84 pnorm(1)-(1-pnorm(1)) #0.68 # 2사분위수의 누적 확률 pnorm(2)-(1-pnorm(2)) pnorm(1.96)-(1-pnorm(1.96)) pnorm(3)-(1-pnorm(3)) pnorm(1)-pnorm(-1) dnorm(0,0,1) # 분위수 -> 확률 dnorm(1,0,1) # 문제 # 평균이 100이고 표준편차가 10인 정규분포에서 50이 나올 확률은? z = (50-100)/10 # z점수 -> 표준정규분포(평균이 0, 표준편차 1인 경우) dnorm(z,0,1) #1.48672e-06 dnorm(50,100,10) # 평균이 90이고 표준편차가 12인 정규분포에서 98이 나올 확률은? # 한 지점에서 평균이 90이고 12일 때의 누적확률, 이것을 표준화(z점수) 했을 때의 # 누적확률은 같은데 한 지점에서 출현확률은 다르다. # 매핑 # 상태가 같은 상태:한 지점에서의 확률값 => 표준정규분포 확률 값 # 범위값에 해당하는 함수 # z점수: 표준정규분포 매핑 # z = (98-90)/12 # 분위수를 구하는 값 => 확률 dnorm # dnorm(z,0,1) #0.319448 -> 잘못된 방법..? dnorm(98,90,12) qnorm(0.319448) plot(seq(-3.2,3.2,length=50),dnorm(seq(-3,3,length=50),0,1),type="l",xlab="",ylab="",ylim=c(0,0.5)) text(x=0,y=0.45,labels=expression("99.7%가 표준편차 3배수 안에 3"~sigma)) arrows(x0=c(-2,2),y0=c(0.45,0.45),x1=c(-3,3),y1=c(0.45,0.45)) segments(x0=c(-2,2),y0=c(-1,-1),x1=c(-2,2),y1=c(0.4,0.4)) text(x=0,y=0.3,labels=expression("95%가 표준편차 2배수 안에 2"~sigma)) arrows(x0=c(-1.5,1.5),y0=c(0.3,0.3),x1=c(-2,2),y1=c(0.3,0.3)) segments(x0=c(-1,1),y0=c(-1,-1),x1=c(-1,1),y1=c(0.25,0.25)) text(x=0,y=0.15,labels=expression("68%가 표준편차 1배수 안에 1"~sigma),cex=0.9) # 확률의 범위 표현 curve(dnorm(x,0,1),-4,4,xlab="z",ylab="f(z)") z=seq(0,4,0.02) lines(z,dnorm(z),type="h",col="grey") # 문제 # 비행기의 평균 비행시간은 120시간이라고 하자. 비행시간이 정규분포이고 # 표준편차가 30시간이라고 가정하고, 9대 표본으로 추출할 때와 36대를 표본으로 추출할 때, # 표본 평균과 분산을 비교해보시오 # 표본이 많아지면->모집단 # 표본이 많아지면 분산이 커짐 # 표본이 많아지면 표준오차는 줄어든다. (t9 = rnorm(9,120,30)) (t36 = rnorm(36,120,30)) mean(t9) var(t9) mean(t36) var(t36) stderr <- function(x) sd(x,na.rm=T)/sqrt(length(na.omit(x))) # 표준오차 stderr(t9) stderr(t36) # 문제 # X~(300,50^2)인 정규분포에서 P(X >= 370)일 확률을 구해라 (z점수 사용) # 평균이 300이고 분산이 50^2 1-pnorm(370,300,50) # 분포를 입력 1-pnorm((370-300)/50) # z점수 (x <- rnorm(1000,0,1)) curve(dnorm(x,0,1),-4,4,xlab="z",ylab="f(z)") z=seq(((370-300)/50),4,0.01) lines(z,dnorm(z),type="h",col="grey") points(((370-300)/50),dnorm(((370-300)/50))) # 브랜드의 백열전구의 수명이 1500시간의 평균값과 75시간의 표준편차로 정규적으로 분포되어있다. # 백열전구가 1410시간보다 덜 오래갈 확률은? z=(1410-1500)/75 pnorm(z) # 백열전구가 1563과 1648시간 사이의 오래갈 확률은? z1=(1563-1500)/75 z2=(1648-1500)/75 pnorm(z2)-pnorm(z1) curve(dnorm(x,0,1),-4,4,xlab="z",ylab="f(z)") z=seq(-4,(1401-1500)/75,0.01) lines(z,dnorm(z),type="h",col="grey") curve(dnorm(x,0,1),-4,4,xlab="z",ylab="f(z)") z=seq((1563-1500)/75,(1648-1500)/75,0.01) #z=seq(0,4,0.02) lines(z,dnorm(z),type="h",col="grey") # 우리나라에서 사육하고 있는 생후 18개월 이상 된 황소 무게는 평균이 500kg이고 # 표준편차가 50kg인 정규분포라 한다. 이제 우량 한우를 집중 육성관리하기 위해 # 무게가 무거운 순서대로 5%에 해당하는 황소를 선발하고자 한다. # 그렇다면 몇 kg 이상인 황소를 선발해야 하는가? qnorm(1-0.05, 500, 50) # 582.2427kg 이상이 황소를 선정 # z점수를 구하는 함수(scale) (x<-matrix(1:10, ncol=2)) scale(x) # 평균, 표준편차(열중심으로 계산):z점수 정규화 (centered.x <- scale(x,scale=F)) # 평균 # 평균비교(정규화):두 집단이 평균이 다르고 표준편차가 다르면 직접 비교가 안됨 # z점수 표준화 후 비교 set.seed(100) h_korean <- rnorm(n=1000,mean=170,sd=10) h_bushman <- rnorm(n=1000,mean=140,sd=8) height <- data.frame(h_korean, h_bushman) # R은 복사에 의해서 전달 # rm(h_korean, h_bushman) -> 데이터 없애는거 par(mfrow=c(1,2)) hist(h_korean, freq=T, main="한국인 키 변도 히스토그램") hist(h_korean, freq=F, main="한국인 키 확률밀도함수 그래프") hist(h_bushman, freq=T, main="부시맨 키 변도 히스토그램") hist(h_bushman, freq=F, main="부시맨 키 확률밀도함수 그래프") heigth <- transform(height, z2.h_korean=(h_korean-mean(h_korean))/sd(h_korean),z2.h_bushman=(h_bushman-mean(h_bushman))/sd(h_bushman)) height <- transform(height,z.h_korean=scale(h_korean),z.h_bushman=scale(h_bushman)) hist(height$z.h_korean,freq=T,main="한국인 표준") hist(height$z.h_bushman,freq=T,main="부시맨 표준") # 문제 # 160cm의 키가 한국인과 부시맨에서 어떤 위치인가 비교해보시오 (160-mean(h_korean))/sd(h_korean) # 역 1배수 (160-mean(h_bushman))/sd(h_bushman) # 정 2.5배수 dnorm(160, mean(h_korean),sd(h_korean)) dnorm(160, mean(h_bushman),sd(h_bushman)) ``` # 3/28 ``` # 카이제곱분포와 자유도 x <- seq(1,10,.1) par(mfrow=c(2,3)) plot(x,dchisq(x,6),type="l") # probability density plot(x,dchisq(x,5),type="l") plot(x,dchisq(x,4),type="l") plot(x,dchisq(x,3),type="l") plot(x,dchisq(x,2),type="l") plot(x,dchisq(x,1),type="l") # factor를 기준으로 한 도수분포표 # 카이제곱분석에서는 인접행렬을 필요로 함함 x = c(1,2,3,4,5,6) y = c(7,8,9,10,11,12) z=10:15 (ta=table(x)) (tb=table(x,y)) # 인접행렬(contingency table), 그래프로 표현하는 방법 (tc=table(x,y,z)) class(tc) #table -> data.frame 아님 ?xtabs (d <- data.frame(x=c("1","2","2","1"), y=c("A","B","A","B"), num=c(3,5,8,7))) class(d) str(d) (xt = xtabs(num ~ y+x,data=d)) # xtabs(도수 ~ 가로+세로) class(xt) margin.table(xt,1) # 1:행방향 합 margin.table(xt,2) # 2:열방향 합 margin.table(xt) #전체합 prop.table(xt,1) # 비율로 표현 prop.table(xt,2) prop.table(xt) Titanic # Class, Sex, Age, Survived class(Titanic) #table margin.table(Titanic,1) # 합계: 선실별 전체 인원수 4가지 margin.table(Titanic,2) # Male, Female margin.table(Titanic,3) # Child, Adult margin.table(Titanic,4) # Np, Yes margin.table(Titanic) barplot(margin.table(Titanic,1)) barplot(margin.table(Titanic,2)) barplot(margin.table(Titanic,3)) barplot(margin.table(Titanic,4)) # Freq:frequency (Titanic.data = xtabs(Freq ~ Survived+Sex, data=Titanic)) # 계층적 인덱스를 이용해서 표현 ftable(Titanic, row.vars = 1:3) ftable(Titanic, row.vars = 2:1, col.vars = "Survived") # xtabs와 ftable의 차이가 무엇인지 확인하시오 # 성별별로 생존여부를 확인하시오 ftable(Titanic, row.vars="Sex",col.vars = "Survived") xtabs(Freq ~ Sex+Survived, data=Titanic) # 객실과 나이별로 데이터를 확인하시오 ftable(Titanic, row.vars="Age",col.vars = "Class") xtabs(Freq ~ Class+Age, data=Titanic) #install.packages("gmodels") # table + chisq(카이제곱분석) library(MASS) library(gmodels) data(mtcars) CrossTable(mtcars$vs,mtcars$am) CrossTable(mtcars$vs,mtcars$am, expected=T) 0.161+0.236+0.207+0.303 # 1 일반횟수 # 2 카이로제곱 # 3 행을 기준으로 비율값 # 4 열을 기준으로 비율 값 # 5 전체를 기준으로 비율 값 19/32 * 18 # 10.6875 13/32 * 18 # 7.3125 # 검정통계랑 -> chisq분포도에서 확률값을 확인한다. (12-10.6875)^2 / 10.6875 # 0.1611842 (6-7.3125)^2 / 7.3125 # 0.2355769 library(MASS) CrossTable(infert$education, infert$induced, expected = T) fisher.test(infert$education, infert$induced, alternative = "two.sided") getwd() "C:\\Users\\401-26\\Documents\\R" data <- read.csv("C:/Users/401-26/Documents/R/cleanDescriptive_1.csv", fileEncoding = "cp949", header=TRUE,) str(data) class(data) # 부모학력수준(level2)과 자녀대학진학여부(pass)에 관련이 있는지 검정하시오(독립성 검정) # 가설설정 # 귀무가설:상호 관련성이 없다. # 대립가설:상호 관련성이 있다. # 유의수준:0.05 x <- data$level2 y <- data$pass result <- data.frame(Level=x, Pass=y) sapply(result, function(x) sum(is.na(x))) #result <- na.omit(result) #sapply(result, function(x) sum(is.na(x))) CrossTable(x=result$Level, y=result$Pass, expected = T) chisq.test(x=result$Level, y=result$Pass) # p-value=0.2507 # 유의수준:0.05 # p-value > 0.05 그러므로 귀무가설을 기각할 수 없다. # 따라서 부모학력수준과 자녀대학진학여부는 상호 관련성이 없다. # 카이제곱분석: 적합성, 독립성, 동질성 검정 # 사이다별로 선호도 차이가 있는가? data <- textConnection( "사이다종류 관측도수 1 14 2 32 3 17 4 9 5 18") (x <- read.table(data, header=T)) # 귀무가설:선호도 차이가 없다. # 대립가설:선호도 차이가 있다. # 유의수준: 0.05 (14 + 32 + 17 + 9 + 18)/5 # 기대도수: 18 res <- (14-18)^2 + (32-18)^2 + (17-18)^2 + (9-18)^2 + (18-18)^2 (res <- res/18) 1 - pchisq(16.333,4) chisq.test(x$관측도수) # x-squared=16.333, df=4, p-value=0.002603 # p-value < 0.05 귀무가설을 기각할 수 있다. # 따라서 사이다 제품별로 선호도 차이가 있다. library(gmodels) library(MASS) str(survey) # 문제 # 성별과 운동량은 상호 관련성이 있는가 # 귀무가설은 관련성이 없다. # 대립가설은 관련성이 있다. # 유의수준: 0.05 sum(is.na(survey)) sapply(survey, function(x) sum(is.na(x))) survey <- na.omit(survey) data(survey) (xt=xtabs(~ Sex+Exer, data=survey)) # formula를 이용해서서 (xt=table(survey$Sex, survey$Exer)) # 기본 (xt=ftable(survey$Sex, survey$Exer)) # 데이터를 요약할 때 편함 CrossTable(survey$Sex, survey$Exer, expected=T) chisq.test(xt) # Chi^2 =4.158519 d.f=2 p=0.1250228 # 검정통계량 chi^2 = 4.158519 # p-value = 0.1250228 # p-value > 0.05 귀무가설을 기각할 수 없다. # 따라서 성별과 운동량은 관련이 없다. # 박수치는 손과 글씨 쓰는 손의 독립성 검정을 하시오 # 귀무가설: 상호 독립적이다. # 대립가설(연구가설): 상호 상관성이 있다. # 유의수준: 0.05 chisq.test(xtabs(~W.Hnd + Clap, data=survey)) fisher.test(xtabs(~W.Hnd + Clap, data=survey)) # p-value > 0.05 귀무가설을 기각할 수 없다. # 성별과 운동량은 상호 관련성이 없다. # p-value=0.007371 # 귀무가설을 기각하고 대립가설을 채택한다. # 글을 쓰는 손과 박수치는 손은 상호 상관성이 있다. # 동질성 검정 # 왼손잡이와 오른손잡이의 비율이 0.3:0.7이다. # 이때 표본이 이런 경향을 보이는가 chisq.test(table(survey$W.Hnd),p=c(.3,.7)) # 귀무가설은 비율이 동일하다. # 대립가설은 비율이 동일하지 않다. # 감정통계량: X-squared=41.796, df = 1, # p-value=1.013e-10 # p-value < 0.05 귀무가설을 기각하고 대립가설을 채택한다. # 비율이 0.3:0.7이 아니다. ``` ``` # 다음 데이터에 대하여 카이제곱식으로 계산해서 p-value값을 구하고 브랜드별 선호도 조사를 검정하시오 data <- textConnection( "브랜드종류 관측도수 BrandA 18 BrandB 23 BrandC 19") (x <- read.table(data, header=T)) # 가설설정 # 계산식으로 하시오 freq.mean = sum(x$관측도수)/length(x$관측도수) #freq.mean=(18+19+23)/3 res <- (18-freq.mean)^2/freq.mean + (23-freq.mean)^2/freq.mean + (19-freq.mean)/freq.mean p.value <- 1-pchisq(res, length(x$관측도수)-1) p.value < 0.05 # 귀무가설을 기각할 수 없다. # 브랜드별 선호도는 차이가 없다. getwd() "C:\\Users\\401-26\\Documents\\R" data <- read.csv("C:/Users/401-26/Documents/R/homogenity.csv", fileEncoding = "cp949", header=TRUE,) str(data) class(data) # method:방법1, 방법2, 방법3 # survey 만족도: 1:매우만족, 2:만족, 3:보통, 4:불만족, 5:매우불만족 # 문제 # 데이터를 변환하고 교차표를 생성한 다음 교육방법별 동질성 검정하고 해석하시오 # 보는 관점으로 생각하면 # 독립성검정하고 같은 방법을 하는데 기대값으로 계산한 검정통계량의 p-value가 유의 수준보다 작지않으면 같다고 본다. data <- subset(data, !is.na(survey),c(method,survey)) nrow(data) data$method2[data$method==1] <- "방법1" data$method2[data$method==2] <- "방법2" data$method2[data$method==3] <- "방법3" data$survey2 = sapply(data$survey, switch, "1"="매우만족","2"="만족","3"="보통","4"="불만족","5"="매우불만족") data$method2 = sapply(data$method, switch, "1"="방법1","2"="방법2","3"="방법3") # 교차표 table(data$method2, data$survey2) xtabs(~method2 + survey2, data=data) ftable(data, row.vars="method2", col.vars="survey") # 동질성 검정 # sigma (관측도수 - 기대도수)^2 chisq.test(data$method2, data$survey2) CrossTable(data$method2, data$survey2, expected=T) # 검정 통계량: X-squared=1.5026, df=4, # p-value=0.8262 > 0.05 # 귀무가설을 기각하지 못한다. # 교육방법에 따른 만족도에는 차이가 없다. # 평균분석(t.test) : 두 집단 분석 # t분포:정규분포보다 완만한 분포 # 정규분포테스트(shapiro.test) # 모수분석(t.test()) 비모수분석(wilcox.test()) # dependent/independent(연구자만 동일 표본을 사용했는지 알 수 있다) pairs # 등분산테스트(var.test()):var.equal=T,F # anova test : 두 집단 이상 # 분산분석(F=군집간분석/군집내분산) 검정통계량 # 자유도에 의해서 영향을 받음 # 정규분포 테스트 # 모수분석(aov()) 비모수분석(kruskal.test) # 등분산테스트(모수적:bartlett.test, 비모수:fliger.test) # 선형회귀처럼 사후분석 : TurkeyHSD() # 단일표본 평균분석 x = c(65,78,88,55,48,95,66,57,79,81) t.test(x, mu=75) # 귀무가설: 두 집단 평균이 같다. # 대립가설: 두 집단 평균이 같지 않다. # 폐암 발생률 평균이 20이라는 가설을 검증하라 (lung <- runif(20,min=19,max=20)) t.test(lung, mu=20, conf.level = .99) # p-value = 5.963e-09 # 귀무가설을 기각하고 대립가설 채택 ``` ``` str(sleep) ?sleep # 두 집단에 대한 평균검정 attach(sleep) par(mfrow=c(1,1)) plot(extra ~ group, data=sleep) var.test(extra[group==1], extra[group==2], data=sleep) # 등분산테스트 # 귀무가설:두 집단의 분산차가 없다. => 등분산이다. # 대립가설:두 집단의 분산차가 있다. with(sleep, t.test(extra[group==1], extra[group==2], var.equal = T)) detach(sleep) # p-value = 0.07919 # 귀무가설을 기각할 수 없고, 두 집단 수면차이는 없다. ?ToothGrowth my_data <- ToothGrowth str(ToothGrowth) res.ftest <- var.test(len ~ supp, data=my_data) res.ftest res.ftest$estimate res.ftest$p.value < 0.05 # 등분산이다. shapiro.test(rnorm(1000)) # p-value = 0.5606 # 귀무가설 채택 set.seed(450) x <- runif(50,min=2,max=4) shapiro.test(x) # p-value = 0.003259 # 대립가설 채택(유의수준(0.05)보다 작은경우) # 귀무가설:정규분포이다. # 대립가설:정규분포가 아니다. # 양측검정과 단측검정 x <- rnorm(1000,5.0,0.5) t.test(x, mu=5.2, alter="two.side",conf.level = 0.95) (result <- t.test(x, mu=5.2, alter="greater", conf.level = 0.95)) (result <- t.test(x, mu=5.2, alter="less", conf.level = 0.95)) ``` ``` # 문제 # iris의 Sepal.Width와 Sepal.Length의 등분산성테스트를 하시오 x <- rnorm(1000,5.0,0.5) var.test(x,x) # 귀무가설은 등분산이다. # 대립가설은 등분산이 아니다. with(iris, var.test(Sepal.width, Sepal.Length)) # p-value = 3.595e-14 등분산이 아니다. # 새로운 당뇨병 치료제를 개발한 제약사의 예 # 치료에 지대한 영향을 주는 외부요인을 통제하기 위해 10명의 당뇨병 환자를 선별하여 1달 동안 #'위약(placebo)'을 투여한 기간의 혈당당(x1)과 '신약(new medicine)'을 투여한 1달 기간 동안의 # 혈당 수치 (x2)를 측정하여 짝을 이루어 혈당 차이를 유의수준 5%에서 비교하시오 # 표본이 동일집단 # 유의미한 차이가 있는가, 차이가 있다면 약효가 있다. x1 <- c(51.4, 58.0, 45.5, 55.5, 52.3, 50.9, 52.7, 50.3, 53.8, 53.1) x2 <- c(50.1, 51.5, 45.9, 53.1, 51.8, 50.3, 52.0, 49.9, 52.5, 53.0) # 정규성테스트 # 등분산테스트 shapiro.test(x1) # 정규성 shapiro.test(x2) # 비정규성 var.test(x1,x2) # p-value = 0.273 등분산이다. t.test(x1,x2, alternative=c('greater'), paired=T, var.equal=T, conf.level=0.95) # 검정통계량 : t = 2.153, df =9, # p-value = 0.02987 # 귀무가설은 평균이 동일하다. # 귀무가설을 기각하지 못함 -> 약효는 없다. ?wilcox.test # 순위테스트 -> 중위수(median)가 같다. wilcox.test(x1, x2, alternative = c("greater"), paired = T, conf.int=F, exact=F, # TIE가 있어 정확한 P값을 계산할 수 없다. conf.level = 0.95) # 검정통계량 V=52.5, # p-value=0.006201 # 평균분석의 귀무가설을 기각하고 대립가설을 채택한다. # 당뇨병치료제는 약효가 있다. wilcox.test(x1, x2, paired = T, var.equal=T, exact=F, conf.level = 0.95) # var.equal=T를 추가한 경우 # 검정통계량 V=52.5, p-value = 0.0124 -> 값이 달라짐 # 평균분석의 귀무가설을 기각하고 대립가설을 채택한다. # anova test ?InsectSprays # 모기약 data(InsectSprays) attach(InsectSprays) str(InsectSprays) InsectSprays xtabs(InsectSprays) # A,B,C,D,E,F 6개의 제품 aov.out = aov(count ~ spray, data=InsectSprays) summary(aov.out) # 회귀분석한 결과와 비슷 TukeyHSD(aov.out) ``` ``` # 평균분석 # 문제 # Cars93 데이터를 이용해 생산국이 USA vs non-USA(Origin) 두 그룹간의 # price 평균이 차이가 있는지를 테스트하시오, 95% 신뢰구간에서 구하시오 ?Cars93 library(MASS) library(tidyverse) library(dplyr) data(Cars93) str(Cars93) (usa_cars <- Cars93 %>% select(Price, Origin) %>% filter(Origin=="USA")) (non_usa_cars <- Cars93 %>% select(Price, Origin) %>% filter(Origin=="non-USA")) (usa_price_mean <- mean(usa_cars$Price)) (non_usa_price_mean <- mean(non_usa_cars$Price)) usa_cars <- subset(Cars93, Origin == "USA") non_usa_cars <- subset(Cars93, Origin == "non-USA") t.test(usa_cars$Price, non_usa_cars$Price, conf.level = 0.95) # 검정통계량: t = -0.95449, df = 77.667, # p-value = 0.3428 p.value < 0.05 # 귀무가설을 기각할 수 없고, 생산국에 따른 price 평균 차이는 없다. ``` # 3/29 ``` # 주성분 분석은 10개의 변수가 입력 -> 10개의 고유치와 10x10 고유벡터행렬이 생성 # 10x10 상관계수행렬이나 공분산행렬이 생성 # 고유값 분해 # 10개의 고유치(순서대로 정렬이 되어서 출력)와 10x10 고유벡터행렬이 생성 # 공분산 (data1 = 1:5) (data2 = 2:6) sum((data1 - mean(data1))*(data2 - mean(data2))) / (length(data1) - 1) cov(1:5, 2:6) cor(1:5, 2:6) #1 -> 완전 정상관계가 있다. # 주성분 분석은 공분산행렬이나 상관계수행렬을 이용 # 상관계수를 구하는 방법: pearson, spearman, kendall => 이 중에서 가장 작은 것을 사용용 (m <- matrix(c(1:10, (1:10)^2),ncol=2)) cor(m, method="pearson") # 모수적 방법(t.test): 데이터가 정규분포 cor(m, method="spearman") # 비모수적 방법(kruskal.test): 정규분포가 아님 -> rank에 의해 cor(m, method="kendall") # 문제 # 다음 두 데이터에 대하여 공분산과 상관계수를 구하시오 (a <- c(1:10)) (b <- (1:10)^2) # 공분산 계산 cov(a,b) sum((a-mean(a))*(b-mean(b)))/(length(a)-1) # 상관계수 계산 cor(a,b) ( sum((a-mean(a))*(b-mean(b)))/(length(a)-1) ) / ( sd(a)*sd(b) ) # 상관분석(cor.test()사용) # 귀무가설은 서로 독립이다.(상관이 없다.) # 대립가설은 서로 상관이 있다. cor.test(c(1,2,3,4,5),c(1,2,3,4,5),method='pearson') cor.test(c(1,2,3,4,5),c(1,2,3,4,5),method='spearman') cor.test(c(1,2,3,4,5),c(1,2,3,4,5),method='kendall') # p-value < 2.2e-16 : 귀무가설을 기각하고 대립가설을 채택한다. (0.05보다 작음음) # 문제 # iris 데이터에 대해서 전체 변수에 대한 상관계수를 구해보시오 # 상관분석이나 회귀분석은 연속형 데이터에서만 가능 library(dplyr) str(iris) cor(iris[,1:4]) iris %>% select(-(Species)) %>% cor() symnum(cor(iris[,1:4])) d <- data.frame(x1=rnorm(10), x2=rnorm(10), x3=rnorm(10)) (M <- cor(d)) #install.packages("corrplot") library(corrplot) corrplot(M,method="shade") # 상관관계 분석 product <- read.csv("C:/Users/401-26/Documents/R/product2.csv", fileEncoding = "cp949", header=TRUE,) head(product) # 5점 척도 summary(product) sd(product$"제품_친밀도");sd(product$"제품_적절성");sd(product$"제품_만족도") cor(product$"제품_친밀도",product$"제품_적절성",product$"제품_만족도") cor(product, method="pearson") #install.packages("corrgram") library(corrgram) corrgram(product) corrgram(product, upper.panel=panel.conf) #install.packages("PerformanceAnalytics") library(PerformanceAnalytics) chart.Correlation(product, histogram=,pch="+") # 문제 # airquality 데이터에 대하여 상관계수를 구해보시오 # 순서: 1 분석 대상 변수만 선택, 2 na는 제거, 3 데이터의 상관정도를 확인, 4 상관계수행렬을 시각화 str(airquality) airquality_1 <- airquality[,c(1:4)] head(airquality_1) cor(airquality_1) (sum(is.na(airquality_1))) supply(x,function(x), sum(is.na(x))) (airquality_2 <- na.omit(airquality_1)) (airquality_cor <- cor(airquality_2)) par(mfrow=c(1,1)) library(corrplot) plot(airquality_2) corrplot(airquality_cor, method="circle") corrplot(airquality_cor, method="square") # ellipse, number, shade, pie (col <- colorRampPalette(c("darkblue","white","darkorange"))(20)) heatmap(x = airquality_cor, col=col, symm=T) library(PerformanceAnalytics) chart.Correlation(airquality_cor, histogram=T, pch=19) cor.test(airquality_2$Wind, airquality_2$Temp, method="spearman", exact=F) # p-value = 2.426e-08 귀무가설을 기각하고 대립가설 채택 # 바람과 온도는 상호 상관관계가 있다. library(corrgram) corrgram(airquality_cor, order=T, lower.panel=panel.shade, upper.panel = panel.pie, text.panel = panel.txt, main="공기질상관계수") ``` ``` # 상관계수행렬과 고유값 분해 a <- c(4.0,4.2,3.9,4.3,4.1) b <- c(2.0,2.1,2.0,2.1,2.2) c <- c(0.60,0.59,0.58,0.62,0.63) (mat <- matrix(c(a,b,c), nrow=5, byrow=F)) # 열로 생성 (avr <- colMeans(mat)) (acha <- a - avr[1]) (bcha <- b - avr[2]) (ccha <- c - avr[3]) (aa_var <- sum(acha*acha)/(length(a)-1)) (ab_var <- sum(acha*bcha)/(length(a)-1)) (ac_var <- sum(acha*ccha)/(length(a)-1)) cov(mat) (aa_cor <- aa_var/(sd(a)*sd(a))) (ab_cor <- ab_var/(sd(a)*sd(b))) (ac_cor <- ac_var/(sd(a)*sd(c))) cor(mat) # 상관계수행렬을 이용해서 고유값 분해 # 변수의 개수:3개 => 상관계수행렬은 3x3행렬로 나타남(정방행렬) (eigenresult <- eigen(cor(mat))) # 고유값 분해 # 크기와 방향으로 분리 eigenresult$values # 고유치(크기) eigenresult$vectors # 고유벡터(방향) # 내적이 0이면 직교다. -> 고유벡터는 정직교하는 축이라고 할 수 있다./고유치는 정렬되어 있다.(큰 값이 주성분이다.) eigenresult$vectors[,1] %*% eigenresult$vectors[,2] #-5.551115e-17 => 0으로 본다. eigenresult$vectors[,2] %*% eigenresult$vectors[,3] #5.551115e-17 => 0으로 본다. eigenresult$vectors[,1] %*% eigenresult$vectors[,3] #0 eigenresult$vectors[1,] %*% eigenresult$vectors[2,] #6.938894e-18 => 0으로본다. eigenresult$values # 정직교하는 축에 데이터를 재표현하면 noise를 제거한 상태 # 주성분을 선택하면 noise가 있는 변수가 제거된다. # 데이터 압축(차원축소)을 하는 효과가 있다. # 이미지 640*640 이미지: 이미지 하나당 409600변수 # 100개만 주성분으로 취하면 원래 이미지로 복원이 가능하다. # 새롭게 생성된 축은 모든 변수가 참여를 해서 만들어진 축 # 이때 시각화를 위해서 선택되는 축은 3개의 변수를 대표하는 축이다. # 그래서 변수만 선택해서 시각화하는 것과는 차이가 있다. # weather 데이터를 이용해서 상관계수행렬을 구하고 고유값 분해를 실시한 다음 # 그 결과가 단위행렬이 되는지를 확인해보시오 weather <- read.csv("C:/Users/401-26/Documents/R/weather.csv", fileEncoding = "cp949", header=TRUE,) str(weather) head(weather) nrow(weather) # 질적 전처리 sum(is.na(weather)) sapply(weather, function(x) sum(is.na(x))) (weather <- na.omit(weather)) summary(weather) weather$key=1:nrow(weather) library(dplyr) outdata=select_if(weather, is.numeric) str(outdata) (res_dir <- cor(outdata)) # 정방행렬이면서 대칭행렬인 상관계수행렬이 구해짐 (res_eig <- eigen(res_dir)) # 고유값 분해 # 숫자변수 고유치:11, 고유벡터:11x11 res_eig$values # 고유치의 크기 해당하는 축방향으로의 분산의 크기 # 통계분석에서는 분산이 크면 중요한 변수(주성분 분석) dim(res_eig$vectors) #11x11 length(res_eig$values)#11 heatmap(res_eig$vectors) res_eig$vectors %*% t(res_eig$vectors) # 대각선만 1, 나머지는 0이 되는 행렬 나옴 -> 단위행렬 # 행렬 %*% 역행렬 = 단위행렬 # 행렬의 항등원 == 단위행렬 # 정방행렬이면서 대칭행렬인 행렬은 전치행렬이 역행렬이다. res_eig$vectors[,1] %*% res_eig$vectors[,2] # 고유벡터는 정직교하는 축벡터의 합이다. # 내적하면 0 나옴 res_eig$vectors[,1] %*% res_eig$vectors[,3] ``` ``` # 주성분 분석 str(USArrests) head(USArrests) pairs(USArrests, panel=panel.smooth, main="USArrests data") # 주성분 선정에 있어서 변수의 크기가 주성분을 결정하면 안되기 때문에 scale prcomp(USArrests, scale=T) # 주성분분석 함수 = 상관계수행렬 + 고유값분해 # 주성분 분석을 하게 되면 정직교하는 축이 결정되는데 이 축은 # 각 변수들이 이 축에 미치는 영향을 고려해서 만들어진 축이다. # 주성분의 의미를 알 수 없어서 PC1이라고 이름을 줌 # 변수의 참여정도에 따라서 전문가가 의미 이름으로 재명명해야 함 plot(prcomp(USArrests)) (fit <- prcomp(USArrests, scale=T)) # princomp보다 낫다. plot(fit, type="lines") names(fit) fit$sdev fit$rotation fit$center summary(prcomp(USArrests, scale=T)) # elbow # 변수가 작으면 bias가 커지는 과소적합문제 발생 # 변수가 많으면 bias는 작아지지만 variance가 커지는 과대적합문제가 발생:차원축소 # 차원축소를 하는 방법: PCA(주성분 분석), FA US.prin <- princomp(USArrests, cor=T) screeplot(US.prin, npcs=4, type="lines") #install.packages("FactoMineR") #install.packages("factoextra") library(FactoMineR) library(factoextra) data(iris) res.pca <- prcomp(iris[,-5], scale=T) get_eig(res.pca) fviz_eig(res.pca, addlabels = T, ylim=c(0,85)) str(mtcars) head(mtcars, 10) fit <- princomp(mtcars, cor=T) # 상관계수로 하는 경우 summary(fit) loadings(fit) plot(fit, type="lines") fit$scores fit$scores[,1:3] fit$scores[,1:4] biplot(fit) # 주성분을 선택, 다중공선성(독립변수가 상호영향을 주면) # 문제 # 다음 데이터를 가지고 주성분 분석을 하여 데이터가 주성분에 미치는 영향을 확인해보시오 # scree plot으로 주성분 확인하기(주성분을 몇개로 할 것인가) # loadings를 이용하여 주성분에 기여하는 정도 확인하기(새로운 축으로 만들어진 주성분 축에 각 변수가 기여한 점) # biplot을 이용하여 주성분에 기여하는 정도를 그래프로 확인하기(주성분의 성분을 확인-새로운 이름으로 결정) (x <- 1:10) (y <- x+runif(10,min=-.5,max=.5)) (z <- x+y+runif(10,min=-10,max=10)) (data <- data.frame(x,y,z)) # 10x3 (pr <- prcomp(data, scale=T)) summary(pr) pr$sdev names(pr) pr$rotation # 고유값 분해 했을 때의 고유벡터 # 3x3 # 회전행렬 # princomp에서는 loadings/ prcomp에서는 rotation pr$x # PC1,PC2,PC3 고유벡터의 축으로 표현된 데이터 biplot(pr) # x값은 원래 데이터가 변환이 되었다(변환된 기준은 새로운 축) # FA(factor analysis): 요인 분석 v1 <- c(1,1,1,1,1,1,1,1,1,1,3,3,3,3,3,4,5,6) v2 <- c(1,2,1,1,1,1,2,1,2,1,3,4,3,3,3,4,6,5) v3 <- c(3,3,3,3,3,1,1,1,1,1,1,1,1,1,1,5,4,6) v4 <- c(3,3,4,3,3,1,1,2,1,1,1,1,2,1,1,5,6,4) v5 <- c(1,1,1,1,1,3,3,3,3,3,1,1,1,1,1,6,4,5) v6 <- c(1,1,1,2,1,3,3,3,4,3,1,1,1,2,1,6,5,4) (med.data <- cbind(v1,v2,v3,v4,v5,v6)) head(med.data) summary(med.data) cor(med.data) #install.packages("psych") #install.packages("GPArotation") library(psych) library(GPArotation) med.factor <- principal(med.data, rotate='none') names(med.factor) med.factor$values op <- par(mfrow=c(1,1)) plot(med.factor$values, type="b") ?factanal # varimax 알고리즘(소스분석) # 요인분석을 하게되면 축이 다시 결정됨(회전해야됨:rotation -> 그 축으로 투영을 해야함) (stats.fact <- factanal(med.data, factors = 3, rotation = "oblimin",scores="Bartlett")) # 요인분석 names(stats.fact) stats.fact$scores # 재표현된 데이터(3개의 축으로) # 3개의 새로운 축에 대한 명명식() stats.fact$loadings # 설명(요인에 기여한 변수들의 역할) stats.fact$factors # 새로운 축 개수 stats.fact$criteria # 최적화의 결과 stats.fact$correlation # 원 데이터에 대한 상관계수 stats.fact$rotmat # 변환을 위한 직교축(회전행렬) stats.fact$converged #True stats.fact$dof #0 stats.fact$method #mle 선형회귀 stats.fact$call # 호출된 함수내용 stats.fact$n.obs ``` ``` dust <- read.csv("C:/Users/401-26/Documents/R/재비산먼지.csv", fileEncoding = "cp949", header=TRUE,) head(dust) str(dust) summary(dust) sd(product$"제품_친밀도");sd(product$"제품_적절성");sd(product$"제품_만족도") cor(product$"제품_친밀도",product$"제품_적절성",product$"제품_만족도") cor(product, method="pearson") sd(dust$"기온");sd(dust$"습도");sd(dust$"재비산먼지평균농도") a <- dust$"기온" b <- dust$"재비산먼지평균농도" ( sum((a-mean(a))*(b-mean(b)))/(length(a)-1) ) / ( sd(a)*sd(b) ) # 상관관계계수 cor.test(a,b,method='pearson') ``` # 3/30 ``` # 잔차제곱합 == PRESS잔차 # PRESS/잔차의 표준편차 => 표준화잔차 # PRESS/잔차의 표준오차 => student 잔차 # 회귀분석의 전제조건 # 정규성:가설검정으로 구성(귀무가설은 기울기가 0이다. 대립가설은 기울기가 0이 아니다.) # 선형성:일반화(비선형->과적합) # 독립성:자기상관성, 다중공선성이 없어야 한다. # 등분산성:잔차가 어느지점에서나 분산이 일치해야 한다.=>잔차가 커기면=>예측이 불확실 # 모델을 평가 # 분산은 카이제곱분포를 따른다.(회귀제곱합이나 잔차제곱합은 분산과 의미가 같다.) # 모델이 유용한가:F통계량(평균비->F분포를 따른다) # 모델이 설명력이 있는가:R^2 # 계수가 유의미한가:계수의 t통계량(표준오차로 나눈 것) # 잘 fitting하고 있는가:잔차의 그래프로 확인 # 다향회귀 y <- c(1,2,3,8,5,3,10,7) x1 <- c(2,4,5,6,8,10,11,13) x2 <- c(1,2,2,4,4,4,6,6) opar=par(mfrow=c(1,3)) plot(x1,y) plot(x2,y) plot(x1,x2) # linear model summary(lm(y~x1)) summary(lm(y~x2)) summary(lm(y~x1+x2)) (Intercept) -0.1041 x1 -1.1566 x2 3.7267 # 기울기가 의미하는 것 -> 종속변수에 대한 영향력 (yhat=-1.1566 * x1 + 3.7267 * x2 + -0.1041) (PRESS <- (y-yhat)^2) # 잔차의 표준편차 : sd(PRESS) sd(PRESS) PRESS/sd(PRESS) # 표준화 잔차 # student 잔차: str(women) head(women) fit=lm(weight~height, data=women) summary(fit) (women.fitted_weight = -87.52+3.45*women$height) residuals <- women$weight - women.fitted_weight min(residuals) max(residuals) quantile(residuals,0.75) quantile(residuals,0.25) # IQR 값을 구해보시오 IQR <- quantile(residuals, 0.75) - quantile(residuals,0.25) plot(weight~height, data=women) abline(fit, col="blue") # 회귀분석의 정당성을 부여하는 것(상관계수) cor(women$weight, women$height) #0.995 -> 정당하다 # 모델이 통계적으로 유의미 한가: F검정통계량:1433, 13DF, p-value:1.091e-14: 유의미하다 # 계수의 유의미는 어떤가: t검정통계량:37.85, p-value:1.09e-14 : 유의미하다 # 모형의 설명력: Adjusted R-squared:0.9903:수정된 결정계수: 우수하다 # 데이터에 적용된 선형회귀는 적절한가: cor:1.0으로 적절하다다ㅡㅅㅊㅁ # cars데이터에 대하여 선형회귀를 하고 모델 평가를 해보시오 data(cars) # F통계량 MSR(회귀)/MSE(잔차) : 1.13 # 자유도로 나누어주어야 하는데 # SSR + SSE = SST head(cars) # 속도와 제동거리(제동거리를 예측) res.lm <- lm(dist~speed, data=cars) summary(res.lm) cor(cars$dist, cars$speed) #0.8068949 # 모델이 통계적으로 유의미 한가: F검정통계량:89.57, p-value:1.49e-12: 유의미하다 # 계수의 유의미는 어떤가: t검정통계량:9.464, p-value:1.49e-12 : 유의미하다 # 모형의 설명력: Adjusted R-squared:0.6438 # 데이터에 적용된 선형회귀는 적절한가: 0.8068949: 적절하다. yhat <- 17.5791 + 3.9324 * cars$speed plot(speed~dist, data=cars) abline(res.lm) names(res.lm) res.lm$soefficients (cof=coef(res.lm)) # 계수 (dist_pred = cof["(Intercept)"]+ cof["speed"]*cars$speed) plot(cars$speed, cars$dist, xlab="속도", ylab="정지거리",xlim=c(0,30),ylim=c(0,125)) abline(coef(res.lm), col=2) res.lm$fitted.values # 모델에서 사용한 데이터에 대하여 예측값 cars$dist - res.lm$fitted.values # 잔차 fitted(res.lm)[1:4] res.lm$residuals residuals(res.lm)[1:4] RSS <- sum((fitted(res.lm)-cars$dist)^2) # 잔차 제곱합 SSR,RSS(residuals sum of square) var(residuals(res.lm)) sd(residuals(res.lm)) (standard_res <- RSS/sd(residuals(res.lm))) # 표준화 잔차 (standard_res <- RSS/sd(residuals(res.lm))/sqrt(length(cars$dist))) # 표준오차 -> student 잔차 # 정규성, 독립성, 등분산성, 선형선 => 잔차를 통해서 확인 par(mfrow=c(1,1)) plot(res.lm) # Fitted values: 예측치(speed -> dist)와 Residuals과의 관계 # Normal Q-Q # Scale-Location: 분포의 일정함 정도를 확인 # Residuals(잔차) vs Leverage(원래 종속변수): 값의 차이가 있음 # Cook Distance -> 이상치의 존재여부를 확인하기 위해서 시각화 # Leverage는 종속변수가 이상치인 경우 # 이상치: 하나의 변수에서, 두개의 변수를 대상으로 하니까. # 이상치 종속변수, 이상치 독립변수 par(mfrow=c(2,2)) plot(res.lm) par(mfrow=c(1,1)) install.packages("lmtest") library(lmtest) bptest(cars$speed ~ cars$dist) # 등분산성을 테스트 # 귀무가설: 등분산이다. # 대립가설: 등분산이 아니다. # 자기상관성 dwtest(res.lm) # 2이면 자기 상관성이 없다. # 0~2 중 특히 0이면 높은 양의 자기 상관관계 # 4에 가까우면 음의 자기 상관관계 predict(res.lm, newdata=data.frame(speed=10)) #21.74499 제동거리(점추정) predict(res.lm, newdata=data.frame(speed=10), interval = "confidence") # (구간추정) # 문제 # 4.0, 25.0, .21일 때의 제동거리를 예측하시오 predict(res.lm, newdata=data.frame(speed=c(4.0, 25.0, .21)), interval = "confidence") # 문제 # MASS 패키지에 Boston 데이터에 대하여 회귀모델을 구축해보시오 library(MASS) data("Boston") str(Boston) ?Boston # medv가 종속변수, 13개 독립변수 -> 총 14개 # 정규화: 독립변수가 종속변수에 미치는 영향을 균등하게 하기 위해서 # 정규화를 안하면 가중치를 부여한 효과가 있다. # 독립변수의 척도가 크면 중요한 변수가 됨 library(lmtest) class(scale(Boston)) boston_df <- as.data.frame(scale(Boston)) head(boston_df,3) nrow(Boston) dim(Boston) # 506x14 # 분석을 할 때는 모델을 만드는 것과 테스트를 하는 데이터를 나눔 # 과적합이 벌어지는 이유: 데이터에 최적화 되는 것이 모델 # 학습 데이터/ 테스트 데이터 분리 set.seed(123) idx <- sample(1:nrow(boston_df),300) trainDF <- boston_df[idx,] testDF <- boston_df[-idx,] # 제외하고 dim(trainDF) # 300x14 dim(testDF) # 206x14 # dist: 범주형 데이터 제외 formula <- medv ~ crim+zn+indus+nox+rm+age+dis+rad+tax+ptratio+black+lstat lm_model <- lm(formula=formula, data=trainDF) lm_model lm_model$fitted.values[1:5] # 예측값 lm_model$residuals[1:5] # 잔차 summary(lm_model) # 예측 pred <- predict(lm_model, testDF) pred2 <- predict(lm_model, trainDF) # 회귀의 평가:연속형 데이터(비율, 구간) # 상관계수를 이용해서 평가 # 학습에 참여하지 않은 테스트로 평가 cor(pred, testDF$medv) #0.84 accuracy # 학습데이터로 평가 cor(pred2, trainDF$medv) #0.86 -> 학습데이터에 모델이 과적합되었다. # 모델은 일반화(어떤 데이터에도 강건하게 예측) dwtest(lm_model) # 귀무가설: 자기상관성이 없다. # 대립가설: 자기상관성이 있다. # p-value = 0.4636 귀무가설을 기각할 수 없으므로 자기상관성이 없다. # 다중공선성 검사 library(car) sqrt(vif(lm_model)) > 2 # nox, tax, dis, rad cor(trainDF[c("nox","tax","dis","rad")]) # tax:rad, nox:dis boston_cor <- cor(trainDF) library(corrplot) corrplot(boston_cor, method="square") # 제거 tax, nox formula <- medv ~ crim+zn+indus+nox+rm+age+dis+rad+tax+ptratio+black+lstat lm_model <- lm(formula = formula, data=trainDF) pred <- predict(lm_model, testDF) cor(pred, testDF$medv) # scale된 예측값을 원래의 값으로 복원 pred[1] * sd(Boston$medv) + mean(Boston$medv) ``` ``` wine <- read.csv("C:/Users/401-26/Documents/R/whitewines.csv", fileEncoding = "cp949", header=TRUE,) str(wine) # 와인의 질:0~10등급 unique(wine$quality) # 범주형 hist(wine$quality) idx <- 1:3700 training_wine <- wine_df[idx,] testing_wine <- wine_df[-idx,] dim(training_wine) dim(testing_wine) # 독립변수에 범주형 데이터가 있으면(3개 이상의 범주형 변수가 있으면 NA) # factor형으로 변경 lm_model <- lm(quality ~. , data=training_wine) pred2 <- predict(lm_model, testing_wine) pred3 <- predict(lm_model, training_wine) ``` ``` # mtcars 데이터에 대하여 선형회귀 모델을 생성하시오 str(mtcars) # 종속변수:mpg(연비) head(mtcars) sapply(mtcars, function(x) sum(is.na(x))) library(psych) pairs.panels(mtcars) mydata <- mtcars mydata$vs <- ifelse(mydata$vs == 0, "브이","직선") mydata$am <- ifelse(mydata$am == 0, "자동","수동") mydata$cyl <- as.factor(mydata$cyl) mydata$vs <- as.factor(mydata$vs) mydata$am <- as.factor(mydata$am) mydata$gear <- as.factor(mydata$gear) mydata$carb <- as.factor(mydata$carb) unique(mydata$cyl) # levels: 4,6,8 unique(mydata$vs) # 브이 직선 unique(mydata$am) # 수동 자동 unique(mydata$gear) # 4 3 5 unique(mydata$carb) # 4 1 2 3 6 8 범주형은 원형인코딩 colnames(mydata) # 1 2 3 4 5 6 7 8 9 10 11 12 13 # "mpg" "cyl" "disp" "hp" "drat" "wt" "qsec" "vs" "am" "gear" "carb" "hp_wt" "mpg_wt" # scale을 결정 mt_sc <- as.data.frame(scale(mydata[-c(2,8,9,10,11)])) mydata <- cbind(mt_sc, mydata[c(2,8,9,10,11)]) head(mydata) # 데이터 분할 idx <- sample(1:nrow(mydata),0.7*nrow(mydata)) train_data = mydata[idx,] test_data = mydata[-idx,] dim(train_data) dim(test_data) # formula를 생성 (vari <- colnames(mydata)) (form <- as.formula(paste(vari[1],paste(vari[-1],collapse="+"),sep="~"))) fit <- lm(form, data=mydata) summary(fit) coefficients(fit) residuals(fit) layout(matrix(c(1,2,3,4),2,2)) plot(fit) pred <- predict(fit, test_data[-1]) pred2 <- predict(fit, train_data[-1]) cor(pred, test_data$mpg) cor(pred2, train_data$mpg) # 과소적합(변수를 추가해서 더 학습) step(fit, direction = "both") # "backward" 후진선택법 or "forward" 전진선택법 # both 단계적 선택법 form <- as.formula("mpg ~ cyl+hp+wt+am") fit <- lm(form, data=mydata) pred <- predict(fit, test_data[-1]) pred2 <- predict(fit, train_data[-1]) cor(pred, test_data$mpg) cor(pred2, train_data$mpg) ``` # 3/31 ``` # dicision tree(의사결정트리) # binary tree # rule base로 작동 # 외삽데이터는 예측불가/회귀분석(기울기, 절편: 회귀선) # 과적합문제(root -> termina node까지 가장 깊은 곳까지를 depth) # 데이터변수와 입력순서를 변경하면 결과가 바뀌는 현상 # 단점을 극복하기 위해서 RandomForest(데이터는 랜덤하게 선택하고 모델의 숲을 만들어서 다수결의 원리) # 100개 중 50개는 동전의 앞면, 50개는 동전의 뒷면 -(50/100) * log2(50/100) # 0.5 (p(x)) library(rpart) (result <- sample(1:nrow(iris),nrow(iris)*0.7)) # 105 train <- iris[result,] test <- iris[-result,] dim(train);dim(test) unique(iris$Species) head(iris) # 분류문제 formula <- Species~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width model <- rpart(formula = formula, data=train) # 불순도를 이용해서 정보이득이 있으면 분할 model pred <- predict(model,test) # 범주별로 확률 -> labeling pred cpred <- ifelse(pred[,1] >= 0.5,"setosa",ifelse(pred[,2] >= 0.5,"versicolor","virginica")) (tb=table(cpred, test$Species)) # 정분류율 오분류율 # 정분류율 sum(diag(tb))/nrow(test) #93% 예측성공 plot(model) text(model, use.n=T, cex=0.6) post(model, use.n=T, file="") #install.packages("rpart.plot") #install.packages("rattle") library(rpart.plot) library('rattle') prp(model) rpart.plot(model) fancyRpartPlot(model) # 통계학에서는 분산이 큰 변수가 중요한 변수다. # ctree: 범주형 데이터를 제외하고 선형분석실행, 입력변수의 레벨이 31개로 제한, # 스스로 가지치기를 해 가지치기가 불필요 # C5.0: 불순도 척도가 Entropy, 범주형 데이터에 대하여 범주 수 만큼 분리가 발발생 # rpart: 지니(지니계수) 인덱스를 사용 #install.packages("party") library(party) result <- sample(1:nrow(iris),nrow(iris)*0.7) table(result) train <- iris[result,] test <- iris[-result,] head(train) formula <- Species ~ Sepal.Length+Sepal.Width+Petal.Length+Petal.Width iris_ctree <- ctree(formula, data=train) plot(iris_ctree, type="simple") plot(iris_ctree) names(train) result <- subset(train, Petal.Length > 1.9 & Petal.Width > 1.7) dim(result) table(result$Species) pred <- predict(iris_ctree, test) pred tr <- table(pred, test$Species) paste(round(sum(diag(tr)) / nrow(test)*100), '%') ``` ``` # 유방암 진단결과 데이터 셋 wdbc <- read.csv("C:/Users/401-26/Documents/R/wdbc_data.csv", fileEncoding = "cp949", header=TRUE,stringsAsFactors = F) str(wdbc) # 569x32 # diagnosis 진단: unique(wdbc$diagnosis) # 'B':양성,'M':악성 (30개 필드의 조합) # DT를 이용해서 모델을 구성해보시오(rpart) # 데이터: train:7, test:3 wdbc <- wdbc[-1] wdbc$diagnosis <- factor(wdbc$diagnosis, levels=c("B","M")) # DT: 데이터에 대해서 강건(범주형/연속형) # 분산이 크면 중요한 변수, 영향력을 동일 normalize <- function(x){ return((x-min(x))/ (max(x)-min(x))) } wdbc_x <- as.data.frame(lapply(wdbc[2:31],normalize)) summary(wdbc_x) # 7:3 wdbc_df <- data.frame(wdbc$diagnosis, wdbc_x) # 569x31 종속변수는 범주형, 독립변수는 정규화 dim(wdbc_df) idx=sample(nrow(wdbc_df), 0.7*nrow(wdbc_df)) wdbc_train = wdbc_df[idx,] wdbc_test = wdbc_df[-idx,] dim(wdbc_train) dim(wdbc_test) names(wdbc_train) # 모델학습 model2 <- rpart(wdbc.diagnosis ~., data=wdbc_train) pred2 <- predict(model2, wdbc_test, type="class") prp(model2) rpart.plot(model2) fancyRpartPlot(model2) # 평가 ( 평가매트릭-혼동행렬) (tr=table(pred2,wdbc_test$wdbc.diagnosis)) paste("정분류율",round(sum(diag(tr))/nrow(wdbc_test),3)*100,'%') paste("오분류율", round(sum(tr[1,2],tr[2,1])/nrow(wdbc_test),3)*100,'%') paste("정밀도",round(tr[1,1]/sum(tr[1,]),3)*100) paste("민감도", round(tr[1,1]/sum(tr[,1]),3)*100) paste("특이도", round(tr[2,2]/sum(tr[,2]),3)*100) # 자동화 패키지 인caret 패키지를 적용(mordern 모델학습) install.packages("dplyr") install.packages("caret") install.packages("recipes", type = 'binary') library(caret) set.seed(300) m <- train(wdbc.diagnosis ~., data=wdbc_train, method="C5.0") # entropy를 이용한 tree분석 p <- predict(m,wdbc_test) tr <- table(p, wdbc_test$wdbc.diagnosis) paste("정분류율",round(sum(diag(tr))/nrow(wdbc_test),3)*100,'%') paste("오분류율", round(sum(tr[1,2],tr[2,1])/nrow(wdbc_test),3)*100,'%') paste("정밀도",round(tr[1,1]/sum(tr[1,]),3)*100) paste("민감도", round(tr[1,1]/sum(tr[,1]),3)*100) paste("특이도", round(tr[2,2]/sum(tr[,2]),3)*100) head(predict(m,wdbc,type="raw")) head(predict(m,wdbc,type="prob")) # cross validation 교차 검증 ctrl <- trainControl(method="cv",number=10, selectionFunction="oneSE") # parameter 조합 expand.grid(height=seq(60,80,5), weight=seq(100,300,50), sex=c("Male","Female")) # 경우의 수를 생성 grid <- expand.grid(.model="tree", .trials=c(1,5,10,15,20,25,30,35), .window="FALSE") m <- train(wdbc.diagnosis ~., data=wdbc_train, method="C5.0", metric="Kappa",trControl=ctrl, tuneGrid=grid) m p <- predict(m,wdbc_test) table(p,wdbc_test$wdbc.diagnosis) m$bestTune m$results ``` # 4/3 ``` library(proxy) #install.packages('kernlab') library(kernlab) (x <- matrix(1:16, nrow=4)) # 4x4 sqrt(sum((x[1,]-x[4,])^2)) # 요소끼리 연산 sqrt(sum((x[,1]-x[4,])^2)) # distance (dist <- dist(x, method="euclidean")) # 두 개의 벡터간의 문제: 4x4에서 3x3(n-1,n-1 행렬) (dist <- dist(x, method="minkowski", p=2)) # 2승을 사용하면 euclidean 됨. (dist <- dist(x, method="binary")) # 2진수로 표현한 다음 (다른 비트의 수-거리) (dist <- dist(x, method="manhattan")) # 직선값대로 계산 #(dist <- dist(x, method="mahalanobis")) # 평균에서 표준편차 몇 배수의 거리값에 있는가 # 거리 행렬도 정방행렬이면서 대칭행렬 # 탐색적 군집분석(계층적 군집분석), 확인적 군집분석 # 군집 수의 계수를 결정하기 위해서 탐색적으로 실시하는 군집분석 idx = sample(1:dim(iris)[1],40) (irisSample = iris[idx,]) species <- irisSample$Species # 백업 irisSample$Species = NULL # 지도->비지도학습습 # dist는 행렬 거리값 (dist(irisSample)) hc = hclust(dist(irisSample),method="ave") # hierachical(계층적) plot(hc, hang=-1,labels=iris$Species[idx]) rect.hclust(hc,k=4) (groups = cutree(hc,k=3)) table(groups, species) clusters <- hclust(dist(iris[,3:4])) plot(clusters) (clusterCut <- cutree(clusters, 3)) # 범주화 table(clusterCut, iris$Species) # 핵심개념 # 거리값: 두개 간의 거리 -> 군집간에는 여러개의 데이터가 있는데 무엇을 기준 # single linkage: 두 군집의 가장 가까운 두 점 사이의 거리 # complete linkage: 두 군집을 가장 먼 두 점 사이의 거리 # average linkage: 각 거리의 평균값으로 # ward: 군집 내 편차들의 제곱합을 고려해서 거리값을 잼 clusters <- hclust(dist(iris[,3:4]), method='average') plot(clusters) library(ggplot2) (clusterCut <- cutree(clusters,3)) table(clusterCut, iris$Species) ggplot(iris, aes(Petal.Length, Petal.Width,color=Species))+geom_point(alpha=0.4, size=3.5)+geom_point(col=clusterCut)+scale_color_manual(values=c('black','red','green')) # kmeans -> 중심값 # 확인적 클러스터링 # 중심값과 k값 # 랜덤하게 중심값을 결정, k값을 결정 -> 각 중심과 데이터의 거리값을 계산 -> 그룹을 지정 # -> 중심값을 재계산 -> 각 중심과 거리값을 계산 # wss가 작아질 때까지(그룹내 응집도 최대화, 그룹간 분리도 최대화) set.seed(100) # random으로 결정되는 값을 일정하게 해서 결과가 동일일 iris_back <- iris iris_back$Species <- NULL # 비지도 학습 # 거리를 기준으로 할 때는 반드시 정규화해야 함 iris.scaled <- scale(iris_back[,-5]) # scale을 하지 않는다면 변수의 사이즈 크기가 거리값 (kmeans.result <- kmeans(iris.scaled,3,trace = T)) table(iris$Species, kmeans.result$cluster) # 산점도에 컬러값을 줌 # 4개의 변수 중에 2개의 변수만 고려해서 시각화 plot(iris.scaled[,c("Sepal.Length","Sepal.Width")],xlim=c(-2,2),ylim=c(-2,2),col=kmeans.result$cluster) points(kmeans.result$centers[,c("Sepal.Length","Sepal.Width")],col=1:3,pch=10,cex=10) iris_back$Species <- iris$Species iris_back$Species <- as.factor(iris_back$Species) iris_back$Species <- kmeans.result$cluster # 종속변수 #sum(iris_back$Species == iris_back$pred) names(kmeans.result) kmeans.result$cluster kmeans.result$centers # 클러스터의 중심값 kmeans.result$withinss # 그룹내 거리의 합 kmeans.result$betweenss # 그룹간 거리의 합 kmeans.result$tot.withinss # 모든 그룹내 거리의 합 sum(kmeans.result$withinss) # 3개의 그룹 합 = kmeans.result$betweenss + kmeans.result$tot.withinss kmeans.result$totss kmeans.result$iter # 2번회전 kmeans.result$size km.res <- kmeans(iris.scaled, 3, nstart=10) library(factoextra) # 4가지 요소를 고려한 시각화 fviz_cluster(km.res, iris[,-5],ellipse.type="norm") # 클러스터 개수 구하기 library(caret) # 자동화 패키지 set.seed(123) # sampling 대신에 createDataPartition으로 데이터 분할 inTrain <- createDataPartition(y=iris$Species, p=0.7, list=F) training <- iris[inTrain,] testing <- iris[-inTrain,] training.data <- scale(training[-5]) summary(training.data) # 중심값을 재계산하면서서 iris.kmeans <- kmeans(training.data[,-5],centers = 3,iter.max = 10000) training$cluster <- as.factor(iris.kmeans$cluster) # 그룹번호 qplot(Petal.Width, Petal.length, color=cluster,data=training) table(training$Species, training$cluster) #install.packages('NbClust') library(NbClust) # 그룹을 조절하면서 kmeans -> 판단근거는 거리값의 합 nc <- NbClust(training.data, min.nc = 2, max.nc = 15, method='kmeans') par(mfrow=c(1,1)) barplot(table(nc$Best.n[1,]), xlab='클러스터수',ylab='기준',main='클러스터 수 확인') # wss: within sum of squares dim(training.data) wssplot <- function(data, nc=15, seed=1234){ wss <- (nrow(data)-1)*sum(apply(data,2,var)) # 분산 for(i in 2:nc){ # 클러스터 개수 set.seed(seed) wss[i] <- sum(kmeans(data,centers = i)$withinss) # 그룹내 거리합 kmeans(data,centers=i)$withinss } plot(1:nc, wss, type='b', xlab='클러스터수',ylab='SS') } wssplot wssplot(training.data) # 중1 신체검사결과 -> 군집분석을 라벨링 해보시오 body <- read.csv('C:/Users/401-26/Documents/R/bodycheck.csv', header=T, fileEncoding='cp949') str(body) names(body) unique(body$안경유무) # 범주형 -> 숫자 df<-body[,-1] df <- na.omit(df) # 열이름 변경 colnames(df) <- c('handpw','height','weight','galss') df <- scale(df) class(df) df <- as.data.frame(df) k2 <- kmeans(df, centers = 2, nstart=25) # 중심점 선택: random한 set을 25번 해서 좋은 결과 사용해라 library(factoextra) fviz_cluster(k2,data=df) df$cluster <- as.factor(k2$cluster) plot(df$handpw,df$weight,col=df$cluster) points(k2$centers[,c('handpw','weight')],col=c(1,2),pch=8,cex=5) k2$centers ``` ``` # 문제 # diamonds 데이터에 대해 price, carat, depth, table 변수 데이터만을 대상으로 # 탐색적 및 확인적 군집분석을 실시하시오 # 전체 53940데이터 중 1000개만을 샘플링해서 작업하시오 # 조건1 탐색적 군집분석을 실시하고 군집수를 결정하시오 # 조건2 가격에 영향을 미치는 요인에 대해서 설명해보시오 # caret과 price의 관계를 군집으로 구분해서 시각화하시오 str(diamonds) (t <- sample(1:nrow(diamonds),1000)) length(t) # 1차원 test <- diamonds[t,]; dim(test) # 1000x10 mydia <- test[c("price","carat","depth","table")] dim(mydia) # 1000x4 # hclust 거리행렬이 입력 class(dist(mydia)) dim(dist(mydia)) # 1000x1000 기본데이터 타입이 아님 # 그룹과 그룹의 거리: single, complete(가장 먼거리), average, ward # 탐색적 군집분석 result <- hclust(dist(mydia),method="average") plot(result, hang=-1) fviz_dend(res, rect=T, cex=0.5, k_colors=c("#00AFBB","#2E9FDF","#E7B800","#FC4E07")) (res <- hcut(mydia, k=4, stand=T)) mydia <- scale(mydia) str(mydia) class(mydia) mydia <- as.data.frame(mydia) result2 <- kmeans(mydia,4) # 군집분석: 행과의 거리 result2 result2$centers unique(result2$cluster) length(result2$cluster) # 모든 거리의 합 = 군집내 거리+군집간 거리 result2$withinss # 군집내 거리 result2$tot.withinss # 군집내 거리 다 더한거 -> WSS(within sum of squares) result2$betweenss # 군집간 거리 # 원본데이터 + 라벨을 결정 => 예측문제가 분류문제로 바뀜 mydia$cluster <- result$cluster head(mydia) plot(mydia[,-5]) cor(mydia[,-5],method="pearson") library(corrgram) corrgram(mydia[,-5]) plot(mydia$carat, mydia$price) #plot(mydia$carat, mydia$price, col=mydia$cluster) points(result2$centers[,c("carat","price")],col=c(3,1,2),pch=8,cex=5) ``` ``` # 청소년 성향에 대한 설문 조사한 데이터 # sns 계정을 가진 4개의 기본정보와 36개의 관심분야 # 청소년 성장 코미디: # 범죄성향, 운동성향, 외모지향, 무기력, 브레임 teens <- read.csv('C:/Users/401-26/Documents/R/snsdata.csv', header=T, fileEncoding='cp949') # 30000x40 종속변수가 없음:비지도학습 set.seed(100) str(teens) head(teens) # gender age 결측치 발생 sum(is.na(teens$gradyear)) sum(is.na(teens$freinds)) sum(is.na(teens)) colSums(apply(teens, 2, is.na)) table(teens$gender) summary(teens$age) teens$age <- ifelse(teens$age >= 13 & teens$age < 20, teens$age, NA) summary(teens$age) mean(teens$age, na.rm=T) # 졸업년도 별로 나이의 평균 (avg <- ave(teens$age, teens$gradyear, FUN=function(x) mean(x, na.rm=T))) teens$age <- ifelse(is.na(teens$age), avg, teens$age) # 3항연산자 teens$age[1:6] summary(teens$age) interests <- teens[5:40] # 1,2,3,4는 기본정보라서 제외 summary(interests) interests_n <- data.frame(lapply(interests, scale)) summary(interests_n) teen_clusters <- kmeans(interests_n, 5) teen_clusters$size table(teen_clusters$cluster) # 그룹의 명명식 => 클러스터를 대표하는 것이 중심값 teen_clusters$centers # 중심값을 그룹성향을 네이밍한다. teens$cluster <- teen_clusters$cluster # 클러스터 종속변수 aggregate(data=teens, age~cluster,mean) # 클러스터 별 나이의 평균을 구하시오 qplot(cluster, age,color=gender,data=teens) res <- aggregate(data=teens,gender=='F'~cluster, mean) plot(res, type='b',lty=2) aggregate(data=teens,softball+volleyball+hair+dress ~gender=='F',mean) ``` ``` # DTW(dynamin tim wraping) # 연속형 수치를 비교 (신호): 지연, 노이즈 install.packages("dtwclust") library(dtwclust) data(uciCT) str(CharTraj) class(CharTraj)# 열별로 데이터사이즈가 상이 # dim(CharTraj) # 선형보간법: 사이즈를 일치시키기 위해 없는 부분은 보간법으로 채워서 계산 series <- reinterpolate(CharTraj, new.length = max(lengths(CharTraj))) (series <- zscore(series)) # NaN => 0으로 채워가면서 pc.dtwlb <- tsclust(series,k=20L, distance="dtw_lb",centroid="pam", seed=3247, trace=T, # 출력 여부 control=partitional_control(pam.precompute=F), # 미리계산 args=tsclust_args(dist=list(window.size=20L))) plot(pc.dtwlb) ``` # 4/4 ``` # 연관분석(Association analysis): 장바구니 분석 #install.packages("arules") library(arules) # apriori를 구현한 패키지 data("Adult") # transactions 거래행렬(범주형데이터에 대해서 행렬로 구성) # ngCMatrix: sparse matrix(희소행렬): text mining # 마트의 품목 많음 -> 장바구니(희소) # 희소행렬 표현하는 방법(안접행렬(공간낭비), 그래프 표현(메모리는 절약->처리시간은 증가)) # @가 나오면 S4 class # S3는 함수를 이용해서 클래스를 생성(R은 함수(1급객체)를 이용하는 프로그래밍) str(Adult) rules <- apriori(Adult) # 6137개의 룰 발생 rules.sub <- subset(rules, subset=rhs %in% "marital-status=Never-married" & lift >2) rules.sub # 58개 룰로 바뀜 inspect(rules.sub) rules.sub <- subset(rules, subset=lhs %ain% c('age=Young','workclass=Private')) # 76개 룰 inspect(rules.sub) # 중요한 rule을 생성하고자 할 때 parameter를 이용한 제약(support, confidence를 이용한 제약) rules <- apriori(Adult, parameter= list(supp=0.5, conf=0.8, minlen=2)) # 80 rules #install.packages('arulesViz') library(arulesViz) plot(rules, method='graph', engine='interactive') # in(섞어도 됨),ain(전체가),pin(partial) white.lhs <- subset(rules, lhs %in% 'race=White') length(white.lhs) inspect(white.lhs) plot(white.lhs, method='graph') # 문제 # race=White이고 native-country=United-State인 것이 lhs에 나타나는 rule 조합을 선택해보시오 whiteORus.lhs <- subset(rules, lhs %in% 'race=White', lhs %in% 'native-country=United-States') inspect(whiteORus.lhs) # 24 # rules에 대해 support와 confidence를 이용해서 내림차순 정렬한 다음 상위 3개만 출력해보시오 inspect(head(rules)) ?sort # 오버로딩 inspect(head(sort(rules, decreasing = T, by= c('support','confidence')),3)) adult <- as(Adult, 'data.frame') str(adult) head(adult) ar <- apriori(Adult, parameter=list(supp=0.1, conf=0.8)) ar2 <- apriori(Adult, parameter=list(supp=0.2, conf=0.95)) ar6 <- apriori(Adult, parameter=list(supp=0.5, conf=0.80)) inspect(head(ar6)) ?apriori rules <- apriori(Adult, parameter=list(support=0.3, confidence=0.50), appearance = list(rhs='relationship=Husband'),conrtol=list(sort=-1)) inspect(rules) # 문제 # 1. 최소support=0.5 최소confidence=0.8을 지정하여 연관규칙 생성 Adult.ap <- apriori(Adult,list(support=0.5,confidence=0.8)) # 2. 연관분석결과를 연관어 네트워크 형태로 시각화 plot(Adult.ap, method='graph') # 'interactive' # 3. lhs가 white(백인)인 규칙만 subset으로 작성하고 연관어 시각화 Adult.lhs <- subset(Adult.ap, subset=lhs %in% 'race=White') inspect(head(sort(Adult.lhs, decreasing=T, by='lift'),10)) plot(Adult.lhs) plot(Adult.lhs, method='graph') data('Groceries') str(Groceries) summary(Groceries) size(head(Groceries)) LIST(head(Groceries,3)) # list # 연관분석인 경우 itemFrequencyPlot(Groceries, topN=15) itemFrequency(Groceries[,1:5]) a_matrix <- matrix(c( 1,1,0,0,0, 0,0,1,0,0, 0,0,0,1,1, 1,1,1,0,0, 0,1,0,0,0 ), ncol=5, byrow=T) dimnames(a_matrix) <- list(paste('Tr',c(1:5),sep=''), c('milk','bread','butter','beer','diapers')) a_matrix image(a_matrix) trans2 <- as(a_matrix, 'transactions') # 희소행렬처리 포맷으로 변환 trans2 rule <- apriori(trans2, parameter=list(supp=0.01, conf=0.01)) # 16개 룰 inspect(rule) ``` 팀과제 ``` # 수작업(지지도, 신뢰도, 향상도 구해라): 맥주->고기 / 라면,우유->맥주 getwd() data <- read.table('./R/tran.txt') data # 맥주->고기 지지도=1/6 신뢰도=1/2 향상도=3/4 # 라면,우유->맥주 지지도=1/6 신뢰도=1/3 향상도=1 # 데이터를 데이터베이스에 입력하고 로딩하여 처리하시오(연관분석:데이터베이스) library(dplyr) library(dbplyr) library(DBI) conn <- DBI::dbConnect(RMySQL::MySQL(), dbname='tran', host='192.168.41.183', user='root', password='acorn1234') dbListTables(conn) a <- dbGetQuery(conn, 'SELECT * FROM tran') class(a) a # 연관분석을 통해 적당한 지지도와 신뢰도를 결정하시오 # 결과를 해석하고 향상도를 이용해서 정렬한 다음 출력하시오 a <- as.matrix(a) image(a) trans2 <- as(a, "transactions") trans2 # db_df <- as.data.frame(a) class(db_df) dbWriteTable(conn, "db_df", db_df) db_df rule2 <- apriori(trans2, parameter = list(supp=0.1, conf=0.1)) inspect(rule2) plot(rule2, method="graph") ``` 선생님 풀이 ``` ######로딩 filepath <- 'tran.txt' (data <- read.table(filepath, header=F, sep='\t')) class(data) length(data) # 데이터가 길이가 다른 경우(list) datalist <- list() for(i in 1:nrow(data)){ vector <- unlist(strsplit(data[i,1],',')) datalist[i] <- list(vector) } datalist # 범주 (items <- unique(unlist(datalist))) # "라면""맥주""우유""고기""과일" # 6x5행렬 matrix matrix <- matrix(0, nrow=nrow(data),ncol=length(items)) dim(matrix) datalist[2] # 확인 for(i in 1:nrow(data)){ for(datai in datalist[i]){ matrix[i, match(datai, items)] <- 1 } } df <- data.frame(matrix, row.names=NULL) colnames(df) <- items print(df) library(dplyr) library(dbplyr) library(DBI) conn <- DBI::dbConnect(RMySQL::MySQL(), dbname='tran', host='192.168.41.183', user='root', password='acorn1234') dbListTables(conn) for(i in 1:nrow(df)){ values <- paste0("(",paste(df[i,],collapse=", "),")") query <- paste0("INSERT INTO tran VALUES ", values) dbExecute(conn,query) } res <- dbSendQuery(conn,"SELECT * FROM transacton_table") df <- dbFetch(res) df df_tran <- as(df, "transactions") rule <- apriori(df_tran, parameter = list(supp=0.6, conf=0.3, minlen=2, maxlen=3)) inspect(rule) plot(rule, method="graph") ``` ``` # data("AdultUCI") data("AdultUCI") str(AdultUCI) # not logical or factor -> 너무 많은 조합이 발생 basket_rules <- apriori(AdultUCI, parameter=list(sup=0.08, conf=0.5, target="rules")) count(unique(AdultUCI['age'])) #74 min(AdultUCI['age']) #17 추가되는 데이터에도 적용가능해야 함 max(AdultUCI['age']) #90 # c(15,25,45,65,100) # ['age']:열선택, [['age']]:요소 AdultUCI[['age']] <- ordered(cut(AdultUCI[['age']],c(15,25,45,65,100)), labels=c('young','middle','senior','old')) AdultUCI["fnlwgt"] <-NULL AdultUCI["education-num"] <- NULL AdultUCI[["capital-gain"]] <- ordered(cut(AdultUCI[["capital-gain"]], c(-Inf, 0, median(AdultUCI[["capital-gain"]][AdultUCI[["capital-gain"]]>0]),Inf)), labels = c("None", "Low", "High")) AdultUCI[["capital-loss"]] <- ordered(cut(AdultUCI[["capital-loss"]], c(-Inf, 0, median(AdultUCI[["capital-loss"]][AdultUCI[["capital-loss"]]>0]),Inf)), labels = c("None", "Low", "High")) Adult_new <- as(AdultUCI,"transactions") basket_rules <- apriori(Adult_new, parameter=list(sup=0.08, conf=0.5, target='rules')) inspect(basket_rules[1:10]) p <- inspectDT(basket_rules) htmlwidgets::saveWidget(p,"arules_2.html",selfcontained = F) browseURL('arules_2.html') ``` ``` # 문제 # 1. 가족관계 및 교육수준의 소득과의 연관성을 확인하시오 str(AdultUCI) summary(AdultUCI) rules <- apriori(cbind(AdultUCI["relationship"], AdultUCI["education"], AdultUCI["income"]), parameter = list(supp=0.01, conf=0.01)) rules = subset(rules, rhs %pin% 'income') inspect(rules) # 2. 주당 일하는 시간과 소득과의 관계를 확인해보시오 str(AdultUCI) summary(AdultUCI) rules2 <- apriori(AdultUCI[c("hours-per-week", "income")], parameter = list(supp=0.01, conf=0.01)) rules2 = subset(rules2, rhs %pin% 'income') inspect(rules2) # 3. 기타 위의 데이터로부터 자기가 주장하고자 하는 내용을 확인하고 의견을 제시하시오 str(AdultUCI) rules3 <- apriori(AdultUCI[c("age", "hours-per-week")], parameter = list(supp=0.01, conf=0.01)) rules3 = subset(rules3, rhs %pin% 'hours-per-week') inspect(rules3) ``` # 4/5 ``` install.packages('tseries') # 시계열 분석 및 계산 재정학(finance)에 대한 함수들 install.packages('forecast') # 선형이나 시계열을 위한 예측 함수들 install.packages('TTR') # Technical Tranding Rule 기술적 트레이닝 룰(주식) # NA(moving average) 이동평균법(5,12,20,60,120) # SMA(simple) 단순이동평균 # WMA(weighted) 가중이동평균(1,2,3,4,5) 가중치(영향력) # EMA(exponentially) 지수이동평균법: 이자-복리계산 library(tseries) library(forecast) library(TTR) ts(1:10, frequency = 4,start=c(1959,2)) # 2사분위수수 print(ts(1:10,frequency = 7,start=c(12,2)),calendar=T) cumsum(1+round(rnorm(100),2)) gnp <- ts(cumsum(1+round(rnorm(100),2)),start=c(1954,7),frequency = 12) class(gnp) # 'ts' 시계열 데이터 포맷 plot(gnp) # kings <- scan('./R/kings.dat') kings kingstimeseries <- ts(kings) kingstimeseries plot(kingstimeseries) # SMA(simple moving average) king_sma3 <- SMA(kingstimeseries, n=3) plot.ts(king_sma3) king_sma8 <- SMA(kingstimeseries, n=8) plot.ts(king_sma3) # kpss.test: Kwiatkowski-Phillips-Schmidt-Shin # 귀무가설: 정상성이다. # 대립가설: 비정상성이다. kpss.test(king_sma8) # Stationarity 안정성테스트: 비정상성 # 차분 king_diff <- diff(kingstimeseries,differences=1) plot.ts(king_diff) kpss.test(king_diff) # 정상성을 가진 데이터 acf(king_diff, lag.max = 20) # 첫번째 lag에서 초과 acf(king_diff, lag.max=20, plot=F) # partial auto corelation function pacf(king_diff, lag.max=20) # lag1,lag2,lag3에 초과 pacf(king_diff, lag.max=20, plot=F) # (p,d,q): AR, 차분회수, MA계수 auto.arima(kings)# 추천: ARIMA(0,1,1) 비정상성을 가진 데이터에도 적용 kingstimeseriesarima <- arima(kingstimeseries, order=c(0,1,1)) # 예측 모델 kingstimeseriesarima # h: Number of periods for forecasting kingstimeseriesforecasts <- forecast(kingstimeseriesarima, h=5) # 예측 함수(predict 아님) kingstimeseriesforecasts plot(kingstimeseriesforecasts) # 문제 x=w=rnorm(300) for(t in 2:300) x[t] = x[t-1]+w[t] # 위의 데이터에 대해 시계열 데이터로 변환해보고 데이터가 정상성을 띄고 있는지 확인해보시오 (x.ts = ts(x)) kpss.test(x.ts) # p-value=0.01: 비정상성 kpss.test(diff(x.ts)) # 정상성 tsdisplay(diff(x.ts)) # 차분한 데이터 # ACF, PACF 데이터의 독립성(시차 데이터가 서로 영향을 미치는가) # 점선 안에 있으면 서로 독립적이다. # 점선을 벗어나면 독립적이지 않은 경우가 있다. # ACF 치수 11에서 독립적이지 않다. # PACF: 10,11이 독립적이지 않다. ## 정상성 확인 library(tseries) class(AirPassengers) (apts <- ts(AirPassengers, frequency = 12)) plot(apts) # adf.test 귀무가설: 비안정성, 대립가설: stationary 정상성 adf.test(apts, alternative = 'stationary',k=0) # 안정성 adf.test(diff(log(AirPassengers)),alternative = 'stationary',k=0) kpss.test(apts) auto.arima(AirPassengers) # ARIMA(2,1,1)(0,1,0) # sesonal 계절적 변동성을 가진 데이터 (fit <- arima(log(AirPassengers),c(2,1,1),seasonal=list(order=c(0,1,0),period=12))) pred <- predict(fit, n.ahead=10*12) # 데이터에 로그를 취해 역으로 exp를 취해줘야 복원 # 지수 exp: e값 지수가 취하는 일정비율(2.718) - 복리계산 ts.plot(AirPassengers, 2.718^pred$pred, log='y', lty=c(1,3)) par(mfrow=c(1,1)) air.model <- Arima(window(AirPassengers, end=1956+11/12),order=c(2,1,1), seasonal=list(order=c(0,1,0),period=12),lambda=0) plot(forecast(air.model, h=48)) lines(AirPassengers) # 문제 # WWWusage의 60기 까지만의 데이터를 분석하고 향후 10기에 대해 예측해보시오 # (분석을 해보니 log가 상황이 더 좋아서 변환해 봄) # 비정상성 -> 정상성으로 변환 ?WWWusage class(WWWusage) # 분당 동접수 WWWusage length(WWWusage) #100 plot(WWWusage) auto.arima(WWWusage) # kpss.test: Kwiatkowski-Phillips-Schmidt-Shin # 귀무가설: 정상성이다. # 대립가설: 비정상성이다. kpss.test(WWWusage) #0.05377: stationarity -> 정상성 # acf, pacf 자기상관성 acf(WWWusage) pacf(WWWusage) kpss.test(log(WWWusage)) plot(log(WWWusage)) acf(log(WWWusage)) pacf(log(WWWusage)) auto.arima(log(WWWusage)) log(WWWusage) %>% Arima(order=c(1,1,1)) %>% forecast(h=20) %>% autoplot data(gas) str(gas) plot(gas) plot.ts(gas) # 분해 시계열 ?stl # Seasonal Decomposition of Time Series by Loess 비선형회귀에 의한 plot(stl(gas,s.window = 'periodic')) kpss.test(gas) # p-value = 0.01 -> 비정상성 acf(gas) pacf(gas) gas.diff <- diff(log(gas)) # 로그.차분 => 정상성을 가진 데이터로 변환 adf.test(gas.diff, alternative = 'stationary',k=0) # stationary kpss.test(gas.diff) # 안정적 acf(gas.diff) pacf(gas.diff) plot(gas.diff) gas.ari <- auto.arima(gas.diff) summary(gas.ari) # ARIMA(0,0,1)(0,1,1) fit <- arima(window(log(gas),end=1980),c(0,0,1), seasonal = list(order=c(0,1,1),period=12)) forres <- forecast(fit, h=100) forres$method forres$residuals # 실제값-예측한 값의 차이 library(scales) library(ggplot2) autoplot(forres$residuals)+xlab('날짜')+ylab('잔차')+ggtitle('모델잔차') ggAcf(forres$residuals)+ggtitle('잔차의 ACF') # 자기상관성 # 귀무가설: 데이터들이 독립적으로 분포(바기 상관성이 0) # 대립가설은 자기상관성이 0이 아니다. Box.test(forres$residuals, lag=10, fitdf=0) # 자기 상관성이 0이 아니다. Box.test(forres$residuals, lag=2, fitdf=0) Box.test(forres$residuals, lag=7, fitdf=0) plot(forecast(fit, h=100)) plot(2.718^forres$fitted) # text mining library(stringr) # 선택, 갯수 # str_extract 처음 패턴만 분리 str_extract('abcd12aaa33','[0-9]{2}') # 정규 표현식(패턴비교) str_extract_all('abcd12aaa33','[0-9]{2}') # if 비교문 사용하는 것이 아닌 패턴을 이용해서 선택 str <- 'hongkildong35lee45kang55안창호25' (result <- str_extract_all(str,'[a-z]{3}')) (result <- str_extract_all(str,'[0-9]{2}')) (result <- str_extract_all(str,'[가-히]{3}')) nchar(str) length(str) head(USArrests) (states=rownames(USArrests)) grep(pattern='k',x=states) grep(pattern='k',x=states, value=T, ignore.case = T) grep(pattern='[ww]',x=states, value=T) money='$money' sub(pattern = '\\$', replacement = '', x=money) gsub('\\d','_',"the dandelion war 2010") # d는 숫자 D는 숫자가 아닌 것 gsub('\\w','_',"the dandelion war 2010") # w는 문자 W는 문자가 아닌 것 gsub(pattern='[[:blank:]]',replacement = '','the dandelion war 2010') data <- c("12.578.69486") # 위 데이터에 .을 없애보시오 gsub(pattern='\\D','',data) gsub('\\.','',data) lorem='대한민국 만세, 모두다 하나 ' str_sub(lorem,start=1,end=5) bad_text=c('This','example ','has serveral ') str_trim(bad_text, side='both') # 공백 제거 str_pad('hola', width=7) str_extract_all('abcd12aaa33',"\\d") #"[1-9]{2} chr6 = '123한글나라456' str_extract(chr6, '[가-흐]+') # 한글만 출력 chr6 ="abdbsdsa 010-1234-5678 dasdasd" pat <- '[0-9]{3}-[0-9]{4}-[0-9]{4}' str_extract(chr6, pat) fruits1 <- c(" 사과 한개","배 두개","바나나 세개 ") str_replace(fruits1,"[개]"," 바구니") shoppinglist <- c("2017-10-19 수입3000원","2017-10-20 수입14500원","2017-10-21 수입2500원") str_extract_all(shoppinglist,"[0-9]{4}원",simplify=T) str_replace_all(shoppinglist,"-","/") str_extract_all(shoppinglist,'[0-9]{4}[-][0-9]{2}[-][0-9]{2}',simplify=T) #install.packages("tm") library(tm) data('crude') # 역문서 빈도 tfidf=function(mat){ tf <- mat id=function(col){sum(!col==0)} idf <- log(nrow(mat)/apply(mat,2,id)) tfidf <- mat for(word in names(idf)){tfidf[,word] <- tf[,word]*idf[word]} return(tfidf) } class(crude) # "VCorpus" "Corpus" ngCMatrix 처럼 tm의 고유포맷(document 20개) inspect(crude[1:3]) (tdm <- TermDocumentMatrix(crude)[1:10,1:20]) # 행이 단어, 열이 문서 # TermDocumentMatrix: 행이 문서, 열이 단어: Sparsity 희소성:94% -> 0(데이터 없다) Docs(tdm) Terms(tdm) tdm <- TermDocumentMatrix(crude, control = list(removePunctuation=T, removeNumbers=T, stopwords=T, wordLengths=c(2,Inf))) findFreqTerms(tdm,lowfreq = 10) findFreqTerms(tdm,2,3) findAssocs(tdm,'dlrs',0.6) findMostFreqTerms(tdm) (termFreauency <- rowSums(as.matrix(tdm))) # 각 단어가 등장하는 회수를 카운팅 (termFreauency <- colSums(as.matrix(tdm))) # 각 문서에 들어있는 단어와 등장횟수 (termFreauency <- subset(termFreauency, termFreauency>=10)) #install.packages('wordcloud') library(wordcloud) m <- as.matrix(tdm) wordFreq <- sort(rowSums(m),decreasing = T) set.seed(375) pal <- brewer.pal(8,"Dark2") wordcloud(words=names(wordFreq),freq=wordFreq, min.freq = 10, random.order = F, rot.per = .1, colors=pal) ```