Resolução analítica de amostras com ou sem reposição após binômio Poisson / Negativo

8

Versão curta

Estou tentando resolver / aproximar analiticamente a probabilidade composta que resulta de desenhos independentes de Poisson e de amostras adicionais com ou sem substituição (eu realmente não me importo com qual). Quero usar a probabilidade com o MCMC (Stan), portanto, preciso da solução apenas até um termo constante. Em última análise, quero modelar um processo em que os desenhos iniciais sejam negativos. distribuição binomial, mas acho que vou conseguir chegar lá com uma solução para o caso Poisson.

É bem possível que a solução não seja viável (não entendo a matemática o suficiente para saber se esse é um problema simples ou muito difícil). Assim, também estou interessado em aproximações, resultados negativos ou intuição por que o problema é provavelmente intratável (por exemplo, comparando com um problema difícil conhecido). Links para artigos / teoremas / truques úteis que me ajudarão a seguir em frente são boas respostas, mesmo que sua conexão com o problema em questão não esteja totalmente resolvida.

Declaração formal

Mais formalmente, primeiro é desenhado independentemente e, em seguida, eu amostro itens aleatoriamente de todo para obter . Ou seja, eu desenho bolas coloridas de uma urna onde a quantidade de bolas de cor é extraída de . Aqui, é assumido conhecido e corrigido e condicionamos em . Tecnicamente, a amostragem é feita sem substituição, mas presumir que a amostragem com substituição não deve ser um grande problema.Y=(y1,...,yN),ynPois(λn)kYZ=(z1,...,zN)knPois(λn)knynk

Eu tentei duas abordagens para resolver a amostragem sem substituição (como esse parecia ser o caso mais fácil devido a alguns termos cancelados), mas fiquei com as duas. A probabilidade de amostragem sem substituição é:

P(Z=(z1,...,zN)|Λ=(λ1,...,λN))=Y;n:ynzn(n=1N(ynzn)(n=1Nynk)n=1NPoisson(yn|λn))P(nynk|Λ)

EDIT: A seção "tentativas de soluções foi removida porque a solução na resposta não se baseia nelas (e é muito melhor)

Martin Modrák
fonte

Respostas:

3

O caso sem substituição

Se você possui variáveis ​​distribuídas independentes de Poissonn

YiPois(λi)

e condição em

j=inYj=K

então

{Yi}Multinom(K,(λij=1nλj))

Assim, você pode encher sua urna com as vezes bolas coloridas como primeiro desenhar o valor do total (que é o ponto de corte distribuído por Poisson pela condição ) e depois encher a urna com bolas de acordo com a distribuição multinomial.nYiKKkK

Esse preenchimento da urna com bolas , de acordo com uma distribuição multinomial, é equivalente a desenhar para cada bola independentemente a cor da distribuição categórica. Em seguida, você pode considerar as primeiras bolas que foram adicionadas à urna como definindo a amostra aleatória (quando essa amostra é extraída sem substituição) e a distribuição para isso é apenas outro vetor distribuído multinomial: Kk{Zi}

{Zi}Multinom(k,(λij=1nλj))

simulação

##### simulating sample process for 3 variables #######


# settings
set.seed(1)
k = 10
lambda = c(4, 4, 4)
trials = 1000000

# observed counts table
Ocounts = array(0, dim = c(k+1, k+1, k+1))

for (i in c(1:trials)) {
  # draw poisson with limit sum(n) >= k
  repeat {
    Y = rpois(3,lambda)
    if (sum(Y) >= k) {break}
  }
  # setup urn
  urn <- c(rep(1, Y[1]), rep(2, Y[2]), rep(3, Y[3]))
  # draw from urn
  s <- sample(urn, k, replace=0)
  Z = c(sum(s==1),sum(s==2),sum(s==3))
  Ocounts[Z[1]+1, Z[2]+1, Z[3]+1] = Ocounts[Z[1]+1, Z[2]+1, Z[3]+1] + 1
}



# comparison
observed = rep(0, 0.5*k*(k+1))
expected = rep(0, 0.5*k*(k+1))   
t = 1

for (z1 in c(0:k)) {
  for (z2 in c(0:(k-z1))) {  
    z3 = k-z1-z2 
    observed[t] = Ocounts[z1+1, z2+1, z3+1]
    expected[t] = trials*dmultinom(c(z1, z2, z3), prob=lambda)
    t = t+1
  }
}

plot(observed,expected)
x2 <- sum((observed-expected)^2/expected)
pvalue <- 1-pchisq(x2, 66-1)

