Harlequin Blue

2013년 6월 24일 월요일

정규분포 관련 성질

정규분포 관련 성질 -> 더 자세한 내용은 [PDF] 다운로드.

표준정규분포의 pdf를 $\phi(x)$로 놓고 스케일링한 것을 $\phi_{\sigma}(x) = \phi(x/\sigma)/\sigma$라고 표기하면 편리하다. 이때 $N(\mu,\sigma^2)$의 pdf를 $\phi_\sigma(x-\mu)$로 쓸 수 있다.

표준정규분포 $\phi$의 도함수들은 $\phi$를 이용하여 표현할 수가 있다. 이 관계를 이용하면 미분 계산을 쉽게할 수 있다.
\begin{align}
\phi'(x) &= -x \phi(x)\\
\phi''(x) &= (x^2-1)\phi(x)
\end{align}

Convolution의 성질을 이용하면 복잡한 계산들을 매우 간단하게 할 수 있다. 우선 convolution은 다음과 같이 정의된다.
\[
(f*g)(x) = \int f(x-y) g(y) dy
\]
정규분포의 pdf에서 다음과 같은 성질을 확인할 수 있다.
\begin{align}
(\phi_{\sigma_1}*\phi_{\sigma_2})(x) = \phi_{\sqrt{\sigma_1^2+\sigma_2^2}}(x)
\end{align}
이 성질을 이용하여 $\int \phi^2(x) dx$, $\int \phi_{\sigma_1}(x)\phi_{\sigma_2}(x) dx$ 등을 쉽게 계산할 수 있다.

2013년 6월 7일 금요일

R 프로그래밍: R에서 C 함수 부르기

이 글에 사용된 코드와 명령은 MacOSX와 Linux에서 확인된 것이다. 특별히 설치되어야 하는 것은 없다. 기본적으로 R과 기본 개발도구가 설치되어 있으면 된다. MacOSX에서는 Xcode를 설치하면 된다. Fortran/C/C++ 무엇으로 짜야 하나 고민이 된다. 다음은 확인된 사실은 아니지만...

R에서 vectorize할수 있는 데까지 다 했는데도 속도를 빠르게 하고 싶다면 Fortran/C/C++을 생각해야 한다.
Fortran이 C보다 우월한가? 논쟁만 부르고 결론은 나지 않는 질문이다. 그럴 수도 있고 아닐 수도 있다.
R에서는 공식적(?)으로 Fortran 77만 지원한다. Fortran 95로 작성한 코드를 R에서 부르는 방법이 없는 건 아니지만... 어떤 사람은 Fortran 77은 절대로 쓰지 말아야한다고도 한다. Fortran 95도 있고 최근 2010년에 표준화된 흔히 부르는 이름으로 "Fortran 2008"도 있다.

.Call

R 패키지를 개발할 때나 아니면 그냥 함수를 짤 때 C/C++/Fortran으로 짠 코드를 불러서 쓸 수 있다. R 함수 중 .C()는 컴파일된 C 코드를 호출할 때, .Fortran()는 컴파일된 Fortran 77 코드를 호출할 때 쓸 수 있다. 이것보다 좀더 유연하고 성능도 나은 것이 .Call()이다. 이것은 C/C++ 코드에 R 객체를 넘겨줄 때 사용된다. .External()도 마찬가지인데 .Call()은 고정된 수의 인자를 넘기고 .External()은 argument list로 넘긴다는 데에 차이가 있다. .Call()을 통해 사용할 C 코드를 작성할 때

R.h, Rinternals.h 헤더를 반드시 include한다.
R 객체를 받을 때 SEXP, REALSXP, INTSXP, STRSXP 타입을 이용한다.
R 객체의 속성은 getAttrib()로 뽑아낼 수 있다.
R 객체 안에 들어있는 데이터에 대한 포인터는 REAL(), INTEGER() 등을 이용하여 얻을 수 있다.
R 객체를 만들 때는 allocVector(), allocMatrix()를 이용하여 SEXP 타입 객체로 만든다.
이때 항상 PROTECT()와 UNPROTECT()를 이용한다. 안전한 메모리 확보를 위한 장치이다.

예를 들어 다음과 같은 함수를 작성하고 myfun.c 파일에 저장하였다고 하자.

#include <R.h>
#include <Rinternals.h>

SEXP say(SEXP x) {
    return x
}

컴파일은 다음과 같이 한다.

$ R CMD SHLIB myfun.c

컴파일하면 만들어지는 myfun.so를 R에서 불러 사용한다. R에서 다음과 같이 로딩한다.

> dyn.load("myfun.so")

로딩이 되었으면 다음처럼 실행할 수 있다.

> .Call("say", 3)
[1] 3

예제: 평균 계산하기

이제 C로 평균을 구하는 함수를 다음과 같이 작성하고 mystat.c 파일로 저장하였다고 하자.

// mystat.c
#include <R.h>
#include <Rinternals.h>

SEXP mean(SEXP x) {
  double *vec = REAL(x);
  int n = length(x);
  SEXP val = PROTECT(allocVector(REALSXP, 1));

  double sum = 0;
  for(int i=0; i < n; i++) sum += vec[i];
  REAL(val)[0] = sum / n;
  UNPROTECT(1);
  return val;
}

컴파일을 하면 mystat.so 파일이 생긴다.

$ R CMD SHLIB mystat.c

R에서 이것을 불러쓰기 위해 다음처럼 wrapper를 작성하자.

mymean <- function(x) {
    if(!is.loaded("mystat")) dyn.load("mystat.so")
    .Call("mean", as.numeric(x))
}

이제 R에서 mymean()이라는 함수를 사용하면 된다.

속도 비교

R에는 물론 평균을 계산해주는 mean() 함수가 있다. 만약 R에서 for으로 평균을 구하면 어떨까? iMac 3.06GHz Intel Core 2 Duo에서 1e7개의 난수를 더한 후 평균을 구하는 데에 다음처럼 오랜 시간이 걸린다.

> n <- 1e7
> x <- runif(n)
> system.time( {
+     total <- 0
+     for(i in 1:n) total <- total + x[i]
+     total / n
+ })
   user  system elapsed 
  8.436   0.076   9.165

물론 R의 mean() 함수를 이용하면 빠르게 계산된다.

> system.time( mean(x) )
   user  system elapsed 
  0.079   0.001   0.086

방금 C로 짠 함수의 경우에는:

> system.time( mymean(x) )
   user  system elapsed 
  0.020   0.000   0.023

R의 내장 mean()보다 빠른 속도를 보여주고 있는데 데이터에 따라서 다를 수 있으니 항상 성능이 좋다고는 못한다. 어쨌든 R에서 for문 돌리는 것에 비하면 빠르다.

행렬

R에서 넘어온 행렬은 컬럼을 쌓은 벡터로 받아서 사용한다. $p\times q$ 행렬 $A$와 $q\times r$ 행렬 $B$의 행렬곱은 다음처럼 계산할 수 있다.

#include <R.h>
#include <Rinternals.h>

/**
 * @SmatA Sp by Sq matrix
 * @SmatB Sq by Sr matrix
 **/
SEXP prod(SEXP SmatA, SEXP SmatB, SEXP Sp, SEXP Sq, SEXP Sr) {

  int p = INTEGER(Sp)[0];
  int q = INTEGER(Sq)[0];
  int r = INTEGER(Sr)[0];

  // R matrix as a vector stacked column by column
  double* a = REAL(SmatA);
  double* b = REAL(SmatB);
  SEXP ans = PROTECT(allocVector(REALSXP, p*r));
  double sum;
  int i, j, k;
  for(i=0; i < p; i++) {
    for(k=0; k < r; k++) {
      sum = 0;
      for(j=0; j < q; j++) {
        sum += a[i + j*p] * b[j + k*q];
      }
      REAL(ans)[i + k*p] = sum;
    }
  }
  UNPROTECT(1);
  return ans;
}

2013년 6월 1일 토요일

R에서 parallel 패키지 이용하여 멀티코어 계산하기

R 3.0 기준으로 얘기하자면 parallel 패키지는 멀티코어 이용도 가능하게 해준다. multicore 패키지는 더이상 사용되지 않는다. 다음처럼 패키지를 로딩한다.

> require(parallel)

Linux, MacOSX에서 간단하게 mclappy와 같은 함수를 이용하는 예를 들면 다음과 같다. iMac 3.06GHz Intel Core 2 Duo에서 돌린 결과이다.

> system.time( out <- sapply(1:100000, function(.) mean(rnorm(100))))
   user  system elapsed
  2.916   0.039   2.959
> system.time( out <- mclapply(mc.cores=2, 1:100000, function(.) mean(rnorm(100))))
   user  system elapsed
  1.487   0.073   1.603

MS Windows에서는 mclapply가 제대로 작동하지 않는다. 대신 clusterApply를 쓸 수 있다. Intel Core i3 3.20GHz에 Windows 7이 설치된 PC에서 돌린 결과이다.

