Estimer le taux auquel l'écart-type évolue avec une variable indépendante

J'ai une expérience dans laquelle je prends des mesures d'une variable normalement distribuée $Y$,

Cependant, des expériences antérieures ont fourni des preuves que l'écart-type est une fonction affine d'une variable indépendante , c'est-à-dire$\sigma$$X$

Je souhaite estimer les paramètres et par échantillonnage à valeurs multiples de . De plus, en raison des limites de l'expérience, je ne peux prendre qu'un nombre limité (environ 30 à 40) d'échantillons de , et je préférerais échantillonner à plusieurs valeurs de pour des raisons expérimentales non liées. Compte tenu de ces contraintes, quelles méthodes sont disponibles pour estimer et ?b Y X Y X a $a$$b$$Y$$X$$Y$$X$$a$$b$

Description de l'expérience

Il s'agit d'informations supplémentaires, si vous souhaitez savoir pourquoi je pose la question ci-dessus. Mon expérience mesure la perception spatiale auditive et visuelle. J'ai une configuration d'expérience dans laquelle je peux présenter soit des cibles visuelles ou auditives de différents endroits, et sujets indiquent l'emplacement perçu de la cible, . La vision * et l'audition deviennent moins précises avec l'augmentation de l'excentricité (c'est-à-dire en augmentant ), que je modélise comme ci-dessus. En fin de compte, je voudrais estimer etY | X | σ a $X$$Y$$|X|$$\sigma$$a$$b$pour la vision et l'audition, je connais donc la précision de chaque sens à travers une gamme d'emplacements dans l'espace. Ces estimations seront utilisées pour prédire la pondération relative des cibles visuelles et auditives lorsqu'elles sont présentées simultanément (similaire à la théorie de l'intégration multisensorielle présentée ici: http://www.ncbi.nlm.nih.gov/pubmed/12868643 ).

* Je sais que ce modèle est imprécis pour la vision lorsque l'on compare le fovéal à l'espace extrafovéal, mais mes mesures sont limitées uniquement à l'espace extrafovéal, où il s'agit d'une approximation décente.

Dans un cas comme le vôtre, où vous avez un modèle génératif relativement simple mais "non standard" pour lequel vous aimeriez estimer des paramètres, ma première pensée serait d'utiliser un programme d'inférence bayésien comme Stan . La description que vous avez donnée se traduirait très clairement en un modèle Stan.

Un exemple de code R, utilisant RStan (l'interface R de 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)

Vous obtiendrez une sortie qui ressemble à ceci (bien que vos nombres aléatoires soient probablement différents des miens):

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

Le modèle a bien convergé (Rhat = 1), et la taille réelle de l'échantillon (n_eff) est raisonnablement grande dans tous les cas, donc sur le plan technique, le modèle se comporte bien. Les meilleures estimations de , b et μ (dans la colonne moyenne) sont également assez proches de ce qui a été fourni.$a$$b$$\mu$

Vous ne pouvez pas vous attendre à des formules fermées, mais vous pouvez toujours noter la fonction de vraisemblance et la maximiser numériquement. Votre modèle est Alors la fonction loglik vraisemblance (à part un terme ne dépendant pas des paramètres) devient l ( μ , a , b ) = - ∑ ln

En R, on peut faire

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

Simulez ensuite quelques données:

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

Ensuite, faites la fonction loglikelihood:

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

Optimisez-le ensuite:

> 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