Estimar a taxa na qual o desvio padrão é escalonado com uma variável independente

11

Eu tenho um experimento no qual estou fazendo medições de uma variável normalmente distribuída ,Y

YN(μ,σ)

Entretanto, experimentos anteriores forneceram algumas evidências de que o desvio padrão é uma função afim de uma variável independente , ou seja,XσX

σ=a|X|+b

YN(μ,a|X|+b)

Gostaria para estimar os parâmetros de e por amostragem em vários valores de . Além disso, devido às limitações do experimento, só posso coletar um número limitado (aproximadamente 30-40) de amostras de e preferiria amostrar em vários valores de por razões experimentais não relacionadas. Dadas estas limitações, o que métodos estão disponíveis para estimar e ?b Y X Y X a babYXYXab

Descrição da experiência

Esta é uma informação extra, se você estiver interessado em saber por que estou fazendo a pergunta acima. Meu experimento mede a percepção espacial auditiva e visual. I têm uma configuração da experiência em que pode apresentar auditivo ou alvos visuais de diferentes locais, , e sujeitos indicam a localização percebida do alvo, . Tanto a visão * quanto a audição ficam menos precisas com o aumento da excentricidade (ou seja, o aumento de ), que eu modelo como acima. Finalmente, eu gostaria de estimar eY | X | σ a bXY|X|σabpara visão e audição, conheço a precisão de cada sentido em vários locais no espaço. Essas estimativas serão usadas para prever a ponderação relativa dos alvos visuais e auditivos quando apresentados simultaneamente (semelhante à teoria da integração multissensorial apresentada aqui: http://www.ncbi.nlm.nih.gov/pubmed/12868643 ).

* Eu sei que esse modelo é impreciso para a visão ao comparar o espaço foveal ao extrafoveal, mas minhas medidas são restritas apenas ao espaço extrafoveal, onde essa é uma aproximação decente.

Adam Bosen
fonte
2
Problema interessante. É provável que as melhores soluções levem em conta os motivos pelos quais você está fazendo esse experimento. Quais são seus objetivos finais? Predição? Estimativa de , e / ou ? Quanto mais você puder nos contar sobre o objetivo, melhores serão as respostas. a σμaσ
whuber
Como o SD não pode ser negativo, é improvável que seja uma função linear de X. Sua sugestão, a | X |, exige uma forma de V mais estreita ou mais estreita com um mínimo em X = 0, o que me parece uma possibilidade pouco natural. . Tem certeza de que isso está certo?
gung - Restabelece Monica
Bom ponto @gung, eu simplificara demais o meu problema. Seria mais realista dizer que é uma função afim de. Vou editar minha pergunta. | X |σ|X|
quer
@whuber A razão para querer isso está um pouco envolvida, mas vou pensar em como explicar o experimento e adicionar mais alguns detalhes à minha pergunta em breve.
Adam Bosen
1
Você tem um bom motivo, a priori, para acreditar que X = 0 representa o DP mínimo, e que f (| X |) é monotônico?
gung - Restabelece Monica

Respostas:

2

Em um caso como o seu, no qual você tem um modelo generativo relativamente simples, mas "não-padrão", para o qual gostaria de estimar parâmetros, meu primeiro pensamento seria usar um programa de inferência bayesiano como Stan . A descrição que você forneceu se traduziria muito claramente para um modelo de Stan.

Alguns exemplos de código R, usando RStan (a interface 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)

Você obterá resultados parecidos com este (embora seus números aleatórios provavelmente sejam diferentes dos meus):

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

O modelo convergiu bem (Rhat = 1), e o tamanho efetivo da amostra (n_eff) é razoavelmente grande em todos os casos, portanto, em um nível técnico, o modelo é bem-comportado. As melhores estimativas de , e (na coluna média) são também bastante próximo ao que foi fornecido.b μabμ

Martin O'Leary
fonte
Ah, eu gosto disso! Eu nunca tinha ouvido falar de Stan antes, obrigado pela referência. Inicialmente, esperava uma solução analítica, mas, dada a falta de respostas, duvido que exista uma. Estou inclinado a acreditar que sua resposta é a melhor abordagem para esse problema.
Adam Bosen
Não me chocaria completamente se existisse uma solução analítica, mas certamente ficaria um pouco surpreso. A força de usar algo como Stan é que é muito fácil fazer alterações em seu modelo - uma solução analítica provavelmente seria muito fortemente restringida ao modelo, conforme indicado.
Martin O'Leary
2

Você não pode esperar fórmulas fechadas, mas ainda pode anotar a função de probabilidade e maximizá-la numericamente. Seu modelo é Então a função de probabilidade de log (exceto um termo que não depende dos parâmetros) se torna e isso é fácil de programar e entregue a um otimizador numérico.

YN(μ,a|x|+b)
l(μ,a,b)=ln(a|xi|+b)12(yiμa|xi|+b)2

Em R, podemos fazer

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

Em seguida, simule alguns dados:

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

Em seguida, faça a função loglikelihood:

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

Em seguida, otimize-o:

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