> system.time( out <- sapply(1:100000, function(.) mean(rnorm(100))) )
   user  system elapsed 
   2.66    0.00    2.67 
> cl <- makeCluster(4)
> system.time( out <- clusterApply(cl, 1:100000, function(.) mean(rnorm(100))) )
   user  system elapsed 
  12.85    8.23   21.94 
> stopCluster(cl)

이해가 가지 않는 결과이다. 코어 하나만 쓰는 것보다 코어 4개를 모두 쓰면 더 느려진다. 작업관리자에서 CPU 사용현황을 모니터링하면 멀티코어로 실행되는 것은 분명히 맞다. 이것은 core에 작업을 할당하는 문제 때문이다. 위 실험의 경우 100000개의 작업을 4개의 코어에 분배하여야한다. 그리고 하나의 작업은 매우 연산량이 매우 작다. 그 분배 처리 작업 때문에 하나의 코어에서 10만개 작업을 순차적으로 할당하는 것보다 느리게 되는 것이다. 코어가 4개이므로 10만번 시뮬레이션을 25000번의 시뮬레이션 4개로 나누어 할당하면 속도가 빨라진다.

> simul <- function(n) sapply(1:n, mean(rnorm(100)))
> system.time( out <- sapply(rep(25000, 4), simul) )
   user  system elapsed 
   2.60    0.01    2.68 
> cl <- makeCluster(4)
> system.time( out <-clusterApply(cl, rep(25000, 4), simul) )
   user  system elapsed 
   0.00    0.02    1.04

2012년 10월 20일 토요일

조건부 확률의 조건부 확률. posterior predictive distribution.

조건부 확률

조건부 확률을 보통 $A$, $B$ 기호를 가지고 이야기하지만 여기서는 좀 헤깔릴 수가 있으니 $Y$, $A$로 표시해 보자. $Y$는 관심 사건이고 $A$는 주어진 조건에 해당하는 사건이다. 조건부 확률의 기본적인 정의에 따라 쓰면 다음과 같다.

\[
P(Y|A) = \frac{P(Y, A)}{P(A)}
\]

그리고 전확률 법칙(Law of Total Probability)라는 것도 있다. 표본 공간이 $B_1,B_2,\cdots,B_n$으로 분할(partition)되어있을 때
\[
P(Y) = \sum_{i=1}^n P(Y|B_i)P(B_i)
\]
이 성립한다는 법칙이다. 이것은 $P(Y)$를 구하기가 힘들 때 $P(Y|B_i)$들은 구하기 쉽고 각 $P(B_i)$들을 이미 알고 있는 경우에 유용하게 쓸 수 있는 법칙이다.

그런데, 이런 경우를 생각해보자. $P(Y|A)$를 구해야하는데 이것을 구하기가 힘들다. 그런데, $B_i$로 이 문제를 분할하여 구하는 것은 쉽다고 하자. 조건부 확률의 조건부 확률을 생각해보는 것이다. 다음처럼 구할 수 있다.
\[
P(Y|A) = \sum_{i=1}^n P(Y|B_i,A)P(B_i|A)
\]
이 식을 바로 생각하려고 하면 헤깔릴 때가 있다. $Y$를 $A$라는 조건 아래서 생각하고 있는 데 거기에 $B_i$라는 조건을 추가로 고려하면 어떻게 될까, 하는 것인데. 두번째 항을 $P(A|B_i)$인 것으로 착각을 할때가 있다. 식을 따져보면 당연한 것이기는 하다.
\[
\sum_{i=1}^n P(Y|B_i,A)P(B_i|A)
= \sum_{i=1}^n \frac{P(Y,B_i,A)}{P(B_i,A)}\frac{P(B_i,A)}{P(A)}
= \sum_{i=1}^n \frac{P(Y,B_i,A)}{P(A)}
= \frac{P(Y,A)}{P(A)}
= P(Y|A)
\]

Posterior Predictive Distribution

미지의 평균과 분산이 $\theta =(\mu, \sigma^2)$인 대상에 대해 $y=(y_1,\cdots,y_n)$을 관측했다. $\tilde y$는 다음 번 관측 대상의 값이다. 이때 $\tilde y$의 분포를 posterior predictive distribution이라고 한다. 다음처럼 구한다.

\begin{align}
p(\tilde y|y) &= \int p(\tilde y, \theta | y) d\theta\\
&= \int p(\tilde y|\theta, y)p(\theta|y) d\theta\\
&= \int p(\tilde y|\theta)p(\theta | y) d\theta.\\
\end{align}
마지막 행이 성립하는 것은 $\tilde y$와 $y$가 given $\theta$에서 조건부 독립이기 때문이다.

예제

여자아이가 태어날 확률을 $\theta$라고 하자. 서로 독립적인 $n$명의 신생아를 관찰하였을 때 여자아이의 수는 이항분포 $Binomial(n, \theta)$를 따른다. 사전분포로는 균일분포 $\theta \sim Unif(0,1)$을 사용하자. 1000명의 신생아를 관찰하였더니 485명이 여자아이였다. 1001번째 신생아가 여자아이일 확률은?

$n$명 중 여자아이의 수 $y$라고 하자. likelihood는
\[
p(y|\theta) = \binom{n}{y} \theta^y(1-\theta)^{n-y}
\]
piror distribution는
\[
\pi(\theta) = 1, \qquad 0 < \theta < 1
\]
posterior distribution은
\begin{align}
p(\theta|y) &\propto p(y|\theta)\pi(\theta)\\
&\propto \theta^y (1-\theta)^{n-y}
\end{align}
형태이므로 $Beta(y+1, n-y +1)$ 분포임을 알 수 있다.

이제 새로 태어날 아이에 대한 관측을 $\tilde y$라고 하면 여자아이일 확률은
\begin{align}
p(\tilde y = 1| y) &= \int_0^1 p(\tilde y = 1 | \theta,y)p(\theta|y) d\theta\\
&= \int_0^1 p(\tilde y = 1 | \theta)p(\theta|y) d\theta \\
&= \int_0^1 \theta p(\theta|y) d\theta\\
&= E(\theta|y) & (\theta|y \sim Beta(y+1, n-y+1))\\
&= \frac{y+1}{n+2} & (\because X\sim Beta(\alpha, \beta),\; EX = \alpha/(\alpha+ \beta))\\
\end{align}

따라서, 새로 태어날 아이가 여자아이일 확률을 486/1002.

2012년 2월 7일 화요일

베이즈통계 책

Berger, J. O. (1985). Statistical Decision Theory and Bayesian Analysis. Springer- Verlag.
Ghosh, J.K. and Ramamoorthi (2002). Bayesian Nonparametrics. Springer.
Lange, K. (1998). Numerical Analysis for Statisticians. Springer.
Robert, C. P. and Casella, G. (1999), Monte Carlo Statistical Methods. Springer.
Norris, J.R. (1997). Markov Chains. Cambridge.
Chen. M.-H., Shao, Q.-M. and Ibrahim, J.G. (2000). Monte Carlo Methods in Bayesian Computation. Springer.
Gelman, A., Carlin, J.B., Stern, H.S. and Rubin, D.B. (2003). Bayesian Data Analysis. Chapman & Hall/CRC.
Ghosh, J.K., Delampady, M. and Samanta, T. (2006). An Introduction to Bayesian Analysis. Springer.
Congdon, P. (2003). Applied Bayesian Modelling. Wiley.
Congdon, P. (2005). Bayesian Models for Categorical Data. Wiley.
Nyhoff, L. R. and Leestma, S. C. (1996). Introduction to Fortran 90 for Engineers and Scientists. Prentince Hall.

선형모형 책

[석사] Alvin C. Rencher and G. Bruce Schaaalje. Linear Models in Statistics (2nd ed.). 2008. Wiley. google books
[대학원] Youngjo Lee, John A. Nelder, Yudi Pawitan. Generalized Linear Models with Random Effects: Unified Analysis via h-likelihood. Chapman & Hall/CRC. 2006. [google books]

확률론, 확률과정 교재

[학부 확률과정] Rick Durrett. Essentials of Stochastic Processes. Springer. 2010. [google books]
George F. Lawler. Introduction to Stochastic Processes (2nd ed.). Chapman & Hall/CRC. 2006. [google books]
Geoffrey Grimmett and David Stirzaker. Probability and Random Process (3rd ed.). Oxford Unisersity Press. 2001. [google books]
Olav Kallenberg. Foundations of Modern Probability. Springer. 2002. [google books]
Sidney I. Resnick. Adventures in Stochastic Processes. Birkhauser. 1992. [google books]
Zdislaw Brzezniak and Tomasz Zastawniak. Basic Stochastic Processes. Springer. 2000. [google books]

교호 작용이 있는 모형에서 다중 비교

R에서 multcomp 패키지의 glht를 이용하면 다중비교(multiple comparison)가 가능하다. 더 자세한 설명은 multcomp 매뉴얼을 참조하자.

