Posso reconstruir uma distribuição normal do tamanho da amostra e dos valores mínimo e máximo? Eu posso usar o ponto médio para proxy da média

14

Eu sei que isso pode ser um pouco complicado, estatisticamente, mas esse é o meu problema.

Eu tenho muitos dados de intervalo, ou seja, o tamanho mínimo, máximo e amostral de uma variável. Para alguns desses dados, também tenho uma média, mas não muitos. Quero comparar esses intervalos entre si para quantificar a variabilidade de cada intervalo e também para comparar as médias. Tenho um bom motivo para supor que a distribuição é simétrica em torno da média e que os dados terão uma distribuição gaussiana. Por esse motivo, acho que posso justificar o uso do ponto médio da distribuição como proxy da média, quando ela estiver ausente.

O que eu quero fazer é reconstruir uma distribuição para cada intervalo e usá-la para fornecer um desvio ou erro padrão para essa distribuição. As únicas informações que tenho são os máximos e mínimos observados em uma amostra e o ponto médio como proxy da média.

Dessa forma, espero poder calcular as médias ponderadas de cada grupo e também calcular o coeficiente de variação de cada grupo, com base nos dados de faixa que tenho e em minhas suposições (de distribuição simétrica e normal).

Eu pretendo usar o R ​​para fazer isso, então qualquer ajuda de código também seria apreciada.

