R 로 함수를 짜보자

뜬금없이 왜 이런 글을 쓰냐하면…

대학원 수업 중에 괜히 있는 함수 직접 짜보라는 수업이 있어서 귀찮게 항상 함수를 짜는데 왜 하는 것인지도 잘 모르겠고… 뭐 막 어려운 것도 아닌 주제에 시간은 꽤 들고… 블로그에 올리면 누군가에겐 훗날 이 과목을 들을 사람들이라든지 도움이 되지 않을까 싶어서이다. 교수님께 들킬까봐 문제라든가 데이터 돌려보는 건 올리지 않고 내가 짠 함수만 올릴 것이다. 아 그리고 과제 하다가 새삼 새로 알게 된 것도 올려볼까 한다.

0. read.fwf

fwf 는 f ixed w idth f ormatted 를 말한다. 도대체 csv 놔두고 왜 쓰는지 모를, 개인적으로 진짜 개 싫어하는 서타일인데 fixed width data의 예를 들면 아래와 같다.

NAME                STATE     TELEPHONE  

John Smith          WA        418-Y11-4111

Mary Hartford       CA        319-Z19-4341

Evan Nolan          IL        219-532-c301

read.table 못쓴다. 딱 봐도 그래보이지… 그래서 read.fwf 를 쓴다. 아유 그럼 다행이네요 왜 싫어하시나요? 그것은 왜냐하면 제가 width를 직접 세어서 입력해줘야 하기 때문입니다. 와 거지같아

예를 들어 저 위의 데이터에서 Mary Hartford 는 알파벳이 12개에 중간에 스페이스 하나 있으니까 width 가 25다 그리고 CA 는 width 가 4다 근데 그렇다고 width를 c(25,4,13) 이라고 설정하고 데이터 읽으면 첫 column 은 제대로 나올지 몰라도 두번째 column 은 빈칸이 나와버리고 그 다음 column에 CA 가 나오고 전화번호는 중간에 잘려버린다 왜냐하면 column들 사이사이의 space 를 알아서 지우지 않고 걔네들도 데이터라고 생각하기 때문이다 진짜 다시 생각해도 개빡쳐

이렇다보니 중간 스페이스 개수도 세어서 새로운 column으로 간주한 다음 나중에 지우든지, 아니면 앞뒤로 흡수시켜버리든지 해야한다. 한학기 내내 fwf 데이터 쓰는 거 같은데 나는 울란다. 더 쉬운 방법 있으면 댓글로 알려주세요 제발…

1. 평균과 sample covariance matrix 와 sample correlation matrix

…를 한번에 계산할 수 있는 함수를 만들었다.

corfunc = function(data,cols) {
  custom_m = c()
  for (i in cols) {
    val_m = sum(data[i])/dim(data[i])[1]
    custom_m = c(custom_m, val_m)
  }

  automat = as.matrix(data[cols])
  meanmat = matrix(rep(custom_m,times=dim(data)[1]),byrow = TRUE, nrow = dim(data)[1])
  X = automat - meanmat

  custom_s = t(X) %*% X / (dim(X)[1]-1)
  sigmat = matrix(rep(diag(custom_s),times = dim(X)[2]),nrow = dim(X)[2])
  sqrtsig = sqrt(sigmat)
  custom_r = custom_s/(sqrtsig * t(sqrtsig))

  output = list(custom_m, custom_s, custom_r)
  return(output)
}

input 으로 dataframe 과 vector 을 받는다. vector 은 사용할 column name 이나 column number 로 이루어진 vector 를 말한다.

1. columnwise mean (열별 평균)을 구해서 custom_m vector 에 차례대로 저장한다.
2. sample covariance matrix 계산을 위해 dataframe 을 matrix 로 만든다. 구해놨던 평균들도 matrix 로 만들어서 계산할 때 쓴다. 계산된 sample covariance matrix 는 custom_s 로 저장한다.
3. custom_s 를 이용해 custom_r 을 계산한다.
4. custom_m, custom_s, custom_r 을 하나의 리스트로 묶고 리스트를 return 한다.

이렇게 만든 함수 결과물이랑 R의 내장함수 (mean(),cov(),cor()) 결과를 all.equal() 로 비교했더니 같다고 나왔다. (== 로 비교하면 당연히 FALSE 로 나온다. 무리수 계산 오차때문에)

T test for correlation

Null hypothesis : $\rho = 0$

$\alpha=0.05$

T test statistics 의 식은 아래와 같다.

$T = \frac{r \sqrt{n-2}}{\sqrt{1-r^2}}$

그냥 똑같이 짜면 된다.

hyptest = function(r, n) {
  T = r * sqrt(n-2) / sqrt(1-r^2)
  t_0975 = qt(p=0.975,df = n-2)
  print(paste("Test Statistic:",T))
  if (abs(T)>t_0975) {
    print("Reject Null Hypothesis")
  } else { print("Do not reject Null Hypothesis")}

}

1. t test statistic 계산하여 프린트
2. t test statistic 의 절댓값과 $t_{0.975}$ 값을 비교 (test statistics 가 크면 null reject)
3. 결과 프린트

Fisher’s Confidence interval, Ruben’s confidence interval

correlation 에 대한 confidence interval 들을 구하는 방법들이다. 수식은 여기를 참고하자.

fisherci = function(r,n) {
  UB = tanh(atanh(r)+qnorm(0.975)/sqrt(n-3))
  LB = tanh(atanh(r)-qnorm(0.975)/sqrt(n-3))
  return(c(LB,UB))
}

rubenci = function(r,n) {
  u = qnorm(0.975)
  r_star = r/sqrt(1-r^2)
  a = 2*n - 3 - u^2
  b = r_star * sqrt((2*n - 3)*(2*n -5))
  c = (2*n -5 -u^2) * r_star^2 - 2*u^2
  y_1 = (b-sqrt(b^2-a*c))/a
  y_2 = (b+sqrt(b^2-a*c))/a
  UB = y_2 / sqrt(1+y_2^2)
  LB = y_1 / sqrt(1+y_1^2)
  return(c(LB,UB))
}

1. input 으로 correlation matrix 랑 observation 개수를 받는다.
2. output 으로 lower bound 와 upper bound 가 저장된다.