음 시계열을 요청하신 분이 계셨는데... 정말 안타깝게도 제가 아직 통계학과 3학년이라 다음학기에 들을 예정이라서요...
아는데까지만 알려드릴게요 죄송합니다 ㅠㅠ
오늘은 행렬로 문을 열어보겠습니다.
1. 행렬 만드는 법.
R로 행렬을 만들어 봅시다. 간단한 재료와 우주의 도움만 있으면 됩니다.
행렬을 만들때는 반드시 행렬의 이름을 지정해야 합니다. A라는 이름으로 저장하겠습니다.
행렬은 2,-5,4 1,-2,1 1,-4,6 으로 3x3 행렬입니다.
행렬을 만드셨다면 이걸 써먹을 줄 아셔야겠죠. 누구나 쉽게 볼 수 있는 방정식을 풀어보겠습니다.
2. 행렬을 이용한 방정식 풀기
위에서 만든 행렬식은 2X-5Y+4Z, X-2Y+Z, X-4Y+6Z 로 바꿔서 볼 수 있습니다.
이걸 이용해서 식을 세워보면
2X-5Y+4Z = -3
X-2Y+Z = 5
X-4Y+6Z = 10 로 가정하고 풀어보겠습니다
어때요 참 쉽죠? 참고로 푸는 방법은 Ax=B, 따라서 x=solve(A)[A의 역행렬]*B 입니다.
%*%라는 수식은 행렬끼리의 곱셈을 나타낼때 쓰는 방법입니다. 저거 안하면 곱셈 안되요.
3. 데이터로 행렬 만들어오기
데이터를 불러오는 작업은 1번에서 했었죠. 이번엔 그 데이터를 행렬로 재구성하는 방법을 해볼껍니다. 이걸 왜하냐구요? 4번에 연관이 되거든요.
1번에서 불러온 P060.txt 데이터를 a로 저장하고 a1이란 이름으로 X1 ~ X6까지 뽑아보겠습니다.
자 위에서 보시면 a에 P060.txt를 불러왔구요 a<-read.table("P060.txt",header=TRUE)
a1이라는 이름으로 a의 2~7까지를 matrix 형식으로 저장했습니다. a1<-as.matrix(a[,2:7])
보시면 아시겠지만 위의 명령어 as.matrix(a[,2:7])은 matrix(행렬)로 a의 2열부터 7열까지 행 구분 없이 전부 뽑아온다. 라는 뜻입니다.
역시 참 쉽죠?
4. 행렬 붙이기
뽑아온 이유가 궁금하시다구요? 간단합니다. X행렬을 만들기 위해서죠.
회귀분석은 단순히 lm 함수만 쓰면 되지만 내가 세운 회귀식이 정말 맞는지를 확인하기 위해서는 여러가지 절차가 필요합니다.
그리고 그 기준이 되는 행렬이 바로 X행렬이죠.
첫 열이 1, 두번째 열부터 X값이 되어 1, X1, X2 ... X6의 순서로 배열되어야 합니다.
이를 위해서 X1 ~ X6까지를 뽑고 앞에 1로 만들어진 nx1행렬을 붙이는 방법입니다.
a2라는 행렬을 만듭니다. 이 행렬은 1을 30개(a1의 열의 갯수) 가지고 있습니다.
그리고 이를 cbind(a2,a1)이라는 함수를 통해 열을 늘리는 방법으로 합칩니다. 참고로 a2가 뒤로 가면 순서가 바뀌니 주의하세요.
위의 사진은 외전인데 rbind라는 명령어로 행을 늘리는 방법입니다. 물론 거의 안쓰니 그냥 알아두시면 좋아요
쉽죠?
행렬 나누기는 제가 위의 3번에서 보여드렸으니 생략하겠습니다.
이번 시간은 여기까지구요. 다음 시간에는 여러 확률분포들을 들고 돌아오겠습니다.
p.s) 참고로 혹시나 R 코드에 관심 있으신 분들을 위해 제공하는 행렬로 방정식 푸는 코드
방정식은 2번에 써있는 그대로입니다.
A <- matrix(c(2,-5,4,1,-2,1,1,-4,6),byrow=T,nrow=3,ncol=3) #A행렬입니다.
b <- matrix(c(-3,5,10),nrow=3,ncol=1) #B행렬입니다.
p <- nrow(A) #p라는 변수를 잡았습니다. A의 열을 나타내는 숫자인데요 자세한건 뒤에...
eps <- 0.001 #작다고 판명할 수 있는 수 입니다. 이것보다 작으면 0으로 취급할겁니다.
U.pls <- cbind(A,b) #위에서 배우셨던 행렬 합치기입니다. A를 먼저 그리고 B를 뒤에 열로 붙였습니다.
for (i in 1:p){ #i라는 변수를 1에서 p까지 적용합니다. p는 위에 보시면 A의 열입니다.
if((abs(U.pls[i,i]) < eps) & (i < p)){
temp.v <- U.pls[i+1,]
U.pls[i+1,] <- U.pls[i,]
U.pls[i,] <- temp.v
}
#U.pls(A와 B를 합친 행렬)의 i,i(대각 행렬) 이 0일 경우(eps보다 작으면 0으로 취급) 계산하지 않고 밑으로 내립니다. 계산 방법은 선형대수를 들으시면 아시게 됩니다.
for (j in (i+1):(p+1)) U.pls[i,j] <- U.pls[i,j]/U.pls[i,i]
U.pls[i,i] <- 1
if (i < p) {
for (k in (i+1):p) U.pls[k,] <- U.pls[k,] - U.pls[k,i]/U.pls[i,i]*U.pls[i,]
}
}
U.pls
#여기까지 해서 rref(Reduced Row Echelon Form)이 종료됩니다.
x <- rep(0,p)
for (i in p:1){
if (i < p){
temp <- 0
for (j in (i+1):p) temp <- temp + U.pls[i,j]*x[j]
x[i] <- U.pls[i,p+1] - temp
}
else x[i] <- U.pls[i,p+1]
}
x
#이제 x를 누르시면 2번의 결과와 같은 124,75,31 이 표시됩니다.
어때요 참 쉽죠?
글쓰기