green_thinlake
fonte
2
Fiquei me perguntando por que você diz que tem dados para valores mínimos e máximos e máximos; depois, você terá informações apenas sobre o mínimo e o máximo esperados. Qual é - observado ou esperado?
Scortchi - Reinstate Monica
Desculpe, esse é o meu erro. Os dados máximos e mínimos são observados (medidos a partir de objetos da vida real). Eu alterei o post.
22414 Greenlight (

Respostas:

11

A função de distribuição cumulativa conjunta para o mínimo de e máximo de x ( n ) para uma amostra de n de uma distribuição gaussiana com μ médio e desvio padrão σ éx(1)x(n)nμσ

F(x(1),x(n);μ,σ)=Pr(X(1)<x(1),X(n)<x(n))=Pr(X(n)<x(n))Pr(X(1)>x(1),X(n)<x(n)=Φ(x(n)μσ)n[Φ(x(n)μσ)Φ(x(1)μσ)]n

onde é o CDF gaussiano padrão. A diferenciação em relação a e fornece a função de densidade de probabilidade conjuntax ( 1 ) x ( n )Φ()x(1)x(n)

f(x(1),x(n);μ,σ)=n(n-1)[Φ(x(n)-μσ)-Φ(x(1)-μσ)]n-2ϕ(x(n)-μσ)ϕ(x(1)-μσ)1σ2

onde é o PDF gaussiano padrão. Tomar o log e soltar os termos que não contêm parâmetros fornece a função de probabilidade de logϕ()

(μ,σ;x(1),x(n))=(n2)log[Φ(x(n)μσ)Φ(x(1)μσ)]+logϕ(x(n)μσ)+logϕ(x(1)μσ)2logσ

Isso não parece muito tratável, mas é fácil ver que ele é maximizado, seja qual for o valor de , definindo , ou seja, o ponto médio - o primeiro termo é maximizado quando o argumento de um CDF é negativo do argumento do outro; os segundo e terceiro termos representam a probabilidade conjunta de duas variáveis ​​normais independentes.μ = μ = x ( n ) + x ( 1 )σμ=μ^=x(n)+x(1)2

Substituindo na probabilidade do log e escrevendo fornece r=x(n)-X(1)(σ;x(1),x(n), μ )=(n-2)log[1-2Φ( - rμ^r=x(n)x(1)

(σ;x(1),x(n),μ^)=(n2)log[12Φ(r2σ)]r24σ22logσ

Esta expressão deve ser maximizada numericamente (por exemplo, com optimizeo statpacote de R ) para encontrar . (Acontece que , em que é uma constante dependendo apenas de talvez alguém mais matematicamente hábil do que eu poderia mostrar o porquê.) σ =k(n)Rknσ^σ^=k(n)rkn

As estimativas não são úteis sem uma medida de precisão. As informações de Fisher observadas podem ser avaliadas numericamente (por exemplo, com hessiano numDerivpacote R ) e usadas para calcular erros padrão aproximados:

I(σ)=-2(σ; μ )

Eu(μ)=-2(μ;σ^)(μ)2|μ=μ^
Eu(σ)=-2(σ;μ^)(σ)2|σ=σ^

Seria interessante comparar a probabilidade e as estimativas do método dos momentos para em termos de viés (o MLE é consistente?), Variância e erro do quadrado médio. Há também a questão da estimativa para os grupos em que a média da amostra é conhecida além do mínimo e do máximo.σ

Scortchi - Restabelecer Monica
fonte
1
+1. Adicionando constante a para a log-probabilidade não irá alterar a localização do seu máximo, mas converte-lo em função da e , onde o valor de que maximiza é alguns função . Equivalentemente, como você afirma. Em outras palavras, a quantidade relevante para trabalhar é a razão do desvio padrão para a faixa (observada), ou igualmente bem sua recíproca - que está intimamente relacionada à faixa Studentizada . 2registro(r)σ/rnσ/rnk(n)σ^=k(n)r
whuber
@ whuber: Obrigado! Parece óbvio em retrospectiva. Vou incorporar isso na resposta.
Scortchi - Restabelece Monica
1

Você precisa relacionar o intervalo ao desvio / variância padrão. seja a média, o desvio padrão e seja o intervalo. Então, para a distribuição normal, temos que % da massa de probabilidade está dentro de 3 desvios padrão da média. Como regra prática, isso significa que, com uma probabilidade muito alta,μσR=x(n)-x(1)99,7

μ+3σx(n)
e

μ-3σx(1)

Subtraindo o segundo do primeiro, obtemos

6σx(n)-x(1)=R
(a propósito, é daí que vem a metodologia de garantia de qualidade "six-sigma" na indústria). Em seguida, você pode obter uma estimativa para o desvio padrão por onde a barra indica médias. É quando você assume que todas as subamostras são da mesma distribuição (você escreveu sobre ter intervalos esperados ). Se cada amostra é um normal diferente, com média e variância diferentes, é possível usar a fórmula para cada amostra, mas a incerteza / imprecisão possível no valor estimado do desvio padrão será muito maior.
σ^=16(x¯(n)-x¯(1))

Ter um valor para a média e para o desvio padrão caracteriza completamente a distribuição normal.

Alecos Papadopoulos
fonte
3
Isso não é uma aproximação aproximada para pequeno nem um resultado assintótico para grande . nn
Scortchi - Restabelece Monica
1
@Stortchi Bem, eu não disse que é uma boa estimativa, mas acredito que é sempre bom ter soluções facilmente implementadas, mesmo que sejam muito difíceis, para ter uma noção quantitativa do problema em mãos, além de outras. abordagens sofisticadas e eficientes, como por exemplo a descrita na outra resposta a essa pergunta.
Alecos Papadopoulos
Eu não insistiria com "a expectativa do intervalo amostral acaba por ser cerca de 6 vezes o desvio padrão para valores de de 200 a 1000". Mas estou sentindo falta de algo sutil em sua derivação, ou não funcionaria tão bem quanto justificar dividir o intervalo por qualquer número? n
Scortchi - Restabelece Monica
@ Scortchi Bem, o espírito da abordagem é "se esperamos que quase todas as realizações caiam dentro de 6 sigmas, é razoável esperar que as realizações extremas estejam próximas da fronteira" - isso é tudo o que realmente existe. Talvez eu estou muito acostumado a operar sob uma informação extremamente incompleta, e obrigado a dizer algo quantitativa sobre isso ... :)
Alecos Papadopoulos
4
Eu poderia responder que ainda mais observações ficariam dentro de da média, dando uma estimativa melhor . Não devo, porque é um absurdo. Qualquer número acima de será uma estimativa aproximada para algum valor de . 10σσ^=R101,13n
Scortchi - Reinstate Monica
1

É simples obter a função de distribuição do máximo da distribuição normal (consulte "P.max.norm" no código). A partir dele (com algum cálculo), você pode obter a função quantil (consulte "Q.max.norm").

Usando "Q.max.norm" e "Q.min.norm", é possível obter a mediana do intervalo relacionado a N. Usando a ideia apresentada por Alecos Papadopoulos (na resposta anterior), é possível calcular sd.

Tente o seguinte:

N = 100000    # the size of the sample

# Probability function given q and N
P.max.norm <- function(q, N=1, mean=0, sd=1){
    pnorm(q,mean,sd)^N
} 
# Quantile functions given p and N
Q.max.norm <- function(p, N=1, mean=0, sd=1){
    qnorm(p^(1/N),mean,sd)
} 
Q.min.norm <- function(p, N=1, mean=0, sd=1){
    mean-(Q.max.norm(p, N=N, mean=mean, sd=sd)-mean)
} 

### lets test it (takes some time)
Q.max.norm(0.5, N=N)  # The median on the maximum
Q.min.norm(0.5, N=N)  # The median on the minimum

iter = 100
median(replicate(iter, max(rnorm(N))))
median(replicate(iter, min(rnorm(N))))
# it is quite OK

### Lets try to get estimations
true_mean = -3
true_sd = 2
N = 100000

x = rnorm(N, true_mean, true_sd)  # simulation
x.vec = range(x)                  # observations

# estimation
est_mean = mean(x.vec)
est_sd = diff(x.vec)/(Q.max.norm(0.5, N=N)-Q.min.norm(0.5, N=N))

c(true_mean, true_sd)
c(est_mean, est_sd)

# Quite good, but only for large N
# -3  2
# -3.252606  1.981593
Vyga
fonte
2
E(R)=σ-1-(1-Φ(x))n-Φ(x)ndx=σd2(n)RΦ()d2nn