http://cran.r-project.org/web/packages/multcomp/index.html

기본

R에 기본으로 들어있는 warpbreaks 데이터세터를 예제로 한다. 이 데이터는 다음과 같이 연속형 반응변수 breaks와 범주형 설명변수 wool, tension으로 구성되어있다.

> summary(warpbreaks)
     breaks      wool   tension
 Min.   :10.00   A:27   L:18   
 1st Qu.:18.25   B:27   M:18   
 Median :26.00          H:18   
 Mean   :28.15                 
 3rd Qu.:34.00                 
 Max.   :70.00

다음과 같은 선형모형을 적합하였다고 하자.

> mod <- lm(breaks ~ wool + tension, data=warpbreaks)

그리고 나서 tension의 세 수준 L, M, H 사이의 다중 비교를 할 필요가 있다고 하자. 물론 우선 multcomp패키지를 로딩해야 하고 glht() (General Linear Hypothesis Test)를 이용한다.

> require(multcomp)
> mod.mc <- glht(mod, linfct=mcp(tension="Tukey"))

결과는 다음과 같이 summary() 함수를 이용해 볼 수 있다. plot(mod.mc)으로 그래프를 그려 볼 수도 있다.

> summary(mod.mc)

  Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = breaks ~ wool + tension, data = warpbreaks)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)   
M - L == 0  -10.000      3.872  -2.582   0.0336 * 
H - L == 0  -14.722      3.872  -3.802   0.0011 **
H - M == 0   -4.722      3.872  -1.219   0.4474   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported -- single-step method)

교호작용이 있는 경우 다중비교

교호작용(interaction)이 있는 모형이 필요하다고 하자. 예를 들면 다음과 같다.

> mod1 <- lm(breaks ~ wool * tension, data=warpbreaks)

이 경우에 모든 가능한 조합은 6가지가 있다. 다음과 같이 6개의 셀을 평균을 비교하는 것이다.

> xtabs(~wool+tension, warpbreaks)
    tension
wool L M H
   A 9 9 9
   B 9 9 9

이때 우리가 다중비교를 원하는 대상이 wool의 동일한 수준 내에서 tension을 비교하는 것이라고 하자. 우선 모형을 만들때 다음과 같이 interaction() 함수를 이용해 새로운 변수를 만들고 절편은 없이 하는 것이 생각하기 편리하다.

tw <- with(warpbreaks, interaction(tension, wool))
mod1 <- lm(breaks ~ tw-1, warpbreaks)

이렇게 했을 때 회귀계수는 각 셀의 평균이 된다. 다음과 같다.

> coef(mod1)
   twL.A    twM.A    twH.A    twL.B    twM.B    twH.B 
44.55556 24.00000 24.55556 28.22222 28.77778 18.77778

이것과 순서를 맞추어 다음처럼 contrast matrix를 만들어주자.

contr <- rbind("A:M-L" = c(-1,1,0,0,0,0),
               "A:H-L" = c(-1,0,1,0,0,0),
               "A:H-M" = c(0,-1,1,0,0,0),
               "B:M-L" = c(0,0,0,-1,1,0),
               "B:H-L" = c(0,0,0,-1,0,1),
               "B:H-M" = c(0,0,0,0,-1,1))

첫행을 보면 "A:M-L"이라고 이름을 붙였는데 wool=A에서 tension=M과 tension=L을 비교하겠다는 것이다. 위에 적합한 모형의 회귀계수 중 두번째 것에서 첫번째 것을 빼준 것이 이에 해당한다. 따라서 c(-1,1,0,0,0,0)이 된다. 이제 다음과 같은 결과를 얻을 수 있다.

> summary(glht(mod1, linfct=contr))

  Simultaneous Tests for General Linear Hypotheses

Fit: lm(formula = breaks ~ tw - 1, data = warpbreaks)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)   
A:M-L == 0 -20.5556     5.1573  -3.986  0.00133 **
A:H-L == 0 -20.0000     5.1573  -3.878  0.00185 **
A:H-M == 0   0.5556     5.1573   0.108  0.99996   
B:M-L == 0   0.5556     5.1573   0.108  0.99996   
B:H-L == 0  -9.4444     5.1573  -1.831  0.30801   
B:H-M == 0 -10.0000     5.1573  -1.939  0.25536   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported -- single-step method)

2012년 2월 1일 수요일

Sparse Matrix

Sparse Matrix는 어디서나 중요한 문제가 된다.

MATLAB/Octave

우선 MATLAB에서다음과 같은 sparse matrix format을 사용할 수 있다. 세 컬럼으로 이루어진 텍스트 파일이다. row index, col index, value로 구성한다. 예를 들어, 다음과 같다.

이 텍스트 파일을 mat.mtl라는 이름으로 저장했다고하자. load를 이용해 파일을 읽고 spconvert()를 이용해 sparse matrix로 바꾼다.

octave-3.4.0:38> load mat.mtl;
octave-3.4.0:39> x = spconvert(mat);

이렇게 만들어진 sparse matrix는 다음과 같다.

octave-3.4.0:40> x
x =

Compressed Column Sparse (rows = 5, cols = 5, nnz = 8 [32%])

  (1, 1) ->  1
  (2, 2) ->  2
  (5, 2) ->  8
  (2, 3) ->  5
  (3, 3) ->  3
  (4, 3) ->  6
  (4, 4) ->  4
  (4, 5) ->  7

일반적인 형렬 형태로 보면 다음과 같다.

octave-3.4.0:41> full(x)
ans =

   1   0   0   0   0
   0   2   5   0   0
   0   0   3   0   0
   0   0   6   4   7
   0   8   0   0   0

2012년 1월 14일 토요일

사전분포

베이지안 통계에서 주관적 사전분포를 정하는 방법이 문제가 된다.

공액 사전 분포 (conjugate prior)
비공액 사전 분포 (non-conjugate prior)
무정보적 사전 분포 (noninformative prior)
제프리스 사전 분포 (Jeffreys prior)

베이지안에서 많이 사용하는 분포

Univariate	pdf	Multivariate	pdf
Binomial$(n,p)$	$\binom{n}{x}p^x(1-p)^{n-x}$	Multinomial$(n,p_1,\cdots,p_k)$	$\frac{n!}{x_1!\cdots x_k!}p_1^{x_1}\cdots p_k^{x_k}$, $\sum x_i = n$
Normal$(\mu,\sigma^2)$	$(2\pi\sigma^2)^{-1/2}\exp\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\}$	Multivariate Normal$(\mu, \Sigma)$	$(2\pi)^{-\frac{k}{2}}\|\Sigma\|^{-\frac{1}{2}}\exp\left\{(x-\mu)'\Sigma^{-1}(x-\mu)\right\}$
Gamma$(\alpha,\beta)$	$\frac{x^{\alpha-1}e^{-x/\beta}}{\Gamma(\alpha)\beta^\alpha}$	Wishart$(n,\Sigma)$
Beta$(\alpha, \beta)$	$\frac{x^{\alpha-1}(1-x)^{\beta-1}}{\text{B}(\alpha,\beta)}$	Dirichlet$(\alpha_1,\cdots,\alpha_{k+1})$

Exponential$(\beta)$: Gamma$(\alpha=1, \beta)$
Chi-squared$(k)$: Gamma$(\alpha=k/2, \beta=1/2)$
Uniform$(0,1)$: Beta$(\alpha=1, \beta=1)$

역감마분포 Inverse-gamma distribution

$1/X$가 Gamma 분포를 따르면 $X$는 Inverse-gamma 분포를 따른다고 한다.
\[
f(x) = \frac{\beta^\alpha)}{\Gamma(\alpha)} x^{-(\alpha+1)} e^{-\beta/x}
\]

디리클레 분포 Dirichlet distribution

디리클레 분포의 pdf:
\[
f(x_1,\cdots,x_{k};\alpha_1,\cdots,\alpha_{k+1}) = \frac{\Gamma(\sum_{i=1}^{k+1}\alpha_i)}{\prod_{i=1}^{k+1}\Gamma(\alpha_i)}\prod_{i=1}^{k+1} x_i^{\alpha_i -1 }
\]

위샤트 분포 Wishart distribution

다변량정규분포를 따르는 $n$개의 iid 확률벡터 $X_1,\cdots,X_n \stackrel{iid}{\sim} N(0,\Sigma)$가 있을 때
\[
\sum_{i=1}^n X_iX_i^T \sim \text{Wishart}(n,\Sigma)
\]
이다.

역위샤트 분포 Inverse-Wishart distribution

2011년 12월 20일 화요일

[Octave] Simple Linear Regression

Octave/MATLAB에서 단순선형회귀

모의실험: 난수생성

단순선형회귀의 간단한 예를 모의생성해 보자.
\[
Y_i = 0.5 + 2 X_{1i} + \epsilon_i, \qquad \epsilon_i \stackrel{iid}{\sim} N(0,1)
\]
와 같은 모형을 생각해 보자.

