-1

I have found the right result of an Algorithm. Can I implement this code in Java?

import numpy as np
from scipy.linalg import eig 
transition_mat = np.matrix([
    [0.8,  0.15,0.05 ],\
    [0.075,0.85,0.075],\
    [0.05, 0.15,0.8  ]])

S, U = eig(transition_mat.T)
stationary = np.array(U[:np.where(np.abs(S - 1.) < 1e-8)[0][0]].flat)
print stationary
print np.sum(stationary)
stationary = stationary / np.sum(stationary)

print stationary

I have implemented this code in Java but the result is wrong

Matrix A = new Matrix(transition);
        A = A.transpose();
        Matrix x = new Matrix(N, 1, 1.0 / N); // initial guess for eigenvector
        for (int i = 0; i < 50; i++) {
            x = A.times(x);
            x = x.times(1.0 / x.norm1());       // rescale
        }
        
        // compute by finding eigenvector corresponding to eigenvalue = 1
        EigenvalueDecomposition eig = new EigenvalueDecomposition(A);
        Matrix V = eig.getV();
        double[] real = eig.getRealEigenvalues();
        for (int i = 0; i < N; i++) {
            if (Math.abs(real[i] - 1.0) < 1E-8) {
                x = V.getMatrix(0, N-1, i, i);
                x = x.times(1.0 / x.norm1());
                System.out.println("Stationary distribution using eigenvector:");
                x.print(9, 6);
            }
        }
LW001
  • 2,452
  • 6
  • 27
  • 36
Fascal
  • 1
  • 1

1 Answers1

0

From a quick look the code seems correct, provided that

  1. your matrix transition is the same as in the python example. From your result I assume
  2. the method Matrix.times does a vector multiplication and not a entry-wise multiplication
  3. the number of 50 iterations should be enough for this example. In cases where your 2nd eigenvalue is also close to one this number will be much larger
  4. the function EigenvalueDecomposition.getV returns the eigenvectors in the orientation you expect (row / column)

Note that a valid transition matrix (the matrix of jumping probabilities) has row or column sum of 1 (a matter of definition) since the chances of staying or switching must be 100%. This implies the existance of an eigenvectors of one. From your use of initial transpose of the matrix, I assume you need row-sum of one.

Since you seem to not have an eigenvalue of 1, the most logical error is a typo in your matrix.

jhp
  • 519
  • 3
  • 13