sexta-feira, 27 de março de 2015

[P6] Brincando com Cadeias de Markov


     A modelagem com Cadeias de Markov está na sua simplicidade: apenas estados e transições e frequências de mudança entre os estados. O simples é facilmente entendível, necessário em muitos aspectos e fundamental. As seguintes frases representam bem este sentimento, aplicável a diversas esferas:

(1) "The KISS (Keep it simple, stupid) Principle."
(2) "Less is more."
(3) "Perfection is achieved, not when there is nothing more to add, but when there is nothing left to take away." - Antoine de Saint-Exupery
(4) "If you can't explain it to a six year old, you don't understand it yourself." – Albert Einstein
(5) "Simplicity is the ultimate sophistication." - Creditada a Leonardo da Vinci, apesar de não existir nenhum local indicando o mesmo.

     A seguir, utilizando estes princípios, é mostrado um estudo sobre as Cadeias de Markov, idealizadas por Andrey Andreyevitch Markov, circa 1905.
     O conceito chave de Cadeias de Markov deve ser obrigatoriamente simples de entender e modelar. Qualquer pessoa com conhecimentos básicos de Álgebra Linear (sim, será necessário...), abstração e criatividade pode utilizar as premissas simples de Cadeias de Markov para modelar e trabalhar com sistemas mais complexos e intrincados, realizando diversas análises e entendendo sua operação fundamental.
     A principal ideia por trás da modelagem Markoviana é considerar um sistema como tendo um espaço de estados bem definido com transições entre os estados e a passagem entre os estados não tem memória (memoryless property), ou seja, para chegar em um estado, não importa os estados visitados anteriormente, apenas os estados possíveis a partir do estado atual (vamos esquecer destes detalhes por enquanto e focar na descrição e solução de um problema).

Os humores do Prof. Steigerschülz

     O Prof. Steigerschülz (lê-se SH-TAY-GUER-SHI-L-TZ), trabalha em uma renomada instituição de ensino superior brasileira e é uma pessoa extremamente simples, inclusive quanto aos seus humores. O professor possui apenas três condições de humor: 'Feliz' (F), 'Irritado' (I) e 'Pensativo' (P). Um dos seus alunos resolveu documentar as mudanças de humor do Prof. Steigerschülz e chegou às seguintes conclusões:

(i) quando o professor está feliz, 30% das vezes ele fica irritado e 40% ele entra em um modo introspectivo de amplo pensamento (o professor chama este estado de profundo PENSO, e não gosta de ser importunado). No resto das vezes (30%) ele permanece feliz;
(ii) quando o professor está irritado, apenas uma piada faz com ele volte a ficar apenas feliz, em 50% das vezes;
(iii) quando está pensativo, pode ficar feliz (em 70% das vezes) ou irritado (em 30% das vezes);

     O problema é que, apesar de saber as vezes que o professor transita entre os estados, o aluno tem a seguinte dúvida: "chegando para trabalhar em um dia qualquer, dadas as chances verificadas acima, qual é a probabilidade do professor estar feliz, irritado ou pensativo?". Para responder esta pergunta, pode-se modelar o problema com uma Cadeia de Markov, mapeando os estados e as transições possíveis para este problema em particular, como segue:
Estados: F, I e P, respectivamente para Feliz, Irritado e Pensativo.
Transições: matriz M
     O modelo é mostrado na figura a seguir (observe que o modelo possui probabilidades nas transições entre os estados):
     Para criar a Matriz M (também chamada de Matriz de Transição ou Matriz Estocástica) é necessário criar uma matriz quadrada para representar os estados que o Prof. Steigerschülz assume:

M=
  F I P
