/********************************************************************
 *  simulation of a Markov chain
 *  @author U. Pantle, 17.05.06
 *******************************************************************/

import java.util.Random;
import java.text.DecimalFormat;

public class Prog2_2b{


    /* output formate */
    static DecimalFormat df = new DecimalFormat ("#0.000");
    static public String out(double y){
	return df.format(y);
    }

/*******************************************************************
 * Simulation by recursion formula; see lecture notes formulas (17)-(19) 
 * 
 * @return state of MC after transition
 *
 * @param x state before transition
 * @param p transition probabilities
 * @param z auxiliary random variable
 ********************************************************************/
    private static int iterate(int i, double[][] p, double z){
	int j;   // state after transition

	if(z <= p[i][0]+p[i][1]){
	    if(z <= p[i][0])
		j = 0;	
	    else
		j=1;
	}
	else{
	   j = 2;
	}
	return j;
    }

/************************ main ******************************/
    public static void main(String args[]) {

	/* initial distribution */
	double[] a = new double[]{1.0/3.0, 1.0/3.0, 1.0/3.0};

	/* transition probabilities */
	double[][] p = new double[3][];   
	p[0] = new double[]{ 0.8, 0.2, 0.0};
	p[1] = new double[]{ 0.4, 0.4, 0.2};
	p[2] = new double[]{ 0.2, 0.6, 0.2};
   
	int n = Integer.parseInt(args[0]);  // number of transitions
	int[] state = new int[n+1];

	Random uni = new Random();

	/* do it several times... just for fun */
	double[] freq = new double[3];
	int sample = Integer.parseInt(args[1]);
	for(int m=0; m< sample; m++){
	    /* generate initial state */
	    double z0 = uni.nextDouble();
	    
	    if(z0 <= a[0]+a[1]){
		if(z0 <= a[0])
		    state[0] = 0;	
		else
		    state[0] = 1;
	    }
	    else{
		state[0] = 2;
	    }
	    
	    /* simulation by recursive representation */
	    for (int k = 1; k <= n; k++){    // ... n transitions
		state[k] = iterate(state[k-1], p, uni.nextDouble());
	    }
	    freq[state[n-1]]++;
	}    

	for(int m=0; m<3; m++){
	    freq[m] /= (double)sample;
	}

	System.out.print("Relative frequencies of final state after n = "+n+" transitions: ("+out(freq[0])+", "+out(freq[1])+", "+out(freq[2])+")\n");
	}
	
}