As Figurinhas dos Vingadores são Honestas? (parte 2)

No último post do blog, eu comprei vários pacotinhos de figurinhas dos Vingadores, para tentar responder à pergunta: as figurinhas são honestas? Após fazer dois testes de hipótese, eu concluí que as figurinhas possuem distribuição uniforme, então parece que sim. Mas tem uma pegadinha: para serem honestas, garantir que são uniformes não é o suficiente, ainda tem outro teste que precisamos fazer!

A Distribuição dos Pacotinhos

Para fazer a análise das figurinhas, eu comprei 106 pacotinhos na banca. Mas enquanto eu abria os pacotinhos, notei uma propriedade curiosa: todos os pacotinhos tinham sempre 4 figurinhas diferentes. Mas olha só, se os pacotinhos não possuem repetidas entre si, então as figurinhas podem até ser uniformes, mas não são independentes! Cada pacotinho tem quatro figurinhas, e se eu tirei a figurinha 42 em um pacotinho, então eu sei que as outras três não são a figurinha 42.

Mas pode ter sido só o acaso né? Vai ver eu dei sorte de comprar exatamente 106 pacotinhos que não tinham repetidas. Mas felizmente nós sabemos como testar se foi sorte ou não, basta fazer o teste do ##\chi^2##, que vimos no último post!

Acho que o jeito mais fácil de testar é contar quantas figurinhas diferentes tem dentro do pacotinho. Quantos pacotinhos possíveis existem com 4 figurinhas diferentes? Esse é fácil, o álbum tem 180 figurinhas, então são ##180\choose 4## pacotinhos. Com 3 diferentes precisa de mais cuidado, existem ##180\choose 3## combinações de figurinhas, mas você precisa contar qual dessas três diferentes está repetida (os casos são AABC, ABBC, ABCC). Então para 3 diferentes são ##3 {180\choose 4}## pacotinhos.

Com 2 diferentes também tem três casos: AAAB, AABB, ABBB, logo temos ##3{180\choose 2}## pacotinhos. Para uma única figurinha diferente, são apenas ##180\choose 1## pacotinhos. Agora podemos montar a tabela com quantos pacotinhos eu deveria ter visto, e quantos pacotinhos eu vi:

	1 diferente	2 diferentes	3 diferentes	4 diferentes
Total de pacotinhos	##180\choose 1##	##3{180\choose 2}##	##3{180\choose 3}##	##180\choose 4##
Total (numérico)	180	48330	2867580	42296805
Total (normalizado)	3.98e-6	1.07e-3	0.634	0.936
Esperado (total*106)	422e-6	0.113	6.72	99.163
Observado	0	0	0	106

Opa, nas duas primeiras categorias ##np\lt 5##, então eu não tenho amostras suficientes para usar o ##\chi^2##! Um jeito de consertar isso sem comprar mais pacotinhos é juntando as três primeiras categorias em uma só:

	Menos de 4 diferentes	4 diferentes
Esperado	6.837	99.163
Observado	0	106

Com isso conseguimos calcular o ##\chi^2## pela definição, chegando no valor de ##\chi^2=7.308##. Qual é o limite para rejeitar a hipótese? Se estivermos trabalhando com ##\alpha=5\%##, e um grau de liberdade, então o limite é:

InverseCDF[ChiSquareDistribution[1], 1 - 0.05]
>> 3.84146

Como ##7.308\gt 3.84146##, então temos que rejeitar a hipótese, e concluir que as figurinhas dos Vingadores não são independentes, e, como consequência, não são honestas, apesar de serem uniformes!

Felizmente para o colecionador, elas são desonestas mas estão jogando do seu lado. Se os pacotinhos diminuem a chance de ter repetidas, então o total esperado de pacotinhos que você precisa comprar é menor que no caso puramente honesto.

Dá para calcular o quanto você economiza com as figurinhas desonestas? Sure, vou mostrar dois métodos para fazer isso:

Cadeias de Markov

Esse primeiro método é mais complicado, mas vou começar por ele porque é mais genérico, e se aplica a mais casos. Imagine que eu já tenho 50 figurinhas no meu álbum, e eu abro mais um pacotinho. Se eu der azar, vão sair quatro figurinhas que eu já tenho. Se eu der sorte, pode sair até quatro que eu ainda não tinha. 

O método mais natural para modelar esse problema é com uma Cadeia de Markov: você considera o estado como sendo quantas figurinhas você tem no álbum, e as transições possíveis são quantas figurinhas você ainda não tinha no pacotinho.

Para montar a matriz de transição, nós precisamos da probabilidade de ir de um estado ##n## até um estado ##m##. Como cada pacotinho tem 4 figurinhas diferentes, então o ##m## precisa ser necessariamente ##n##, ##n+1##, ..., ##n+4##, e outras transições não são possíveis. 

