# R code to calculate the SPECTRAL ENVELOPE for a categorical time series # *** Must source mvspec.R first *** # The data set is a column of INTEGERS representing the categories u = factor(scan("/mydata/bnrf1ebv.dat")) # first, input the data as factors and then x = model.matrix(~u-1)[,1:3] # make an indicator matrix #x = x[1:1000,] # select subsequence if desired Var = var(x) # var-cov matrix xspec = mvspec(x, spans=c(7,7)) # spectral matrices are an array called fxx fxxr = Re(xspec\$fxx) # fxxr is real(fxx) # #----- compute Q = Var^-1/2 -----# ev = eigen(Var) Q = ev\$vectors%*%diag(1/sqrt(ev\$values))%*%t(ev\$vectors) # #--- compute spec env and scale vectors ---# num = xspec\$n.used # sample size used for FFT nfreq = length(xspec\$freq) # number of freqs used specenv = matrix(0,nfreq,1) # initialize the spec envelope beta=matrix(0,nfreq,3) # initialize the scale vectors for (k in 1:nfreq){ ev = eigen(2*Q%*%fxxr[,,k]%*%Q/num) # get evalues of normalized spectral matrix at freq k/n specenv[k] = ev\$values[1] # spec env at freq k/n is max evalue b = Q%*%ev\$vectors[,1] # beta at freq k/n beta[k,] = b/sqrt(sum(b^2)) # helps to normalize beta } # #--- output and graphics ---# frequency = (0:(nfreq-1))/num plot(frequency, 100*specenv, type="l", ylab="Spectral Envelope (%)") ## add significance threshold to plot ## m=xspec\$kernel\$m etainv=sqrt(sum(xspec\$kernel[-m:m]^2)) thresh=100*(2/num)*exp(qnorm(.9999)*etainv)*matrix(1,nfreq,1) lines(frequency,thresh, lty="dashed", col="blue") #-- details --# output = cbind(frequency, specenv, beta) colnames(output)=c("freq","specenv","A", "C", "G") # output # uncomment to write to screen # write.table(output, file="output.txt", quote=FALSE, row.names=FALSE) # uncomment to write to file