F 0,3000 0,3000 0,4000
I 0,5000 0,5000 0,0000
P 0,7000 0,3000 0,0000

     A ideia é simular o dia-a-dia do professor e tentar encontrar um regime onde não existem mais mudanças (não é garantido que isto será conseguido, pois as cadeias devem ser ergódicas), ou seja, estimando seu comportamento no infinito. Para fazer isso, podem ser usados Métodos Diretos (Método de Gauss, fatoração LU, etc) ou Iterativos (Método da Potência, Subespaços de Krylov - Arnoldi, GMRES - etc) de solução (olhar a referência PHILIPPE 1989).
     Nesta postagem será utilizada a técnica de elevar uma matriz à sua enésima potência (URGENTE: revisar multiplicação entre matrizes de mesma ordem), como mostrado a seguir (estão sendo mostrados os valores arredondando para 4 casas após a vírgula, para simplificação do método - utilize o máximo de precisão sempre).

M2=
  F I P
F 0,5200 0,3600 0,1200
I 0,4000 0,4000 0,2000
P 0,3600 0,3600 0,2800
 
M3=
  F I P
F 0,4576 0,3744 0,1680
I 0,4400 0,3760 0,1840
P 0,4320 0,3744 0,1936





M4=
  F I P
F 0,4467 0,3750 0,1783
I 0,4463 0,3750 0,1787
P 0,4461 0,3750 0,1789
 
M5=
  F I P
F 0,4464 0,3750 0,1786
I 0,4464 0,3750 0,1786
P 0,4464 0,3750 0,1786

     Observe que ao elevar a matriz à quinta potência (M5), os valores das linhas já não se alteram mais. Este é o resultado (corresponde ao eigenvector, ou autovetor, da matriz em questão), e demonstra as probabilidades de permanência de cada estado do professor:
     Isto quer dizer que o Professor Steigerschülz está Feliz mais ou menos 45% (44,64% ou 0,4464) do seu dia (ou mês, ou ano). Mas, existe uma probabilidade não irrisória de 38% (37,50% ou 0,3750) do professor estar Irritado e, tecnicamente, não deveria ser importunado. Por fim, em 18% (17,86% ou 0,1786) do tempo, ele está Pensativo, ou seja, uma baixa probabilidade do professor estar em um estado de pensamento puro.

Questões relativas à solução e ao desempenho do processo

     De uma maneira relativamente simples, foi construído um modelo que foi resolvido e teve seus resultados obtidos. Entretanto, por mais interessante que seja a abordagem de elevar uma matriz à sua enésima potência, ela é computacionalmente custosa, pois consome muitos recursos (muitas multiplicações e somas cada vez). Isso é uma pena pois a convergência deste modelo em particular foi atingida logo na quinta potência, evidenciando o fato que este processo encontra resultados rapidamente, mas é custoso em termos de recursos dispendidos.
     Cabe ressaltar que a matriz dos modelos desta postagem possuem baixa ordem, ou seja, matrizes com poucas células, por exemplo, 3x3 (3 linhas e 3 colunas). Pense que o modelo do Professor Steigerschülz poderia considerar toda a gama das suas emoções, aumentando a ordem da matriz consideravelmente.
     Logo, é necessário que se pense em uma maneira alternativa de calcular estas probabilidades sem gastar muitos recursos computacionais. A técnica mostrada aqui consiste em multiplicar um vetor de probabilidades por uma matriz de transição (soma das linhas igual a 1). O vetor de probabilidades inicial pode ser arbitrário (ou seja, pode ser representado por quaisquer valores diferentes de zero) e uma matriz de transição M (que já tínhamos antes). Ao realizar este processo iterativamente, usando o vetor como entrada e calcular quando este não mais varia, chega-se às probabilidades de permanência de cada estado, ou seja, o nosso objetivo.

     Então, vamos utilizar a mesma matriz M definida anteriormente, mas desta vez, vamos criar um vetor de probabilidades inicial V, como mostrado abaixo.