octave-3.4.0:171> x1 = [1:10]';
octave-3.4.0:171> n = length(x1);
octave-3.4.0:173> x = [ones(n,1), x1]
x =

    1    1
    1    2
    1    3
    1    4
    1    5
    1    6
    1    7
    1    8
    1    9
    1   10

octave-3.4.0:174> beta = [0.5; 2]
beta =

   0.50000
   2.00000

다음과 같이 randn() 함수를 이용하여 표준정규분포에서 난수를 생성할 수 있다. normrnd() 함수를 이용할 수도 있다. 여기에서는 동일한 실험 결과를 얻기 위하여 seed를 지정하여 randn으로 난수를 생성하였다.

octave-3.4.0:175> randn("seed", 12345);
octave-3.4.0:176> y = x * beta + randn(n,1)
y =

    3.9025
    3.6416
    6.3307
    8.7463
   11.6325
   12.6323
   14.5467
   16.9846
   18.5778
   21.1140

그래프는 다음과 같이 그릴 수 있다.

octave-3.4.0:177> plot(x1, y, "o")

회귀계수의 추정: 최소제곱법

최소제곱법(Least Square Method)를 이용한 회귀계수 추청치는 $b_1 = S_{XY}/S_{XX}$, $b_0 = \overline Y - b_1 \overline X$라는 공식을 이용하여 구할 수 있다.

octave-3.4.0:50> Sxy = sum((x1 - mean(x1)) .* (y - mean(y)))
Sxy =  165.56
octave-3.4.0:51> Sxx = sum(((x1 - mean(x1)).^2))
Sxx =  82.500
octave-3.4.0:52> b1 = Sxy / Sxx
ans =  2.0068
octave-3.4.0:53> b0 = mean(y) - b1 * mean(x1)
b0 =  0.77332

적합된 직선 $y= 0.77332 + 2.0068x$를 원데이터와 함께 그리면 다음과 같다.

octave-3.4.0:59> plot(x1, y, "o");
octave-3.4.0:60> hold on;
octave-3.4.0:61> plot(x1, b0 + b1 * x1);

