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!

A Marcha dos Lemmings

Lemmings são pequenos roedores que vivem na tundra ártica. Eles possuem pelos longos e macios, rabo curto, são herbívoros e explodem.

Sim, eles explodem! Os lemmings possuem a capacidade de ligar um contador luminoso sobre suas cabeças. A cada passo que dão, o contador decrementa, até que finalmente explodem quando o contador chega a zero.

O que você faria se soubesse que vai explodir em seis passos? Parece que não tem muitos passeios possíveis com exatamente seis passos, mas a intuição é meio falha aqui. Na verdade, o número de passeios possíveis é exponencial no número de passos! Como fazer um algoritmo que calcula o número de caminhos possíveis que o lemming pode fazer antes de explodir?

As Regras do Jogo

Para modelar o problema, vamos supor que os passos que lemming pode fazer são descritos por um grafo orientado. Por exemplo, no grafo abaixo o lemming começa no nó A. Do A ele pode ir para B ou C, de B ele pode ir para C, e de C ele pode ir para A.

Se o nosso lemming pode dar seis passos, então existem sete caminhos diferentes que ele pode fazer:

ABCABCA

ABCACAB

ABCACAC

ACABCAB

ACABCAC

ACACABC

ACACACA

Logo ##p(6)=7##, e podemos até tabelar ##p(n)## para todos os ##n## de 0 a 6:

$$p(0..6)=1,2,2,3,4,5,7$$
Dado ##n##, como fazer um algoritmo que calcule ##p(n)##?

Eu vou mostrar três algoritmos possíveis usando a linguagem do Wolfram Mathematica. Para isso, o grafo precisa ser descrito de alguma maneira que os algoritmos entendam. Eu vou usar uma matriz de adjacências: para um grafo com ##k## nós, a matriz é ##k\times k## e o elemento na linha ##i## da coluna ##j## é 1 se existe uma ligação saindo do nó ##j## e indo em direção ao nó ##i##, e 0 em caso contrário. Para o nosso exemplo, a matriz é:

$$M=\left[\begin{matrix}0 && 0 && 1 \\ 1 && 0 && 0 \\ 1 && 1 && 0\end{matrix}\right]$$

Força Bruta

Sempre que eu tenho que fazer um algoritmo do zero, prefiro começar com a força bruta. É fácil e rápido de implementar, e embora não seja eficiente, serve para conferir o resultado dos algoritmos mais avançados.

Nesse caso, para implementar a força bruta você passa o grafo, o nó em que o lemming está, quantos passos faltam para terminar, e o resto é uma recursão simples. Como toda recursão, não esqueça de fazer o caso base! Para esse problema, o caso base é ##p(0)=1##, ou seja, tem um único caminho possível sem dar passo nenhum, que é não sair do lugar.

No fim, a solução em força bruta fica assim:

count[graph_, pos_, 0] := 1
count[graph_, pos_, size_] := Sum[
    If[graph[[i, pos]] == 1, count[graph, i, size - 1], 0],
    {i, 1, Length[graph]}]

Vamos conferir o resultado:

graph := {{0, 0, 1}, {1, 0, 0}, {1, 1, 0}}
Print[Table[count[graph, 1, n], {n, 0, 6}]]

{1,2,2,3,4,5,7}

Perfeito! Mas qual é a complexidade desse algoritmo? Bem, ele precisa visitar cada um dos caminhos possíveis, e o número total de caminhos é exponencial, logo a complexidade vai ser do tipo ##O(c^n)##, onde a constante ##c## depende do grafo. Ou seja, você não quer usar esse algoritmo quando o ##n## é grande.

Multiplicação de Matrizes

Como melhorar então? Tem uma maneira calcular ##p(n)## mais rápido usando multiplicação de matrizes. É mais fácil mostrar como funciona por exemplos: no estado inicial, nós temos um lemming no nó A, então vamos montar um vetor correspondente: ##V=[1 \;0\; 0]^T##. Se multiplicarmos ##M## por ##V##, qual o resultado?

$$M V=\left[\begin{matrix}0 && 0 && 1 \\ 1 && 0 && 0 \\ 1 && 1 && 0\end{matrix}\right] \left[\begin{matrix}1\\0\\0\end{matrix}\right]=\left[\begin{matrix}0\\1\\1\end{matrix}\right]$$

