Memprediksi varian data heteroskedastik

Saya mencoba melakukan regresi pada data heteroskedastik di mana saya mencoba untuk memprediksi varians kesalahan serta nilai rata-rata dalam hal model linier. Sesuatu seperti ini:

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

Dengan kata lain, data terdiri dari pengukuran berulang pada berbagai nilai dan . Saya menganggap pengukuran ini terdiri dari nilai rata-rata "benar" yang merupakan fungsi linier dari dan , dengan aditif Gaussian noise yang standar deviasinya (atau varians, Saya belum memutuskan) juga tergantung secara linear pada . (Saya bisa membiarkan dependensi yang lebih rumit pada dan - tidak ada motivasi teoritis yang kuat untuk bentuk linear - tapi saya lebih suka tidak terlalu rumit pada tahap ini.)$y(x,t)$$x$$t$$\bar{y}(x,t)$$x$$t$$\xi(x,t)$$x,t$$x$$t$

Saya tahu istilah pencarian di sini adalah "heteroskedastisitas," tetapi semua yang saya dapat temukan sejauh ini adalah diskusi tentang cara mengurangi / menghapusnya untuk memprediksi , tetapi tidak ada yang mencoba memprediksi dalam hal variabel independen. Saya ingin memperkirakan dan d dengan interval kepercayaan (atau setara Bayesian), dan jika ada cara mudah untuk melakukannya di SPSS, semakin baik! Apa yang harus saya lakukan? Terima kasih.$\bar{y}$ $\sigma$$y_0, a, b, \sigma_0, c$$d$

Saya pikir masalah pertama Anda adalah bahwa bukan lagi distribusi normal, dan bagaimana data perlu ditransformasikan menjadi homoscedastic tergantung pada apa σ ( x , t ) itu. Misalnya, jika σ ( x , t ) = a x + b t , maka kesalahannya adalah tipe proporsional dan logaritma data y harus diambil sebelum regresi, atau, regresi disesuaikan dari kuadrat terkecil kuadrat (OLS) menjadi berbobot. kuadrat terkecil dengan $N\left(0,\sigma\left(x,t\right)\right)$$\sigma\left(x,t\right)$$\sigma\left(x,t\right)= ax+bt$ berat (yang mengubah regresi untuk meminimalkan kesalahan jenis proporsional). Demikian pula, jika σ ( x , t ) = e a x + b t , salah satu harus mengambil logaritma logaritma dan kemunduran itu.$1/y^2$$\sigma\left(x,t\right)= e^{a x+b t}$

Saya pikir alasan mengapa prediksi jenis kesalahan kurang tercakup adalah bahwa orang pertama melakukan regresi lama (mengeluh, biasanya kuadrat terkecil, OLS). Dan dari plot residual, yaitu, , seseorang mengamati bentuk residu, dan satu plot histogram frekuensi data, dan melihatnya. Kemudian, jika residu adalah bukaan sinar kipas ke kanan, seseorang mencoba pemodelan data proporsional, jika histogramnya tampak seperti peluruhan eksponensial, ia mungkin mencoba membalas, 1 / y , dan seterusnya untuk akar kuadrat, kuadrat, eksponensial , mengambil eksponensial-y.$model-y$$1/y$

Nah, itu hanya cerita pendek. Versi yang lebih panjang mencakup lebih banyak jenis regresi termasuk regresi median Theil, regresi bivariat Deming, dan regresi untuk meminimalkan kesalahan kesalahan-masalah yang tidak memiliki hubungan baik-of-kurva-fit dengan kesalahan yang diperbanyak yang diminimalkan. Yang terakhir adalah tipuan, tapi lihat inisebagai contoh. Sehingga itu membuat perbedaan besar apa jawaban yang seseorang coba dapatkan Biasanya, jika seseorang ingin membangun hubungan antara variabel, OLS rutin bukanlah metode pilihan, dan regresi Theil akan menjadi perbaikan cepat dan kotor pada itu. OLS hanya meminimalkan dalam arah y, sehingga kemiringannya terlalu dangkal, dan intersepnya terlalu besar untuk menentukan apa hubungan mendasar antara variabel-variabel tersebut. Untuk mengatakan ini dengan cara lain, OLS memberikan estimasi kesalahan paling kecil jika diberikan x, itu tidak memberikan perkiraan bagaimana x berubah dengan y. Ketika nilai-r sangat tinggi (0.99999+) membuat sedikit perbedaan apa yang digunakan regresi dan OLS dalam y kira-kira sama dengan OLS di x, tetapi, ketika nilai-r rendah, OLS dalam y sangat berbeda dari OLS dalam x.

Singkatnya, banyak tergantung pada apa alasan yang memotivasi melakukan analisis regresi di tempat pertama. Itu menentukan metode numerik yang diperlukan. Setelah pilihan itu dibuat, residu kemudian memiliki struktur yang terkait dengan tujuan regresi, dan perlu dianalisis dalam konteks yang lebih besar.

Perintah ekstensi STATS BREUSCH PAGAN dapat menguji residu untuk heteroskedastisitas dan memperkirakannya sebagai fungsi dari beberapa atau semua regressor.

Pendekatan umum untuk masalah seperti ini adalah untuk memaksimalkan kemungkinan (diatur) dari data Anda.

Dalam kasus Anda, log-kemungkinan akan terlihat seperti mana ϕ ( x

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

Anda dapat mengkode ekspresi ini menjadi fungsi dalam paket statistik favorit Anda (saya lebih suka Python, R atau Stata, karena saya tidak pernah melakukan pemrograman dalam SPSS). Kemudian Anda dapat memberi makan ke sebuah optimizer numerik, yang akan memperkirakan nilai optimal θ parameter Anda θ = ( y 0 , a , b , σ 0 , c , d ) .$\hat{\theta}$$\theta=(y_0, a, b, \sigma_0, c, d)$

Jika Anda memerlukan interval kepercayaan, pengoptimal ini juga dapat memperkirakan matriks Hessian dari θ (turunan kedua) di sekitar yang optimal. Teori estimasi kemungkinan maksimum mengatakan bahwa untuk besar n matriks kovarians dari θ dapat diperkirakan sebagai H - 1 .$H$$\theta$$n$$\hat{\theta}$$H^{-1}$

Berikut ini contoh kode dengan 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)

Perhatikan bahwa rumusan masalah Anda dapat menghasilkan negatif , dan saya harus mempertahankan diri dari itu dengan penggantian brute force terlalu kecil σ dengan 10 - 10 .$\sigma$$\sigma$$10^{-10}$

Hasil (estimasi parameter dan kesalahan standarnya) yang dihasilkan oleh kode adalah:

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

You can see that estimates are close to their true values, which confirms correctness of this simulation.