Intervalos de tolerância não paramétricos para variáveis ​​discretas

8

Suponha que você tenha um monte de gente avaliando o quanto eles gostaram de um filme em uma escala discreta de 1 a 10, e você gostaria de um intervalo [ l , u ] tal que com (pelo menos) 95% de confiança, (pelo menos) 90 % de todas as pessoas que assistem ao filme classificá-lo-á não inferior a 1 e superior a u . [ l , u ] é então um intervalo de tolerância (bilateral) com 95% de confiança e 90% de cobertura. (Para ser claro, 95% de confiança implica que, se você repetir esse procedimento várias vezes, 95% dos intervalos produzidos receberão pelo menos 90% de cobertura populacional.) É claro que geralmente queremos que [ l , u ] seja o mais estreito possível. possível enquanto ainda atendemos aos nossos requisitos.

Eu já vi vários métodos não paramétricos para construir intervalos de tolerância para variáveis ​​aleatórias contínuas. Também vi métodos para construir intervalos de tolerância para variáveis ​​binomiais e de Poisson. (O pacote R toleranceimplementa vários desses métodos; Young, 2010.) Mas e as variáveis ​​discretas quando a distribuição é desconhecida? Esse é geralmente o caso de escalas de classificação como a do meu exemplo, e supondo que uma distribuição binomial não pareça segura porque dados reais em escala de classificação geralmente exibem estranheza, como a multimodalidade.

Faria sentido recorrer aos métodos não paramétricos para variáveis ​​contínuas? Como alternativa, o que dizer de um método de Monte Carlo, como gerar 1.000 réplicas de autoinicialização da amostra e encontrar um intervalo que captura pelo menos 90% da amostra em pelo menos 950 das réplicas?

Young, DS (2010). tolerance: Um pacote R para estimar intervalos de tolerância. Journal of Statistical Software, 36 (5), 1–39. Recuperado em http://www.jstatsoft.org/v36/i05

Kodiologist
fonte
você quer dizer binomial ou multinomial? multinomial permitiria comportamento multimodal?
seanv507
Eu quero dizer binomial. No caso de uma escala de classificação, por exemplo, você definiria o número de ensaios de Bernoulli com o número de pontos da escala. Intervalos entre as categorias de uma distribuição multinomial não faria muito sentido, eu acho, uma vez que as categorias são desordenadas.
Kodiologist
@Kodiologist, sua variável de resultado é uma "escala discreta de 1 a 10", mas isso significa que é uma resposta multinomial ordenada. (Ou que não estou recebendo alguma coisa?)
Jim
@ Jim "Multinomial ordenado" é um pouco de oxímoro. Em uma distribuição multinomial, a ordem das categorias é arbitrária.
Kodiologist

Respostas:

1

A variável de interesse é distribuída multinomialmente com probabilidades de classe (célula): . Além disso, as classes são dotadas de uma ordem natural.p1,p2,...,p10

Primeira tentativa: menor "intervalo preditivo" contendo90%

p     = [p1, ..., p10] # empirical proportions summing to 1
l     = 1
u     = length(p)
cover = 0.9

pmass = sum(p)

while (pmass - p[l] >= cover) OR (pmass - p[u] >= cover)
    if p[l] <= p[u]
       pmass = pmass - p[l]
       l     = l + 1  
    else # p[l] > p[u]
       pmass = pmass - p[u]
       u     = u - 1
    end        
end

Uma medida não-paramétrica da incerteza (por exemplo, variância, confiança) nas estimativas quantitativas de poderia realmente ser obtida por métodos padrão de autoinicialização .l,u

Segunda abordagem: "pesquisa de inicialização" direta

Abaixo, forneço código Matlab executável que aborda a questão diretamente de uma perspectiva de autoinicialização (o código não é idealmente vetorizado).

%% set DGP parameters:
p = [0.35, 0.8, 3.5, 2.2, 0.3, 2.9, 4.3, 2.1, 0.4, 0.2];
p = p./sum(p); % true probabilities

ncat = numel(p);
cats = 1:ncat;

% draw a sample:
rng(1703) % set seed
nsamp = 10^3; 
samp  = datasample(1:10, nsamp, 'Weights', p, 'Replace', true);

Verifique se isso faz sentido.

