Прогнозирование дисперсии гетероскедастических данных

Я пытаюсь сделать регрессию на гетероскедастических данных, где я пытаюсь предсказать отклонения ошибки, а также средние значения в терминах линейной модели. Что-то вроде этого:

$$\begin{align}\\ y\left(x,t\right) &= \bar{y}\left(x,t\right)+\xi\left(x,t\right),\\ \xi\left(x,t\right) &\sim N\left(0,\sigma\left(x,t\right)\right),\\ \bar{y}\left(x,t\right) &= y_{0}+ax+bt,\\ \sigma\left(x,t\right) &= \sigma_{0}+cx+dt. \end{align}$$

Словом, данные состоят из повторных измерений при различных значениях и . Я предполагаю, что эти измерения состоят из «истинного» среднего значения которое является линейной функцией и , с аддитивным гауссовым шумом , стандартное отклонение которого (или дисперсия, Я не решил) также линейно зависит от . (Я мог бы допустить более сложные зависимости от и - нет сильной теоретической мотивации для линейной формы - но я бы не стал слишком усложнять вещи на этом этапе.)$y(x,t)$$x$$t$$\bar{y}(x,t)$$x$$t$$\xi(x,t)$$x,t$$x$$t$

Я знаю, что поисковым термином здесь является «гетероскедастичность», но все, что я смог найти, это дискуссии о том, как уменьшить / удалить его, чтобы лучше предсказать , но ничего с точки зрения попыток предсказать в терминах независимых переменных. Я хотел бы оценить и с доверительными интервалами (или байесовскими эквивалентами), и если есть простой способ сделать это в SPSS, тем лучше! Что я должен делать? Благодарю.$\bar{y}$ $\sigma$$y_0, a, b, \sigma_0, c$$d$

Я думаю, что ваша первая проблема заключается в том, что больше не является нормальным распределением, и то, как данные должны быть преобразованы в гомоскедастические, зависит именно от того, что такое σ ( x , t ) . Например, если σ ( x , t ) = a x + b t , то ошибка является пропорциональной, и логарифм y-данных должен быть взят перед регрессией, или регрессия, скорректированная с обычного наименьших квадратов (OLS) до взвешенного наименьших квадратов с $N\left(0,\sigma\left(x,t\right)\right)$$\sigma\left(x,t\right)$$\sigma\left(x,t\right)= ax+bt$ вес (который изменяет регрессию до минимизированной пропорциональной погрешности типа). Точно так же, если σ ( x , t ) = e a x + b t , нужно будет взять логарифм логарифма и регрессировать.$1/y^2$$\sigma\left(x,t\right)= e^{a x+b t}$

Я думаю, что причина, по которой предсказание типов ошибок плохо освещено, заключается в том, что сначала нужно выполнить любую старую регрессию (стон, обычно обычные наименьшие квадраты, OLS). А из остаточного графика, т. Е. , можно наблюдать остаточную форму, а другой - гистограмму частоты данных и смотреть на это. Затем, если остатки представляют собой веерные балки, открывающиеся справа, можно попробовать пропорциональное моделирование данных, если гистограмма выглядит как экспоненциальный спад, можно попробовать возвратно-поступательное движение, 1 / y и т. Д. Для квадратных корней, возведения в квадрат, возведения в степень , принимая экспоненциально-у.$model-y$$1/y$

Теперь это только небольшой рассказ. Более длинная версия включает в себя гораздо больше типов регрессии, в том числе регрессию медианы Тейла, двумерную регрессию Деминга и регрессию для минимизации ошибки некорректных задач, которые не имеют особой зависимости соответствия качества от минимизируемой распространяемой ошибки. Это последнее является громадина, но, увидеть этоВ качестве примера. Так что это имеет большое значение, какие ответы каждый пытается получить. Как правило, если кто-то хочет установить связь между переменными, рутинная OLS не является предпочтительным методом, и регрессия Тейла может быть быстрым и грязным улучшением этого. OLS минимизируется только в направлении y, поэтому наклон слишком мал, а пересечение слишком велико, чтобы установить, каковы основные отношения между переменными. Иными словами, OLS дает наименьшую оценку ошибки для y при заданном x, но не дает оценку того, как x изменяется с y. Когда значения r очень высоки (0.99999+), не имеет большого значения, какая регрессия используется, и OLS в y приблизительно такой же, как OLS в x, но, когда значения r низкие, OLS в y очень отличается от МЖС в х.

Таким образом, многое зависит от того, какие именно мотивы в первую очередь мотивировали проведение регрессионного анализа. Это диктует необходимые численные методы. После того, как этот выбор сделан, остатки тогда имеют структуру, которая связана с целью регрессии, и должны быть проанализированы в этом более широком контексте.

Команда расширения STATS BREUSCH PAGAN может одновременно проверять невязки на гетероскедастичность и оценивать ее как функцию некоторых или всех регрессоров.

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

$$LL(y_0, a, b, \sigma_0, c, d) = \sum_{i=1}^n \log \phi(y_i, y_0 + a x_i + b t_i, \sigma_0 + c x_i + d t_i)$$ $$\phi(x, \mu, \sigma) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$$

$\hat{\theta}$$\theta=(y_0, a, b, \sigma_0, c, d)$

$H$$\theta$$n$$\hat{\theta}$$H^{-1}$

Вот пример кода в Python:

import scipy
import numpy as np

# generate toy data for the problem
np.random.seed(1) # fix random seed
n = 1000 # fix problem size
x = np.random.normal(size=n)
t = np.random.normal(size=n)
mean = 1 + x * 2 + t * 3
std = 4 + x * 0.5 + t * 0.6
y = np.random.normal(size=n, loc=mean, scale=std)

# create negative log likelihood
def neg_log_lik(theta):
    est_mean = theta[0] + x * theta[1] + t * theta[2]
    est_std = np.maximum(theta[3] + x * theta[4] + t * theta[5], 1e-10)
    return -sum(scipy.stats.norm.logpdf(y, loc=est_mean, scale=est_std))

# maximize
initial = np.array([0,0,0,1,0,0])
result = scipy.optimize.minimize(neg_log_lik, initial)
# extract point estimation
param = result.x
print(param)
# extract standard error for confidence intervals
std_error = np.sqrt(np.diag(result.hess_inv))
print(std_error)

$\sigma$$\sigma$$10^{-10}$

Результат (оценки параметров и их стандартные ошибки), произведенные кодом:

[ 0.8724218   1.75510897  2.87661843  3.88917283  0.63696726  0.5788625 ]
[ 0.15073344  0.07351353  0.09515104  0.08086239  0.08422978  0.0853192 ]

Вы можете видеть, что оценки близки к своим истинным значениям, что подтверждает правильность этого моделирования.