Schätzen Sie die Rate, mit der die Standardabweichung mit einer unabhängigen Variablen skaliert

Ich habe ein Experiment, in dem ich Messungen einer normalverteilten Variablen .$Y$

Frühere Experimente haben jedoch einige Beweise dafür geliefert, dass die Standardabweichung eine affine Funktion einer unabhängigen Variablen X ist , d. H.$\sigma$$X$

Ich möchte die Parameter und b durch Abtasten von Y bei mehreren Werten von X schätzen . Außerdem kann ich aufgrund von experimentellen Einschränkungen nur eine begrenzte (ungefähr 30-40) Anzahl von Proben von Y entnehmen und würde es aus nicht verwandten experimentellen Gründen vorziehen, bei mehreren Werten von X zu probieren . Welche Methoden stehen angesichts dieser Einschränkungen zur Schätzung von a und b zur Verfügung ?$a$$b$$Y$$X$$Y$$X$$a$$b$

Versuchsbeschreibung

Dies sind zusätzliche Informationen, wenn Sie daran interessiert sind, warum ich die obige Frage stelle. Mein Experiment misst die auditive und visuelle räumliche Wahrnehmung. Ich habe einen Versuchsaufbau, in dem ich entweder akustische oder visuelle Ziele von verschiedenen Orten , und die Probanden geben den wahrgenommenen Ort des Ziels Y an . Sowohl das Sehen * als auch das Hören werden mit zunehmender Exzentrizität (dh zunehmendem | X | ) ungenauer , was ich oben als σ modelliere . Letztendlich möchte ich a und b schätzen$X$$Y$$|X|$$\sigma$$a$$b$Ich kenne also die Präzision jedes Sinnes an verschiedenen Orten im Raum. Diese Schätzungen werden verwendet, um die relative Gewichtung von visuellen und auditorischen Zielen bei gleichzeitiger Darstellung vorherzusagen (ähnlich der hier vorgestellten Theorie der multisensorischen Integration: http://www.ncbi.nlm.nih.gov/pubmed/12868643 ).

* Ich weiß, dass dieses Modell beim Vergleich des fovealen mit dem extrafovealen Raum für das Sehen ungenau ist, aber meine Messungen beschränken sich ausschließlich auf den extrafovealen Raum, wo dies eine anständige Annäherung ist.

In einem Fall wie Ihrem, in dem Sie ein relativ einfaches, aber "nicht standardmäßiges" generatives Modell haben, für das Sie Parameter schätzen möchten, wäre mein erster Gedanke, ein Bayes'sches Inferenzprogramm wie Stan zu verwenden . Die Beschreibung, die Sie gegeben haben, würde sich sehr sauber auf ein Stan-Modell übertragen lassen.

Einige Beispiele für R-Code mit RStan (der R-Schnittstelle zu 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)

Sie erhalten eine Ausgabe, die ungefähr so ​​aussieht (obwohl Ihre Zufallszahlen wahrscheinlich von meinen abweichen werden):

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

Das Modell ist gut konvergiert (Rhat = 1), und die effektive Stichprobengröße (n_eff) ist in allen Fällen relativ groß, sodass sich das Modell auf technischer Ebene gut verhält. Die besten Schätzungen von , b und μ (in der mittleren Spalte) liegen ebenfalls ziemlich nahe an dem, was bereitgestellt wurde.$a$$b$$\mu$

In R können wir tun

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

Dann simulieren Sie einige Daten:

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

Dann machen Sie die Loglikelihood-Funktion:

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

Dann optimieren Sie es:

> 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