Vetor de probabilidades inicial - V
1,0000 0,0000 0,0000

     Vamos utilizar este vetor V para multiplicar a matriz iteradas vezes (onde busca-se um vetor estacionário de probabilidades, ou seja, um vetor que ao ser multiplicado por uma matriz resulta no próprio vetor), usando o vetor como entrada de dados a cada iteração e observando um resíduo entre cada iteração do método (vamos supor 0,0002, ou seja, quando a diferença de cada elemento do vetor da iteração i menos o vetor da iteração (i-1) for menor que 0,0002 - neste caso, pois normalmente o resíduo é da ordem de 0,000000001). Nota: a função ABS retorna o valor absoluto de um parâmetro.
Iteração Vetor de probabilidades - V Resíduo Total
ABS(Vi,j - Vi-1,j)
0 1,0000 0,0000 0,0000 1,4000
1 0,3000 0,3000 0,4000 0,5600
2 0,5200 0,3600 0,1200 0,2000
3 0,4200 0,3720 0,2080 0,0800
4 0,4576 0,3744 0,1680 0,0310
5 0,4421 0,3749 0,1830 0,0124
6 0,4482 0,3750 0,1768 0,0050
7 0,4457 0,3750 0,1793 0,0020
8 0,4467 0,3750 0,1783 0,0008
9 0,4463 0,3750 0,1787 0,0004
10 0,4465 0,3750 0,1785 0,0002
11 0,4464 0,3750 0,1786 -

     Com um processo um pouco diferente, mas computacionalmente menos custoso, chega-se aos mesmos resultados de antes: o vetor de probabilidades final contém as probabilidades de permanência. Foram necessárias 11 iterações de multiplicações entre um vetor e uma matriz. Como a matriz é pequena, o ganho não é muito aparente, mas, à medida que a matriz aumenta, os custos começam a ficar bem mais interessantes.

Tipos de Cadeias de Markov

      As Cadeias de Markov vistas no exemplo do Professor Steigerschülz são conhecidas por serem definidas utilizando-se uma escala de tempo discreta (ou Discrete Time Markov Chains, DTMC). Isso quer dizer que o sistema evolui conforme passos discretos bem definidos (por isso que é dito que o tempo passa de forma discreta). O exemplo anterior modela um sistema que evolui conforme probabilidade em cada estado assumido pelo professor, entretanto, uma outra forma de modelagem Markoviana é possível, utilizando-se uma escala de tempo dita contínua (ou Continuous Time Markov Chains, CTMC), onde as transições entre os estados mapeiam taxas e a passagem do tempo segue premissas da distribuição exponencial (que é memoryless).
      Em termos de solução numérica, o procedimento é bem parecido, exceto que a matriz não é mais de probabilidades e sim de taxas e a diagonal principal da matriz contém a soma das taxas de saída de cada estado de forma negativada, ou seja, a linha soma zero. Estas matrizes especiais são chamadas pelo nome pomposo de Gerador Infinitesimal Q.

Matriz Exemplo
       
       
       
       

      Para resolver esta matriz, será necessário utilizar o conceito de Matriz Identidade. Uma matriz identidade possui valores iguais a 1 na diagonal principal e 0 nos demais elementos, conforme mostrado a seguir:

Matriz Identidade I3 (ordem 3)
1 0 0
0 1 0
0 0 1

      A ideia é modelar um sistema com taxas entre os estados, criar uma matriz que mapeie estas transições, modificar a diagonal principal para conter a soma das taxas negativadas e depois aplicar um mecanismo de solução numérica.
      Considere o seguinte problema: um aluno possui três estados, trabalhando (T), estudando (E) e surfando na internet (I). Observou-se que este aluno em questão visita os estados conforme taxas ou frequências (estes dados foram obtidos a partir de observações ou logs de entrada em laboratórios de estudo/Internet, etc), que aconteceram conforme uma unidade de tempo (horas, minutos ou dias).

     A figura a seguir explica o modelo (observe que existem taxas nas transições entre os estados). Pode-se pensar, por exemplo, que o aluno em questão, em uma hora (ou outra unidade de tempo qualquer), vai de Estudando para Internet com uma frequência de 30 vezes (também pode-se considerar tempo em cada estado, uma vez que Tempo=1/Frequência).
     As taxas com que este aluno visita os estados estão definidas na matriz a seguir.

