모의실험: 난수생성
단순선형회귀의 간단한 예를 모의생성해 보자.\[
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.00683Octave에서는 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
댓글 없음:
댓글 쓰기