Convergência de um processo de Markov

10

Desafio

Dada uma matriz estocástica esquerda ou direita, em que o limite em que x se aproxima do infinito da matriz à potência de x se aproxima de uma matriz com todos os valores finitos, retorne a matriz à qual a matriz converge. Basicamente, você deseja continuar multiplicando a matriz por si só até que o resultado não seja mais alterado.

Casos de teste

[[7/10, 4/10], [3/10, 6/10]] -> [[4/7, 4/7], [3/7, 3/7]]
[[2/5, 4/5], [3/5, 1/5]] -> [[4/7, 4/7], [3/7, 3/7]]
[[1/2, 1/2], [1/2, 1/2]] -> [[1/2, 1/2], [1/2, 1/2]]
[[1/3, 2/3], [2/3, 1/3]] -> [[1/2, 1/2], [1/2, 1/2]]
[[1/10, 2/10, 3/10], [4/10, 5/10, 6/10], [5/10, 3/10, 1/10]] -> [[27/130, 27/130, 27/130], [66/130, 66/130, 66/130], [37/130, 37/130, 37/130]]
[[1/7, 2/7, 4/7], [2/7, 4/7, 1/7], [4/7, 1/7, 2/7]] -> [[1/3, 1/3, 1/3], [1/3, 1/3, 1/3], [1/3, 1/3, 1/3]]

Regras

  • As brechas padrão se aplicam
  • Você pode escolher se deseja uma matriz estocástica direita ou esquerda
  • Você pode usar qualquer tipo de número razoável, como valores flutuantes, racionais, decimais de precisão infinita, etc.
  • Isso é , portanto, o envio mais curto em bytes para cada idioma é declarado vencedor em seu idioma. Nenhuma resposta será aceita
HyperNeutrino
fonte
@FryAmTheEggman parece que alguns comentários anteriores foram excluídos, portanto isso pode ser redundante, mas matrizes redutíveis e periódicas já são excluídas por "Dada uma matriz estocástica esquerda ou direita, em que o limite em que x se aproxima do infinito da matriz ao poder de x se aproxima de uma matriz com todos os valores finitos ", que li como dizendo que a entrada é garantida para convergir para uma solução única. (ou seja, a matriz de entrada é garantida para ser ergodic.)
Nathaniel
@ Nathaniel Isso não é bem verdade, como se a cadeia fosse redutível, você pode obter um resultado (como para a matriz de identidade) que atenda ao que você disse, mas a cadeia descrita não é irredutível e, portanto, não é garantido que a entrada seja ser ergódico (já que não será recorrente positivo). Garantir a ergodicidade é o que o OP deseja, e acho que eles têm isso agora, graças à restrição adicional de todos os valores de linha serem idênticos. Se você conhece uma maneira melhor de restringir a entrada sem precisar adicionar uma explicação das cadeias de Markov, tenho certeza que o HyperNeutrino apreciaria! :)
FryAmTheEggman
11
@FryAmTheEggman ah, você está certo, desculpe. Eu estava pensando em iteração de poder em vez de elevar a matriz a um poder. (Por "solução única", eu quis dizer "uma que seja independente do ponto de partida do processo de iteração", mas que não seja relevante aqui.) Concordo que a condição 'todas as linhas são idênticas' faz o trabalho. Suponho que o OP poderia apenas dizer "a cadeia de Markov é garantida como ergódica", o que satisfaria pessoas como nós, que provavelmente se preocuparão com isso!
Nathaniel
Na verdade, se B é uma solução para BA = B , o mesmo acontece com cB para qualquer constante escalar c . Portanto, uma solução diferente de zero não pode ser estritamente única, a menos que você determine a escala. (Exigir que B seja estocástico faria isso.) Além disso, obviamente, se as linhas ou as colunas de B iguais são dependentes de A ser estocástico esquerdo ou direito.
Ilmari Karonen 8/0318
2
Para quem nunca aprendeu sobre matrizes durante a aula de matemática no ensino médio e como multiplicá-las: mathsisfun.com/algebra/matrix-multiplying.html . Eu tive que procurá-lo para entender o que estava sendo solicitado .. Talvez seja de conhecimento comum em outros lugares na Terra ..: S
Kevin Cruijssen

Respostas:

7

R ,  44  43 bytes

function(m){X=m
while(any(X-(X=X%*%m)))0
X}

Experimente online!