Ou seja, um lemming que está em A vai parar em B ou C. Continuando:

$$M V=\left[\begin{matrix}0 && 0 && 1 \\ 1 && 0 && 0 \\ 1 && 1 && 0\end{matrix}\right] \left[\begin{matrix}0\\1\\1\end{matrix}\right]=\left[\begin{matrix}1\\0\\1\end{matrix}\right]$$
Dado um lemming em B ou C, o resultado é um lemming em A ou C. Se o vetor de entrada mostra quantos lemmings tem em cada nó, multiplicar pela matriz M vai mostrar onde eles podem estar depois de dar um passo.

Agora é só estender o raciocínio. Se multiplicar por M correponde a um passo, então multiplicar ##n## vezes por M é o equivalente a ##n## passos. E se você começou o processo com um vetor unitário, então a soma do elementos de V é o número de caminhos possíveis partindo do nó inicial!

Uma implementação possível é a abaixo:

count[graph_, pos_, size_] :=
  Total[MatrixPower[graph, size] . UnitVector[Length[graph], pos]]

Conferindo:

graph := {{0, 0, 1}, {1, 0, 0}, {1, 1, 0}}
Print[Table[count[graph, 1, n], {n, 0, 6}]]

{1,2,2,3,4,5,7}

Bateu certinho novamente! E qual é a complexidade? Se você implementou a multiplicação e a exponenciação de matrizes do jeito mais simples, vai ser da ordem de ##O(k^3 n)##, que certamente é melhor que exponencial mas ainda não é boa o suficiente para usar na prática.

No entanto, você pode melhorar isso com implementações melhores! A multiplicação de matrizes pode usar o algoritmo de Strassen, e a exponenciação pode usar o método binário, aí o total fica ##O(k^{2.8} \log_2 n)##, bem melhor; e se a matriz for esparsa dá para ficar ainda mais rápido.

Função Geradora

Mas tem um método que supera o anterior, que é usar funções geradoras! Nós começamos tudo com um grafo orientado, e todo grafo orientado é equivalente a um autômato finito. Por sua vez, todo autômato finito é equivalente a uma gramática regular, então deve ter um jeito de descrever os caminhos do lemming usando uma gramática. Para o exemplo dado, a gramática fica assim:

$$\begin{align*}\alpha &\to A \;|\; A \beta \;|\; A \gamma \\ \beta &\to B \;|\;  B \gamma \\ \gamma &\to C \;|\; C \alpha\end{align*}$$
Todos os caminhos possíveis produzem strings que são aceitas por essa gramática, e a gramática não aceita nenhuma string que não seja um caminho possível.

O truque agora é usar combinatória analítica, que para o problema de enumerar strings de uma gramática é bem simples: basta trocar cada token por um ##z##, e cada OR por uma adição. A gramática vira um sistema de equações:

$$\begin{align*}\alpha &= z + z\beta +z \gamma \\ \beta &= z+z \gamma \\ \gamma &= z+z \alpha\end{align*}$$
Resolvendo para ##\alpha##, chegamos na função geradora:

$$\alpha=\frac{z + 2 z^2 + z^3}{1 - z^2 - z^3}$$
E para que serve a função geradora? Olhe só o que acontece quando abrimos a função em série de potências:

$$\frac{z + 2 z^2 + z^3}{1 - z^2 - z^3}=z+2z^2+2z^3+3z^4+4z^5+5z^6+7z^7+\dots$$

Ahá! Os coeficientes de ##z## são exatamente o número de caminhos no grafo! Para calcular ##p(6)##, basta olhar o coeficiente de ##z^7## (você precisa somar um porque a função geradora está contando o número de tokens na strings, e nós queremos o número de passos).

Podemos ir além agora. Quando a função geradora é uma função racional, você sempre consegue inverter e conseguir uma fórmula explícita para ##p(n)##:

$$p(n)=0.957\times 1.325^n+0.426\times 0.869^n\cos(1.469+ 2.438 n)$$
Essa fórmula é exata (ou seria, se eu não tivesse truncado os números para caber na tela). Note que o resultado sempre vai ser um inteiro, e que a parte oscilante da fórmula sempre vai ser menor que 0.4, então dá para simplificar mais ainda:

