Distribuição Preditiva Frequentista de uma variável de Cauchy

7

Não consegui encontrar isso na literatura, mas isso provavelmente significa que estou procurando no lugar errado. Eu estou procurando encontrar a distribuição preditiva freqüentista, supondo que ela exista, para uma variável Cauchy unidimensional e n-dimensional.

O problema com a versão n-dimensional é que não há nada como uma matriz covariável; em vez disso, existe apenas um parâmetro de escala que torna os erros hiper-circulares. Eu pude ver isso interferindo na existência de um valor crucial.

EDITAR

Eu estou procurando prever partir de um conjunto de observações extraído de uma distribuição de Cauchy com centro e escala \ sigma, ou prever y_ {i + 1} a partir de alguma equação y = mx + b, onde x é extraído de uma distribuição de Cauchy como acima. Pode ser um vetor ou multidimensional, mas estou tentando determinar as propriedades relativas da previsão Bayesiana versus Frequentista. Meus dados são extraídos de um Cauchy truncado ou um Cauchy, dependendo de qual conjunto.xi+1x1xiμσ,yi+1y=mx+b,x

Um intervalo de previsão funcionará, pois apenas definirei o intervalo para 100%.

Dave Harris
fonte
11
Quando diz Cauchy variado , você quer dizer um regressor na equação de regressão , e procura a confiança intervalo da estimativa OLS para ? Não sei se entendi a pergunta. Xj,iYi=j=1JXj,iβj+εiβj
precisa
11
Existem trabalhos, como este de 2008 , sobre estimativa com distribuições t multivariadas. O Cauchy multivariado é um caso especial do t multivariado. Como um aparte, isso permite uma estrutura de correlação totalmente flexível. Isso ajuda ou não corresponde à sua pergunta?
Eric_kernfeld
11
Você está tentando fazer isso stats.stackexchange.com/questions/16349 para uma distribuição multivariada com média zero?
Sextus Empiricus
11
@eric_kernfeld Eu tenho que ler com atenção, mas sim, é isso, exceto que eu quero saber como encontrar sua densidade preditiva usando métodos Frequentist.
Dave Harris
11
parece que você está tentando estimar os parâmetros de distribuição Cauchy de . Isso está certo? xi
Aksakal

Respostas:

2

A solução geral para o seu problema é a estimativa de máxima verossimilhança (MLE) dos seus parâmetros . Uma vez que eles são obtidos comoθθ^, você os substitui no seu pdf pelos parâmetros desconhecidos, ou seja, estima o pdf da sua variável aleatória como f^(xi)=f(xi|θ^). Isso permite que você construa a distribuição preditiva de sua variável aleatória Cauchy.

Para o caso univariado , este artigo é um excelente recurso . Para o Cauchy univariado com centroμ e escala σ, um tem um formulário fechado se você tiver 34observações. Se você temn>4 observações, o MLE existe. Se você temnobservações, você terá que resolver duas equações que são facilmente derivadas definindo a primeira derivada da probabilidade logarítmica como zero; veja aqui a forma exata. (Na notação deles,x0=μ e σ=γ.) A solução desse problema numericamente tem uma implementação na linguagem R, veja aqui .

Para o caso multivariado , tudo o que você precisa observar é que a distribuição multivariada de Cauchy é simplesmente uma distribuição multivariadat-distribuição em que o parâmetro grau de liberdade está definido como 1, como já foi apontado nos comentários. Para os multivariadost, você pode fazer inferência do MLE como explicado excelentemente nesta resposta , que se baseia no artigo que eric_kernfeld apontou. Não encontrei uma implementação pronta para esse algoritmo, mas como você verá quando examinar a resposta fornecida no post, será realmente fácil implementá-la.

Diferença da previsão bayesiana : No cenário bayesiano, você colocaria um prévio nos parâmetrosμ e σ, modelando sua incerteza sobre eles como uma variável aleatória. Assim, você obterá distribuições posteriores para ambos os parâmetros, que indicam a certeza relativa que você tem sobre eles, dados seus dados. Se você tem a parte posteriorq(μ,σ|x1,,xn), você obtém sua distribuição preditiva como f(x|μ,σ)q(μ,σ|x1,,xn)dμdσ, integrando sua incerteza. Por outro lado, a configuração do MLE fornecerá estimativas pontuais deμ e σque você conecte no formulário funcional do seu pdf. Equivalentemente, você poderia dizer que o MLE leva a um posterior com massa pontual1 na tupla (μ^,σ^) e 0probabilidade em qualquer outro valor. Assim, você ignora toda a incerteza de parâmetro nesse caso e depende do fato de queθ^ é assintoticamente equivalente a θ, significa que f^(x)f(x) (uniformemente x)