Apenas continua se multiplicando até encontrar uma matriz fixa. Aparentemente X!=(X=X%*%m), a comparação é reatribuída X, o que é legal.

Obrigado ao @Vlo por remover um byte; apesar de riscado 44 ainda é regular 44.

Giuseppe
fonte
11
Eu me pergunto por function(m){ while(any(m!=(m=m%*%m)))0 m}que não funciona. Imprecisões numéricas que impedem que a condição de terminação seja acionada?
CodesInChaos
@CodesInChaos provavelmente é falta de precisão. Mudar para uma biblioteca de precisão arbitrária também não ajuda - eles atingem o tempo limite (Rmpfr) ou falham (gmp) da mesma maneira, embora eu provavelmente esteja fazendo algo errado.
Giuseppe
@ Giuseppe Acho que a abordagem sugerida é repetida ao quadrado até que não mude mais? (Eu não sei ler R) #
3120 user202729 16:00
@ user202729 sim, é. R usa números de ponto flutuante de 64 bits e sei que os erros se propagam rapidamente.
Giuseppe
Eu acho que esse algoritmo é instável. A geléia também tem o mesmo problema. (TODO provar isso e encontrar uma alternativa)
user202729
5

Octave ,45 42. 35 bytes

@(A)([v,~]=eigs(A,1))/sum(v)*any(A)

Experimente online!

Economizou 3 bytes graças a Giuseppe e mais 7 graças a Luis Mendo!

Isso usa que o vetor próprio correspondente ao valor próprio 1 (também o valor máximo máximo) é o vetor da coluna que é repetido para cada valor da matriz limitante. Temos que normalizar o vetor para ter a soma 1 para que seja estocástico, depois basta repeti-lo para preencher a matriz. Não sou muito versado em golfe no Octave, mas não consegui encontrar uma maneira funcional de fazer a multiplicação repetida, e um programa completo parece que sempre será mais longo do que isso.

Podemos usar any(A)uma vez que, a partir das restrições, sabemos que a matriz deve descrever uma cadeia de Markov irredutível e, portanto, cada estado deve ser acessível a partir dos outros estados. Portanto, pelo menos um valor em cada coluna deve ser diferente de zero.

FryAmTheEggman
fonte
Como eigssempre retorna o vetor próprio correspondente 1? Minha memória das cadeias de markov é um pouco confusa.
Giuseppe
42 bytes
Giuseppe
@ Giuseppe Como a matriz é estocástica, e algumas outras coisas, seu valor próprio máximo é 1 e eigsretorna a partir do maior valor próprio. Além disso, obrigado pelo golfe!
FryAmTheEggman
Ah, certo. As cadeias de Markov estão no meu próximo exame, mas como é para atuários, alguns detalhes estão faltando.
Giuseppe
3

APL (Dyalog) , 18 6 bytes

12 bytes salvos graças a @ H.PWiz

+.×⍨⍣≡

Experimente online!

Uriel
fonte
+.×⍨⍣≡por 6 bytes. Isto é, praça até que nada muda
H.PWiz
@ H.PWiz Eu acredito que sim. Eu não tenho idéia do por que eu não tinha feito isso em primeiro lugar. Obrigado!
Uriel
3

k / q, 10 bytes

k / q porque o programa é idêntico nos dois idiomas. O código é realmente apenas k / q idiomático.

{$[x]/[x]}

Explicação

{..}está fora da sintaxe lambda, com xum parâmetro implícito

$[x] tem $ como operador de multiplicação de matriz binária, fornecendo apenas um parâmetro cria um operador unário que é multiplicado pela matriz de Markov

/[x] aplica a multiplicação esquerda até a convergência, começando com o próprio x.

ulucs
fonte
3

C (gcc) , 207 192 190 181 176 bytes + 2 bytes de flag-lm

  • Economizou quinze dezessete e vinte e dois bytes graças ao ceilingcat .
  • Salva nove bytes; removendo return A;.
float*f(A,l,k,j)float*A;{for(float B[l*l],S,M=0,m=1;fabs(m-M)>1e-7;wmemcpy(A,B,l*l))for(M=m,m=0,k=l*l;k--;)for(S=j=0;j<l;)m=fmax(m,fdim(A[k],B[k]=S+=A[k/l*l+j]*A[k%l+j++*l]));}

Experimente online!