$$p(n)=\left\lfloor 0.5+0.957\times 1.325^n\right\rfloor$$
Basta pegar a parte exponencial e arredondar para o inteiro mais próximo. Podemos observar duas coisas aqui. Primeiro, agora nós temos um algoritmo onde o número de operações independe do valor de ##n##, então efetivamente esse algoritmo é ##O(1)##. Na prática você vai implementar usando bignum, aí não fica ##O(1)## de verdade; mas os outros algoritmos mostrados também vão ficar mais lentos na mesma proporção.

Segundo, nós provamos que aquele primeiro algoritmo era de fato ##O(c^n)##! Conseguimos até calcular quem é o tal do ##c##: é aproximadamente 1.325; ou, se você quiser o valor exato com radicais:

$$c=\frac{\sqrt[3]{9+\sqrt{69}}+\sqrt[3]{9-\sqrt{69}}}{\sqrt[3]{18}}$$
Para implementar esse método computacionalmente, é só reescrever aquele sistema de equações em forma matricial:

$$\left[\begin{matrix}\alpha\\ \beta \\ \gamma\end{matrix}\right] = z \left[\begin{matrix}0 && 1 && 1\\0 && 0 && 1\\ 1 &&0 && 0\end{matrix}\right] \left[\begin{matrix}\alpha\\ \beta \\ \gamma\end{matrix}\right]  + z\left[\begin{matrix}1\\ 1\\1\end{matrix}\right] $$

A matriz ali no meio é a transposta da matriz de adjacências! Agora já podemos implementar:

count[graph_, pos_, size_] := With[{n = Length[graph]},
  SeriesCoefficient[LinearSolve[
      IdentityMatrix[n] - z Transpose[graph],
      ConstantArray[z, n]][[pos]],
    {z, 0, size + 1}]]

Como esperado, a resposta é correta:

graph := {{0, 0, 1}, {1, 0, 0}, {1, 1, 0}}
Print[Table[count[graph, 1, n], {n, 0, 6}]]

{1,2,2,3,4,5,7}

Um Problema em Aberto

Nós chegamos em um algoritmo que é ##O(1)##, dá para otimizar mais? Não! Esse é um limite teórico.

Quer dizer, existem algoritmos que são ##O(0)##. Por exemplo, "crie um algoritmo que, dado uma circunferência e um diâmetro, calcule a razão entre os dois". Esse é algoritmo é ##O(0)## porque a resposta é sempre ##\pi##, independe do input. Mas, se a resposta não é constante, ##O(1)## é o melhor algoritmo possível. (A não ser que você crie um algoritmo que dê a resposta antes do input chegar. Nesse caso me ensine como faz, porque você inventou a viagem no tempo).

Mas existe uma pequena variação no problema dos lemmings que o torna muito, muito mais difícil. Suponha que agora o lemming não pode passar duas vezes pelo mesmo nó. A força bruta para essa variação ainda é ##O(c^n)##. Mas ninguém bolou um algoritmo que seja melhor que a força bruta nesse caso!

Quando o grafo é finito, você pode usar a força bruta para enumerar todas as soluções, já que eventualmente o grafo acaba. Mas se o grafo for infinito, aí danou-se. Para esse caso, tem gente procurando a solução desde a década de 50. Já são setenta anos procurando o algoritmo e ninguém achou nada para um grafo genérico, nem mesmo um assintótico! O melhor resultado até hoje é um de 2011 onde acharam o ##c## para um grafo específico, o lattice hexagonal. Para todos os outros grafos, o melhor que temos são heurísticas e aproximações numéricas.

Se você estiver sem nada o que fazer no próximo fim de semana, essa pode ser uma boa maneira de ocupar o tempo!

O Objetivo Final da Matemática

Quem acompanhou as contas na caixa azul do post sobre o Pokèmon Twitch deve ter notado que eu pulei um pedaço. Em certo ponto da demonstração, eu magicamente tirei do chapéu uma fórmula nada óbvia:

E como foi que eu cheguei nessa fórmula? Simples, eu calculei no Wolfram Alpha:

"Ricbit, você trapaceou! Que roubalheira! Um matemático de verdade jamais faria isso!"

Mas será que é trapaça mesmo? Nesse post, eu vou mostrar como se resolve essa conta sem usar o Wolfram Alpha; e, se você acha que é trapaça, talvez acabe mudando de idéia :)

A progressão hipergeométrica