Matriz de Taxas T
  Trabalhando  Estudando      Internet    
Trabalhando   10 20
 Estudando  5   30
Internet 7 8  

     Para transformar esta Matriz de Taxas T em um Gerador Infinitesimal Q (onde a diagonal contém a correção das taxas do modelo) precisamos somar os valores das taxas da linha e colocar este número na diagonal, de forma negativada, conforme a matriz a seguir:

Gerador Infinitesimal Q
  Trabalhando  Estudando      Internet    
Trabalhando -30 10 20
 Estudando  5 -35 30
Internet 7 8 -15

     A ideia agora é criar uma nova matriz, onde a soma das linhas é igual a 1, usando a matriz identidade. Então, segue-se a seguinte operação: cada célula desta nova matriz contém o valor da célula correspondente da Matriz Identidade menos o valor da célula da matriz sendo calculada dividido pelo maior valor da matriz (em módulo). O maior valor (em módulo) da matriz certamente encontra-se na diagonal. Neste caso, o maior valor (em módulo), ou MAXij=-35 encontra-se na linha 2 e coluna 2.

Matriz de Transição M
0,1429 0,2857 0,5714
0,1429 0,0000 0,8571
0,2000 0,2286 0,5714

     Por exemplo, vamos calcular juntos a célula M00 que corresponde à célula da primeira linha e primeira coluna: equivale a pegar o valor da célula da matriz identidade (que é igual a 1), menos o valor da célula (-30) dividido por (-35).
Isso é igual a 1-(-30/-35)=1-30/35=1-0,8571=0,1429.
     Observe que esta nova matriz é uma matriz onde a soma das linhas é igual a 1, ou seja, voltamos ao problema anterior das Cadeias de Markov do tipo DTMC, onde a solução foi descrita anteriormente. Então, usa-se um método numérico de solução, ou seja, eleva-se a matriz à uma potência ou pode-se utilizar a multiplicação de um vetor de probabilidades (que guarda o resultado) por uma matriz de transição. A seguir é mostrado a operação de elevar a matriz à uma potência:

M2=
  T E I
T 0,1755 0,1714 0,6531
E 0,1918 0,2367 0,5714
I 0,1755 0,1878 0,6367
 
M3=
  T E I
T 0,1783 0,1933 0,6284
E 0,1794 0,1962 0,6244
I 0,1786 0,1941 0,6273
 
M4=
  T E I
T 0,1787 0,1944 0,6270
E 0,1787 0,1944 0,6270
I 0,1787 0,1944 0,6270

     Este é o resultado, e demonstra as probabilidades de permanência de cada estado do professor:

Probabilidades de permanência (resultado)
Trabalhando 17,87%
Estudando 19,44%
Internet 62,70%
     Os resultados mostram que o aluno em questão encontra-se aproximadamente 63% do seu tempo surfando na Internet e o resto do tempo mais ou menos dividido entre Trabalhando (18%) e Estudando (20%). Esta análise não era clara antes de modelar o problema em uma Cadeia de Markov e solucionar o problema obtendo as probabilidades de permanência de cada estado. Isso demonstra o poder e a simplicidade do formalismo.
     Pode-se utilizar a multiplicação do vetor de probabilidades pela matriz de transição para resolver o problema, assim como poderia ter sido utilizado um método direto ou outro método iterativo.