psamp = mean(bsxfun(@eq, samp', cats)); % sample probabilities
bar([p(:), psamp(:)])

insira a descrição da imagem aqui

Execute a simulação de autoinicialização.

%% bootstrap simulation:
rng(240947)

nboots = 2*10^3;
cover  = 0.9;
conf   = 0.95;    

tic
Pmat = nan(nboots, ncat, ncat);
for b = 1:nboots

    boot  = datasample(samp, nsamp, 'Replace', true); % draw bootstrap sample    
    pboot = mean(bsxfun(@eq, boot', cats));      

    for l = 1:ncat
        for u = l:ncat
            Pmat(b, l, u) = sum(pboot(l:u));   
        end
    end

end
toc % Elapsed time is 0.442703 seconds.

O filtro de cada autoinicialização replica os intervalos, , que contêm pelo menos massa de probabilidade e calcula uma estimativa de confiança (freqüentista) desses intervalos.90 %[l,u]90%

conf_mat = squeeze(mean(Pmat >= cover, 1))

     0         0         0         0         0         0         0    1.0000    1.0000    1.0000
     0         0         0         0         0         0         0    1.0000    1.0000    1.0000
     0         0         0         0         0         0         0    0.3360    0.9770    1.0000
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0

Selecione aqueles que satisfazem a confiança desejada.

[L, U] = find(conf_mat >= conf);
[L, U]

 1     8
 2     8
 1     9
 2     9
 3     9
 1    10
 2    10
 3    10

Convencendo-se de que o método de inicialização acima é válido

As amostras de bootstrap destinam-se a substituir algo que gostaríamos de ter, mas não o fazem, ou seja: novos e independentes resultados da verdadeira população subjacente (abreviação: novos dados).

No exemplo que eu dei, conhecemos o processo de geração de dados (DGP), portanto, poderíamos "trapacear" e substituir as linhas de código referentes às novas amostras de bootstrap por novos e independentes sorteios do DGP real.

newsamp = datasample(cats, nsamp, 'Weights', p, 'Replace', true);
pnew    = mean(bsxfun(@eq, newsamp', cats));

Em seguida, podemos validar a abordagem de bootstrap comparando-a com o ideal. Abaixo estão os resultados.

A matriz de confiança de dados novos e independentes extrai:

     0         0         0         0         0         0         0    1.0000    1.0000    1.0000
     0         0         0         0         0         0         0    1.0000    1.0000    1.0000
     0         0         0         0         0         0         0    0.4075    0.9925    1.0000
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0
     0         0         0         0         0         0         0         0         0         0

Os limites inferior e superior de confiança correspondente a :95%

 1     8
 2     8
 1     9
 2     9
 3     9
 1    10
 2    10
 3    10

Concluímos que as matrizes de confiança concordam estreitamente e que os limites são idênticos ... Validar a abordagem de autoinicialização.

Jim
fonte
2
Intervalos de tolerância e confiança são coisas diferentes. Na verdade, o que você descreveu não é um intervalo de confiança, mas um intervalo de previsão, que é outro tipo distinto de intervalo.
Kodiologist 14/09/16
1
Sua edição parece ser uma implementação do que eu quis dizer quando escrevi "um método de Monte Carlo, como gerar 1.000 repetições de autoinicialização da amostra e encontrar um intervalo que capture pelo menos 90% da amostra em pelo menos 950 das replicadas". Embora intuitivo, não tenho certeza se isso realmente funciona ou faz sentido, e foi por isso que criei essa pergunta.
Kodiologist 15/09/16
@ Kodiologist A resposta agora contém uma seção validando a abordagem de inicialização. Obviamente, isso poderia ser levado adiante, por exemplo, aninhado em loops sobre o tamanho da amostra e as probabilidades de classe.
Jim Jim
Mostrar que o método de autoinicialização está correto para esse problema em sua generalidade total significa mostrar que ele possui a confiança e a cobertura corretas, independentemente da distribuição anterior dos parâmetros (é, afinal, um método freqüentista). Para isso, acho que a simulação não seria suficiente; você precisaria de prova matemática. Mas você foi extraordinariamente persistente e mostrou que o bootstrap funciona pelo menos algumas vezes, então você merece um A por esforço.
Kodiologist 16/09/16