Bem, isso é a menos que no caso exótico em que n é par e n/2 das suas observações têm valor x1 enquanto a outra metade valoriza x2, o que acontece com probabilidade zero porque a distribuição de Cauchy é contínua.

Jeremias K
fonte
Jeremias. Você acha que existem possibilidades para incorporar a incerteza sobreθ^no intervalo de previsão? E como construímos uma distribuição preditiva a partir do pdf de uma distribuição cauchy multidimensional?
Sextus Empiricus
Se você assumir que o xi são sorteios aleatórios de uma variável aleatória Cauchy com parâmetros desconhecidos, a inserção direta dos parâmetros estimados na forma funcional fornece a distribuição preditiva dos próximos sorteios de xi.
precisa
Quanto à incorporação de incerteza de parâmetro, se você quiser fazer isso, terá que seguir o caminho bayesiano. Observe que, como produto secundário, o parâmetro posterior da inferência bayesiana concentra-se no MLE por meio do teorema de Bernstein Mises assintoticamente.
precisa
1

Pode-se usar um método de Monte Carlo para obter estimativas empíricas para as relações entre os x1....xi e o intervalo de previsão para xi+n.

Motivação: Se estimarmos o intervalo de previsão com base nos quartis / CDF de uma distribuição que segue das estimativas de probabilidade máxima (ou outro tipo de estimativa de parâmetros), subestimamos o tamanho do intervalo. Efetivamente, na prática, o pontoxi+n cairá fora da faixa com mais frequência do que o previsto.

A figura abaixo demonstra o quanto subestimamos o tamanho do intervalo, expressando quantas vezes mais uma nova medição xiestá fora do intervalo preditivo com base nas estimativas de parâmetros. (com base em cálculos com 2000 repetições para a previsão)

Por exemplo, se usarmos um intervalo de previsão de 99% (portanto, esperamos 1% de erros), obteremos 5 vezes mais erros se o tamanho da amostra for 3.

Esse tipo de cálculo pode ser usado para estabelecer relações empíricas sobre como podemos corrigir o intervalo, assim como os cálculos mostram que, para grandes n a diferença se torna menor (e em algum momento pode-se considerar irrelevante).

diferença entre a estimativa do MLE e o intervalo de confiança efetivo

set.seed(1)

# likelihood calculation
like<-function(par, x){
  scale = abs(par[2])
  pos   = par[1]
  n <- length(x)
  like <- -n*log(scale*pi) - sum(log(1+((x-pos)/scale)^2))
  -like
}

# obtain effective predictive failure rate rate
tryf <- function(pos, scale, perc, n) {

  # random distribution
  draw <- rcauchy(n, pos, scale)

  # estimating distribution parameters based on median and interquartile range
  first_est <- c(median(draw), 0.5*IQR(draw))

  # estimating distribution parameters based on likelihood
  out <- optim(par=first_est, like, method='CG', x=draw)
  # making scale parameter positive (we used an absolute valuer in the optim function)
  out$par[2] <- abs(out$par[2])

  # calculate predictive interval
  ql <- qcauchy(perc/2, out$par[1], out$par[2])
  qh <- qcauchy(1-perc/2, out$par[1], out$par[2])

  # calculate effective percentage outside predicted predictive interval
  pl <- pcauchy(ql, pos, scale)
  ph <- pcauchy(qh, pos, scale)
  error <- pl+1-ph
  error
}

# obtain mean of predictive interval in 2000 runs
meanf <- function(pos,scale,perc,n) {
  trueval <- sapply(1:2000,FUN <- function(x) tryf(pos,scale,perc,n))
  mean(trueval)
}


#################### generate image

# x-axis chosen desired interval percentage
percentages <- 0.2/1.2^c(0:30)

# desired sample sizes n
ns <- c(3,4,5,6,7,8,9,10,20,30)