Como replicar estes resultados de forma fácil

     É possível replicar tudo que foi feito nesta postagem utilizando-se uma ferramenta específica para este propósito (Matlab, Maple, Matematica ou resolvedor próprio de Cadeias de Markov) ou então utilizando-se o MS-Excel ou qualquer planilha eletrônica que forneça operações com matrizes.
     Aqui, a ferramenta a ser utilizada é o MS-Excel, e a ideia é utilizar a função MATRIZ.MULT().
     Para este exemplo, estou assumindo que você já está trabalhando com a matriz na forma DTMC (soma das linhas igual a 1). Depois de criar uma matriz, vamos supor, nas células D3:F5, vamos selecionar nove células em branco (em outro local da planilha), vamos supor, D7:F9, apertar a tecla F2 e inserir a fórmula =MATRIZ.MULT(D3:F5;D3:F5) e então, magicamente, apertar a sequência de teclas CTRL-SHIFT-ENTER. Estas teclas realizam a operação da multiplicação de matrizes no MS-Excel. Faça isso até a matriz não variar mais nas linhas. Você chegará aos resultados apresentados aqui. O mesmo procedimento pode ser feito para realizar a operação da multiplicação do vetor pela matriz.
OBSERVAÇÃO: no MS-Excel, você *deve* selecionar as células em branco - onde serão calculados os novos valores da matriz - e então apertar F2, inserir a fórmula e teclar CTRL-SHIFT-ENTER. Se estes passos não forem seguidos, não funcionará!

Implementando seu próprio resolvedor

     Implementar um resolvedor de Cadeias de Markov é simples. Basta criar um formato de entrada para receber matrizes descritas com modelos Markovianos do tipo DTMC ou CTMC e entrar outras informações secundárias (mas não obrigatórias), tais como a descrição dos estados. Por exemplo, a seguir, é mostrado um esboço de formato de entrada (arquivo chamado markov.txt):
