https://github.com/cran/fields
Raw File
Tip revision: ce722edae3c1b9e1af2985ce3500b11058facf0e authored by Doug Nychka on 24 August 2006, 01:46:17 UTC
version 3.04
Tip revision: ce722ed
Krig.updateY.R
"Krig.updateY" <-
function(out, Y, verbose = FALSE, yM=NA )
{
	#given new Y values but keeping everything else the same finds the 
	#new u vector and pure error SS associated with the Kriging estimate
	# the steps are
	# 1) standardize if neccesary
	# 2) find means, in the case of replicates
	# 3) based on the decomposition, multiply a weighted version of yM
	#    with a large matrix extracted from teh Krig object out. 
	#
	# The out object may be large. This function is written so that out is # #not changed with the hope that it is not copied locally  in this  #function .
	# All of the output is accumulated in the list out2
	#STEP 1
	#
	# transform Y by mean and sd if needed
	#
	if(out$correlation.model) {
		Y <- (Y - predict(out$mean.obj, out$x))/predict(out$sd.obj,
			out$x)
		if(verbose)
			print(Y)
	}
	#
	#STEP 2
        if( is.na( yM[1])){
          out2 <- Krig.ynew(out, Y)}
        else{
          out2<- list( yM= yM, 
          shat.rep = NA,shat.pure.error = NA, pure.ss = NA)}

	if(verbose) {
		print(out2)
	}
	#
	#STEP3
	#
	# Note how matrices are grabbed from the Krig object
	#
	if(verbose) cat("Type of decomposition", out$decomp, fill = TRUE)

	if(out$decomp == "DR") {
		#
		#
		u <- t(out$matrices$G) %*% t(out$matrices$X) %*% (out$weightsM *
			out2$yM)
		#
		# find the pure error sums of sqaures.  
		#
		temp <- out$matrices$X %*% out$matrices$G %*% u
		temp <- sum(out$weightsM * (out2$yM - temp)^2)
		out2$pure.ss <- temp + out2$pure.ss
		if(verbose) {
			cat("pure.ss", fill = TRUE)
			print(temp)
			print(out2$pure.ss)
		}
	}
	#####
	##### end DR decomposition block 
	#####
	####
	#### begin WBW decomposition block
	####
	if(out$decomp == "WBW") {
		#### decomposition of Q2TKQ2
		u <- c(rep(0, out$nt), t(out$matrices$V) %*% qr.q2ty(out$
			matrices$qr.T, sqrt(out$weightsM) * out2$yM))
		if(verbose)
			cat("u", u, fill = TRUE)
		#
		# pure error in this case from 1way ANOVA 
		#
		if(verbose) {
			cat("pure.ss", fill = TRUE)
			print(out2$pure.ss)
		}
	}
	#####
	##### end WBW block
	#####
	out2$u <- u
	out2
}
back to top