Qual a chance de ir do estado ##n## para ##n+k##? Nós sabemos que as quatro figurinhas do pacotinho são diferentes entre si, então temos ##k## figurinhas novas e ##4-k## figurinhas que já tínhamos. A quantidade de pacotinhos que fazem essa transição é igual a ##{180-n\choose k}{n\choose 4-k}##, já que precisamos escolher ##k## figurinhas dentre que as não temos, e ##4-k## dentre as que temos. O total de pacotinhos é ##180\choose 4##, então a probabilidade total do estado ##n## para ##n+k## é ##{180-n\choose k}{n\choose 4-k}/{180\choose 4}## . Com isso podemos montar a matriz de transição ##M##:

$$ M= \begin{bmatrix} 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & \dots & 0 &  0 & 0 \\ 0 & 0 & 0 & 0 & 0.02 & 0.97 & 0 & 0 & \dots & 0 &  0 & 0 \\ 0 & 0 & 0 & 0 & \approx0 & 0.04 & 0.95 & 0 & \dots & 0 &  0 & 0 \\ 0 & 0 & 0 & 0 & \approx0 & \approx0 & 0.06 & 0.93 & \dots & 0 &  0 & 0 \\ \vdots &\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots  \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots & 0.95 &  0.04 & \approx0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots & 0 &  0.97 & 0.02 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots & 0 &  0 & 1 \\ \end{bmatrix}$$

Pela matriz podemos conferir a intuição: se você não tem nenhuma figurinha, o primeiro pacotinho vai garantidamente te dar 4 figurinhas novas, então o estado 0->4 na matriz vale 1, e todos os outros nessa linha são zero. Por outro lado, se está faltando só uma figurinha para completar, comprar um pacotinho novo dá um chance de 2% de conseguir a faltante. E, claro, se você tem todas as figurinhas, então a chance de você não conseguir mais nenhuma é 100%.

O que queremos calcular é quanto economizamos devido às figurinhas desonestas. Para isso, vamos calcular qual o valor esperado do número de pacotinhos que precisamos comprar para completar o álbum. Tendo a matriz de transição em mãos, isso é equivalente ao número médio de transições para ir do estado 0 (álbum vazio) ao estado 180 (álbum completo). E como calcular isso? Usando um algoritmo secreto!

Really, um algoritmo secreto. Não sei o motivo, mas você não acha esse algoritmo em lugar nenhum; não tem na wikipedia e nem nos livros de processos estocásticos que eu tenho aqui em casa. Então aproveitem para anotar porque é bem útil:

1. Dada a matriz de transição ##M##, calcule a matriz reduzida ##R##, retirando da matriz original todas as linhas e colunas associadas aos estados finais. (No nosso caso só a última linha e coluna são finais).
2. Calcule a matriz fundamental ##N## do processo, através da fórmula ##N=(I-R)^{-1}##.
3. O valor esperado do número de transições é a soma dos elementos da primeira linha da matriz ##N##.

Implementando o algoritmo acima com um scriptinho rápido no Mathematica, concluímos que o valor esperado do número de pacotinhos que você precisa comprar é ##257.98##. Como cada pacotinho custa R$ 1.20, então você precisa de R$309.58 em média para completar o álbum.

E quanto você gastaria com figurinhas 100% honestas? Esse cálculo eu já fiz aqui no blog anteriormente, se as figurinhas são uniformes e independentes então esse é o Coupon Collector Problem: para ##n## figurinhas a solução é ##n H(n)##, onde ##H(n)## é o número harmônico de ##n##. Essa fórmula dá o total de figurinhas, como elas são independentes, para o total de pacotinhos é só dividir por 4, chegando no total de ##259.78## pacotinhos, ou R$ 311.74, uma economia de 2 reais! (Eu falei que economizava, só não falei que era muito).

Embora o algoritmo secreto seja bacana e funcione em vários casos, ele é meio lento. Você precisa inverter uma matriz 181x181, o que não é muito fácil. O fato dela ser triangular e esparsa ajuda, mas tem um método melhor, que calcula o valor em O(n)!

Linearidade da Expectativa

Qual é o teorema mais surpreendente que você conhece? Só em probabilidade tem vários que são difíceis de acreditar à primeira vista, como o Monty Hall Problem e o Paradoxo do Aniversário. Mas um que sempre me espanta é a linearidade da expectativa (o valor esperado da soma é a soma dos valores esperados):

$$E[a+b]=E[a]+E[b]$$

Para mim a parte curiosa é que ele funciona mesmo se as variáveis ##a## e ##b## não forem independentes! Para o nosso problema isso ajuda bastante: o valor esperado para conseguir 180 figurinhas claramente é correlacionado com o valor esperado para conseguir 179 figurinhas, mas eu posso usar relações lineares entre eles sem medo.

Para essa análise eu acho mais fácil pensar invertido. Se estiver faltando 0 figurinhas para completar o album, quantos pacotinhos eu preciso? Nenhum né, ele já está completo. E se estiver faltando só uma figurinha? Aí depende do que sair no pacotinho, se sair a figurinha faltante, um pacotinho basta, se não sair, então você adiciona um pacotinho no seu valor esperado e repete o procedimento:

$$\begin{align*}\\E[1]&=p_0(1+E[0]) + p_1(1+E[1])\\ (1-p_1)E[1] &= p_1 + p_0(1+E[0]) \\ E[1] &=\frac{1}{1-p_1}\left(p_1+p_0(1+E[0])\right)\end{align*}$$

