Harlequin Blue

2012년 2월 7일 화요일

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

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일 수요일

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

블로그 보관함