Perkirakan Tingkat Deviasi Timbangan Standar Dengan Variabel Independen

Saya memiliki percobaan di mana saya melakukan pengukuran variabel terdistribusi normal $Y$,

Namun, percobaan sebelumnya telah memberikan beberapa bukti bahwa standar deviasi $\sigma$ adalah fungsi affine dari variabel independen $X$ , yaitu

Saya ingin memperkirakan parameter $a$ dan $b$ dengan sampling $Y$ di beberapa nilai dari $X$ . Selain itu, karena keterbatasan percobaan, saya hanya dapat mengambil jumlah sampel terbatas (sekitar 30-40) $Y$, dan lebih suka sampel pada beberapa nilai $X$ untuk alasan eksperimental yang tidak terkait. Dengan adanya kendala-kendala ini, metode apa yang tersedia untuk memperkirakan $a$ dan $b$ ?

Deskripsi Eksperimen

Ini adalah informasi tambahan, jika Anda tertarik mengapa saya mengajukan pertanyaan di atas. Eksperimen saya mengukur persepsi spasial pendengaran dan visual. Saya memiliki setup percobaan di mana saya bisa hadir baik auditori atau target visual dari lokasi yang berbeda, , dan mata pelajaran menunjukkan lokasi yang dirasakan dari target, Y . Baik visi * dan audisi menjadi kurang tepat dengan meningkatnya eksentrisitas (yaitu peningkatan | X | ), yang saya modelkan sebagai σ di atas. Pada akhirnya, saya ingin memperkirakan a dan $X$$Y$$|X|$$\sigma$$a$$b$untuk kedua visi dan audisi, jadi saya tahu ketepatan setiap indera di berbagai lokasi di ruang angkasa. Perkiraan ini akan digunakan untuk memprediksi bobot relatif dari target visual dan pendengaran ketika disajikan secara bersamaan (mirip dengan teori integrasi multisensor yang disajikan di sini: http://www.ncbi.nlm.nih.gov/pubmed/12868643 ).

* Saya tahu bahwa model ini tidak akurat untuk penglihatan ketika membandingkan ruang foveal dengan ruang ekstrafoveal, tetapi pengukuran saya dibatasi hanya untuk ruang ekstrafoveal, di mana ini merupakan perkiraan yang layak.

Dalam kasus seperti milik Anda, di mana Anda memiliki model generatif yang relatif sederhana, tetapi "tidak standar" yang ingin Anda perkirakan parameternya, pemikiran pertama saya adalah menggunakan program inferensi Bayesian seperti Stan . Deskripsi yang Anda berikan akan menerjemahkan dengan sangat bersih ke model Stan.

Beberapa contoh kode R, menggunakan RStan (antarmuka R ke Stan).

library(rstan)

model_code <- "
data {
    int<lower=0> n; // number of observations
    real y[n];
    real x[n];
}
parameters {
    real mu; // I've assumed mu is to be fit.
             // Move this to the data section if you know the value of mu.
    real<lower=0> a;
    real<lower=0> b;
}
transformed parameters {
    real sigma[n];
    for (i in 1:n) {
        sigma[i] <- a + b * fabs(x[i]);
    }
}
model {
    y ~ normal(mu, sigma);
}
"

# Let's generate some test data with known parameters

mu <- 0
a <- 2
b <- 1

n <- 30
x <- runif(n, -3, 3)
sigma <- a + b * abs(x)
y <- rnorm(n, mu, sigma)

# And now let's fit our model to those "observations"

fit <- stan(model_code=model_code,
            data=list(n=n, x=x, y=y))

print(fit, pars=c("a", "b", "mu"), digits=1)

Anda akan mendapatkan hasil yang terlihat seperti ini (walaupun angka acak Anda mungkin akan berbeda dengan milik saya):

Inference for Stan model: model_code.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

   mean se_mean  sd 2.5%  25% 50% 75% 97.5% n_eff Rhat
a   2.3       0 0.7  1.2  1.8 2.2 2.8   3.9  1091    1
b   0.9       0 0.5  0.1  0.6 0.9 1.2   1.9  1194    1
mu  0.1       0 0.6 -1.1 -0.3 0.1 0.5   1.4  1262    1

Samples were drawn using NUTS(diag_e) at Thu Jan 22 14:26:16 2015.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

$a$$b$$\mu$

Anda tidak bisa mengharapkan formula tertutup, tetapi Anda masih bisa menuliskan fungsi kemungkinan dan memaksimalkannya secara numerik. Model Anda adalah Kemudian fungsi kemungkinan log (terlepas dari istilah yang tidak tergantung pada parameter) menjadi dan itu mudah diprogram dan berikan kepada pengoptimal angka.

Di R, bisa kita lakukan

make_lik  <-  function(x,y){
    x  <-  abs(x)
    function(par) {
        mu <- par[1];a  <-  par[2];  b <-  par[3]
        axpb <-  a*x+b
        -sum(log(axpb)) -0.5*sum( ((y-mu)/axpb)^2 )
    }
}

Kemudian simulasikan beberapa data:

> x <-  rep(c(2,4,6,8),10)
> x
 [1] 2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8 2 4 6 8 2 4
[39] 6 8
> a <- 1
> b<-  3
> sigma <-  a*x+b
> mu  <-  10
> y  <-  rnorm(40,mu, sd=sigma)

Kemudian buat fungsi kemungkinan loglikel:

> lik <-  make_lik(x,y)
> lik(c(10,1,3))
[1] -99.53438

Kemudian optimalkan:

> optim(c(9.5,1.2,3.1),fn=function(par)-lik(par))
$par
[1] 9.275943 1.043019 2.392660

$value
[1] 99.12962

$counts
function gradient 
     136       NA 

$convergence
[1] 0

$message
NULL