Olha só, eu sei que ##E[0]=0##; para calcular ##E[1]## eu preciso só do ##E[0]##; para o ##E[2]## eu preciso do ##E[1]## e ##E[0]##; e assim por diante. Como cada pacotinho tem no máximo 4 figurinhas que eu não tenho, então eu só preciso olhar quatro valores para trás na recorrência. Portanto, esse algoritmo é linear!

Preenchendo as probabilidades corretas, o algoritmo completo fica assim:

$$\begin{align*} \\ E[0] &= 0 \\ E[n] &= \frac{1}{1-{{180-n}\choose 4}/{180\choose 4}}\left(\frac{{180-n}\choose 4}{180\choose 4}+\sum_{1\le k \le 4}\frac{{n\choose k}{{180-n}\choose{4-k}}}{180\choose 4}E[n-k]\right) \\ \end{align*} $$

Essa implementação realmente é bem mais rápida que o anterior. A versão com matrizes no Mathematica leva 1s para rodar, enquanto que um scriptinho python com o algoritmo linear leva só 55ms!

Agora que temos os pacotinhos corretamente modelados, finalmente vamos conseguir bolar uma estratégia para minimizar o custo de completar o álbum. Mas vai ficar para o próximo post porque esse já ficou comprido demais :)

As Figurinhas dos Vingadores são Honestas? (parte 1)

O Brasil é um país muito estranho. Uma das inúmeras bizarrices brasileiras é que temos tabus temporários. Por exemplo, homem se vestir de mulher é tabu (vai acabar com a família, fim dos tempos, etc), mas no carnaval pode! Da mesma maneira, adulto comprando álbum de figurinhas é tabu (precisa crescer, largar os brinquedos, etc.), mas em época de Copa pode!

Sempre que chega época de Copa eu fico espantado com a dedicação que meus amigos tem para completar o álbum com os jogadores. E pensando um pouco no fenômeno, me veio a dúvida: quanto dinheiro você precisa para terminar um álbum desses?

Resolvi fazer um experimento para tentar calcular esse valor. Mas, ao invés de usar o álbum da Copa, escolhi completar o álbum dos Vingadores, já que é um universo que eu tenho mais afinidade. Os valores originais são os seguintes:

• O álbum em si custa R$ 6.90, e você precisa de 180 figurinhas para completar. 
• Cada pacotinho com 4 figurinhas custa R$ 1.20.
• Existe um kit com 15 pacotinhos que você pode comprar com desconto por R$ 13.90. Mas é difícil de achar: eu só encontrei na loja Geek na av. Paulista.
• Por fim, você sempre pode pedir figurinhas avulsas no site da editora (nesse caso você pode comprar até 40 figurinhas com preço individual de R$ 0.30, mais o frete de R$ 10.00).

Se você for uma pessoa sem imaginação, de cara já tem uma estratégia fácil para completar o álbum: é só pedir tudo pelo correio. Nesse caso, basta comprar o álbum e mais 180 figurinhas em 5 fretes, num total de R$110.90.

Tem como gastar menos que isso? Para descobrir, temos que analisar as figurinhas!

A Distribuição das Figurinhas

Para descobrir a distribuição das figurinhas, eu comprei ##t=106## pacotinhos (ou seja, ##n=106\times 4=424## figurinhas), e anotei o conteúdo de cada um deles (os dados originais estão no google docs). A primeira coisa que fiz foi traçar o histograma para ver o jeitão dos dados:

Pelo histograma dá para ver que faltam nove figurinhas para completar meu álbum (são os nove buracos no histograma). Por outro lado tem três figurinhas que repetiram 6 vezes, o que indica que possivelmente esse álbum não é honesto, já que algumas figurinhas são mais fáceis que outras.

Ou... talvez não? A intuição humana é notoriamente falha em intuir probabilidades, de repente essa impressão de que algumas figurinhas são mais fáceis é falsa. O jeito correto de decidir isso é fazendo um teste de hipótese. O teste em si é simples, mas eu vou fazer em câmera lenta caso alguém nunca tenha visto.

Para começar, precisamos de uma hipótese nula ##H_0## e um nível de significância ##\alpha##.

Nossa hipótese nula ##H_0## é que as figurinhas tem distribuição uniforme e são independentes. O teste vai nos falar se esse hipótese precisa ser rejeitada. O nível de significância nós podemos escolher, mas escolher o ##\alpha## é uma arte: você precisa levar em conta o número de amostras, e as taxas aceitáveis de falsos positivos e falsos negativos. Usualmente, ##\alpha##=5% funciona bem.

Qual o significado desses parâmetros? Fica simples de entender com um exemplo: suponha que todas as 424 figurinhas que eu comprei vieram repetidas. Em uma distribuição é uniforme, isso não é impossível, mas é extremamente improvável: a chance é de 1 em ##180^{424}\sim 1.7\times 10^{956}##. Esse valor é maior que o número de átomos no universo (muito maior!), então a chance de ter acontecido ao acaso é quase nula.