resultados

> # results from chi-sq test
> x2
[1] 75.49286
> pvalue
[1] 0.1754805

Binomial negativo

Os argumentos funcionariam da mesma forma no caso de uma distribuição binomial negativa que resulta ( sob certas condições ) em uma distribuição multinomial Dirichlet.

Abaixo está um exemplo de simulação

##### simulating sample process for 3 variables #######

# dirichlet multinomial for vectors of size 3
ddirmultinom =  function(x1,x2,x3,p1,p2,p3) {
  (factorial(x1+x2+x3)*gamma(p1+p2+p3)/gamma(x1+x2+x3+p1+p2+p3))/
  (factorial(x1)*gamma(p1)/gamma(x1+p1))/
  (factorial(x2)*gamma(p2)/gamma(x2+p2))/
  (factorial(x3)*gamma(p3)/gamma(x3+p3))
}

# settings
set.seed(1)
k = 10
theta = 1
lambda = c(4,4,4)
trials = 1000000

# calculating negative binomials pars
means = lambda
vars = lambda*(1+theta)

ps = (vars-means)/(vars)
rs = means^2/(vars-means)


# observed counts table
Ocounts = array(0, dim = c(k+1,k+1,k+1))

for (i in c(1:trials)) {
  # draw poisson with limit sum(n) >= k
  repeat {
    Y = rnbinom(3,rs,ps)
    if (sum(Y) >= k) {break}
  }
  # setup urn
  urn <- c(rep(1,Y[1]),rep(2,Y[2]),rep(3,Y[3]))
  # draw from urn
  s <- sample(urn,k,replace=0)
  Z = c(sum(s==1),sum(s==2),sum(s==3))
  Ocounts[Z[1]+1,Z[2]+1,Z[3]+1] = Ocounts[Z[1]+1,Z[2]+1,Z[3]+1] + 1
}



# comparison
observed = rep(0,0.5*k*(k+1))
expected = rep(0,0.5*k*(k+1))   
t = 1

for (z1 in c(0:k)) {
  for (z2 in c(0:(k-z1))) {  
    z3 = k-z1-z2 
    observed[t]=Ocounts[z1+1,z2+1,z3+1]
    expected[t]=trials*ddirmultinom(z1,z2,z3,lambda[1]/theta,lambda[2]/theta,lambda[3]/theta)
    t = t+1
  }
}

plot(observed,expected)
x2 <- sum((observed-expected)^2/expected)
pvalue <- 1-pchisq(x2,66-1)

# results from chi-sq test
x2
pvalue
Sextus Empiricus
fonte
1
Obrigado pela resposta. Eu já tentei essa abordagem e tenho dois problemas: a) não vejo por que o condicionamento na soma é o mesmo que reamostragem (matematicamente) eb) meus dados (reconhecidamente limitados) simulados da reamostragem Poisson + não se encaixavam bem com distribuição multinomial. Você poderia explicar mais sobre por que o condicionamento e a reamostragem seriam os mesmos?
Martin Modrák 03/08/1918
Oh, em uma leitura mais aprofundada, acho que entendi o seu ponto. Talvez eu tenha cometido um erro estúpido com as simulações, vá verificar.
Martin Modrák
Eu uso dois passos. 1) encher a urna com bolas coloridas desenhando de distribuições independentes de Poisson, é equivalente a encher a urna desenhando um total da distribuição Poisson e depois o de uma distribuição multinomial. 2) desenhar um subconjunto de bolas de uma urna cheia de bolas coloridas cujos números são definidos por uma distribuição multinomial, equivale à distribuição multinomial do anoter (intuição, você pode ver a urna como sendo preenchida pela cor de cada bola que um empate a distribuição categórica). YiYiYi=KYi
Sextus Empiricus
Sim, houve um erro no meu código de simulação :-) Obrigado pela ajuda! Você está correto e eu também entendo a etapa do seu raciocínio que eu não vi anteriormente. Estou tentando entender por que a mesma relação não vale para neg. distribuição multinomial binomial e Dirichlet. (supondo que as variáveis ​​RN tenham fatores Fano constantes, condicionando sua soma em obter DM) - porque "amostragem sem substituição de multinomial" é multinomial, mas "amostragem sem substituição de DM" não é DM (esse raciocínio está correto na sua opinião? )
Martin Modrák
1
Acontece que eu também estraguei a segunda simulação. Parece funcionar também para o caso do DM.
Martin Modrák