Konvergenz eines Markov-Prozesses

10

Herausforderung

Bei einer links- oder rechtsstochastischen Matrix, bei der sich die Grenze, wenn sich x der Unendlichkeit der Matrix nähert, der Potenz von x einer Matrix mit allen endlichen Werten nähert, die Matrix zurückgeben, zu der die Matrix konvergiert. Grundsätzlich möchten Sie die Matrix so lange mit sich selbst multiplizieren, bis sich das Ergebnis nicht mehr ändert.

Testfälle

[[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]]

Regeln

  • Standardschlupflöcher gelten
  • Sie können wählen, ob Sie eine rechts- oder eine linksstochastische Matrix wünschen
  • Sie können jeden vernünftigen Zahlentyp verwenden, z. B. Gleitkommazahlen, Rationals, Dezimalstellen mit unendlicher Genauigkeit usw.
  • Dies ist , daher wird die kürzeste Übermittlung in Bytes für jede Sprache zum Gewinner ihrer Sprache erklärt. Es wird keine Antwort akzeptiert
HyperNeutrino
quelle
@FryAmTheEggman Es scheint, dass einige frühere Kommentare gelöscht wurden, so dass dies möglicherweise redundant ist, aber reduzierbare und periodische Matrizen werden bereits ausgeschlossen durch "Bei einer links- oder rechtsstochastischen Matrix, bei der sich die Grenze als x der Unendlichkeit der Matrix zur Potenz nähert of x nähert sich einer Matrix mit allen endlichen Werten ", was ich so gelesen habe, dass die Eingabe garantiert zu einer einzigartigen Lösung konvergiert. (dh die Eingabematrix ist garantiert ergodisch.)
Nathaniel
@ Nathaniel Das ist nicht ganz richtig, als ob die Kette reduzierbar wäre, können Sie ein Ergebnis (wie für die Identitätsmatrix) erzielen, das Ihren Aussagen entspricht, aber die darin beschriebene Kette ist nicht irreduzibel und daher wäre die Eingabe nicht garantiert ergodisch sein (da es nicht immer wieder positiv sein wird). Die Gewährleistung der Ergodizität ist das, was das OP wünscht, und ich denke, das haben sie jetzt, dank der zusätzlichen Einschränkung, dass alle Zeilenwerte identisch sind. Wenn Sie einen besseren Weg kennen, um die Eingabe einzuschränken, ohne eine Erklärung der Markov-Ketten hinzufügen zu müssen, würde HyperNeutrino das sicher zu schätzen wissen! :)
FryAmTheEggman
1
@FryAmTheEggman ah, du hast recht, sorry. Ich dachte eher an Power-Iteration, als die Matrix zu einer Power zu machen. (Mit "eindeutige Lösung" meinte ich "eine, die vom Startpunkt des Iterationsprozesses unabhängig ist", aber das ist hier nicht relevant.) Ich stimme zu, dass die Bedingung "Alle Zeilen sind identisch" den Job erledigt. Ich nehme an, das OP könnte einfach sagen "Die Markov-Kette ist garantiert ergodisch", was Leute wie uns zufriedenstellen würde, die sich wahrscheinlich darüber Sorgen machen werden!
Nathaniel
Wenn B eine Lösung für BA = B ist , gilt dies auch für cB für jede Skalarkonstante c . Eine Lösung ungleich Null kann also nicht eindeutig sein, es sei denn, Sie legen die Skala irgendwie fest. (Das Erfordernis, dass B stochastisch ist, würde dies tun.) Es hängt natürlich auch davon ab, ob A links oder rechts stochastisch ist , ob die Zeilen oder Spalten von B gleich sind .
Ilmari Karonen
2
Für alle anderen, die während des Mathematikunterrichts in der High School noch nie etwas über Matrizen gelernt haben und wie man sie multipliziert: mathsisfun.com/algebra/matrix-multiplying.html . Ich musste nachschlagen, um zu verstehen, was gefragt wurde. Vielleicht ist es anderswo auf der Erde allgemein bekannt. S
Kevin Cruijssen

Antworten:

7

R ,  44  43 Bytes

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

Probieren Sie es online aus!

Multipliziert einfach weiter, bis eine feste Matrix gefunden wird. Anscheinend X!=(X=X%*%m)führt der Vergleich dann eine Neuzuweisung durch X, das ist also ordentlich.

Vielen Dank an @Vlo für das Rasieren eines Bytes. obwohl durchgestrichen 44 ist immer noch regulär 44.

Giuseppe
quelle
1
Ich frage mich, warum function(m){ while(any(m!=(m=m%*%m)))0 m}nicht funktioniert. Numerische Ungenauigkeiten, die verhindern, dass die Beendigungsbedingung ausgelöst wird?
CodesInChaos
@CodesInChaos höchstwahrscheinlich ist es ein Mangel an Präzision. Das Wechseln zu einer Bibliothek mit beliebiger Genauigkeit hilft auch nicht - sie haben entweder eine Zeitüberschreitung (Rmpfr) oder schlagen auf die gleiche Weise fehl (gmp), obwohl ich wahrscheinlich etwas falsch mache.
Giuseppe
@ Giuseppe Ich denke, der vorgeschlagene Ansatz wird wiederholt quadriert, bis sich nichts mehr ändert? (Ich kann R nicht lesen)
user202729
@ user202729 Ja, das ist es. R verwendet 64-Bit-Gleitkommazahlen und ich weiß, dass sich Fehler ziemlich schnell ausbreiten.
Giuseppe
Ich denke, dass der Algorithmus instabil ist. Gelee hat das gleiche Problem auch. (TODO beweisen Sie dies und finden Sie eine Alternative)
user202729
5