Se eu tirasse 424 figurinhas repetidas, eu teria bastante convicção de que as figurinhas não são uniformes. Quando fazemos um teste de hipótese, é essa ideia que estamos tentando capturar: os dados que coletamos são tão raros que não poderiam ter acontecido na prática? A chance de terem acontecido ao acaso é menor que ##\alpha##=5%? Se for, então talvez a hipótese ##H_0## seja falsa.

Existem vários testes de hipótese possíveis. Como queremos checar se nossas amostras seguem uma distribuição dada (uniforme), então o teste mais natural a se fazer é o teste do ##\chi^2## (pronuncia-se quiquadrado), que serve justamente para verificar se duas distribuições são iguais.

Você começa o teste do ##\chi^2## separando suas amostras em categorias. No nosso caso, o mais natural seria fazer cada figurinha ser uma categoria, mas isso não vai funcionar. O motivo é que o teste do ##\chi^2##, estritamente falando, só funciona com infinitas amostras. Como é impossível coletar infinitas amostras, nós assumimos que o resultado não muda muito de "infinitas amostras" para "muitas amostras".

Mas quantas amostras são "muitas amostras"? Em geral, a regra usada é que ##np>5##, onde ##n## é o número de amostras e ##p## é a menor probabilidade de uma amostra estar em uma categoria. Vamos fazer as contas: se cada figurinha é uma amostra, e as 180 figurinhas são uniformes, então ##p=1/180##. Como tenho ##n=424## figurinhas, então ##np=424/180\sim 2.35##. Oops, 2.35 é menor que 5, então 180 categorias é muito.

Vamos reduzir então. Eu vou escolher ##c=5## categorias, porque 180 divide certinho por 5. Cada categoria vai ter ##180/5=36## figurinhas, logo ##p=1/36## e ##np=424/36\sim 11.77##, agora sim! O histograma com cinco categorias é assim:

Agora nós vamos tabelar nossas amostras. A linha Observado mostra quantas amostras nós coletamos em cada categoria, a linha Esperado mostra quantas amostras uma distribuição ideal teria (no caso uniforme, todas são iguais a ##n/c=424/5=84.8##):

	Categoria 1	Categoria 2	Categoria 3	Categoria 4	Categoria 5
Observado	75	90	80	95	84
Esperado	84.8	84.8	84.8	84.8	84.8

Podemos calcular o ##\chi^2## agora. A intuição do teste é simples, vamos calcular qual o quadrado da diferença entre o observado e o esperado, e normalizar. Se esse número for muito grande, então tem algo errado. Calculando pela fórmula, temos:

$$\chi^2=\sum\frac{(Y_i-E_i)^2}{E_i}\sim 2.957$$
Esse valor, 2.957 é grande ou é pequeno? Para decidir isso precisamos primeiro calcular quantos graus de liberdade tem esse modelo. O número de graus de liberdade é o número de parâmetros que estamos estimando. Como nossa tabela tem 5 categorias, então o número de graus de liberdade é 4.

Opa! Se são 5 categorias, não tinha que ser 5 graus de liberdade? É verdade que estamos tentando estimar ##Y_1##, ##Y_2##, ##Y_3##, ##Y_4## e ##Y_5##. Mas nós sabemos que a soma de todos eles é ##n##, então se tivermos os quatro primeiros, o quinto é só calcular!

Em tempos idos esse era o momento em que você precisava sacar do bolso uma tabela de ##\chi^2## (como curiosidade, coloquei online a tabela que eu usava na Poli). Mas hoje em dia você pode ir direto no Mathematica para descobrir o valor de ##\chi^2## com ##\alpha##=5% e quatro graus de liberdade:

InverseCDF[ChiSquareDistribution[4], 1 - 0.05]
>> 9.48773

Como 2.957 é menor que 9.488, então não podemos rejeitar a hipótese de que as figurinhas são uniformes e independentes, e concluímos que sim, as figurinhas provavelmente são honestas.

Alguns de vocês irão chiar, eu sei. "Mas veja, algumas figurinhas saíram repetidas seis vezes, outras não apareceram nenhuma vez! Não tem como isso ser uniforme, claramente tem figurinhas que são mais fáceis!".

Bem, se você ainda estiver na dúvida, você pode comprar mais pacotinhos para deixar o teste mais robusto. Ao invés disso, eu vou usar só as figurinhas que eu tenho para ir direto na raiz da dúvida: afinal de contas, é normal tirar muitas figurinhas repetidas?

A Análise das Repetidas

A intuição nos diz que, em uma distribuição uniforme, todas as figurinhas tem que aparecer mais ou menos com a mesma frequência. Só que esse é mais um caso onde a intuição é falha! Quem lê o blog há mais tempo deve lembrar de um post sobre o paradoxo do aniversário, que diz exatamente isso: mesmo em distribuições uniformes, a chance de colisões é muito alta.

Podemos usar esse fato para tentar um novo teste de hipótese. Dessa vez, vamos contar quantas repetidas nós temos para cada figurinha, e separá-las por número de repetições. O histograma da quantidade de figurinhas repetidas fica assim:

O histograma foi traçado da quantidade de figurinhas repetidas:

Quantidade de repetições	0	1	2	3	4	5	6
Figurinhas	9	41	50	48	24	5	3

O álbum tem 180 figurinhas, então a soma da segunda linha vale 180. Do total de figurinhas, nove delas não apareceram, três apareceram 6 vezes, e a maioria apareceu entre 1 e 3 vezes.

O formato do histograma parece familiar né? Se nós soubéssemos qual a é a fórmula exata para o número de figurinhas repetidas em uma distribuição uniforme, poderíamos usar o teste do ##\chi^2## para ver se a hipótese ainda se mantém. Só que agora não tem jeito, para descobrir a fórmula precisamos fazer as contas. Pule a caixa azul se você tem medo de combinatória:

O truque para usar combinatória analítica numerada (labelled) é sempre tentar achar um isomorfismo entre o seu problema e uma lista de inteiros. Para o nosso caso, imagine o seguinte experimento: temos ##n## bolinhas numeradas de ##1## a ##n##, e temos também ##m## caixinhas, cada caixa representando uma figurinha. A bolinha representa a ordem em que eu comprei a figurinha. Se eu atirar as bolinhas em caixas aleatórias, tenho um problema isomórfico a comprar figurinhas aleatórias.

Por exemplo: se a bolinha 5 cair na caixa 32, então significa que a quinta figurinha que eu comprei foi a de número 32. Se no fim do experimento a caixa 32 tiver seis bolinhas, então eu tenho a figurinha 32 repetida seis vezes.

No nosso experimento, as caixinhas são diferentes umas das outras, então vamos usar o operador ##\text{SEQ}_m(z)## para modelar as caixinhas. A ordem das bolinhas dentro da caixinha não importa (tanto faz se dissermos que a caixa 32 tem as bolinhas 3, 5, 9, 10 e 15; ou se tiver as bolinhas 9, 5, 3, 15, e 10, é a mesma coisa). Então dentro das caixinhas, usaremos o operador ##\text{SET}(z)##. O experimento completo é ##\text{SEQ}_m(\text{SET}(z))##.

Digamos que estamos interessados em saber quantas caixinhas vão terminar com ##s## bolinhas. Um jeito de medir isso é colocar um marcador nas caixinhas com ##s## elementos. Para isso, nós tiramos as das caixinhas os conjuntos de ##s## elementos, e colocamos eles de volta com o marcador, digamos que o marcador seja ##u##. Então o experimento para medir ##s## repetidas é ##\text{SEQ}_m(\text{SET}(z)-\text{SET}_s(z)+u \text{SET}_s(z))##.

Os experimentos vão mapear para funções geradoras exponenciais de duas variáveis. Das tabelas de combinatória analítica:
$$\begin{align*}
\text{SEQ}_m(z) &= z^m \\
\text{SET}(z) &= e^z \\
\text{SET}_s(z) &= \frac{z^s}{s!}\\
\text{SEQ}_m(\text{SET}(z)-\text{SET}_s(z)+u \text{SET}_s(z)) &=
\left(e^z-\frac{z^s}{s!}+u \frac{z^s}{s!}\right)^m \\
&= \left(e^z+ (u-1) \frac{z^s}{s!}\right)^m
\end{align*}$$
Se abrirmos essa função geradora em série, então o coeficiente do termo ##z^n u^k## indica quantas combinações possíveis de ##n## bolinhas vão ter ##k## caixinhas com ##s## bolinhas repetidas (como a função geradora é exponencial, você precisa multiplicar o coeficiente por ##n!## para ter o valor real). Se chamarmos esse coeficiente de ##R_s(n, k)##, então temos:
$$\sum_{n}\sum_{k}R_s(n,k) z^n u^k = \left(e^z+ (u-1) \frac{z^s}{s!}\right)^m $$
Vamos listar com cuidado o que temos. ##n! R_s(n,k)## é número de combinações possíveis com ##k## caixinhas contendo exatamente ##s## bolinhas. Logo, de todas as combinações possíveis de ##n## bolinhas, existem ##\sum k n! R_s(n,k)## caixinhas com ##s## bolinhas.

Por outro lado, o número total de combinações de ##n## bolinhas em ##m## caixinhas é ##m^n##. Como cada combinação dessas tem ##m## caixinhas, então o número total de caixinhas que podem ser preenchidas em todas as combinações é ##m^{n+1}##.

Juntando as duas contagens, a probabilidade de um caixinha tem ##s## bolinhas em ##n## arremessos é:
$$p_s(n)=\sum_{k>0}\frac{k \;n!}{m^{n+1}}R_s(n, k)$$
Mas olha só! Nós podemos conseguir essa somatória derivando a função geradora em relação a ##u##, e depois setando ##u=1##.
$$\left.\frac{\partial}{\partial u}\sum_{n}\sum_{k}R_s(n,k) z^n u^k \right|_{u=1}=\sum_n \sum_{k}k R_s(n,k) z^n$$
Agora é só fazer as contas:
$$\begin{align*}p_s(n)
&=\sum_k \frac{n!}{m^{n+1}}k R_s(n, k) \\
&=\frac{n!}{m^{n+1}} [z^n] \left.\frac{\partial}{\partial u}\sum_{n}\sum_{k}R_s(n,k) z^n u^k \right|_{u=1}\\
&=\frac{n!}{m^{n+1}} [z^n] \left.\frac{\partial}{\partial u}\left(e^z+ (u-1) \frac{z^s}{s!}\right)^m  \right|_{u=1}\\
&=\frac{n!}{m^{n+1}} [z^n] \left. m\left(e^z+ (u-1) \frac{z^s}{s!}\right)^{m-1} \frac{z^s}{s!} \right|_{u=1}\\
&=\frac{n!}{m^{n+1}} [z^n]  m\left(e^z\right)^{m-1} \frac{z^s}{s!} \\
&=\frac{n!}{s!\;m^n} [z^n]z^s e^{z(m-1)}\\
\end{align*}$$

Abrindo a exponencial em série:

$$\begin{align*}p_s(n)
&=\frac{n!}{s!\;m^{n}} [z^n]z^s\sum_{k\ge 0}\frac{(z(m-1))^k}{k!} \\ &=\frac{n!}{s!\;m^{n}} [z^n]\sum_{k\ge 0}\frac{z^{k+s}(m-1)^k}{k!} \\ &=\frac{n!}{s!\;m^{n}} [z^n]\sum_{p\ge s}\frac{z^p(m-1)^{p-s}}{(p-s)!} \\ &=\frac{n!}{s!\;m^{n}}\; \frac{(m-1)^{n-s}}{(n-s)!} \\ &=\frac{n!}{s!(n-s)!}\; \frac{(m-1)^{n-s}}{m^{n}} \\ &={{n}\choose{s}}\frac{(m-1)^{n-s}}{m^{n}} \\ \end{align*}$$
Essa é a solução exata, em forma fechada, que queríamos! Mas tem uma massagem que podemos fazer que vai ajudar a entender por que essa fórmula é assim. Em geral ##n## é bem maior que ##m##, então vamos reescrever ##m## de modo que ##n=\lambda m##, e calcular o limite quando ##n## é muito grande: $$\begin{align*}\lim_{n\to\infty}p_s(n) &=\lim_{n\to\infty}\frac{n!}{s!(n-s)!}\; \frac{(n/\lambda-1)^{n-s}}{(n/\lambda)^{n}} \\ &=\lim_{n\to\infty}\frac{n!}{s!(n-s)!}\; \left(\frac{n-\lambda}{\lambda}\right)^{n-s} \left(\frac{\lambda}{n}\right)^{n}\\ &=\lim_{n\to\infty}\frac{n!}{s!(n-s)!}\; \lambda^{s} \left(\frac{n-\lambda}{n}\right)^{n-s}\left(\frac{1}{n}\right)^{s} \\ &=\lim_{n\to\infty}\frac{n!\;\lambda^{s}}{s!(n-s)!\;n^{s}}\; \left(1-\frac{\lambda}{n}\right)^{n-s}\\ &=\lim_{n\to\infty}\left(\frac{\lambda^{s}}{s!}\right) \left(\frac{n!}{(n-s)!\;n^s}\right) \left(1-\frac{\lambda}{n}\right)^{n} \left(1-\frac{\lambda}{n}\right)^{-s} \\ \end{align*}$$ O limite do produto é o produto dos limites, desde que todos existam e sejam finitos, então vamos conferir um por um. O primeiro termo é constante, easy. No segundo, o numerador é um polinômio de grau ##n##, o denominador é um polinômio de grau ##n##, então a fração é constante e vale 1, já que os coeficientes iniciais ambos valem 1. O terceiro é a definição de ##e^{-\lambda}##. No quarto a segunda fração vai para zero, então o total tende a 1. Combinando tudo: $$\lim_{n\to\infty}p_s(n)= e^{-\lambda}\frac{\lambda^s}{s!}$$ Ahá, é uma distribuição de Poisson!

Não foi à toa que aquele histograma das figurinhas repetidas parecia familiar, elas seguem uma distribuição de Poisson! Em retrospecto até faz sentido. Se a distribuição das figurinhas é uniforme, então em um longo período de tempo, cada figurinha deve ter uma taxa de aparição constante (igual a ##\lambda=n/m##). Se a taxa é constante ao longo do tempo, então a quantidade de figurinhas em um pequeno intervalo tem que seguir a distribuição de Poisson.

Agora falta pouco para conseguir usar o teste do ##\chi^2## nas repetidas. Precisamos primeiro separar as amostras em categorias. O natural seria usar uma categoria por quantidade de repetidas, mas em teoria, um cara muito azarado poderia ter 424 figurinhas repetidas, e nesse caso o produto ##np## seria muito, muito menor que 5. Vamos então aglutinar todas as categorias de 6 repetidas em diante em um grupo só. A tabela fica assim:

Quantidade de repetições	0	1	2	3	4	5	6 ou mais
Figurinhas, medido	9	41	50	48	24	5	3
Figurinhas, teórico	17.1	40.2	47.4	37.2	21.9	10.3	5.9

Agora podemos calcular o ##\chi^2## medido e o teórico, lembrando que temos 6 graus de liberdade e ##\alpha## valendo 5%: $$\chi^2=\sum\frac{(Y_i-E_i)^2}{E_i}\sim 11.5275$$

InverseCDF[ChiSquareDistribution[6], 1 - 0.05]
>> 12.59

Conclusão: como 11.53 é menor que 12.59, então novamente não podemos rejeitar a hipótese de que as figurinhas são uniformes e independentes, e concluímos que sim, as figurinhas provavelmente são honestas. É verdade que tem muitas repetidas, mas o número de repetidas está dentro do esperado.

Quanto gastar para completar?

No pŕoximo post, vamos finalmente calcular qual o minimo de dinheiro que você precisa gastar para completar o álbum. Stay tuned!

Matemática x Tanques de Guerra

Quando eu era criança pequena lá na década de 80, eu tinha um brinquedo de tanque de guerra. Ele chamava Panzer, e era fabricado pela brinquedos Mimo. Para a época era bem divertido: tinha controle remoto, e, além de andar, ainda disparava mísseis!

Mas eu não percebia a ironia do brinquedo. O Panzer não era só um brinquedo, ele existiu de verdade. Na vida real, o Panzer era uma máquina nazista de destruição. A série Panzer de tanques alemães foi responsável pela morte direta de incontáveis aliados. Ver um Panzer surgindo no horizonte era como ver a morte vindo em sua direção. Era um mimo esse Panzer.

Por outro lado, o Panzer era forte, mas não era invencível. Os aliados tinham armas capazes de deter o Panzer, o problema era descobrir quantas dessas armas eles deveriam levar para o front. Poucas semanas antes do Dia D, os alemães começaram a usar o novo modelo Panzer V, e os aliados não tinham como estimar quantas armas levar sem saber quantos Panzer V os inimigos fabricaram.

Foi então que os aliados usaram a arma definitiva contra o Panzer: a Matemática! Usando um truque muito esperto, eles estimaram que os alemães estavam fabricando 270 Panzers por mês. Depois que a guerra acabou, eles foram nos arquivos alemães conferir os números, e na verdade estavam fazendo 276 Panzers por mês. A estimativa foi muito boa!

A Quantidade de Tanques

O truque dos aliados se baseava em um erro dos alemães e em uma hipótese estatística. O erro dos alemães era que os tanques tinham número de série sequencial. Quando saíam da fabrica, eram numerados como 1, 2, 3, e assim por diante. Ao fazer isso, eles deixaram uma brecha para os hackers aliados.

Suponha que você capturou um Panzer alemão, e o serial dele era 14. Você consegue estimar quantos tanques eles produziram baseado só nesse número? No caso geral não, mas você pode fazer uma hipótese simplificadora, que é supor que o tanque que você capturou é aleatório (ou seja, os alemães estavam distribuindo igualmente os tanques por toda a Europa).

Com essa hipótese a intuição é clara. Vamos supor, hipoteticamente, que os alemães fizeram 200 tanques. Qual chance de você capturar um tanque de serial menor ou igual a 14 nessa situação? É 14/200 né? E se os alemães fizeram só 20 tanques? Aí a chance é 14/20, dez vezes maior que 14/200. Com a hipótese de tanque aleatório, é mais provável que os alemães tenham feito 20 tanques do que 200 tanques.

Você pode formalizar esse argumento usando o teorema de Bayes, e calcular o valor esperado do número de tanques. Vamos fazer as contas na caixa azul:

Suponha que conseguimos capturar k tanques. Dentre esses k tanques, o maior serial encontrado tem o valor m. Com essas informações, vamos calcular qual é o valor esperado do número de tanques produzidos, que vamos chamar de n. Podemos abrir esse valor esperado pela definição:
$$E[n\;|\;m,k] = \sum_{n} n \; p(n\; |\; m,k)$$
O problema então é achar a probabilidade condicional de n, dados m e k. Nós ainda não sabemos quanto vale essa probabilidade, mas podemos abri-la usando o teorema de Bayes:
$$p(n\;|\;m,k)=\frac{p(m,k\;|\;n)\;p(n)}{p(m,k)}$$
Essas três probabilidades nós conseguimos calcular! A primeira delas, p(m,k|n), sai por combinatória. Dado um vetor de tamanho n, de quantas maneiras conseguimos escolher k elementos, de modo que o maior deles seja m?

Um dos elementos precisa necessariamente ser o m, então só precisamos calcular a posição dos outros k-1 restantes. E eles precisam ser menores que m, então só podem estar nas m-1 posições menores que m. Por isso, essa quantidade vale binomial(m-1, k-1). Note, entretanto, que isso só vale quando n>=m, caso contrário a probabilidade é zero (não tem como o máximo ser m se o vetor tem um tamanho menor que m).

Para calcular a probabilidade desejada, temos que achar de quantas maneiras podemos escolher k elementos dentre n, sem considerar a restrição do maior deles ser m. Mas isso é o próprio binomial(n, k). Então, a probabilidade é:
$$p(m,k\;|\;n) = \frac{\displaystyle{m-1\choose k-1}} {\displaystyle{n\choose k}}[n\ge m]$$
Vamos agora para a segunda, p(n). Eu não tenho nenhuma informação sobre quantos tanques foram produzidos, então não sei qual é a probabilidade de ter n tanques. Mas eu posso chutar que a distribuição é uniforme. A chance de ter 15 tanques, ou de ter 500 tanques, a priori, deve ser a mesma. Por isso, vamos fazer p(n) ser uma constante independente de n:
$$p(n)=c$$
A última, p(m,k), é a probabilidade do máximo ser m em k escolhas, independente do valor de n. Para isso, podemos usar a lei da probabilidade total: 
$$p(m,k) = \sum_n p(m,k\;|\;n) \;p(n)$$ 
Essas duas já sabemos calcular! É só substituir:
 $$\begin{align*} p(m,k) &= \sum_n p(m,k\;|\;n) \;p(n) \\ &= \sum_n \frac{{m-1\choose k-1}} {{n\choose k}}[n\ge m]\;c\\ &= \sum_{n\ge m} c {m-1\choose k-1} {n\choose k}^{-1}\\ &= c {m-1 \choose k-1} \sum_{n\ge m} {n\choose k}^{-1}\\ \end{align*} $$
A somatória tem uma forma fechada, que você pode achar com o algoritmo de Gosper:
$$\sum_{n=m}^{\infty}{n\choose k}^{-1}=\frac{m}{k-1}{m\choose k}^{-1}$$ 
Substituindo:
$$\begin{align*} p(m,k) &= c {m-1 \choose k-1} \sum_{n\ge m} {n\choose k}^{-1}\\ &= c {m-1 \choose k-1} \times \frac{m}{k-1}{m\choose k}^{-1} \\ &= \frac{cm}{k-1} {m-1 \choose k-1} {m\choose k}^{-1} \\ \end{align*} $$
Agora podemos voltar e substituir as três probalidades na fórmula de Bayes:
$$ \begin{align*} p(n\;|\;m,k)&=\frac{p(m,k\;|\;n)\;p(n)}{p(m,k)} \\ &= \frac{\displaystyle{m-1\choose k-1}{n\choose k}^{-1}[n\ge m]\times c} {\displaystyle\frac{cm}{k-1}{m-1\choose k-1}{m\choose k}^{-1}}\\ &= \frac{k-1}{m}{m\choose k}{n\choose k}^{-1}[n\ge m] \end{align*}$$
Por fim, podemos substituir na fórmula do valor esperado:
 $$\begin{align*} E[n\;|\;m,k] &= \sum_{n} n \; p(n\; |\; m,k)\\ &= \sum_n n \frac{k-1}{m}{m\choose k}{n\choose k}^{-1}[n\ge m]\\ &= \frac{k-1}{m}{m\choose k}\sum_{n\ge m} n {n\choose k}^{-1}\\ \end{align*}$$
Essa somatória também tem uma forma fechada por Gosper:
$$ \sum_{n=m}^{\infty} n{n\choose k}^{-1}=\frac{m(m-1)}{k-2}{m\choose k}^{-1}$$
Chegamos então na substituição final:
$$\begin{align*} E[n\;|\;m,k] &= \frac{k-1}{m}{m\choose k}\sum_{n\ge m} n {n\choose k}^{-1}\\ &= \frac{k-1}{m}{m\choose k}\times\frac{m(m-1)}{k-2}{m\choose k}^{-1} \\ &= \frac{(k-1)(m-1)}{k-2} \\ \end{align*}$$

A Estimativa na Prática

A fórmula final é super simples, mas ela funciona mesmo na prática? Podemos testá-la fazendo uma simulação numérica. Para isso, eu fiz dez mil simulações. Em cada uma eu sorteio o número de tanques fabricados, e escolho aleatoriamente 3 deles. Aí eu calculo a estimativa pela fórmula, e normalizo o erro. O código está no github, e o histograma resultante é o abaixo. A normalização é tal que o valor 0 é uma estimativa totalmente correta:

A estimativa não ficou muito boa não! Ao invés de ter uma gaussiana em torno do zero, o estimador tem um bias muito forte para os negativos. Felizmente, esse problema só acontece porque capturamos apenas 3 tanques. Se, ao invés disso, tivéssemos capturado 30 tanques, então o histograma seria bem melhor. Compare o histograma de 30 tanques sobreposto ao de 3 tanques:

Bem melhor né? Agora quase todas as estimativas estão bem próximas do zero. Esse método é bastante sensível ao número de tanques capturados. Os aliados sabiam disso, por isso que a estimativa deles foi tão boa (estimaram 270 tanques/mês, quando o valor real era 276). Sabe quantos tanques eles capturaram para fazer essa estimativa? Só dois tanques!

Claro que tem um truque. A estimativa com dois tanques é muito ruim, mas eles notaram que os alemães colocavam serial em todas as peças do tanque. Em especial, cada tanque tinha 48 rodas, cada uma com um serial único. Por isso, eles conseguiram usar a fórmula com k=96, o que deu a precisão alta que queriam. E no fim a matemática ganhou a guerra :)