#   SYMBOLIC COMPUTATION FOR STATISTICAL INFERENCE IN R
#   D.F. Andrews, University of Toronto
#   2001
#
# This file presents examples of the use of Symbolic Computation
# for Statistical Inference.  A full discussion of the theory
# on which the computations are based is found in the book by 
# Jamie Stafford and the author of the same title.
# The file contains examples of computing moments, cumulants,
# unbiased estimates of cumulants, expansions of mle's and
# mean log likelihood rations and their expectations to order n^(-2).
#
# On a Mac cube, execution of this file takes 100 seconds 
# 75 of which are for the last example.
#
# The basic procedures are contained in the file SCSI.R.
#
source("SCSI.R")
#
# OBJECTS, LISTS and AVERAGES
#
#The objects of statistics have forms such as averages,
#moments and cumulants.  All objects are represented here as
#lists. 

#For example, here we make an object representing the product of 3
#averages and print it. The function AV is used to cast the list
#as a product of averages denoted by A.  Here we create A(x)^3 
#
l <- list(list(list("X")),list(list("X")),list(list("X")))
axayaz <- AV(l); axayaz
#
# CALCULATIONS as TRANSFORMATIONS of LISTS;
# EXPECTATION of PRODUCTS of AVERAGES
#
#Objects of one form may be transformed to another form.
#The function for computing the transformed lists to type X
#from type Y is denoted by TRXY.  Such transformations may be used, 
#for example, to compute the expected value of a product of averages.
#
#This computation introduces another form: disjoint averages
#denoted by D. This represent objects of the form 
#sum(X_i Y_j)/n^2 where i !=j in the sum.
#Moments (expectations) are denoted by M.
# TRDA transforms averages to disjoint averages.
# TRMD transforms disjoint averages to moments by
# computing their expected value.
# These transformations can be combined.
#
# Compute the expected value of the product of averages.
#
m3 <- TRMD(TRDA(axayaz)); m3
#
# MOMENTS <-> CUMULANTS
#
#CUM creates a cumulant from a list.
#MOM creates a moment.
#
#Make an object representing the covariance K(X,Y), a 2nd cumulant.
#
k <- CUM(list(list(list("X"),list("Y")))); k
#
#TRMK transforms cumulants to moments.
#TRKM is the inverse transformation.
#
#Express the cumulant in terms of moments.
#
m <- TRMK(k);m
#
#UNBIASED ESTIMATES of CUMULANTS
#
#TRDM transforms moments to disjoint averages which
#are the unbiased estimates of the moments.
#TRAD transforms the disjoint averages to averages.
#These can be combined to compute the unbiased estimate of
#the products of moments and hence of the cumulant.
#
#Compute  unbiased estimate of the cumulant
#in a form that can be computed in R.
#
a <- TRAD(TRDM(m)); a
#
#NUMERICAL EVALUATION OF EXPRESSIONS
#
#EVAL evaluates an expression composed of averages.
#
#Evaluate the unbiased estimate for simulated data.
#
#First we create some simulated data
n <- 20; X <- rnorm(n); Y <- X
# and evaluate the result
EVAL(a)
#check -- the result should be the variance of X.
var(X)
#
#Now we do the same thing for the third cumulant.
#
k3 <- CUM(list(list(list("X"),list("X"),list("X")))); k3
#
#Express it in terms of moments.
#
m3 <- TRMK(k3); m3
#
#Compute its unbiased estimate in a form that can be computed in R.
#
a3 <-TRSIMP(TRSIMP( TRAD(TRDM(m3)))); a3
#
#Note Fisher's calculation is rather less elementary.
#TRSIMP was used to simplify expressions.
#
#Evaluate the expression for the estimate of the cumulant.
#
n <- 20; X <- rnorm(n)
EVAL(a3)
#
#GENERALIZED CUMULANTS; COMPLEMENTARY SET PARTITIONS
#
# Create a generalized cumulant.
#
gk <- CUM(list(list(list("X","Y"),list("Z")))); gk
#
# Express it in terms of simple cumulants.
#
kgk <- TRKM(TRMK(gk)); kgk
#
#Note that this calculation serves to compute the 
#complementary set partition of the list defining the cumulant.
#
#
#SERIES EXPANSIONS
#
#SERAEP creates a series expansion.  
#The length of the series is specified by the first argument.
#The second argument is the object to be expanded.
#Here we expand a simple average.
#The result is expressed in terms of expectations and
#centred averages denoted by Z.
#
l <- list(list(list("X"))); ax <-AV(l); ax
sax <- SERAEP(3,ax); sax
#
#SERadd adds two series
saxplussax <- SERadd(sax,sax); saxplussax
#
#SERtimes computes a product
saxtimessax <- SERtimes(sax,sax); saxtimessax
#
#SERpower computes powers
saxpower4 <- SERpower(sax,4); saxpower4
#
#SERConst creates a series consisting of a constant.
#The length of the series is specified by the second argument.
st4 <- SERConst("t0",4); st4
#
#Much of statistics deals with averages of functions of observations.
#For example, the log-likelihood is n times the average of 
#log(f(x|theta).
#
#The multivariate case is more interesting. Products are denoted by DOT.
#The form the the result is designed for subsequent evaluation.
#It is not pretty.
#
#Here we compute some average functions and roots of equations.
#The derivative of a function is denoted by its second argument.
#
#Suppose f(theta) = psi(theta,x_i), a function of 
# theta with a value that depends on x_i.
# Here we compute the AV(f(AV(X)) i.e. where theta = AV(X)
#First we create a series expansion for AV(X)
#
l <- list(list(list("X"))); ax <- VAV(l);
sax <-SERAEP(2,ax); sax
#
#Sfun expands average functions in power (Taylor) series
#Now we compute the expansion of AV(f(A(X))
fax <- Sfun("f",sax,0,"A"); fax
#
#ROOTS OF ESTIMATING EQUATIONS
#
#MLE and M- estimates are defined as roots of equations.
#Now we compute the root of the equation AV(f(t)) = f(t0)
#The inverse of the expected hessian is denoted by IE.
#
st3 <- SERConst("t0",3)
rootfax <- Sroot("f",st3,0,"A"); rootfax
#
#Check the root
check <- Sfun("f",rootfax,0,"A"); check
#
#The likelihood function, l, is special in that its expected score is 0.
#The MLE is computed.
st3 <- SERConst("t0",3)
rootl13 <- Sroot("l",st3,1,"A"); rootl13
#
#LIKELIHOOD RATIOS
#
#The loglikelihood ratio is computed.
st4 <- SERConst("t0",4); rootl14 <- Sroot("l",st4,1,"A")
lrootl14 <- Sfun("l",rootl14,0,"A"); lrootl14
#
#Now to evaluate for the Gau(mu,sigma) likelihood.  This evaluation
#requires functions to evaluate derivatives etc.  These need
#to be specified for the particular distribution, here gaussian.
#This code is a tad lengthy.
#
#first specify expected values of derivatives of l.
el <- function(x,der=0){
	if(der == 0) {
		e <- -1/2 - log(x[2])
	}
	else{
		if(der == 2){
			e <- matrix(c(-1/(x[2]^2),0,0,-2/(x[2]^2)),2,2)
		}
		else{
			if(der == 3){
				e <-
array(c(0,2/(x[2]^3),2/(x[2]^3),0,2/(x[2]^3),0,0,10/(x[2]^3)),
					dim=c(2,2,2))
			}
			else{
				e <- 0
			}
		}
	}
	return(e)
}
#
#Specify standardized averages of derivatives of l
lx <- function(y,t=t0,der = 0){
	if(der == 0) {
		lxv <- -1/2*((y - t[1])/t[2])^2 - log(t[2]) + 1/2
	}
	else{
		if(der == 1){
			lxv <- array(c((y-t[1])/(t[2]^2),(y-t[1])^2/(t[2])^3 - 1/t[2]),2)
		}
		else{
			if(der == 2){
				lxv <-
matrix(c(0,-2*(y-t[1])/(t[2]^3),-2*(y-t[1])/(t[2]^3),-3*(y-t[1])^2/(t[2]^4)+3/(t[2]^2)),
					2,2)
			}
			else{lxv<-0
			}
		}
	}
	return(lxv)
}
#
zl <- function(t,der=0,dater=200){
	zxl <- (lapply(xdata,lx,der=der,t=t));
	s <- zxl[[1]];
	for(i in 2:length(zxl)){
		s <- s + zxl[[i]]
	}
	s/length(zxl)
}
#
#Specify the data
t0<- c(0,1); nx <- 40; xdata <- rnorm(nx)
#
#Evaluate mle
SEReval(rootl13)
#
#Compare to numerical value
#
t1 <- mean(xdata); t2 <- sqrt(sum((xdata-t1)^2)/length(xdata)); c(t1,t2)
#
#Evaluate mean log likelihood ratio
SEReval(lrootl14)
#
#Compare the series approximation to value computed numerically.
lt12 <- -log(t2) -1/2*(mean((xdata-t1)^2/(t2^2))); lt12
#
#EXPECTED VALUES FOR SERIES: MLE and LIKELIHOOD RATIO
#
#Expected values can be easily computed for reals
#using the methods above.
#To use these methods for the multivariate case
#we must first express the series in Einstein (summation) notation.
#
#Add subscripts to the mle
#
srootl13 <- Ssubscript(rootl13); srootl13
#
#Now compute its expected value
srootl13 <- Ssubscript(rootl13);
er<-SexpectZ(srootl13); er
#
#Compute the expected value of the mean likelihood ratio.
#This takes about 1 minute on a Mac cube.
#(The warnings arise because the first term has no subscripts.)
#
st5 <- SERConst("t0",5); rootl15 <- Sroot("l",st5,1,"A")
lrootl15 <- Sfun("l",rootl15,0,"A"); 
slrootl15 <- Ssubscript(lrootl15);
er<-SexpectZ(slrootl15); er