Oktave ,45 42 35 Bytes

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

Probieren Sie es online aus!

3 Bytes dank Giuseppe und 7 weitere dank Luis Mendo!

Dies verwendet, dass der dem Eigenwert 1 entsprechende Eigenvektor (auch der maximale Eigenwert) der Spaltenvektor ist, der für jeden Wert der Grenzmatrix wiederholt wird. Wir müssen den Vektor normalisieren, um die Summe 1 zu haben, damit er stochastisch ist. Dann wiederholen wir ihn einfach, um die Matrix auszufüllen. Ich bin nicht besonders gut mit dem Golfen von Octave vertraut, aber ich konnte keinen funktionalen Weg finden, um wiederholte Multiplikationen durchzuführen, und ein vollständiges Programm scheint immer länger zu sein.

Wir können verwenden, any(A)da wir aufgrund der Einschränkungen wissen, dass die Matrix eine irreduzible Markov-Kette beschreiben muss und daher jeder Zustand von den anderen Zuständen aus erreichbar sein muss. Daher muss mindestens ein Wert in jeder Spalte ungleich Null sein.

FryAmTheEggman
quelle
Wie wird eigsimmer der entsprechende Eigenvektor zurückgegeben 1? Meine Erinnerung an Markov-Ketten ist etwas verschwommen.
Giuseppe
42 Bytes
Giuseppe
@Giuseppe Da die Matrix stochastisch ist und einige andere Dinge, ist ihr maximaler Eigenwert 1 und eigskehrt ausgehend vom größten Eigenwert zurück. Danke auch für den Golf!
FryAmTheEggman
Ah richtig. Markov-Ketten sind bei meiner nächsten Prüfung, aber da es sich um Aktuare handelt, fehlen einige Details vollständig.
Giuseppe
3

APL (Dyalog) , 18 6 Bytes

12 Bytes dank @ H.PWiz gespeichert

+.×⍨⍣≡

Probieren Sie es online aus!

Uriel
quelle
+.×⍨⍣≡für 6 Bytes. Das heißt, quadratisch, bis sich nichts mehr ändert
H.PWiz
@ H.PWiz Ich glaube es ist. Ich habe keine Ahnung, warum ich es überhaupt nicht getan habe. Vielen Dank!
Uriel
3

k / q, 10 Bytes

k / q, weil das Programm in beiden Sprachen identisch ist. Der Code ist wirklich nur idiomatisch k / q.

{$[x]/[x]}

Erläuterung

{..}ist aus Lambda-Syntax, mit xals impliziter Parameter

$[x] hat $ als Multiplikationsoperator für die binäre Matrix. Wenn nur ein Parameter angegeben wird, wird ein unärer Operator erstellt, der Multiplikationen mit der Markov-Matrix hinterlässt

/[x] wendet die linke Multiplikation bis zur Konvergenz an, beginnend mit x selbst.

ulucs
quelle
3

C (gcc) , 207 192 190 181 176 Bytes + 2 Flagbytes-lm

  • Gespeichert fünfzehn siebzehn zwanzig-zwei Bytes dank ceilingcat .
  • Neun Bytes gespeichert; entfernen 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]));}

Probieren Sie es online aus!

Jonathan Frech
quelle
@ceilingcat Wenn die Compiler-Flag-Bytes gezählt werden, ergibt sich erneut 192. Hat Ihren Vorschlag trotzdem aufgenommen.
Jonathan Frech
@ceilingcat Danke.
Jonathan Frech
2

Python 3 , 75 61 Bytes

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

Probieren Sie es online aus!

In Testfällen treten Float-Ungenauigkeiten auf, sodass die Werte geringfügig von den genauen Brüchen abweichen können.

PS. numpy.allclose()wird verwendet, weil numpy.array_equal()es länger ist und zu Ungenauigkeiten neigt.

-14 Bytes Danke HyperNeutrino;) Oh ja, ich habe den @ -Operator vergessen; P.

Shieru Asakoto
quelle
Verwenden Sie dotanstelle von matmul: D
HyperNeutrino
Nehmen Sie tatsächlich numpy Arrays als Eingabe und führen Siex=n@n
Folgendes aus
59 Bytes
HyperNeutrino
Hinten hinten hinzugefügt, f=weil es rekursiv genannt wird;)
Shieru Asakoto
Oh ja du hast recht :) guter Anruf!
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 dank @ceilingcat .

Auf jeden Fall nicht die richtige Sprache. Verdammte Gleitkommapräzision.

Erläuterung:

Probieren Sie es online aus.

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
quelle
Warum ist die Hauptfunktion so ausführlich?
user202729
@ user202729 Da float/ doublenicht die richtige Gleitkommapräzision haben, java.math.BigDecimalsollte stattdessen verwendet werden. Und statt einfach +-*/, BigDecimals verwenden .add(...), .subtract(...), .multiply(...), .divide(...). Also etwas so einfach wie es 7/10wird new BigDecimal(7).divide(new BigDecimal(10)). Außerdem sind die ,scale,RoundingModein the dividefür Werte mit 'unendlichen' Dezimalwerten (wie 1/3Sein 0.333...) erforderlich . Die Hauptmethode kann natürlich Golf gespielt werden, aber ich habe mich nicht darum gekümmert, als ich gesucht und ersetzt habe, um die Floats in BigDecimals umzuwandeln.
Kevin Cruijssen
@ceilingcat Danke!
Kevin Cruijssen