Para resolver essa somatória sem usar o Wolfram Alpha, nós precisamos relembrar as progressões que aprendemos no colégio. A primeira delas foi a progressão aritmética, onde a diferença entre dois números consecutivos é constante:

   

A segunda delas é a progressão geométrica, onde a razão entre dois números consecutivos é constante:

   

Uma progressão extremamente útil que não ensinam no colégio é a progressão hipergeométrica. Ela parece com a progressão geométrica, porque é definida a partir da razão entre dois elementos consecutivos. Mas agora, ao invés dessa razão ser constante, nós permitimos que a razão seja uma função racional (que é só um jeito chique de dizer que ela é o quociente entre dois polinômios):

 

Naquela somatória que queremos resolver, os termos sendo somados formam uma progressão hipergeométrica. Confira:

   

O numerador é um polinômio em q, o denominador é um polinômio em q. Então confere com a nossa definição de progressão hipergeométrica.

A soma da progressão hipergeométrica

Certamente vocês lembram que as progressões do colégio tinham uma fórmula para a soma de seus termos. A soma dos n primeiros termos da progressão aritmética é:

A soma dos n primeiros termos da progressão geométrica é:

 

A soma dos n primeiros termos da progressão hipergeométrica é... oops! Não existe uma fórmula geral para a soma de termos hipergeométricos.

Antigamente, achar uma fórmula para uma soma hipergeométrica era tão raro que elas até ganhavam nomes especiais. Por exemplo, a igualdade abaixo tem um nome, é o Binômio de Newton:

 

Outras somas hipergeométricas não tem fórmula, mas são tão úteis, que a própria soma ganha um nome e passa a valer como se fosse uma fórmula:

Por muitos séculos esse foi o estado da arte. Alguém acha uma fórmula aqui, alguém nomeia uma soma ali. Mas isso mudou no meio do século 20.

O algoritmo de Gosper

O primeiro progresso na direção de um método geral de resolução de somas hipergeométricas foi feito pela Mary Celine Fasenmyer em 1945. Depois de se formar em matemática, ela se juntou à ordem religiosa Sisters of Mercy (sim, a mesma que inspirou o nome da banda), e lá ela criou o Método da Irmã Celine, um algoritmo que consegue resolver automaticamente alguns tipos de somas hipergeométicas.

Em 1978 o método foi aprimorado por Bill Gosper e se tornou o algoritmo de Gosper (muito embora eu desconfie que, se a Irmã Celine fosse homem, o algoritmo seria de Gosper-Celine). O método de Gosper se baseia em um teorema surpreendente: A solução de uma soma hipergeométrica é outro hipergeométrico, se e somente se ela puder ser escrita na forma abaixo, onde R[j] é uma função racional.

Essa somatória parece esquisita. Se j é o índice da somatória, então o que ele está fazendo no lado direito da equação? E cadê os limites da somatória?

A resposta é que essa é uma somatória indefinida. Funciona do mesmo jeito que no cálculo: lá você tem integrais definidas e integrais indefinidas, aqui nós temos somatórias definidas (com limites), e somatórias indefinidas (sem limites). Para uma somatória indefinida, sempre vale a propriedade abaixo:

 

Essas duas equações são tudo que nós precisamos para usar o algoritmo de Gosper. Vamos à caixa azul calcular a nossa somatória original:

Começamos substituindo uma equação na outra:

Podemos então dividir tudo por h[j]:

 

Opa, a razão entre os h nós já calculamos lá em cima para a nossa somatória!

 

Falta achar R[j]. Eu sei que ele é um quociente de polinômios, mas não sei qual o grau desses polinômios. Eu vou chutar que o grau é 1, se der errado depois eu volto e aumento o grau. Pois bem, se o grau é 1, então R[j] é da forma abaixo:

 

Agora eu vou substituir. Aperte o cinto que lá vem turbulência:

   

O resultado é um polinômio gigante em j que é identicamente nulo. Se ele é identicamente nulo, então cada um dos seus coeficientes precisa ser zero. Com isso podemos montar um sistema de equações:

 

Agora é "só" resolver o sistema! Vou poupá-los de páginas e páginas, o resultado é o seguinte:

   

Se o sistema não tivesse solução, eu teria que voltar lá em cima e aumentar o grau do polinômio (se as contas foram titânicas com grau 1, imagine com grau 2). Mas o processo não é infinito, existem métodos que permitem achar um limite superior para o grau do polinômio. Se você atingir esse limite sem achar nenhuma solução, então você provou que sua somatória não pode ser escrita como um hipergeométrico.

