explode = ts(scan("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/eqexp.dat")) explode = ts(scan("eqexp.dat")) data.combined = ts(matrix(explode, nrow = length(explode)/17, ncol = 17), frequency = 40) colnames(data.combined) = c(paste("eq", 1:8, sep=""),paste("ex", 1:8, sep=""), "NZ") ## Seventeen columns - 8 earthquakes, 8 explosions, 1 event of unknown origin ## Each column contains the 'P' and 'S' signals ## Split each column into these two signals: data.P = ts(data.combined[1:(nrow(data.combined)/2), ], frequency = 40) data.S = ts(data.combined[(nrow(data.combined)/2+1):(nrow(data.combined)), ], frequency = 40) ############ Set some options: ############ L = 5 scaling = T if(scaling == T) { for (j in 1:ncol(data.P)) data.P[,j] = data.P[,j]/(max(data.P[,j])-min(data.P[,j])) for (j in 1:ncol(data.S)) data.S[,j] = data.S[,j]/(max(data.S[,j])-min(data.S[,j])) } ######################################## data.eq.P = data.P[,1:8] data.ex.P = data.P[,9:16] data.NZ.P = data.P[,17] data.eq.S = data.S[,1:8] data.ex.S = data.S[,9:16] data.NZ.S = data.S[,17] ############ Frequency domain discrimination analysis ################### #Compute all 17 spectral density matrices data = cbind(data.NZ.P, data.NZ.S) # SPEC[i,j,k] is the spectrum between the i-th and j-th series at frequency 2k/n': SPEC = array(dim = c(ncol(data),ncol(data),nextn(nrow(data))/2)) for(i in 1:ncol(data)) { for (j in i:ncol(data)) { power = spec.pgram(data[,c(i,j)], kernel("daniell",(L-1)/2), plot = F) SPEC[i,i, ] = power$spec[,1] SPEC[j,j, ] = power$spec[,2] coh.ij = power$coh phase.ij = power$phase SPEC[i,j, ] = sqrt(coh.ij*power$spec[,1]*power$spec[,2])*exp(1i*phase.ij) SPEC[j,i, ] = Conj(SPEC[i,j, ]) }} SPEC.NZ = SPEC SPEC.EQ = array(dim = c(8, dim(SPEC))) SPEC.EX = array(dim = c(8, dim(SPEC))) for(s in 1:8) { # First get the spectra from each training series of earthquakes: data = cbind(data.eq.P[,s], data.eq.S[,s]) # SPEC[i,j,k] is the spectrum between the i-th and j-th series at frequency 2k/n': SPEC = array(dim = c(ncol(data),ncol(data),nextn(nrow(data))/2)) for(i in 1:ncol(data)) { for (j in i:ncol(data)) { power = spec.pgram(data[,c(i,j)], kernel("daniell",(L-1)/2), plot = F) SPEC[i,i, ] = power$spec[,1] SPEC[j,j, ] = power$spec[,2] coh.ij = power$coh phase.ij = power$phase SPEC[i,j, ] = sqrt(coh.ij*power$spec[,1]*power$spec[,2])*exp(1i*phase.ij) SPEC[j,i, ] = Conj(SPEC[i,j, ]) }} SPEC.EQ[s,,,] = SPEC # Now get the spectra from each training series of explosions: data = cbind(data.ex.P[,s], data.ex.S[,s]) # SPEC[i,j,k] is the spectrum between the i-th and j-th series at frequency 2k/n': SPEC = array(dim = c(ncol(data),ncol(data),nextn(nrow(data))/2)) for(i in 1:ncol(data)) { for (j in i:ncol(data)) { power = spec.pgram(data[,c(i,j)], kernel("daniell",(L-1)/2), plot = F) SPEC[i,i, ] = power$spec[,1] SPEC[j,j, ] = power$spec[,2] coh.ij = power$coh phase.ij = power$phase SPEC[i,j, ] = sqrt(coh.ij*power$spec[,1]*power$spec[,2])*exp(1i*phase.ij) SPEC[j,i, ] = Conj(SPEC[i,j, ]) }} SPEC.EX[s,,,] = SPEC } f.NZ = SPEC.NZ f.eq = apply(SPEC.EQ, c(2:4), mean) f.ex = apply(SPEC.EX, c(2:4), mean) I1.all = vector(length = dim(f.eq)[3]) for (k in 1:length(I1.all)) { tr = sum(diag(solve(f.eq[,,k],f.NZ[,,k]))) deter = prod(Re(eigen(f.eq[,,k], only.values = TRUE)$values)) I1.all[k] = tr + log(deter) } I1 = sum(I1.all) I.eq = Re(I1)/nrow(data) I2.all = vector(length = dim(f.ex)[3]) for (k in 1:length(I2.all)) { tr = sum(diag(solve(f.ex[,,k],f.NZ[,,k]))) deter = prod(Re(eigen(f.ex[,,k], only.values = TRUE)$values)) I2.all[k] = tr + log(deter) } I2 = sum(I2.all) I.ex = Re(I2)/1024 cat("KL discrepancies between f = f.NZ and the spectra based on the training samples are", "\n") print(cbind(I.eq, I.ex)) cat("I.eq - I.ex =", I.eq - I.ex, "\n") ############ Frequency domain cross-validation ################### I.eq.minus.I.ex = vector(length=16) for(m in 1:17) { # Compute the spectra from each holdout series: if(m <= 8) data = cbind(data.eq.P[,m], data.eq.S[,m]) if(m > 8 & m <=16) data = cbind(data.ex.P[,m-8], data.ex.S[,m-8]) if(m == 17) data = cbind(data.NZ.P, data.NZ.S) # SPEC[i,j,k] is the spectrum between the i-th and j-th series at frequency 2k/n': SPEC = array(dim = c(ncol(data),ncol(data),nextn(nrow(data))/2)) for(i in 1:ncol(data)) { for (j in i:ncol(data)) { power = spec.pgram(data[,c(i,j)], kernel("daniell",(L-1)/2), plot = F) SPEC[i,i, ] = power$spec[,1] SPEC[j,j, ] = power$spec[,2] coh.ij = power$coh phase.ij = power$phase SPEC[i,j, ] = sqrt(coh.ij*power$spec[,1]*power$spec[,2])*exp(1i*phase.ij) SPEC[j,i, ] = Conj(SPEC[i,j, ]) }} SPEC.HOLDOUT = SPEC SPEC.EQ = array(dim = c(length(setdiff(1:8,m)), dim(SPEC))) SPEC.EX = array(dim = c(length(setdiff(1:8,m-8)), dim(SPEC))) # Compute the spectra from each remaining training series of earthquakes: for(ss in 1:length(setdiff(1:8,m))) { s = setdiff(1:8,m)[ss] data = cbind(data.eq.P[,s], data.eq.S[,s]) # SPEC[i,j,k] is the spectrum between the i-th and j-th series at frequency 2k/n': SPEC = array(dim = c(ncol(data),ncol(data),nextn(nrow(data))/2)) for(i in 1:ncol(data)) { for (j in i:ncol(data)) { power = spec.pgram(data[,c(i,j)], kernel("daniell",(L-1)/2), plot = F) SPEC[i,i, ] = power$spec[,1] SPEC[j,j, ] = power$spec[,2] coh.ij = power$coh phase.ij = power$phase SPEC[i,j, ] = sqrt(coh.ij*power$spec[,1]*power$spec[,2])*exp(1i*phase.ij) SPEC[j,i, ] = Conj(SPEC[i,j, ]) }} SPEC.EQ[ss,,,] = SPEC } # Compute the spectra from each remaining training series of explosions: for(ss in 1:length(setdiff(1:8,m-8))) { s = setdiff(1:8,m-8)[ss] data = cbind(data.ex.P[,s], data.ex.S[,s]) # SPEC[i,j,k] is the spectrum between the i-th and j-th series at frequency 2k/n': SPEC = array(dim = c(ncol(data),ncol(data),nextn(nrow(data))/2)) for(i in 1:ncol(data)) { for (j in i:ncol(data)) { power = spec.pgram(data[,c(i,j)], kernel("daniell",(L-1)/2), plot = F) SPEC[i,i, ] = power$spec[,1] SPEC[j,j, ] = power$spec[,2] coh.ij = power$coh phase.ij = power$phase SPEC[i,j, ] = sqrt(coh.ij*power$spec[,1]*power$spec[,2])*exp(1i*phase.ij) SPEC[j,i, ] = Conj(SPEC[i,j, ]) }} SPEC.EX[ss,,,] = SPEC } f.holdout = SPEC.HOLDOUT f.eq = apply(SPEC.EQ, c(2:4), mean) f.ex = apply(SPEC.EX, c(2:4), mean) I1.all = vector(length = dim(f.eq)[3]) for (k in 1:length(I1.all)) { tr = sum(diag(solve(f.eq[,,k],f.holdout[,,k]))) deter = prod(Re(eigen(f.eq[,,k], only.values = TRUE)$values)) I1.all[k] = tr + log(deter) } I1 = sum(I1.all) I.eq = Re(I1) I2.all = vector(length = dim(f.ex)[3]) for (k in 1:length(I2.all)) { tr = sum(diag(solve(f.ex[,,k],f.holdout[,,k]))) deter = prod(Re(eigen(f.ex[,,k], only.values = TRUE)$values)) I2.all[k] = tr + log(deter) } I2 = sum(I2.all) I.ex = Re(I2) I.eq.minus.I.ex[m] = (I.eq - I.ex)/1024 if (m <=8) cat("m = ", m, "I.ex.minus.I.eq =", -I.eq.minus.I.ex[m], "\n") if (m>8) cat("m = ", m, "I.eq.minus.I.ex =", I.eq.minus.I.ex[m], "\n") } mat = cbind(1:8, -I.eq.minus.I.ex[1:8], 1:8, I.eq.minus.I.ex[9:16]) colnames(mat) = c("EQ", "I.ex-I.eq", "EX", "I.eq-I.ex") cat("Differences in KL discrepancies; negative values indicate misclassifications","\n") print(round(mat,2))