Каков возраст Вселенной?

Как было показано в предыдущем сообщении, статистическая модель представляет собой упрощенное математическое представление процесса, который, как мы полагаем, привел к генерации наблюдаемых значений изучаемой переменной. Статистические модели используются для решения двух типов задач: 1) выявление характера взаимоотношений между переменной-откликом и некоторыми внешними факторами (предикторами) и 2) предсказание значений переменной-отклика на основе известных значений предикторов. В данном сообщении будет представлен пример использования простой линейной регрессии для решения задачи первого типа. Этой задачей будет установление... возраста Вселенной. Но сначала, как обычно, немного теории...

Простая линейная регрессия

Линейная регрессия является одним из наиболее широко используемых методов статистики и сводится к нахождению зависимости следующего вида:

\[y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \dots + \beta_p x_{ip} + \epsilon_i ,\]

где \(y_i\) - это i-тое значение зависимой переменной (также называемой "переменной-окликом"), \(x_{i1}, x_{i2}, \dots x_{ip}\) - соответствующие значения независимых переменных \(x_{1}, x_{2}, \dots x_{p}\) (также называемых "предикторами"), а \(\epsilon_i\) - остатки, т.е. необъясненные моделью расхождения между реально наблюдаемыми и предсказанными значениями зависимой переменной:

\[ \epsilon_i = y_i - (\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \dots + \beta_p x_{ip}) \] 

Предполагается, что остатки модели имеют нормальное распределение со средним значением 0 и некоторым стандартным отклонением \(\sigma\). Предикторы могут быть представлены как количественными, так и качественными переменными. В случае, когда все предикторы являются качественными переменными, речь идет о дисперсионном анализе (подробнее см. здесь, здесь и здесь). При наличии в модели как качественных, так и количественных переменных, говорят о ковариационном анализе (англ. "analysis of covariance", ANCOVA). Модель, в которую входит только один количественный предиктор, называют простой линейной регрессией. 

Оценка параметров линейной регрессии (\(\beta_0\), \(\beta_1\), и т.д.) выполняется по методу наименьших квадратов (МНК) (англ. least squares), который представляет собой частный случай более общего метода максимального правдоподобия (англ. maximum likelihood). Предположим, что мы имеем дело с простой линейной регрессией:

\[y_i = \beta_0 + \beta_1 x_i + \epsilon_i \]

МНК заключается в нахождении таких значений \(\beta_0\) и \(\beta_1\), которые минимизируют сумму квадратов остатков (англ. residual sum of squares) модели, т.е.

\[ RSS = \sum (y_i - (\beta_0 + \beta_1 x_i))^2 \]

Можно показать, что минимальное значение \(RSS\) будет достигаться при следующих значениях коэффициентов модели:

\[ \hat{\beta_1} = \frac{\sum (x_i - \bar{x})(y_i - \bar{y})}{\sum (x_i - \bar{x})^2}, \]

\[ \hat{\beta_0} = \bar{y} - \hat{\beta_1}\bar{x}, \] 

где \( \hat{\beta_1} \) и \( \hat{\beta_0}\) - выборочные оценки соответствующих параметров модели, а \( \bar{y} \) и \( \bar{x}\) - выборочные средние значения переменных \(y\) и \(x\).

Наконец, стандартное отклонение остатков (англ. residual standard error) оценивается как

\[ \sigma = \sqrt{RSS/(n - 2)} \]

Данные по скорости движения галактик

Freedman et al. (2001) опубликовали данные по расстоянию до 24 галактик, а также по скорости удаления этих галактик, полученные при помощи космического телескопа "Хаббл". Данные были собраны в рамках проекта (т.н. Key Project  - "ключевой проект"), целью которого являлось уточнение значения постоянной Хаббла (см. таблицы 4 и 5 в оригинальной статье, а также рисунок ниже). Эта постоянная представляет собой коэффициент в уравнении закона Хаббла, который описывает связь между расстоянием до внегалактического объекта (например, галактики, квазара) и скоростью его удаления, обусловленного расширением Вселенной после Большого взрыва. Этот закон выражается простой линейной зависимостью, которая может быть записана следующим образом:

\[y = \beta x, \]

где \(y\) - относительная скорость движения любых двух галактик, разграниченных в данный момент времени расстоянием \(y\). Постоянная Хаббла, обозначенная здесь как \(\beta\), выражается в км/с на мегапарсек. Величина, обратная постоянной Хаббла, дает приблизительный возраст Вселенной (т.н. "Хаббловский возраст").

Используя данные Freedman et al. (2011), мы можем оценить значение постоянной Хаббла при помощи простой линейной регрессии (обратите внимание: свободный член уравнения здесь приравнен нулю, поскольку в момент, когда Вселенная находилась в состоянии сингулярности, галактик не существовало и они, конечно же, не могли удаляться друг от друга):

\[y_i = \beta x_i + \epsilon_i\]

Оценивание параметров линейной регрессии в R

Готовые для работы данные из статьи Freedman et al. (2001) можно найти в составе таблицы hubble из пакета gamair:

install.packages("gamair")
library(gamair)
data(hubble)
 
str(hubble)
'data.frame': 24 obs. of  3 variables:
 $ Galaxy: Factor w/ 24 levels "IC4182","NGC0300",..: 2 3 4 5 6 8 9 7 10 11 ...
 $ y     : int  133 664 1794 1594 1473 278 714 882 80 772 ...
 $ x     : num  2 9.16 16.14 17.95 21.88 ...

Как было показано в нескольких предыдущих сообщениях (см., например, здесь и здесь), в R параметры линейной модели можно оценить при помощи стандартной функции lm():

M <- lm(y ~ x - 1, data = hubble)

Обратите внимание на -1 в формуле модели - такое обозначение используется для того, чтобы исключить свободный член регрессионной модели (см. объяснение выше).

Посмотрим на результат подгонки модели:

summary(M)
 
Call:
lm(formula = y ~ x - 1, data = hubble)
 
Residuals:
   Min     1Q Median     3Q    Max 
-736.5 -132.5  -19.0  172.2  558.0 
 
Coefficients:
  Estimate Std. Error t value Pr(>|t|)    
x   76.581      3.965   19.32 1.03e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 258.9 on 23 degrees of freedom
Multiple R-squared:  0.9419, Adjusted R-squared:  0.9394 
F-statistic: 373.1 on 1 and 23 DF,  p-value: 1.032e-15

Как видим, оцененное значение постоянной Хаббла составило 76.581 км/с на мегапарсек. Это значение существенно отличается от нуля (см. Р-значение соответствующего t-теста в столбце Pr(>|t|)). На рисунке ниже приведена линия регрессии, описываемая полученным нами уравнением \(y = 76.581x\):

Расчет возраста Вселенной теперь не представляет труда. Один мегапарсек - это \(3.09\times10^{19}\) км. Разделим полученную выше постоянную Хаббла на это значение, чтобы выразить ее в секундах:

hub.const <- 76.581/3.09e19
 
hub.const
[1] 2.47835e-18

Тогда возраст Вселенной, выраженный в секундах, составит:

age <- 1/hub.const
 
age
[1] 4.034943e+17

Выполнив простое преобразование, получим возраст, выраженный в годах:

age/(60^2*24*365)
[1] 12794721567

Другими словами, данные Freedman et al. (2001) указывают на то, что возраст Вселенной составляет ~12.8 миллиардов лет. Следующее сообщение будет посвящено обсуждению того, насколько велика неопределенность в отношении этой величины. Пока же лишь отметим, что по последним оценкам возраст Вселенной составляет 13.798 млрд. лет.