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

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

Y \sim N (μ, σ)

$Y \sim N(\mu,\sigma)$

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

σ = a | X | + b

$\sigma = a|X| + b$

Y \sim N (μ, a | X | + b)

$Y \sim N(\mu,a|X| + b)$

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

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 $\sigma$ powyżej. Docelowo chciałbym oszacować i $a$ $b$ zaró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

Ciekawy problem. Prawdopodobnie najlepsze rozwiązania uwzględnią powody, dla których przeprowadzasz ten eksperyment. Jakie są twoje ostateczne cele? Prognoza? Oszacowanie

i / lub

? Im więcej możesz nam powiedzieć o celu, tym lepsze mogą być odpowiedzi.

μ

$\mu$

a

$a$

σ

$\sigma$

— 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ą

. Zmienię moje pytanie.

σ

$\sigma$

| X |

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

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:

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 , i (w średniej kolumnie) są również dość blisko tego, co zostało przewidziane. $a$ $b$ $\mu$

— 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

Y \sim N. (μ, za | x | + b)

$\newcommand{\dist}{\sim} Y \dist N(\mu, a|x|+b)$

l (μ, za, b) = - \sum \ln (za | x_{ja} | + b) - \frac{1}{2)} \sum {(\frac{y_{ja} - μ}{za | x_{ja} | + b})}^{2)}

$l(\mu, a, b) = -\sum \ln(a|x_i|+b) -\frac12\sum\left(\frac{y_i-\mu}{a|x_i|+b}\right)^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