loop=function(y, coords, cov.fun, cov.par, means=NULL){ n=length(y) if(!n==nrow(coords))stop("The length of y and the number of locations do not match") vars=lapply(list(x=0), FUN=cov.fun, para=cov.par)[[1]] if(length(means)==1) means=rep(means, n) result=matrix(nrow=n, ncol=4) colnames(result)=c("Observed", "Predict", "SD", "Z") for(i in 1:n){ d=as.matrix(dist(coords[-i,])) V=lapply(list(d), FUN=cov.fun, para=cov.par)[[1]] V=solve(V) d=dist(rbind(coords[i,], coords[-i,]))[1:(n-1)] k=lapply(list(d), FUN=cov.fun, para=cov.par)[[1]] m=ifelse(is.null(means),(rep(1, n-1)%*%V%*%y[-i])/sum(V), means[i]) result[i,2]=m+k%*%V%*%(y[-i]-m) result[i,3]=(vars-k%*%V%*%k+(1-sum(V%*%k))^2/(sum(V)))^(1/2) # see page 17 in my note } result[,1]=y result[,4]=(y-result[,2])/result[,3] Z=result[,4] CRPS=mean(result[,3]*(Z*(2*pnorm(Z)-1)+2*pnorm(Z)-(1/pi)^(1/2)) ) LogS=0.5*log(2*pi)+mean(log(result[,3])+0.5*Z^2) return(list(output=result, scores=c(Z=mean(result[,4]), Z2=mean(result[,4]^2), RMSE=sqrt(mean((result[,1]-result[,2])^2)), CRPS=CRPS, LogS=LogS))) } cov.exp=function(x, para){ result=para[2]*exp(-x/para[3]) result[x==0]=result[x==0]+para[1] result } likfit1=likfit(coords = paws[, c("X", "Y")], data = paws\$Soil8, ini.cov.pars=c(13,4.4), nugget=0.3) loop(paws\$Soil8, coords=paws[,c("X", "Y")], cov.fun="cov.exp", cov.par=c(likfit1\$nugget, likfit1\$cov.par)) paws.cv=xvalid(data=paws\$Soil8, coords=paws[,c("X", "Y")],model=likfit1) #To get RMSE sqrt(mean((paws.cv\$error)^2)) # CRPS Z=paws.cv\$std.error mean((paws.cv\$krige.var)^0.5*((Z*(2*pnorm(Z)-1)+2*pnorm(Z)-(1/pi)^(1/2)))) #LogS 0.5*log(2*pi)+mean(log(sqrt(paws.cv\$krige.var))+0.5*Z^2)