La Estimación De La Velocidad A La Que La Desviación Estándar De Las Escalas Con Una Variable Independiente

Question

Tengo un experimento en el que estoy tomando las mediciones de una distribución normal de la variable $Y$,

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

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

$$\sigma = a|X| + b$$

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

Me gustaría estimar los parámetros $a$ $b$ por muestreo $Y$ a varios valores de $X$. Además, debido a experimentar limitaciones que sólo puede tomar un número limitado (aproximadamente 30-40) número de muestras de $Y$, y prefiere muestra varios valores de $X$ para sufragar experimental razones. Dadas estas limitaciones, los métodos disponibles para estimar el $a$$b$?

Experimento Descripción

Esta es la información adicional, si usted está interesado en por qué le estoy preguntando a la pregunta de arriba. Mi experimento medidas auditivo y visual de la percepción espacial. Tengo un experimento de instalación en la que me puedo presentar ya sea visual o auditiva objetivos de diferentes lugares, $X$, y los sujetos indican que la percepción de la localización del objetivo, $Y$. Visión* audición y obtener menos precisa con el aumento de la excentricidad (es decir el aumento de $|X|$), modelo $\sigma$ por encima. En última instancia, me gustaría estimación $a$ $b$ tanto de la visión y la audición, así que yo sé que la precisión del sentido a través de una serie de ubicaciones en el espacio. Estas estimaciones serán utilizados para predecir la ponderación relativa de visual y auditiva de los objetivos cuando se presentan simultáneamente (similar a la teoría de la integración multisensorial que aquí se presenta: http://www.ncbi.nlm.nih.gov/pubmed/12868643).

*Sé que este modelo es inexacta de la visión cuando se comparan foveal a extrafoveal espacio, pero mis medidas son limitados únicamente a extrafoveal espacio, donde se trata de una aproximación decente.

Answer 1

Usted no puede esperar cerrado fórmulas, pero todavía se puede escribir la función de probabilidad y maximizar numéricamente. Su modelo es $$ \newcommand{\dist}{\sim} Y \dist N(\mu,|x|+b) $$ A continuación, el loglikelihoodfunction (aparte de un plazo no depende de los parámetros) se convierte en $$ l(\mu, a, b) = -\sum \ln(a|x_i|+b) -\frac12\sum\left(\frac{y_i-\mu}{un|x_i|+b}\right)^2 $$ y que es fácil de programar y dar numéricos optimizador.

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

A continuación, hacer la loglikelihood función:

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

A continuación, optimizar 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

Answer 2

En un caso como el tuyo, donde se tiene una relativamente simple, pero "no estándar" modelo generativo que te gustaría para estimar los parámetros para, mi primer pensamiento sería el uso de una inferencia Bayesiana programa como Stan. La descripción que has dado traduciría muy limpiamente a un modelo Stan.

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

Usted obtendrá una salida que se ve algo como esto (aunque sus números aleatorios probablemente será diferente a la mía):

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 $a$, $b$ y $\mu$ (en la media de la columna) son también bastante cerca de lo previsto.