Jonathan Frech
fonte
@ceilingcat Contando os bytes do sinalizador do compilador, isso resulta em 192 novamente. Incluiu sua sugestão, no entanto.
precisa
@ceilingcat Obrigado.
Jonathan Frech
2

Python 3 , 75 61 bytes

f=lambda n:n if allclose(n@n,n)else f(n@n)
from numpy import*

Experimente online!

Nos casos de teste, existem imprecisões de flutuação, portanto os valores podem diferir um pouco das frações exatas.

PS. numpy.allclose()é usado porque numpy.array_equal()é mais longo e propenso a flutuações de imprecisões.

-14 bytes Obrigado HyperNeutrino;) Ah, sim, esqueci o operador @; P

Shieru Asakoto
fonte
Use em dotvez de matmul: D
HyperNeutrino
Na verdade, ter matrizes numpy como entrada e fazer x=n@n: P tio.run/...
HyperNeutrino
59 bytes
HyperNeutrino
Adicionado de volta o f=em frente porque é recursivamente chamado;)
Shieru Asakoto
Oh sim, você está certo :) boa ligação!
HyperNeutrino
2

Java 8, 356 339 bytes

import java.math.*;m->{BigDecimal t[][],q;RoundingMode x=null;for(int l=m.length,f=1,i,k;f>0;m=t.clone()){for(t=new BigDecimal[l][l],i=l*l;i-->0;)for(f=k=0;k<l;t[i/l][i%l]=(q!=null?q:q.ZERO).add(m[i/l][k].multiply(m[k++][i%l])))q=t[i/l][i%l];for(;++i<l*l;)f=t[i/l][i%l].setScale(9,x.UP).equals(m[i/l][i%l].setScale(9,x.UP))?f:1;}return m;}

-17 bytes graças a @ceilingcat .

Definitivamente não é o idioma certo .. Maldita precisão de ponto flutuante ..

Explicação:

Experimente online.

import java.math.*;     // Required import for BigDecimal and RoundingMode
m->{                    // Method with BigDecimal-matrix as both parameter and return-type
  BigDecimal t[][],q;   //  Temp BigDecimal-matrix
  RoundingMode x=null;  //  Static RoundingMode value to reduce bytes
  for(int l=m.length,   //  The size of the input-matrix
          f=1,          //  Flag-integer, starting at 1
          i,k;          //  Index-integers
      f>0;              //  Loop as long as the flag-integer is still 1
      m=t.clone()){     //    After every iteration: replace matrix `m` with `t`
    for(t=new BigDecimal[l][l],
                        //   Reset matrix `t`
        i=l*l;i-->0;)   //   Inner loop over the rows and columns
      for(f=k=0;        //    Set the flag-integer to 0
          k<l           //    Inner loop to multiply
          ;             //      After every iteration:
           t[i/l][i%l]=(q!=null?q:q.ZERO).add(
                        //       Sum the value at the current cell in matrix `t` with:
             m[i/l][k]  //        the same row, but column `k` of matrix `m`,
             .multiply(m[k++][i%l])))
                        //        multiplied with the same column, but row `k` of matrix `m`
        q=t[i/l][i%l];  //     Set temp `q` to the value of the current cell of `t`
    for(;++i<l*l;)      //   Loop over the rows and columns again
      f=t[i/l][i%l].setScale(9,x.UP).equals(m[i/l][i%l].setScale(9,x.UP))?
                        //    If any value in matrices `t` and `m` are the same:
         f              //     Leave the flag-integer unchanged
        :               //    Else (they aren't the same):
         1;}            //     Set the flag-integer to 1 again
  return m;}            //  Return the modified input-matrix `m` as our result
Kevin Cruijssen
fonte
Por que a função principal é tão detalhada?
user202729
@ user202729 Porque float/ doublenão tem a precisão de ponto flutuante adequada, java.math.BigDecimaldeve ser utilizada. E, em vez de simplesmente +-*/, BigDecimals usar .add(...), .subtract(...), .multiply(...), .divide(...). Então, algo tão simples quanto 7/10se torna new BigDecimal(7).divide(new BigDecimal(10)). Além disso, o ,scale,RoundingModeno divideé necessário para valores com valores decimais 'infinitos' (como 1/3ser 0.333...). O método principal pode, é claro, ser jogado no golfe, mas eu não me incomodei quando fiz uma pesquisa e substituição para converter os carros alegóricos em BigDecimals.
21718 Kevin Kevin Krijssen
@ceilingcat Thanks!
Kevin Cruijssen