Oszacuj szybkość, z jaką odchylenie standardowe skaluje się z niezależną zmienną

11

Mam eksperyment, w którym wykonuję pomiary normalnie rozłożonej zmiennej Y ,

YN(μ,σ)

Jednak poprzednie eksperymenty dostarczyły pewnych dowodów, że odchylenie standardowe σ jest funkcją afiniczną zmiennej niezależnej X , tj

σ=a|X|+b

YN(μ,a|X|+b)

Ja jak oszacowanie parametrów i B przez próbkowanie Y w wielu wartości X . Dodatkowo, z powodu ograniczeń eksperymentu, mogę pobrać tylko ograniczoną (około 30-40) liczbę próbek Y , i wolałbym pobierać próbki przy kilku wartościach X z niepowiązanych powodów eksperymentalnych. Biorąc pod uwagę te ograniczenia, jakie metody są dostępne do oszacowania a i b ?abYXYXab

Opis eksperymentu

To jest dodatkowa informacja, jeśli jesteś zainteresowany, dlaczego zadaję powyższe pytanie. Mój eksperyment mierzy słuchową i wizualną percepcję przestrzenną. Mogę mieć konfigurację eksperyment, w którym można przedstawić dźwiękowy lub wizualny celów z różnych miejsc, X , i pacjenci wskazuje na postrzeganie lokalizacji docelowej, Y . Zarówno widzenie *, jak i przesłuchanie stają się mniej precyzyjne wraz ze wzrostem mimośrodowości (tj. Wzrostem |X| ), który jako σσ powyżej. Docelowo chciałbym oszacować i Bzabzarówno dla widzenia, jak i przesłuchania, więc znam precyzję każdego zmysłu w różnych miejscach w przestrzeni. Szacunki te zostaną wykorzystane do przewidywania względnej wagi celów wizualnych i słuchowych, gdy zostaną przedstawione jednocześnie (podobnie do teorii integracji multisensorycznej przedstawionej tutaj: http://www.ncbi.nlm.nih.gov/pubmed/12868643 ).

* Wiem, że ten model jest niedokładny dla widzenia podczas porównywania foveal z przestrzenią pozaziemską, ale moje pomiary są ograniczone wyłącznie do przestrzeni pozaziemskiej, gdzie jest to przyzwoite przybliżenie.

Adam Bosen
źródło
2
Ciekawy problem. Prawdopodobnie najlepsze rozwiązania uwzględnią powody, dla których przeprowadzasz ten eksperyment. Jakie są twoje ostateczne cele? Prognoza? Oszacowanie , a i / lub σ ? Im więcej możesz nam powiedzieć o celu, tym lepsze mogą być odpowiedzi. μzaσ
whuber
Ponieważ SD nie może być ujemna, jest mało prawdopodobne, aby była to funkcja liniowa X. Twoja sugestia, a | X |, wymaga węższego lub szerszego kształtu V z minimum przy X = 0, co wydaje mi się raczej nienaturalną możliwością . Jesteś pewien, że to prawda?
gung - Przywróć Monikę
Dobrze, @gung, niewłaściwie uprościłem mój problem. Bardziej realistycznie byłoby powiedzieć, że jest funkcją afiniczną | X | . Zmienię moje pytanie. σ|X|
Adam Bosen
@whuber Powód tego jest trochę zaangażowany, ale zastanowię się, jak wyjaśnić eksperyment i wkrótce dodać więcej szczegółów do mojego pytania.
Adam Bosen
1
Czy masz dobry powód, aby a-priori wierzyć, że X = 0 reprezentuje minimalne SD i że f (| X |) jest monotoniczny?
gung - Przywróć Monikę

Odpowiedzi:

2

W przypadku takim jak twój, w którym masz stosunkowo prosty, ale „niestandardowy” model generatywny, dla którego chcesz oszacować parametry, najpierw pomyślałem o użyciu programu wnioskowania bayesowskiego, takiego jak Stan . Podany opis przełoży się bardzo czysto na model Stana.

Przykładowy kod R, wykorzystujący RStan (interfejs R do 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)

Otrzymasz wynik, który wygląda mniej więcej tak (chociaż twoje liczby losowe prawdopodobnie będą inne niż moje):

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).

Model zszedł dobrze (Rhat = 1), a efektywna wielkość próby (n_eff) jest dość duża we wszystkich przypadkach, więc na poziomie technicznym model jest dobrze zachowany. Najlepsze szacunki , b i ľ (w średniej kolumnie) są również dość blisko tego, co zostało przewidziane.zabμ

Martin O'Leary
źródło
Och, podoba mi się to! Nie słyszałem wcześniej o Stanie, dzięki za referencje. Początkowo liczyłem na rozwiązanie analityczne, ale biorąc pod uwagę brak odpowiedzi, wątpię, by takie istniało. Jestem skłonny wierzyć, że twoja odpowiedź jest najlepszym podejściem do tego problemu.
Adam Bosen
Nie byłoby dla mnie kompletnym szokiem, gdyby istniało rozwiązanie analityczne, ale z pewnością byłbym trochę zaskoczony. Siła używania czegoś takiego jak Stan polega na tym, że bardzo łatwo jest wprowadzać zmiany w modelu - rozwiązanie analityczne prawdopodobnie byłoby bardzo mocno ograniczone do modelu, jak podano.
Martin O'Leary
2

YN.(μ,za|x|+b)
l(μ,za,b)=-ln(za|xja|+b)-12)(yja-μza|xja|+b)2)

W R możemy to zrobić

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

Następnie symuluj niektóre dane:

> 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)

Następnie włącz funkcję loglikelihood:

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

Następnie zoptymalizuj:

> 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
kjetil b halvorsen
źródło