# fmc.r

# functions for messing with (discrete time) finite state Markov chains

# demo at end of file



# Darren J Wilkinson
# d.j.wilkinson@ncl.ac.uk
# http://www.staff.ncl.ac.uk/d.j.wilkinson/

# Last updated: 18/9/02

# first the useful fuctions:

"pfmc" <-
function(P,check=TRUE){

# function to return a probability vector representing a stationary
# distribution of the Markov chain with transition matrix P

	if (check == TRUE) {
		stochcheck(P)
	}
        e_eigen(t(P))$vectors[,1]
        e/sum(e)
}

"pvcheck" <-
function(v,tol=1e-15)
{

# function to check if the given vector v is a probability vector

	if (min(v) < 0) {
		stop("Non-pv: neg elts")
	}
	if (max(v) > 1) {
		stop("Non-pv: elts > 1")
	}
	zero_sum(v)-1
	if (zero*zero > tol) {
		stop("Non-pv: non-unit sum")
	}
}

"rdunif" <-
function(p,check=TRUE)
{

# function to simulate a single discrete random quantity with
# given probability vector p

	if (check == TRUE) {
		pvcheck(p)
	}
        cp_cumsum(p)
        length(cp[cp<runif(1)])+1
}

"rfmc" <-
function(P,x0=1,n=100,check=TRUE)
{

# function to simulate a single sample path for a Markov
# chain with given transition matrix P

	if (check == TRUE) {
		stochcheck(P)
	}
        x_x0
        v_vector("numeric",n)
        for (i in 1:n) {
                x_rdunif(P[x,],check=FALSE)
                v[i]_x
        }
        ts(v)
}

"stochcheck" <-
function(P,tol=1e-12)
{

# function to check if the given P is a stochastic matrix

	d_dim(P)
	if ( d[1] != d[2] ) {
		stop("Non-stoch: not square")
	}
	if ( min(P) < 0 ) {
		stop("Non-stoch: neg elts")
	}
	if ( max(P) > 1 ) {
		stop("Non-stoch: elts > 1")
	}
	zero_(P %*% rep(1,d[1]) - rep(1,d[1]) )
	if ( sum(zero*zero) > tol ) {
		stop("Non-stoch: non-unit row sums")
	}
}




# now a demo of how the above functions can be used...


# first a function to build an example transition matrix:

"buildP" <-
function(n=10000)
{
        P_matrix(0,ncol=20,nrow=20)
        P[row(P)==col(P)]_0.7
        P[row(P)==col(P)-1]_0.2
        P[row(P)==col(P)+1]_0.1
        P[1,1]_0.8
        P[20,20]_0.9
	P
}

# now a function to use the functions and plot the results

"mcdemo" <-
function(Pmat=P,n=1000)
{
        op_par(mfrow=c(3,1))
        sample_rfmc(Pmat,n=n)
        plot(sample,main="Sample trace")
        barplot(tabulate(sample)/n,main="Empirical distribution")
        barplot(pfmc(Pmat),main="Theoretical stationary distribution")
        par(op)
}

# run the demo

mcdemo(buildP())



# that's all...