-
Notifications
You must be signed in to change notification settings - Fork 0
/
hw4.Rmd
182 lines (142 loc) · 8.88 KB
/
hw4.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
---
title: ""
output: html_notebook
---
#### Практическая работа №4
### Метод наименьших квадратов в задаче линейной и нелинейной регрессии
> Глушков Егор Александрович, гр. 20.М04-мм
> Вариант 4
---
1. Промоделировать нелинейную модель $y = f(x, a, b) + \delta$ с несмещённой нормально распределенной ошибкой, дисперсия которой равна $\varepsilon$, считая, $x$ стандартно нормально распределённой случайной величиной.
* $f(x, a, b) = a \sin x + bx^2, \; a=2, b=1, \varepsilon=0.8$
2. Оценить параметры нелинейной модели по методу наименьших квадратов (численно). Применить к модельным данным линейную модель и оценить параметры. Построить на двумерной диаграмме основную и линейную модель. Сравнить невязки для обеих моделей.
3. Для линейной модели выполнить дисперсионный анализ, проверить значимость прогноза и коэффициентов регрессии. Сравнить непосредственные вычисления с результатами встроенной функции.
4. Промоделировать данные для множественной регрессии. Применить функцию $lm$. Ответить на вопросы о значимости коэффициента детерминации, частных коэффициентов регрессии, о коэффициенте корреляции между остатком и независимыми переменными.
---
> Задание линейной и нелинейной моделей регрессии
Нелинейная модель
$$y_i = a \sin x_i + bx_i^2 + \delta_i, \quad \delta_i \sim \mathcal{N}(0, \sigma)$$
Модель одномерной линейной регрессии
$$y_i = \alpha + \beta x_i, \quad \delta_i \sim \mathcal{N}(0, \sigma)$$
Задаем нелинейную и линейную модели и их функции ошибки *L*
```{r}
N <- 100
f <- function(x, ab) ab[1] * sin(x) + ab[2] * x^2
L <- function(X, Y, ab) sum((Y - f(X, ab))^2)
f_lin <- function(x, ab_lin) ab_lin[1] + ab_lin[2] * x
L_lin <- function(X, Y, ab_lin) sum((Y - f_lin(X, ab_lin))^2)
```
Моделируем данные нелинейной модели
```{r}
ab <- c(2, 1)
eps <- 0.8
X <- rnorm(N)
Y <- f(X, ab) + rnorm(N, 0, eps^2)
```
> Оценка параметров линейной модели и наилучший прогноз
$$ \hat{\beta} = \frac{\sum\nolimits_i x_i y_i - n \bar{x} \bar{y}}{{\sum\nolimits_i x_i^2 - n \bar{x}^2}}, \quad \hat{\alpha} = \bar{y} - \hat{\beta} \bar{x}, \quad \hat{y_i} = \hat{\alpha} - \hat{\beta} x_i$$
```{r}
b. <- (sum(X * Y) - N * mean(X) * mean(Y)) / (sum(X * X) - N * mean(X)^2)
a. <- mean(Y) - b. * mean(X)
ab_lin_est <- c(a., b.)
Y. <- f_lin(X, ab_lin_est)
ab_lin_est
```
> Дисперсионный анализ
Источники вариации: общий $Q_T$, обусловленный регрессией $Q_R$, невязка $Q_E$. Коэффициент детерминации $R^2$
```{r}
QT <- sum((Y - mean(Y))^2)
QR <- sum((Y. - mean(Y))^2)
QE <- sum((Y - Y.)^2)
R2 <- QR / QT
c(QT=QT, QR=QR, QE=QE, R2=R2)
```
`r round(R2, 3)*100`% дисперсии объясняется данной моделью -- линией регресии.
Проверим равенство источников вариации [посчитаны верно]:
```{r}
print(paste("(QR + QE == QT):", QR + QE - QT < 1e-9))
```
Дисперсии общая и коэффициентов регрессии
```{r}
xx <- sum((X - mean(X))^2)
S2 <- QE / (N - 2)
S2.a <- S2 * sum(X^2) / (xx * N)
S2.b <- S2 / xx
c(S2, S2.a, S2.b)
```
Статистики для проверки значимости прогноза и коэффициентов регрессии
```{r}
F_stat <- QR / QE * (N - 2)
Pf <- 1 - pf(F_stat, 1, N - 2)
T.a <- ab_lin_est[1] / sqrt(S2.a)
T.b <- ab_lin_est[2] / sqrt(S2.b)
P.a <- 2*(1 - pt(abs(T.a), N - 2))
P.b <- 2*(1 - pt(abs(T.b), N - 2))
P.model <- 1 - pf(F_stat, 1, N - 2)
```
Выведем самостоятельно посчитанные величины:
```{r}
show_stats_coef <- rbind(c(ab_lin_est[1], sqrt(S2.a), T.a, P.a), c(ab_lin_est[2], sqrt(S2.b), T.b, P.b))
colnames(show_stats_coef) <- c("Estimate", "Std.Error","t value", "Pr(>|t|)")
rownames(show_stats_coef) <- c("(Intercept)", "X"); show_stats_coef
show_stats_model <- c(sd=sqrt(S2), R2=R2, F_stats=F_stat, p.value.F=P.model); show_stats_model
```
Проверим при помощи встроенной функции *lm*:
```{r}
summary(lm(Y~X))
```
Как видно, все посчитанные величины совпадают. Вычисления были верны.
Что касается самой модели, то итоговое *p-value* (близкое к нулю) говорит о значимости регрессии, также малые *Pr* говорят о том, что соответствующие этим уровням значимости коэффициенты отличны от нуля и вносят вклад в модель. $R^2$ сообщает, что `r round(R2, 3)*100`% дисперсии объясняется данной моделью -- линией регресии.
> Оценка параметров нелинейной модели по МНК
Параметры нелинейной модели оцениваем методом наименьших квадратов (МНК)
```{r}
NLM <- nlm(function(ab) L(X, Y, ab), ab)
ab. <- NLM$estimate
rbind(ab=ab, ab.=ab.)
```
```{r}
plot(X, Y)
f_ <- function(x) f(x, ab); curve(f_, -3, 3, add=TRUE, col=2)
f_ <- function(x) f(x, ab.); curve(f_, -3, 3, add=TRUE, col=3)
f_ <- function(x) f_lin(x, ab_lin_est); curve(f_, -3, 3, add=TRUE, col=4, lty=2)
legend('bottomright', c('hyp', 'mnk', 'linear'), pch=20, col=2:4)
```
Невязки
```{r}
c(Q.linear=L_lin(X, Y, ab_lin_est), Q.hyp=L(X, Y, ab), Q.mnk=L(X, Y, ab.))
```
Заметим, что линейная модель оценивает достаточно грубо, тогда как нелинейная модель достаточно точно описывает данные. Конечно, меньшей невязкой как суммы квадратов отклонений точек от графика обладает модель, построенная по МНК.
> Множественная регрессия
```{r}
N <- 100
a <- c(2, -1, 5)
eps <- 3
X1 <- rnorm(N, 1, 1.5)
X2 <- rnorm(N, 2, 0.5)
Y <- a[1] * X1 + a[2] * X2 + a[3] + rnorm(N, 0, eps^2)
SLM <- summary(lm(Y ~ X1 + X2)); SLM
```
Согласно *p-value* < 0.05 отвергаем гипотезу о том, что данная модель не зависит от предикторов (т.е. их коэффициенты равны нулю -- отвергаем). Отвергается и незначимость коэффициента при *X1*, а вот *X2*, как следует, не вносит весомый вклад ($p.value = `r round(SLM$coefficients[,4][3], 3)` > 0.05$, то не отвергается факт, коэффициент незначим)
Проверим некоррелированность остатков [практически равны 0]
```{r}
c(cor(X1, SLM$residuals), cor(X1, SLM$residuals))
```
```{r}
op <- par(mfrow = c(2, 2))
plot(X1, SLM$residuals)
title(sub=paste("r", round(cor(SLM$residuals, X1), 3), sep="="))
plot(X2, SLM$residuals)
title(sub=paste("r", round(cor(SLM$residuals, X2), 3), sep="="))
plot(X1, Y)
title(sub=paste("r", round(cor(X1, Y),3), sep="="))
plot(X2, Y)
title(sub=paste("r", round(cor(X2, Y), 3), sep="="))
```
Как видно на верхних двух рисунках, коэффициент корреляции остатков с независимыми переменными равен нулю, тогда как в целом некоторая зависимость (линейная) в данных наблюдается (см. нижние рисунки).
```{r}
sh_t <- shapiro.test(SLM$residuals); sh_t
```
Согласно критерию Шапиро-Уилка, нормальность остатков не отвергается с $p.value = `r round(sh_t$p.value, 3)`$
```{r}
qqnorm(SLM$residuals)
```