Лабораторная работа №6

Задание 1. Коофициент корреляции цены рубля и нефти в американских долларах.

rm(list=ls()) # clears out memory

usd.data="./QUANDL_USDRUB.csv"
oil.data="./QUANDL_OPECOIL.csv"

url_usd="https://www.quandl.com/api/v1/datasets/CURRFX/USDRUB.csv?trim_start=2014-08-01"
download.file(url_usd, usd.data, "wget")

url_oil="https://www.quandl.com/api/v1/datasets/OPEC/ORB.csv?trim_start=2014-08-01"
download.file(url_oil, oil.data, "wget")


usd.data <- read.csv(usd.data, header = TRUE, sep = ",")
oil.data <- read.csv(oil.data, header = TRUE, sep = ",")

merged.data <- merge(usd.data,oil.data, by="Date")
rub_vs_oil = data.frame(DATE=merged.data$Date, 
                        RUB=merged.data$Rate,
                        OIL=merged.data$Value)

Расчитаем коэффициент корреляции по формуле

\[ r_{xy} = \frac{\sum(X - \widetilde{X})(Y - \widetilde{Y})}{\sqrt{\sum(X - \widetilde{X})^2}\sqrt{\sum(Y - \widetilde{Y})^2}} \]

rub_avg = weighted.mean(rub_vs_oil$RUB)
oil_avg = weighted.mean(rub_vs_oil$OIL)
x = rub_vs_oil$RUB - rub_avg
y = rub_vs_oil$OIL - oil_avg
r_xy = (t(x) %*% y) / (sqrt(t(x)%*%x) * sqrt(t(y)%*%y))
r_xy

           [,1]
[1,] -0.9610832

Вывод: выборочне величины имеют явно выраженную линейную зависимость

Задание 2. Регрессионный анализ цены рубля и нефти

Начнем с модели линейной регрессии \(y^* = a + bx\). Вычислим \(b\) и \(a\) по формулам

\[ b = r_{xy} \frac{s_{y}}{s_{x}} \] \[ a = E[Y] - bE[X] \]

b = r_xy * sqrt( t(y)%*%y / t(x)%*%x)
a = weighted.mean(rub_vs_oil$OIL) - b * weighted.mean(rub_vs_oil$RUB)

a = 152.11 b = -1.67

Отобразим полученную линию регресии на графике

Интересно сравнить результаты полученные врчучную с пакетным методом. Видим, что они совпадают.

lm(formula = OIL ~ RUB, data = rub_vs_oil)

Call:
lm(formula = OIL ~ RUB, data = rub_vs_oil)

Coefficients:
(Intercept)          RUB  
    152.111       -1.669  

Посмотрим на гистограмму остатков и вычислим их дисперсию

Дисперсия: 21.07

Рассмотри модель не линейной регрессии \[ y^*=sin(w_1x+w_2) + w_3x + w_4 \]

Nonlinear regression model
  model: OIL ~ sin(w1 * RUB + w2) + w3 * RUB + w4
   data: rub_vs_oil
     w1      w2      w3      w4 
  1.041   1.529  -1.668 152.088 
 residual sum-of-squares: 14183

Number of iterations to convergence: 19 
Achieved convergence tolerance: 1.908e-06

Дисперсия такой модели 19.78, это лишь немного ниже, чем в случае линейной регрессии.

Задание 3. Дисперсионный анализ

Школьник получил в целом 30 оценок по пяти дисциплинам. Выборку постройте сами. Найдите коэффициент детерминации и определите, является уровень значимости дисциплины в общей вариации оценок.

subjects <- gl(5, 6, labels=c('Дисциплина1','Дисциплина2','Дисциплина3','Дисциплина4','Дисциплина5'))
marks <- c(5,4,5,4,5,4,  4,3,4,4,3,3,  5,5,5,4,5,5,  5,4,3,2,5,3,  4,5,3,5,3,4)
ds <- data.frame(subject=subjects, mark=marks)
ds

Нулевая гипотеза \(H_0\): средние значения по дисциплинам равны

Расчитаем средние по дисциплинам

means <- aggregate(mark ~ subject, data=ds, mean)
means

Расчитаем дисперсии по дисциплинам

variances <- aggregate(mark ~ subject, data=ds, var)
variances

Получим оценку дисперсии поппуляции с двух сторон, как внутригрупповую (MSE) и межгрупповую (MSB)

MSE = mean(variances$mark)
MSE

[1] 0.6066667

GM = mean(ds$mark)
MSB = 6 * (t((means$mark - GM)) %*% (means$mark - GM)) / 4

обшее среднее GM = 4.1

MSB = 1.8833333

DFN = 5 - 1
DFD = 30 - 5
MSB/MSE

         [,1]
[1,] 3.104396

Отношение \(\frac{MSB}{MSE}=\) 3.1043956

P-значение: 0.0333, что чуть ниже уровня значимости 0.05

Расчитаем также, как раскладывается общая сумма квадратов на составляющие \[ SSQ_{общая} = SSQ_{дисциплины} + SSQ_{ошибки} \]

ssq_total = (ds$mark - GM) %*% (ds$mark - GM)
ssq_condition = 6 * (t((means$mark - GM)) %*% (means$mark - GM))
ssq_error = (ds[ds$subject=='Дисциплина1',2] - means[1,2]) %*% (ds[ds$subject=='Дисциплина1',2] - means[1,2]) + 
  (ds[ds$subject=='Дисциплина2',2] - means[2,2]) %*% (ds[ds$subject=='Дисциплина2',2] - means[2,2]) +
  (ds[ds$subject=='Дисциплина3',2] - means[3,2]) %*% (ds[ds$subject=='Дисциплина3',2] - means[3,2]) +
  (ds[ds$subject=='Дисциплина4',2] - means[4,2]) %*% (ds[ds$subject=='Дисциплина4',2] - means[4,2]) +
  (ds[ds$subject=='Дисциплина5',2] - means[5,2]) %*% (ds[ds$subject=='Дисциплина5',2] - means[5,2]) 

\[ 22.7 = 7.5 + 15.2 \]

Интересно сравнить результаты наших вычислений с результатами пакетных. Видно, что они совпадают

da <- aov(mark ~ subject, data=ds)
summary(da)

            Df Sum Sq Mean Sq F value Pr(>F)  
subject      4  7.533  1.8833   3.104 0.0333 *
Residuals   25 15.167  0.6067                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Делаем вывод, что покрайней мере по одной из дисциплин средняя успеваемость отличается от других групп. Можно также отметить, размер групп 6 элементов сильно вляет на точность метода в худшую сторону.