#!/usr/bin/python -Wall # ================================================================ # Some Markov-chain hackwork. # John Kerl # kerl.john.r@gmail.com # 2007-07-09 # ================================================================ # This software is released under the terms of the GNU GPL. # Please see LICENSE.txt in the same directory as this file. # Copyright (c) John Kerl 2007 # ================================================================ from __future__ import division # 1/2 = 0.5, not 0. import sackmat_m import random import copy import sys import pmtc_tm # ---------------------------------------------------------------- # Scales the vector so its row sum is 1. def row_normalize_vector(v): sum = sackmat_m.vec_contract(v) if (sum == 0.0): return n = len(v) for i in range(0, n): v[i] /= sum # ---------------------------------------------------------------- # Scales each row of the matrix so the row sums are 1. def row_normalize_matrix(A): [nr, nc] = A.dims() for i in range(0, nr): row_normalize_vector(A[i]) # ---------------------------------------------------------------- # Use the metric induced by the max norm. def are_close(A, B, tol = 1e-4): # Assume same dims n = A.square_dim() for i in range(0, n): for j in range(0, n): d = abs(A[i][j] - B[i][j]) if (d > tol): return [0, d] return [1, tol] # ---------------------------------------------------------------- # Iterates the chain, looking for stability. def is_stable(A, verbose=0): Ap = copy.copy(A) #maxits = 20 maxits = 1000 k = 0 while (k < maxits): Ao = Ap Ap = Ap * A [yn, err] = are_close(Ao, Ap) if (verbose): print "-- k = %d err = %.4e" % (k, err) Ap.printf() print "" if (yn == 1): return [1, Ap] k += 1 return [0, Ap] # ---------------------------------------------------------------- # Let q = 1-p. # p q 0 0 0 0 # 0 p q 0 0 0 # 0 0 p q 0 0 # 0 0 0 p q 0 # 0 0 0 0 p q # q 0 0 0 0 p def circpmat1(n, p): A = sackmat_m.make_zero_matrix(n, n) for i in range(0, n): A[i][i] = p A[i][(i+1)%n] = 1-p return A # ---------------------------------------------------------------- # Let q = 1-p. # 0 q 0 0 0 p # p 0 q 0 0 0 # 0 p 0 q 0 0 # 0 0 p 0 q 0 # 0 0 0 p 0 q # q 0 0 0 p 0 def circpmat2(n, p): A = sackmat_m.make_zero_matrix(n, n) for i in range(0, n): A[i][(i-1)%n] = p A[i][(i+1)%n] = 1-p return A # ---------------------------------------------------------------- # Assign matrix elements using uniform distribution on [0, 1), then # normalize rows. def randprbmat(n): A = sackmat_m.make_zero_matrix(n, n) for i in range(0, n): for j in range(0, n): A[i][j] = random.random() row_normalize_matrix(A) return A # ---------------------------------------------------------------- def pqmat(p, q): return sackmat_m.sackmat([[p, 1-p],[1-q, q]]) # ================================================================ #A = randprbmat(6,6) #A = pmtc_tm.from_cycles([[1,2],[3,4,5]],5).to_permutation_matrix() #A = circpmat1(6, .5) #A = circpmat2(6, .1) p = 0.9 q = 0.8 if (len(sys.argv) == 3): p=float(sys.argv[1]) q=float(sys.argv[2]) A=pqmat(p, q) print "A:" A.printf() print "" #sys.exit(0) [yn, As] = is_stable(A, 1) if (yn == 1): print "stable" As.printf() else: print "unstable" As.printf()