행렬 표현을 이용하면 $b = (X'X)^{-1}X'y$이고 다음과 같이 계산할 수 있다.

octave-3.4.0:206> b = inv(x' * x) * x' * y
ans =

   0.77332
   2.00683

Octave에서는 Ordinary Least Square 문제를 풀 수 있는 ols(Y, X) 함수를 제공한다.

octave-3.4.0:207> ols(y, x)
ans =

   0.77332
   2.00683

적합치와 잔차

반응변수의 적합치(fitted values)는 회귀계수를 구하였으니 당연히

octave-3.4.0:54> yhat = b0 + b1 * x1;

으로 구할 수 있다. 잔차(residuals)는 관측치와 기대치의 차이이므로

octave-3.4.0:55> res = y - yhat;

로 구할 수 있다.

행렬표현을 이용하면 $\hat y = Xb=X(X'X)^{-1}X'y$로 구할 수 있다. 이때 $H=X(X'X)^{-1}X'$을 모자행렬(hat matrix)라고 한다.

octave-3.4.0:52> H = x * inv(x' * x) * x';
octave-3.4.0:52> yhat = H * y
ans =

    2.7802
    4.7870
    6.7938
    8.8007
   10.8075
   12.8143
   14.8212
   16.8280
   18.8348
   20.8417

잔차(residuals)는 $e = (I-H)y$로 구할 수 있다.

octave-3.4.0:53> n
n =  10
octave-3.4.0:54> res = (eye(n) - H) * y
ans =

   1.122329
  -1.145389
  -0.463069
  -0.054352
   0.824982
  -0.182018
  -0.274477
   0.156663
  -0.256991
   0.272321

결정계수

적합된 회귀식이 얼마나 타당한가를 알아보기 위한 것으로 결정계수(coefficient of determination)가 사용된다. 결졍계수는
\[
R^2 = \frac{SSR}{SST} = \frac{SST -SSE}{SST}
\]
이고 여기서 총제곱합(SST), 회귀제곱합(SSR), 오차제곱합(SSE)은 $SST=SSR+SSE$ 관계를 가지며 다음과 같다.
\[
SST = \sum_{i=1}^n (Y_i - \overline Y)^2, \qquad SSR = \sum_{i=1}^n (\hat Y_i - \overline Y)^2, \qquad SSE=\sum_{i=1}^n(Y_i-\hat Y_i)^2
\]
각각을 계산하면 다음과 같다.

octave-3.4.0:68> SST = sum((y - mean(y)).^2)
SST =  336.00
octave-3.4.0:70> yhat = b0 + b1*x1;
octave-3.4.0:71> SSR = sum((yhat - mean(y)).^2)
SSR =  332.26
octave-3.4.0:72> SSE = sum((y - yhat).^2)
SSE =  3.7427

결정계수 $R^2$는 다음과 같다.

octave-3.4.0:73> R2 = SSR / SST
R2 =  0.98886

분산분석

단순선형회귀의 분산분석은 귀무가설 $H_0: \beta_1 =0$ 대 $H_1: \beta_1 \neq 0$을 검정하기 위하여 귀무가설 하에
\[
F_0 = \frac{(SSR/\sigma^2)/1}{(SSE/\sigma^2)/(n-2)} = \frac{MSR}{MSE} \sim F(1, n-2)
\]
이라는 사실을 이용한다.

octave-3.4.0:75> MSR = SSR / 1
MSR =  332.26
octave-3.4.0:76> MSE = SSE / (n-2)
MSE =  0.46784
octave-3.4.0:77> F0 = MSR / MSE
ans =  710.19

유의확률(p-value)를 계산하면 다음과 같이 매우 작으므로 귀무가설을 기각한다.

octave-3.4.0:80> 1 - fcdf(F0, 1, n-2)
ans =  4.2286e-09

또는 $F(1, 8)$ 분포에서 유의수준 0.05에 해당하는 5.3177보다 $F_0=710.19$가 매우 크므로 귀무가설을 기각한다.

octave-3.4.0:81> finv(0.95, 1, n-2)
ans =  5.3177

추론

기울기에 대한 추론

기울기 $\beta_1$의 추정량 $b_1$은 $E(b_1)=\beta_1$이고 $\text{Var}(b_1) = \sigma^2/S_{XX}$인 정규분포를 따른다. 표준화하면
\[
\frac{b_1 - \beta_1}{\sigma/\sqrt{S_{XX}}} \sim N(0,1)
\]
이고 $\sigma^2$을 모르므로 대신에 추정량 $s^2$을 사용하면
\[
\frac{b_1 - \beta_1}{s/\sqrt{S_{XX}}} \sim t(n-2)
\]
이다. 여기서 분모는 표준오차 $\text{SE}(b_1)= s/\sqrt{S_{XX}}$이다.

기울기 추정량 $b_1$의 표준오차는

octave-3.4.0:88> sqrt(MSE/Sxx)
ans =  0.075305

귀무가설 $H_0: \beta_1 = 0$ 대 $H_1: \beta_1 \neq 0$을 검정하기 위하여 귀무가설 하에서 t-통계량을 계산하면

octave-3.4.0:89> t0 = b1 /sqrt(MSE/Sxx)
t0 =  26.649

이다. 유의확률이 매우 작으므로 또는 t0가 매우 크므로 귀무가설은 기각된다.

octave-3.4.0:99> 2 * (1 - tcdf(t0, n-2))
ans =  4.2286e-09
octave-3.4.0:100> tinv(0.975, n-2)
ans =  2.3060

[Octave] Matrix Algebra

Octave/MATLAB에서 행렬대수

기본 조작법

행렬 만들기

꺽쇠괄호로 묶고 행은 세미콜론, 열은 쉼표로 구분하여 입력한다. 열은 쉼표 없이 스페이스로 구분할 수도 있다. 행이나 열의 개수가 맞지 않으면 에러가 난다.

octave-3.4.0:67> A = [1, 2, 3; 4, 5, 6]
A =

   1   2   3
   4   5   6

일정한 간격의 숫자를 자동으로 생성할 때는 start:end 또는 start:step:end 형식의 문법을 이용할 수 있다.

octave-3.4.0:101> 1:10
ans =

    1    2    3    4    5    6    7    8    9   10

octave-3.4.0:102> 1:2:10
ans =

    1    3    5    7    9

0으로 채워진 행렬은 zeros(N, M), 1로 채워진 행렬은 ones(N, M)으로 생성할 수 있다.

octave-3.4.0:152> ones(4,1)
ans =

   1
   1
   1
   1

항등 행렬(identity matrix)는 eye(X)로 생성할 수 있고 대각원소 x를 가지는 대각행렬을 diag(x)로 생성할 수 있다.



octave-3.4.0:156> eye(3)

ans =



Diagonal Matrix



1   0   0

0   1   0

0   0   1



octave-3.4.0:157> diag([1 2 3])

ans =



Diagonal Matrix



1   0   0

0   2   0

0   0   3

인덱싱/슬라이싱

행렬 A에 대해 A(i,j)로 i행, j열의 원소를 가리킬 수 있다.

octave-3.4.0:92> A
A =

   1   2   3
   4   5   6

octave-3.4.0:90> A(2, 2)
ans =  5
octave-3.4.0:91> A(1, 2)
ans =  2

행렬 내의 일부 범위를 지정할 수도 있다.

octave-3.4.0:93> A(1:2, 1:2)
ans =

   1   2
   4   5

합치기: 행/열 묶기

행 벡터 x, y를 묶어 행렬을 만들 수 있다.

octave-3.4.0:116> x = 1:10
x =

    1    2    3    4    5    6    7    8    9   10

octave-3.4.0:117> y = 11:20
y =

   11   12   13   14   15   16   17   18   19   20

octave-3.4.0:118> [x; y]
ans =

    1    2    3    4    5    6    7    8    9   10
   11   12   13   14   15   16   17   18   19   20

열 벡터를 묶어 행렬을 만들 수도 있다.

octave-3.4.0:123> x = 1:3; y = 4:6; z = 7:9
z =

   7   8   9

octave-3.4.0:126> [x' y' z']
ans =

   1   4   7
   2   5   8
   3   6   9

두 행렬을 합칠 수도 있다.

octave-3.4.0:130> A
A =

   1   4   7
   2   5   8
   3   6   9

octave-3.4.0:131> B
B =

   11   12
   13   14
   15   16

octave-3.4.0:132> [A B]
ans =

    1    4    7   11   12
    2    5    8   13   14
    3    6    9   15   16

행렬의 크기

행렬 A에 대해 ndims(A) 명령을 내리면 2가 나온다. 2차원이기 때문이다. 컬럼의 수는 columns(A), 행의 수는 rows(A), 원소의 개수는 numel(A)로 알 수 있다. size(A)는 행과 열의 개수를 동시에 알려주고 length(A)는 행과 열 중 긴 것의 길이를 알려준다.

기본 연산

전치

실수 행렬 A의 전치(transpose)는 A' 또는 A.'로 가능하다.

octave-3.4.0:123> A = [1, 2, 3; 4, 5, 6]
A =

   1   2   3
   4   5   6

octave-3.4.0:124> A'
ans =

   1   4
   2   5
   3   6

참고로, 복소수 행렬 B의 경우에는 B.'는 전치, B'는 전치 켤레 행렬을 뜻한다.

행렬의 곱

두 행렬 A, B의 곱은 A*B로 계산할 수 있다.

octave-3.4.0:125> A = [1, 2, 3; 4, 5, 6]
A =

   1   2   3
   4   5   6

octave-3.4.0:126> B = [1, 1; 1, -1; -1, 1]
B =

   1   1
   1  -1
  -1   1

octave-3.4.0:127> A * B
ans =

   0   2
   3   5

행렬의 곱이 아니라 원소끼리의 곱(element by element multiplication)을 계산하려고 할 때는, A .* B와 같은 연산자를 이용한다.

octave-3.4.0:128> A = [1, 2, 3; 4, 5, 6]
A =

   1   2   3
   4   5   6

octave-3.4.0:129> B = [1, 0, 1; 0, 1, 0]
B =

   1   0   1
   0   1   0

octave-3.4.0:130> A .* B
ans =

   1   0   3
   0   5   0

트레이스

정방행렬 A의 대각원소의 합을 행렬 A의 트레이스(trace)라고 한다. trace(A)로 구할 수 있다.

octave-3.4.0:9> A = [1, 0, -1; 0, 2, 2; -1, 4, 5]
A =

   1   0  -1
   0   2   2
  -1   4   5

octave-3.4.0:10> trace(A)
ans =  8

행렬식

정방행렬 A의 행렬식(determinant)는 det(A) 또는 inverse(A)로 구할 수 있다.

octave-3.4.0:16> A = [1, 2; 3, 4]
A =

   1   2
   3   4

octave-3.4.0:17> det(A)
ans = -2

역행렬

정방행렬 A의 역행렬은 inv(A)로 구할 수 있다.

octave-3.4.0:16> A = [1, 2; 3, 4]
A =

   1   2
   3   4

octave-3.4.0:18> inv(A)
ans =

  -2.00000   1.00000
   1.50000  -0.50000

고유치와 고유벡터

정방행렬 A의 고유치(eigen value)와 고유벡터(eigenvalue)는 eig(A)로 구할 수 있다.

octave-3.4.0:30> A = [1, 4; 9, 1]
A =

   1   4
   9   1

octave-3.4.0:31> eig(A)
ans =

   7
  -5

octave-3.4.0:32> [v, lambda] = eig(A)
v =

   0.55470  -0.55470
   0.83205   0.83205

lambda =

Diagonal Matrix

   7   0
   0  -5

vec

행렬 $A$를 열벡터로 바꾸는 연산은 vec(A)로 할 수 있다.

octave-3.4.0:149> A
A =

   1   2   3
   4   5   7

octave-3.4.0:150> vec(A)
ans =

   1
   4
   2
   5
   3
   7

행렬의 유형

주어진 행렬 A의 유형을 matrix_type(A)를 이용하여 알 수 있다. unknown, full, positive definite, diagonal, permuted diagonal, upper, lower, banded, banded positive definite, singular로 구별해 준다.

octave-3.4.0:33> A = [2, 2, 1; 2, 5, 1; 1, 1, 2]
A =

   2   2   1
   2   5   1
   1   1   2

octave-3.4.0:34> matrix_type(A)
ans = Positive Definite

행렬의 분해

비정칙치 분해

행렬 $A$는
\[
A = U S V'
\]
로 분해될 수 있다. 이때 $U$, $V$는 직교행렬(orthogonal)이고 $S$는 대각행렬이다. $S$의 대각원소는 $A$의 비정칙치(singular values)라고 하며 이것은 $A'A$의 고유치의 양의 제곱근과 같다.

A =

   2   2   1
   2   5   1
   1   1   2

octave-3.4.0:180> [U, S, V] = svd(A)
U =

  -0.44910   0.29313  -0.84403
  -0.84403  -0.44910   0.29313
  -0.29313   0.84403   0.44910

S =

Diagonal Matrix

   6.41147         0         0
         0   1.81521         0
         0         0   0.77332

V =

  -0.44910   0.29313  -0.84403
  -0.84403  -0.44910   0.29313
  -0.29313   0.84403   0.44910

결과를 살펴보자. 우선, $U$와 $V$는 직교행렬이다.

octave-3.4.0:183> U' * U
ans =

   1.0000e+00  -1.9429e-16   1.6653e-16
  -1.9429e-16   1.0000e+00   2.2204e-16
   1.6653e-16   2.2204e-16   1.0000e+00

octave-3.4.0:184> V' * V
ans =

   1.00000  -0.00000   0.00000
  -0.00000   1.00000  -0.00000
   0.00000  -0.00000   1.00000

다음으로 $A=USV'$임을 확인해 보자.

octave-3.4.0:187> U * S * V'
ans =

   2.00000   2.00000   1.00000
   2.00000   5.00000   1.00000
   1.00000   1.00000   2.00000

다음으로 $A'A$의 고유치와 $S$의 대각원소의 제곱이 일치하는지 확인해 보자.

octave-3.4.0:193> eig(A' * A)
ans =

    0.59802
    3.29498
   41.10700

octave-3.4.0:192> diag(S) .** 2
ans =

   41.10700
    3.29498
    0.59802

촐레스키 분해

양정치 행렬(definite positive matrix) $A$는 항상
\[
A = R' R
\]
로 표현할 수 있으며 이때 $R$은 대각원소가 양수인 상삼각행렬(upper triangular matrix)이다. 이것을
촐레스키 분해(Cholesky decomposition)라고 하며 chol(A)로 계산할 수 있다.

octave-3.4.0:43> A
A =

   2   2   1
   2   5   1
   1   1   2

octave-3.4.0:46> matrix_type(A)
ans = Positive Definite
octave-3.4.0:44> R = chol(A)
R =

   1.41421   1.41421   0.70711
   0.00000   1.73205   0.00000
   0.00000   0.00000   1.22474

결과를 확인해 보자. 위에서 보듯이 R은 대각원소가 모두 양수이고 상삼각행렬이다. $R'R$을 계산해 보면 $A$를 얻을 수 있다.

octave-3.4.0:120> R' * R
ans =

   2.0000   2.0000   1.0000
   2.0000   5.0000   1.0000
   1.0000   1.0000   2.0000

QR 분해

임의의 실수 정방행렬 $A$의 경우에 항상
\[
A = QR
\]
로 나타낼 수 있다. 이때 $Q$는 직교행렬(orthogonal matrix)이다. 즉, $Q'Q=I$이다. $R$은 상삼각행렬(upper triangular matrix)이다. 이것을
QR decomposition이라고 하며 qr(A)로 계산할 수 있다.

octave-3.4.0:131> A = [2, 2, 1; 2, 5, 1; 1, 1, 2]
A =

   2   2   1
   2   5   1
   1   1   2

octave-3.4.0:132> [Q, R] = qr(A)
Q =

   0.66667  -0.59628  -0.44721
   0.66667   0.74536   0.00000
   0.33333  -0.29814   0.89443

R =

   3.00000   5.00000   2.00000
   0.00000   2.23607  -0.44721
   0.00000   0.00000   1.34164

결과를 확인해 보자. R은 위 결과에서 보듯이 상삼각행렬이고 다음에서 Q가 직교행렬인 것을 확인할 수 있다. 그리고 $QR =A$를 만족한다.

octave-3.4.0:134> Q' * Q
ans =

   1.00000  -0.00000   0.00000
  -0.00000   1.00000  -0.00000
   0.00000  -0.00000   1.00000

octave-3.4.0:140> Q * R
ans =

   2.00000   2.00000   1.00000
   2.00000   5.00000   1.00000
   1.00000   1.00000   2.00000

QR 분해는 정방행렬이 아닌 경우에도 가능하다.

2011년 12월 16일 금요일

R 속도비교: loop, apply, vectorization

이 실험은 R version 2.14.0 (2011-10-31)에서 이루어진 것이다.

R version 2.14.0 (2011-10-31)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

컴퓨터는 iMac 3.06GHz Intel Core 2 Duo, 메모리는 12GB 1067 MHz DDR3. 운영체제는 Mac OS X Lion 10.7.2이다.

loop를 돌릴 때

우선, 결과값을 저장할 벡터를 미리 크기를 정해서 만들어 놓는 것과 아닌 것에는 엄청난 차이가 난다.

> x <- 1:(10^5)
> y <- c()
> system.time(for(i in 1:length(x)) y[i] <- x[i]^2)
   user  system elapsed 
 13.054  33.554  65.792
> y <- numeric(length(x))
> system.time(for(i in 1:length(x)) y[i] <- x[i]^2)
   user  system elapsed 
  0.312   0.008   0.455

loop를 돌리며 결과를 벡터에 저장하려고 할 때는 반드시 해당 벡터를 먼저 초기화해야 한다.

loop, apply, vectorization

이제 계산할 벡터의 크기를 10배 늘려서 loop, apply, vectorization을 비교해 보자. apply 함수를 사용할 때는 결과값이 저장될 벡터를 미리 초기화할 필요가 없다. vectorization의 속도가 엄청나게 빠른 것을 확인할 수 있다.

> x <- 1:(10^6)
> y <- numeric(length(x))
> system.time(for(i in 1:length(x)) y[i] <- x[i]^2)
   user  system elapsed 
  3.019   0.045   4.533 
> rm(y)
> system.time(y <- sapply(x, function(i) i^2))
   user  system elapsed 
  3.248   0.071   4.605 
> rm(y)
> system.time(y <- x^2)
   user  system elapsed 
  0.006   0.000   0.006

흔히 for loop보다 apply()가 "훨씬" 빠르다고 이야기하는 경우가 있는데 그건 사실이 아니다. 혹시 S-Plus에서는 그렇다는 말이 있는데 사실인지 모르겠다. R에서는 루프도 apply만큼 빠르다. 단 결과를 저장할 벡터를 초기화하는 것을 잊어서는 안 된다. 둘의 속도 차이는 경우에 따라 없을 수도 있고 루프가 빠를 수도 있고 apply가 빠를 수도 있다. 물론, apply()를 쓰는 것이 코드가 더 깔끔하다.

Vectorization의 엄청난 속도

경우에 따라 다르겠지만 원소들의 제곱으로 이루어진 벡터를 만드는 바로 이 문제의 경우, 루프나 apply와 비교했을 때 vectoriation이 500배 이상 빠른 속도를 보이는 것을 확인할 수 있다. 벡터의 크기를 늘려가며 vectorization을 해 보았다. 1천만(10^8)개까지는 큰 문제가 없이 산술적으로 속도가 증가한다. 1억개(10^9)부터는 연산에 걸린 시간은 정상적으로 증가했지만 elapsed가 갑자기 시간이 오래 걸렸는데 메모리 부족 문제인 듯하다.

> x <- 1:(10^6)
> system.time(y <- x^2)
   user  system elapsed 
  0.005   0.000   0.005 
> x <- 1:(10^7)
> system.time(y <- x^2)
   user  system elapsed 
  0.059   0.001   0.065 
> x <- 1:(10^8)
> system.time(y <- x^2)
   user  system elapsed 
  0.575   0.509   1.246 
> x <- 1:(10^9)
> system.time(y <- x^2)
   user  system elapsed 
  6.266   8.735 155.847 
> x <- 1:(10^10)
이하에 에러1:(10^10) : 결과의 벡터가 너무 깁니다

마지막 에러가 난 것은 벡터의 크기가 너무 크기 때문이다. R에서 허용하는 벡터의 최대 크기는 2^31 - 1, 약 2*10^9이다.

`matrix, array의 경우: rowSums()`, `colSums()` 사용하기

vectoriation을 모든 경우에 할 수 있는 것은 아니다. 행렬(matrix)이나 배열(array)에서 행이나 열의 합을 구하는 경우에 루프나 apply보다는 rowSums(), colSums()를 쓰는 것이 빠르다.

> x <- matrix(runif(10^8), ncol=10^4)
> str(x)
 num [1:10000, 1:10000] 0.3492 0.0621 0.2593 0.8765 0.5281 ...
> y <- numeric(10000)
> system.time(for(i in 1:10000) y[i] <- sum(x[i,]))
   user  system elapsed 
  6.836   0.755   9.653 
> rm(y)
> system.time(y <- apply(x, 1, sum))
   user  system elapsed 
  6.597   1.512  10.263 
> rm(y)
> system.time(y <- rowSums(x))
   user  system elapsed 
  0.263   0.003   0.326

list의 경우: `lapply()`

list는 vector가 아니라 list이니까 Vectorizatioin하다는 것은 이름 그대로 말이 안 되지 않는가? 루프나 apply를 써야 하는데, list를 대상으로 할 때는 루프보다 lapply()의 코드가 효율적이라고 알려져 있고 실제로도 조금 빠른 것 같다.

> x <- list()
> for(i in 1:10000) x[[i]] <- runif(10000)
> y <- numeric(10000)
> system.time(for(i in 1:10000) y[i] <- sum(x[[i]]))
   user  system elapsed 
  0.257   0.004   0.390 
> system.time(y <- lapply(x, sum))
   user  system elapsed 
  0.228   0.005   0.281

약간 다르게:

> x <- list()
> for(i in 1:100000) x[[i]] <- runif(10)
> y <- numeric(100000)
> system.time(for(i in 1:100000) y[i] <- sum(x[[i]]))
   user  system elapsed 
  0.393   0.007   0.543 
> system.time(y <- lapply(x, sum))
   user  system elapsed 
  0.147   0.003   0.229

누적 계산이 필요한 경우

> x <- runif(1000000)
> system.time(y <- cumsum(x))
   user  system elapsed 
  0.032   0.006   0.102 
> system.time({y[1] <- x[1]; for(i in 2:length(x)) y[i] <- y[i-1] + x[i]})
   user  system elapsed 
  4.140   0.029   4.330

2011년 12월 13일 화요일

MCMC

Markov Chain Monte Carlo methods

S. Chib and E. Greenberg (1995). Understanding the Metropolis-Hastings Algorithm. The American Statistician, 49-4:327-335.
G. Casella and E. I. George (1992). Explaining the Gibbs Sampler. The American Statistician, 46-3:167-174.
G. Cassella, M. Lavine, and C. P. Robert (2001). Explaining the Perfect Sampler. The American Statistician, 55-4:299-305.

Markov Chain Monte Carlo and Gibbs Sampling by B. Walsh (2004)
MCMC Sampling (slide) by T. Bahadori (2011)

2011년 12월 6일 화요일

MA(q) 모형 적합

시계열. MA(q) 모형. 최소제곱추정법.

MA(1) 모형을 상수항 없이 \[ x_t = \epsilon_t - \theta \epsilon_{t-1}, \quad \epsilon_t \sim \textit{WN}(0, \sigma^2) \] 이라고 하자. \[ \epsilon_t = x_t + \theta \epsilon_{t-1} \] 이므로 최소제곱법으로 추정량 $\hat\theta$를 구하기 위해 \[ S(\theta) = \sum_{t=1}^n \epsilon_t^2 = \sum_{t=1}^n (x_t - \theta\epsilon_{t-1})^2 \] 를 최소로 하는 $\theta$를 찾을 수 있다. 그런데 $\epsilon_t$가 $\theta$에 따라 달라지므로 위 식은 간단히 풀릴 수가 없다. $\epsilon_t$가 $\theta$의 함수라는 사실을 드러내기 위해 $\epsilon_t(\theta)$로 표기하자. 이 함수가 어떤 함수인지 모르지만 테일러 전개를 해서 1차 함수로 근사하면 위 식을 쉽게 풀 수 있을 것이다. $\epsilon_t(\theta)$를 어떤 상수 $\theta^*$에서 테일러 전개를 1차항까지만 하면 다음과 같다. \[ \epsilon_t(\theta) \approx \epsilon_t(\theta^*) + (\theta - \theta^*) \omega_t(\theta^*), \qquad \omega_t(\theta) = \frac{\partial \epsilon_t(\theta)}{\partial \theta} \] 이제 이것을 이용하며 $S(\theta)$를 나타내면 \[ S(\theta) = \sum_{t=1}^n [\epsilon_t(\theta^*) + (\theta-\theta^*)\omega_t(\theta^*)]^2 \] 가 된다. 여기서 $\theta^*$은 어떤 상수이고 따라서 $\theta$에 관한 1차식을 제곱한 것의 합에 불과하므로 $S(\theta)$를 최소로 하는 $\theta$를 쉽게 찾아낼 수 있다. 위 식은 $ S(\beta) = \sum_{t=1}^n (y_t - \beta x_t)^2 $ 과 동일한 형태이고 우리는 이때 $S(\beta)$를 최소로 하는 것이 $\beta = S_{XY}/S_{XX}$라는 것을 잘 알고 있다. 즉, \[ \theta - \theta^* = \frac{-\sum_{t=1}^n \epsilon_t(\theta^*)\omega_t(\theta^*)}{\sum_{t=1}^n[\omega_t(\theta^*)]^2} \] 이다. 이 식에서 $\theta^*$에 어떤 초기값 $\theta_{(0)}$을 넣고 계산할 수 있다. 초기값은 예를 들어 적률이용추정량을 사용할 수 있다. 이렇게 해서 얻어지는 $\theta$를 $\theta_{(1)}$라고 하고 이 과정을 계속 반복할 수 있다. 즉,
\[ \theta_{(j+1)} = \theta_{(j)} +
\frac{-\sum_{t=1}^n \epsilon_t(\theta_{(j)})\omega_t(\theta_{(j)})}{\sum_{t=1}^n[\omega_t(\theta_{(j)})]^2} \]
와 같이 충분하다고 생각될 때까지 반복 계산을 하여 $\hat\theta$를 얻을 수 있다.

이때 $\epsilon_t(\theta^*)$, $\omega_t(\theta^*)$를 어떻게 구하는지가 문제가 된다. 우선 $\epsilon_0(\theta)=0$이라고 가정하고 $\epsilon_t(\theta)$를 살펴보자. \begin{align*} \epsilon_t(\theta) &= x_t + \theta \epsilon_{t-1}(\theta) = x_t + \theta (x_{t-1} + \theta \epsilon_{t-2}) = \cdots \\ &= x_t + \theta x_{t-1} + \theta^2 x_{t-2} + \cdots + \theta^{t-1} x_1 + \theta^t \epsilon_0(\theta)\\ &=\sum_{k=0}^{t-1} \theta^k x_{t-k} \end{align*} 여기에서 $x_i$는 관측치이므로 어떤 상수 $\theta^*$에 대해 $\epsilon(\theta^*)$를 구할 수 있다. 그리고 $\epsilon_t(\theta) = x_t + \theta\epsilon_{t-1}(\theta)$의 양변을 미분하여 얻어지는 식 또한 점화식의 형태가 되어 동일한 방법으로 계산할 수 있다. \begin{align*} \frac{\partial\epsilon_t(\theta)}{\partial \theta} &= \epsilon_{t-1}(\theta) + \theta \frac{\partial \epsilon_{t-1}(\theta)}{\partial\theta} = \cdots \\ &= \epsilon_{t-1}(\theta) + \theta \epsilon_{t-2}(\theta) + \theta^2 \epsilon_{t-3}(\theta) + \cdots + \theta^{t-2}\epsilon_1(\theta) +\theta^{t-1}\frac{\partial\epsilon_0(\theta)}{\partial \theta}\\ &= \sum_{k=0}^{t-2} \theta^k \epsilon_{t-1-k} \end{align*}

2011년 9월 17일 토요일

책: Bootstrap, Resampling, Subsampling, MCMC

책: Bootstrap, Resampling, Subsamplling, MCMC

An introduction to the bootstrap by Efron and Tibshirani. 1993. Chapman & Hall.
Bootstrap Methods and their Application by Davison and Hinkley. 1997. CUP
Resampling Methods: A Practical Guide to Data Analysis (3rd ed.) by Phillip I. Good. 2006. Birkhauser Boston. [ebook]
Permutation, Parametric, and Bootstrap Tests of Hypothesis by Pillip I. Good. 2005. Springer. [ebook]
The Bootstrap and Edgeworth Expansion by Peter Hall. 1995. Springer.
The Jackknife, the bootstrap, and Other Resampling Plans by Bradley Efron. 1987. SIM.
The Jackknife and Bootstrap by Jun Shao. 1995. Springer.
Resampling Methods for Dependent Data by S. N. Lahiri. 2003. Springer
Subsampling by Politis, Romano, and Wolf. 1999. Springer.
Monte Carlo Statistical Methods by Robert and Casella. 2010. Springer.
Introducing Monte Carlo Methods with R by Rober and Casella. 2009. Springer. [ebook]

Read:

Bootstrap Methods and Permutation Tests by Hesterberg et al.
The bootstrap and Markov chain Monte Carlo by Bradley Efron.
An Introduction to the Bootstrap with Application in R by Davison and Kuonen.
Bootstrapping Regression Models by John Fox.

2011년 7월 13일 수요일

R에서 생존분석: Kaplan-Meier 방법

준비

생존 분석 패키지를 다음 명령으로 읽어들인다.

> require(survival)

급성 골수성 백혈병(Acute Myelogenous Leukemia) 데이터가 예제로 들어있다. 작은 데이터이다. 내용을 보면:

> aml
   time status             x
1     9      1    Maintained
2    13      1    Maintained
3    13      0    Maintained
4    18      1    Maintained
5    23      1    Maintained
6    28      0    Maintained
7    31      1    Maintained
8    34      1    Maintained
9    45      0    Maintained
10   48      1    Maintained
11  161      0    Maintained
12    5      1 Nonmaintained
13    5      1 Nonmaintained
14    8      1 Nonmaintained
15    8      1 Nonmaintained
16   12      1 Nonmaintained
17   16      0 Nonmaintained
18   23      1 Nonmaintained
19   27      1 Nonmaintained
20   30      1 Nonmaintained
21   33      1 Nonmaintained
22   43      1 Nonmaintained
23   45      1 Nonmaintained

이 데이터는 단순한 데이터프레임이다. 확인해 보자.

> class(aml)
[1] "data.frame"
> str(aml)
'data.frame': 23 obs. of  3 variables:
 $ time  : num  9 13 13 18 23 28 31 34 45 48 ...
 $ status: num  1 1 0 1 1 0 1 1 0 1 ...
 $ x     : Factor w/ 2 levels "Maintained","Nonmaintained": 1 1 1 1 1 1 1 1 1 1 ...

즉, 단순한 테이블 데이터라는 뜻의다. 아직 특별히 생존 분석을 위한 가공이 이루어지 않은 것이다. 위와 같은 형태로만 데이터가 입력되어 있으면 된다.

Survival Object

생존 분석을 하기 위해 데이터를 Survival Object로 만들어 보자. 이 객체는 survival 패키지에서 제공하는 것이다. Survival Object는 Surv()를 이용하여 만든다. 이 함수는 기본적으로 Surv(time, event)로 사용할 수 있다. 두 변수가 필요하다. time은 생존시간이고 event는 상태(status)를 나타내는 변수이다. 상태는 기본적으로

0=right censored
1=event at time
2=left censored
3=interval censored

로 코딩한다. 지금 주어진 aml데이터는 status 변수에 0/1밖에 없다. 다음과 같이 Survival Object를 만들면 status 값이 0인 데이터 즉 right censored 데이터인 경우에 시간 뒤에 +가 붙어 표시되는 것을 확인할 수 있다.

> (aml.Surv <- with(aml, Surv(time, status)))
 [1]   9   13   13+  18   23   28+  31   34   45+  48  161+   5    5    8    8   12   16+  23   27   30   33   43 
[23]  45

Kaplan-Meier method

Kaplan-Meier 방법으로 중도 절단이 있는 자료(censored data)의 생존 함수를 추정하는 것은 survfit()로 할 수 있다. 앞서 만든 Survival Object를 이용해서

> aml.survfit <- survfit(aml.Surv ~ aml$x)

라고 할 수도 있지만,

> aml.survfit <- survfit(Surv(time, status) ~ x, aml)

라고 하는 것이 보기 좋다. 분석 자료 개요는 다음과 같다.

> aml.survfit
Call: survfit(formula = Surv(time, status) ~ x, data = aml)

                records n.max n.start events median 0.95LCL 0.95UCL
x=Maintained         11    11      11      7     31      18      NA
x=Nonmaintained      12    12      12     11     23       8      NA

분석 결과를 보려면 summary() 함수를 이용하면 된다. 결과는 다음과 같다.

> summary(aml.survfit)
Call: survfit(formula = Surv(time, status) ~ x, data = aml)

                x=Maintained 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    9     11       1    0.909  0.0867       0.7541        1.000
   13     10       1    0.818  0.1163       0.6192        1.000
   18      8       1    0.716  0.1397       0.4884        1.000
   23      7       1    0.614  0.1526       0.3769        0.999
   31      5       1    0.491  0.1642       0.2549        0.946
   34      4       1    0.368  0.1627       0.1549        0.875
   48      2       1    0.184  0.1535       0.0359        0.944

                x=Nonmaintained 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    5     12       2   0.8333  0.1076       0.6470        1.000
    8     10       2   0.6667  0.1361       0.4468        0.995
   12      8       1   0.5833  0.1423       0.3616        0.941
   23      6       1   0.4861  0.1481       0.2675        0.883
   27      5       1   0.3889  0.1470       0.1854        0.816
   30      4       1   0.2917  0.1387       0.1148        0.741
   33      3       1   0.1944  0.1219       0.0569        0.664
   43      2       1   0.0972  0.0919       0.0153        0.620
   45      1       1   0.0000     NaN           NA           NA

생존 곡선 그리기

생존 곡선은 plot() 명령으로 간단히 그릴 수 있다.

> plot(aml.survfit)

위 데이터에서 이미 보았듯이 aml$x를 보면

> levels(aml$x)
[1] "Maintained"    "Nonmaintained"

두 개의 집단으로 이루어져 있다. 두 생존곡선을 다른 형태의 선으로 그려서 구별하고 싶다면

> plot(aml.survfit, lty=1:2)

라고 lty (line type) 옵션을 주는 것만으로도 충분하다. 범례도 달고 싶다면 일반적으로는 다음과 같이 하면 된다.

K <- length(aml.survfit$strata)
LEGEND <- attr(aml.survfit$strata, "names")
plot(aml.survfit, lty=1:K)
legend("topright", LEGEND, lty=1:K)

물론, 여기서 K나 LEGEND를 직접 타이핑해도 된다. 직접 써줄 때는 범례 이름을 바꾸어 쓰는 일이 없도록 주의해야 한다.

plot(aml.survfit, lty=1:2)
legend("topright", c("Maintained", "Nonmaintained"), lty=1:2)

2011년 7월 12일 화요일

생존 분석 기본 개념

생존 분석 기본 개념
생존 분석(survival analysis)이라고 하면 주로 생명체의 죽음 또는 기계 장치의 고장에 관련된 통계 분석이지만 이 분야에 한정된 것은 아니다. 공학 분야에서 reliability analysis, failure time analysis, 경제학 분야에서 duration analysis, transition analysis, 사회학 분야에서 event history analysis 등의 이름으로 불리며 이들 사이에 실질적인 분석 방법의 차이는 없다. 질병의 발생이나 기계의 오작동 뿐만 아니라 지진의 발생, 자동차 사고 발생, 주가의 추락, 등을 분석하는 데에 사용된다[1].
생존 분석을 더 일반적인 개념을 사용하여 'time to event' 자료의 분석이라고도 하는데 사망, 고장, 사고, 폭락 등을 사건(event)이라고 보고 해당 사건이 일어날 때까지 걸린 시간을 데이터로 다룬다는 의미이다. 이제 이런 사건을 생존 분석에서는 사망(death)라고 지칭하는 경우가 많고 좀더 일반적으로는 고장/실패(failure)로 지칭하는 경우가 많다. 하지만 반드시 '사망'이거나 '실패'일 필요는 없다. 무엇이든 우리가 관심을 가지는 사건이면 된다.

생존 시간

우리는 사건이 일어나는 시간에 관심이 있다. 이 사건이 발생하는 시간을 $T$라고 하자. 이것은 확률 변수이다. 생존 분석에서는 이것을 생존 시간(survival time)이라고 한다. 시간이 항상 0에서부터 시작한다고 보면 사건이 일어나는 시점 $T$는 사건이 일어나기까지 걸린 시간이기도 하다.

생존 함수

생존 함수(survival function 또는 suvivor function) $S(t)$는 생존 시간 $T$가 특정 시간 $t$ 보다 클 확률을 가리킨다. $T$의 cdf를 $F(t)$라고 하면 다음 관계가 성립한다. \[ S(t) = P(T > t) = 1 - F(t) \] 생존 함수 $S(t)$는 $t$ 시점에서 살아있을 또는 $t$ 시점까지 살아남을 확률이다. 이때 $t$ 시점까지 살아남는다는 것은 $t$ 시점까지만 살고 바로 그때 죽는다는 뜻이 아니다. $t$ 시점 이후에 죽는다는 뜻이다. 더 일반적으로 말로 표현하면 $t$ 시점까지 관심 사건이 일어나지 않을 확률이다.

위험 함수

위험 함수(hazard function) $h(t)$는 주어진 시간 $t$에서 단위 시간 당 사건 발생율이다. \[ h(t) = \lim_{\Delta t \to 0} \frac{P(t \leq T < t+ \Delta t |T \geq t)}{\Delta t} \] 사건이 일어날 순간 가능성(instantaneous potential)이라고도 할 수 있다. 위험 함수를 조건부 실패율(conditional failure rate)라고도 한다.

함수들 사이의 관계

생존 함수 $S(t) = P(T>t)$와 생존 시간 $T$의 cdf인 $F(t)$사이에는 \[ S(t) = 1 - F(t) \] 관계가 성립하므로 양변을 $t$에 관해 미분하면 \[ \frac{d}{dt}S(t) = - f(t) \] 와 같은 관계를 얻을 수 있다. 여기에서 $f(t)$는 $T$의 pdf이다.

위험 함수 $h(t)$의 정의를 따라 식을 바꾸어 보면 \begin{align*} h(t) &= \lim_{\Delta t \to 0} \frac{P(t \leq T < t+ \Delta t |T \geq t)}{\Delta t}\\ &= \lim_{\Delta t \to 0} \frac{P(t \leq T < t+\Delta t)}{P(T \geq t)}\frac{1}{\Delta t}\\ & = \frac{1}{P(T \geq t)}\lim_{\Delta t \to 0} \frac{P(t\leq T < t+\Delta t)}{\Delta t}\\ & = \frac{1}{P(T \geq t)}f(t)\\ & = \frac{f(t)}{S(t)}\\ & = - \frac{d}{dt}\ln S(t) \end{align*}

누적 위험 함수

누적 위험 함수 $H(t)$는 다음과 같다. \[ H(t) = \int_0^t h(u)du = - \ln S(t) \] 따라서 \[ S(t) = e^{H(t)} \] 이다.

참고

Survival Analysis Using SAS: A Practical Guide

2011년 7월 6일 수요일

R 객체 지향 프로그래밍: S4 기초

R OOP S4 기초

먼저 S3에 대해서는 여기 참조. S4 클래스로 프로그래밍을 하기 위해서 기본적으로 다음 함수들을 시용해야 한다.

setClass()
new()
setGeneric()
setMethod()

클래스

새로운 클래스는 다음처럼 setClass() 함수를 이용하여 정의한다.

setClass("circle", 
  representation(x="numeric", y="numeric", r="numeric"))

이제 new()를 이용하여 인스턴스를 만들 수 있다. 객체의 속성을 S4에서늘 슬롯(slot)이라고 한다. 슬롯은 object@slot 형태로 접근할 수 있다.

> a <- new("circle", x=5, y=10, r=4)
> str(a)
Formal class 'circle' [package ".GlobalEnv"] with 3 slots
  ..@ x: num 5
  ..@ y: num 10
  ..@ r: num 4
> a@x
[1] 5
> a@y
[1] 10
> a@r
[1] 4

한 객체가 특정 클래스에 속하는지 체크할 때는 is()를 이용한다.

> is(a, "circle")
[1] TRUE

메소드

메소드는 setMethod()를 이용하여 정의한다. 기존의 정의된 generic function이 없다면 우선 generic function부터 만들어야 하는데 이때에는 setGeneric()을 이용한다. 예를 들면 다음처럼 한다.

setGeneric("area", function(object) standardGeneric("area"))

setMethod("area", "circle",
          function(object) pi * object@r^2
          )

이제 메소드를 사용해 보면

> area(a)
[1] 50.26548
> area(new("circle", x=3, y=4, r=3))
[1] 28.27433

한 generic function에 어떤 method들이 있는지 알고 싶을 때는 showMethods()를 이용한다.

> showMethods("area")
Function: area (package .GlobalEnv)
object="circle"

참조

http://www.soph.uab.edu/Statgenetics/Events/Rshort/060227-8-s4slides.pdf
John M. Chambers (2008), Software for data analysis: programming with R, Springer.

블로그 보관함