#FELIZ
#IRRITADO
#PENSATIVO
0.3 0.3 0.4
0.5 0.5 0.0
0.7 0.3 0.0
     Para facilitar, observe que a entrada de dados está no formato americano, com pontos ao invés de vírgulas. Isso facilita a leitura do arquivo, pois este já se encontra com os valores necessários.
     Você precisa fazer a consistência dos dados de entrada, por exemplo, duas consistências obrigatórias seriam (1) o número de descrições de estados equivalem ao total de estados informados? (três linhas e três colunas); (2) a soma das linhas fecha 1? Caso estes problemas sejam detectados, você deve informar ao usuário. Outros problemas que talvez aconteçam é a divergência do modelo. Para contornar este problema, você deve implementar o conceito de total de iterações, que o usuário pode passar no arquivo (#iteracoes=1000) ou pela linha de comando (./markovsolver 1000 para o GNU/Linux). A seguir, é mostrado um exemplo onde passam-se outros parâmetros para o resolvedor:
#iteracoes=1000
#residuo=0.0000001
#precisao=2
...
     A saída poderia ser um arquivo ou a tela, contendo os resultados obtidos:
O modelo obteve a convergência em 15 iterações 
e os resultados estão apresentados a seguir:
FELIZ=44,64%
IRRITADO=37,50%
PENSATIVO=17,86%
Tempo de execução: 00m15s user time
Data: 23/03/2015 09:22:05
     Ou então: "Para X iterações, não foi obtida a convergência do modelo".
     Para os modelos CTMC, não é necessário informar os valores da diagonal principal, já que estes podem ser calculados diretamente por você no programa.
#iteracoes=1000
#residuo=0,0000001
#precisao=2
0 10 20
5 0 30
7 8 0
     O ideal é o seu programa identificar se o modelo é DTMC (soma das linhas é igual a 1) ou CTMC (matriz de taxas). Para resolver o modelo, utilize sempre a técnica iterativa do Método da Potência usando a multiplicação de um vetor por uma matriz. Se for uma CTMC, você deve calcular a diagonal principal, negativá-la e então retirar a matriz identidade e dividir cada célula pelo maior valor (em módulo) da matriz. Se isso não for feito, a solução não funcionará.

Considerações especiais

     A pesquisa em Cadeias de Markov evolui constantemente onde pesquisadores trabalham com formalismos ditos estruturados que trabalham com representações de sistemas situadas acima das Cadeias de Markov, ou seja, definem um formato de descrição de sistemas que gera uma Cadeia de Markov (underlying Markov Chain). Estes formatos permitem a fácil decomposição de sistemas e utilizam primitivas básicas de modelagem que auxiliam os usuários a mapear e entender sistemas complexos, além das análises proporcionadas. Exemplos de tais formalismos são as Redes de Petri Estocásticas, as Álgebras de Processos de Sistemas e suas variações, como Redes de Petri Coloridas, Superpostas, Generalizadas, entre outros exemplos.
     Cada formalismo estruturado apresenta e explica sua lista de primitivas de modelagem e uma técnica específica de solução do modelo. Por exemplo, as Álgebras de Processo combinam os estados dos seus inúmeros processos em uma grande Cadeia de Markov onde os usuários podem usar qualquer ferramenta matemática de solução (por exemplo, GNU/Octave, Matlab, Maple, Matematica, etc). Outros formalismos concentram-se em mecanismos que representam a grande Cadeia de Markov de forma eficiente em memória, utilizando por exemplo, Descritores Markovianos.
     Existem muitas pesquisas correntes em formalismos estruturados baseados em ideias de Cadeias de Markov e palavras-chave são: structured stochastic formalism, process algebra e stochastic petri nets.

Epílogo

     Cadeias de Markov são amplamente utilizadas e empregadas em diversos contextos e domínios de aplicação, desde economia até o Algoritmo principal executado pelo Google (chamado de PageRank, mas isso é uma história para outra postagem...).
     Pesquisadores ao redor do mundo estudam Cadeias de Markov e processos Markovianos com aplicação, por exemplo, em Sistemas Paralelos e Distribuídos, biologia marítima, química, física de multidões, e muitas outras aplicações.
     Se você se interessou pelo assunto, procure na Internet outras informações sobre Cadeias de Markov.
     Alguns links interessantes sobre o assunto:
1. Markov Chains @ Khan Academy
2. Material sobre Cadeias de Markov bem explicativo (em inglês)
3. Gerador de textos
4. Markov and You
... ...

Referências

(PHILIPPE et al. 1989) Bernard Philippe, Yousef Saad, William Stewart. Numerical methods in Markov chain modeling. [Research Report] RR-1115, 1989. . Disponível em HAL

Questões interessantes deste projeto
  1. Crie as estruturas de dados para trabalhar com vetores de alta-precisão (double) e matrizes
  2. Faça uma função que testa se uma matriz tem a soma das linhas iguais a 1,0000
  3. Implemente a leitura de um arquivo de entrada (modelo.txt) e a saída em um arquivo contendo os resultados da convergência ou divergência do modelo (mesmo nome de entrada com terminação .out, por exemplo, modelo.out)
  4. Eleve uma matriz à uma potência qualquer (se a soma das linhas for igual a 1,0000)
  5. Eleve uma matriz até que cada linha da matriz não é mais alterada conforme um resíduo igual a 0,00000001, ou informado em um arquivo de entrada (o resíduo é verificado quando a diferença entre os valores das linhas não excede um valor predeterminado, como 0,00000001)
  6. Faça uma função que multiplica, iteradas vezes (até que o vetor não é mais alterado conforme um resíduo - mesmo do item anterior), não deixando que o método ultrapasse X iterações, onde X é informada em um arquivo de entrada (caso ultrapasse, avisar o usuário)
  7. Faça uma função que lê um arquivo contendo uma matriz de um modelo Markoviano (DTMC ou CTMC - você deve decidir o que a matriz é em tempo de execução e converte-la caso necessário) e retorna o vetor de probabilidades de permanência caso exista convergência até 1000 iterações
  8. Faça as funções para um valor de iterações configurável (ITER=1000, mas podendo variar), passando por linha de comando ou por arquivo de entrada
  9. Faça os resultados aparecerem em um novo arquivo gerado pela sua ferramenta
  10. Implemente uma ferramenta visual para abrir modelos e resolver Cadeias de Markov (dica: usar o wxWidgets)

Nenhum comentário:

Postar um comentário