Estime la tasa a la que la desviación estándar escala con una variable independiente

Tengo un experimento en el que estoy tomando medidas de una variable normalmente distribuida ,$Y$

Sin embargo, experimentos anteriores han proporcionado alguna evidencia de que la desviación estándar es una función afín de una variable independiente , es decir$\sigma$$X$

Me gustaría para estimar los parámetros de y mediante el muestreo de en múltiples valores de . Además, debido a las limitaciones del experimento, solo puedo tomar un número limitado (aproximadamente 30-40) de muestras de , y preferiría tomar muestras a varios valores de por razones experimentales no relacionadas. Teniendo en cuenta estas limitaciones, qué métodos están disponibles para estimar y ?$a$$b$$Y$$X$$Y$$X$$a$$b$

Descripción del experimento

Esta es información adicional, si está interesado en por qué estoy haciendo la pregunta anterior. Mi experimento mide la percepción espacial auditiva y visual. Tengo una configuración del experimento en el que puede presentar auditiva o objetivos visuales de diferentes sitios, , y los sujetos indican la localización percibida de la diana, . Tanto la visión * como la audición se vuelven menos precisas al aumentar la excentricidad (es decir, aumentar ), que modelé como arriba. En última instancia, me gustaría estimar y$X$$Y$$|X|$$\sigma$$a$$b$tanto para la visión como para la audición, por lo que sé la precisión de cada sentido en una variedad de ubicaciones en el espacio. Estas estimaciones se utilizarán para predecir la ponderación relativa de los objetivos visuales y auditivos cuando se presenten simultáneamente (similar a la teoría de integración multisensorial presentada aquí: http://www.ncbi.nlm.nih.gov/pubmed/12868643 ).

* Sé que este modelo es inexacto para la visión al comparar el espacio foveal con el extrafoveal, pero mis mediciones se limitan únicamente al espacio extrafoveal, donde esta es una aproximación decente.

En un caso como el suyo, donde tiene un modelo generativo relativamente simple, pero "no estándar" para el que desea estimar los parámetros, mi primer pensamiento sería usar un programa de inferencia bayesiano como Stan . La descripción que ha dado se traduciría muy limpiamente a un modelo Stan.

Algunos ejemplos de código R, utilizando RStan (la interfaz R para 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)

Obtendrá una salida similar a esta (aunque sus números aleatorios probablemente serán diferentes a los míos):

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

El modelo ha convergido bien (Rhat = 1), y el tamaño efectivo de la muestra (n_eff) es razonablemente grande en todos los casos, por lo que a nivel técnico el modelo se comporta bien. Las mejores estimaciones de , y (en la columna media) son también bastante cerca de lo que estaba previsto.$a$$b$$\mu$

No puede esperar fórmulas cerradas, pero aún puede escribir la función de probabilidad y maximizarla numéricamente. Su modelo es Entonces la función loglikelihood (aparte de un término que no depende de los parámetros) se convierte en y eso es fácil de programar y dar a un optimizador numérico.

En R, podemos hacer

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

Luego simule algunos datos:

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

Luego haga que la función loglikelihood funcione:

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

Luego optimízalo:

> 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