Prognozowanie wariancji danych heteroscedastycznych

Próbuję wykonać regresję danych heteroscedastycznych, w których próbuję przewidzieć wariancje błędów, a także wartości średnie w odniesieniu do modelu liniowego. Coś takiego:

$$\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}$$

Słowami, dane składa się z powtarzalnych pomiarów przy różnych wartościach i . Sądzę pomiary te składają się z „prawdziwego” Średnia wartość , który jest liniową funkcją i , w dodatku do szumu gaussowskiego , którego odchylenie standardowe (lub wariancji Nie zdecydowałem) zależy również liniowo od . (Mógłbym pozwolić na bardziej skomplikowane zależności od i - nie ma silnej teoretycznej motywacji dla formy liniowej - ale wolałbym nie komplikować rzeczy na tym etapie).$y(x,t)$$x$$t$$\bar{y}(x,t)$$x$$t$$\xi(x,t)$$x,t$$x$$t$

Wiem, że wyszukiwanym terminem jest „heteroscedastyczność”, ale jak dotąd udało mi się znaleźć dyskusje o tym, jak go zmniejszyć / usunąć, aby lepiej przewidzieć , ale nic w kategoriach próby przewidzenia$\bar{y}$ $\sigma$ pod względem zmiennych niezależnych. Chciałbym szacują, i D z przedziałami ufności (lub ekwiwalentów Bayesa), a jeśli nie jest to łatwy sposób, aby to zrobić w SPSS tym lepiej! Co powinienem zrobić? Dzięki.$y_0, a, b, \sigma_0, c$$d$

Myślę, że twoim pierwszym problemem jest to, że nie jest już rozkładem normalnym, a to, jak dane muszą zostać przekształcone, aby były homoscedastyczne, zależy dokładnie od tego, czym jest σ ( x , t ) . Na przykład, jeśli σ ( x , t ) = a x + b t , to błąd ma charakter proporcjonalny, a logarytm danych y powinien zostać przyjęty przed regresją lub regresja skorygowana ze zwykłych najmniejszych kwadratów (OLS) do ważonej najmniej kwadratów z $N\left(0,\sigma\left(x,t\right)\right)$$\sigma\left(x,t\right)$$\sigma\left(x,t\right)= ax+bt$ waga (która zmienia regresję na zminimalizowany błąd typu proporcjonalnego). Podobnie, jeśli σ ( x , t ) = e a x + b t , należałoby wziąć logarytm logarytmu i go regresować.$1/y^2$$\sigma\left(x,t\right)= e^{a x+b t}$

Myślę, że powodem słabego przewidywania typów błędów jest to, że najpierw wykonuje się jakąkolwiek starą regresję (jęk, zwykle zwykłe najmniejsze kwadraty, OLS). Oraz z linii dostarczającej pozostałości działki, tj , obserwuje się resztkowe kształt i przedstawi się histogram częstotliwości danych i wygląda na to. Następnie, jeśli reszty są wiązką wachlarza otwierającą się w prawo, próbuje się proporcjonalnego modelowania danych, jeśli histogram wygląda jak rozkład wykładniczy, można spróbować odwrotności, 1 / y itd. Itd. Dla pierwiastków kwadratowych, kwadratów, potęgowania , biorąc wykładniczy-y.$model-y$$1/y$

To tylko krótka historia. Dłuższa wersja zawiera o wiele więcej rodzajów regresji, w tym regresję medianową Theila, regresję dwuwymiarową Deminga i regresję w celu minimalizacji błędu źle przedstawionych problemów, które nie mają szczególnego związku dopasowania dopasowania krzywej do pomniejszonego propagowanego błędu. Ten ostatni jest niesamowity, ale zobacz tojako przykład. Tak, że robi to dużą różnicę, jakie odpowiedzi próbuje się uzyskać. Zazwyczaj, jeśli ktoś chce ustalić związek między zmiennymi, rutynowy OLS nie jest metodą z wyboru, a regresja Theila byłaby szybką i nieprzyzwoitą poprawą. OLS minimalizuje się tylko w kierunku y, więc nachylenie jest zbyt płytkie, a przecięcie zbyt duże, aby ustalić, jaka jest podstawowa zależność między zmiennymi. Innymi słowy, OLS podaje oszacowanie najmniejszego błędu ay przy x, nie podaje oszacowania, jak x zmienia się zy. Gdy wartości r są bardzo wysokie (0,99999+), nie ma znaczenia, jaką regresję stosuje się, a OLS w y jest w przybliżeniu taki sam jak OLS w x, ale gdy wartości r są niskie, OLS w y bardzo różni się od OLS w x.

Podsumowując, wiele zależy dokładnie od tego, jakie jest uzasadnienie przeprowadzania analizy regresji w pierwszej kolejności. To dyktuje potrzebne metody numeryczne. Po dokonaniu tego wyboru reszty mają następnie strukturę związaną z celem regresji i muszą być analizowane w tym szerszym kontekście.

Polecenie rozszerzenia STATUS BREUSCH PAGAN może zarówno testować resztki pod kątem heteroscedastyczności, jak i szacować je jako funkcję niektórych lub wszystkich regresorów.

Ogólne podejście do tego rodzaju problemów polega na maksymalizacji (uregulowanego) prawdopodobieństwa danych.

$$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}$

Oto przykładowy kod w Pythonie:

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}$

Wynik (oszacowania parametrów i ich błędy standardowe) wygenerowany przez kod to:

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

Widać, że szacunki są zbliżone do ich prawdziwych wartości, co potwierdza poprawność tej symulacji.