# computations
y <- matrix(rep(percentages, length(ns)), length(percentages))
for (i in which(ns>0)) {
  y[,i] <- sapply(percentages, FUN <- function(x) meanf(0,1,x,ns[i]))
}

# plotting
plot(NULL,
     xlim=c(0.0008,1), ylim=c(0,10),
     log="x",
     xlab="aimed error rate",
     ylab="effective error rate / aimed error rate",
     yaxt="n",xaxt="n",axes=FALSE)
axis(1,las=2,tck=-0.0,cex.axis=1,labels=rep("",2),at=c(0.0008,1),pos=0.0008)
axis(1,las=2,tck=-0.005,cex.axis=1,at=c(0.001*c(1:9),0.01*c(1:9),0.1*c(1:9)),labels=rep("",27),mgp=c(1.5,1,0),pos=0.0008)
axis(1,las=2,tck=-0.01,cex.axis=1,labels=c(0.001,0.01,0.1,1), at=c(0.001,0.01,0.1,1),mgp=c(1.5,1,0),pos=0.000)
#axis(2,las=1,tck=-0.0,cex.axis=1,labels=rep("",2),at=c(0.0008,1),pos=0.0008)
#axis(2,las=1,tck=-0.005,cex.axis=1,at=c(0.001*c(1:9),0.01*c(1:9),0.1*c(1:9)),labels=rep("",27),mgp=c(1.5,1,0),pos=0.0008)
#axis(2,las=1,tck=-0.01,cex.axis=1,labels=c(0.001,0.01,0.1,1), at=c(0.001,0.01,0.1,1),mgp=c(1.5,1,0),pos=0.0008)
axis(2,las=2,tck=-0.01,cex.axis=1,labels=0:15, at=0:15,mgp=c(1.5,1,0),pos=0.0008)


colours <- hsv(c(1:10)/20,1,1-c(1:10)/15)
for (i in which(ns>0)) {
  points(percentages,y[,i]/percentages,pch=21,cex=0.5,col=colours[i],bg=colours[i])
}

legend(x=0.4,y=4.5,pch=21,legend=ns,col=colours,pt.bg=colours,title="sample size")

title("difference between confidence interval and effective confidence interval")


plot(ns,y[31,]/percentages[31],log="")
Sextus Empiricus
fonte
O que o gráfico nos diz que, além de usar um tamanho pequeno de amostra, produzirá uma estimativa ruim de seus parâmetros ao usar mle ? Não vejo como isso invalida o uso do mle, uma vez que as taxas de erro parecem excelentes, mesmo para um tamanho de amostra muito pequeno de 30. Também não tenho certeza de entender qual é a alternativa que você propõe, você se importaria em expandir os métodos computacionais você mencionou no início da sua resposta?
Jeremias K
11
@ JeremiasK Em aplicações práticas, com amostras pequenas, pode-se usar esses cálculos como fatores de correção empiricamente determinados.
Sextus Empiricus
Isso faz sentido! Eu não acho que você mencione isso no post, talvez você deva editá-lo para que as pessoas não precisem ler os comentários #
Jeremias K
@MartijnWeterings até agora você faz mais sentido. O pivôn(μ^μ)σ^$ segue o padrão normal quando o tamanho da amostra chega a cerca de 100, mas percebi que estou além das minhas habilidades para relaxar isso porque, em vez de escolher uma variável, estou escolhendo uma função para a minimização e ainda não o fiz antes.
Dave Harris
@DaveHarris Acredito que meu método não seja tão diferente do de Jeremia, exceto que faço uma expressão (e apenas por uma abordagem experimental em matemática) para o intervalo subestimado que ocorre porque a distribuição f(x,x^0,γ^) é uma versão super dispersa do f(x,x0,γ).
Sextus Empiricus
0

Parece que tudo que você precisa é estimar os parâmetros da distribuição de Cauchy a partir do conjunto de dados xi. Aqui está o que Stephens propõe, não é o MLE, e o autor afirma que esse método é consistente e mais estável que o MLE, embora você deva levar em conta que isso foi escrito no século passado.

insira a descrição da imagem aqui

onde Cauchy é parametrizado da seguinte maneira: insira a descrição da imagem aqui

Depois de ter a distribuição, sua previsão de pontos será α^. Observe que, como não há momentos, você não poderá mostrar que sua previsão é ótima no sentido usual, como minimizar o custo quadrado esperado.

Aksakal
fonte