Voltando à nossa equação original, agora nós sabemos quanto vale a somatória indefinida:

   

Nós queremos a soma de 1 até m, então precisamos converter a soma indefinida em uma soma definida, usando o teorema fundamental do cálculo discreto:

 

Aplicando o teorema, temos:

   

Pronto, esse é o resultado final. Eu acho muito importante que a pessoa resolva o algoritmo de Gosper uma vez na vida, para nunca mais cometer a sandice de repetir o processo. O método é mecânico, você só precisa seguir as regras e não errar as contas. O computador faz isso melhor que você. Use o Wolfram Alpha e economize grafite.

O objetivo final da Matemática

No final do século 19, o Bertrand Russell e o Alfred Whitehead tiveram uma idéia. Ele bolaram um plano para axiomatizar toda a matemática, de uma maneira completa e livre de contradições. Se o plano funcionasse, então toda a Matemática seria reduzida à manipulação de símbolos, que é um processo puramente mecânico e automatizável. O plano começou muito bem, a ponto do Whitehead fazer a seguinte declaração:

"A Civilização avança estendendo o número de operações importantes que podemos fazer sem ter que pensar."

Baseado nessa citação, o Concrete Mathematics do Knuth vai além e conclui:

"O objetivo final da matemática é eliminar a necessidade de todo pensamento inteligente."

É verdade, o objetivo sempre foi esse. Para os babilônios, resolver uma equação de segundo grau era dificílimo, só os grandes mestres conseguiam. Hoje em dia, nós transformamos a equação de segundo grau em uma fórmula que é ensinada para crianças de 14 anos. Você não precisa ser inteligente para usar a fórmula, é só substituir os números e fazer as contas.

Imagine se fosse possível fazer isso com toda a matemática, ao invés de só com equações de segundo grau! Você poderia usar computadores para fazer as contas, no lugar de crianças de 14 anos. A Matemática seria um problema resolvido, e aí você pode se concentrar no que interessa, que são as aplicações. Um físico acabou de descobrir uma propriedade quântica nova, será que ela permite criar um aparelho de teletransporte? Eu monto a equação, jogo no computador, e tcharam! Está aqui a resposta!

Toda a saga por trás do plano do Russell e do Whitehead é contada em quadrinhos no fabuloso Logicomix. Infelizmente, o final da história é que o plano falhou. E falhou espetacularmente. Em 1931, o Gödel provou que nenhum sistema axiomático é capaz de provar todas as sentenças verdadeiras que é capaz de descrever, acabando de uma vez por todas com o sonhado plano da axiomatização completa.

E o século 20 desceu ladeira abaixo. O Church provou que, mesmo se a axiomatização completa tivesse tido sucesso, seria impossível criar um método que prove um teorema qualquer a partir dos axiomas. O Turing provou que é impossível criar um método que analise um programa de computador e diga se ele termina ou se entra em loop infinito. O Matiyasevich provou que é impossível criar um método que resolva todas as equações diofantinas. O Richardson provou que é impossível criar um método que prove todas as identidades trigonométricas. Foi o fim da matemática automatizada... ou não?

Claro que não! Todos esses teoremas provam a impossibilidade do caso geral. Mas você sempre pode automatizar um caso específico. Não existe uma fórmula geral para a soma hipergeométrica, mas o algoritmo de Gosper prova que pelo menos um caso especial pode ser automatizado: o caso da soma hipergeométrica cujo resultado também é um termo hipergeométrico. E esse caso é imensamente útil na prática, especialmente se você estuda computação ou combinatória.

E no fim das contas, usar o Wolfram Alpha para resolver uma somatória é trapacear ou não? Por dentro, tudo que o Wolfram Alpha faz é usar o algoritmo de Gosper. Se você acredita que usar um método computacional é roubar, então por consistência precisa abrir de mão de todos os outros métodos computacionais. Não pode mais resolver a equação de segundo grau usando a fórmula de Bhaskara, não pode mais calcular o máximo divisor comum com o algoritmo de Euclides, não pode nem calcular quanto é 463 vezes 367 sem antes decorar a tabuada do 367.

Para mim isso não faz sentido. E por isso eu uso o Wolfram Alpha sem medo de ser feliz :)