rm(list=ls())
#
# das Kovarianz- und Variogrammmodell (exponentielles Modell)
#
cov.mod <- function(x){        # Kovarianzmodell
  cov.mod <- exp(-abs(x)/5)
}
gamma.mod <- function(x) {     # Variogrammmodell
  gamma.mod <- cov.mod(0) - cov.mod(x)
}
#
# Koordinaten des zu schätzenden Punktes:
#
y.coords <- c(3,4)
#
# Punkte des "Messgitters"
#
p.x <- c(2,4,0)
p.y <- c(6,4,2)
p.coords <- matrix(c(p.x,p.y),nrow=2,byrow=T)
#
# Aufbau der Krigematrix und der rechten Seite des LGS
#
D <- matrix(1:16,ncol=4)
dist <- p.coords[,1]-p.coords[,1]; D[1,1] <- gamma.mod(sqrt(dist%*%dist))
dist <- p.coords[,1]-p.coords[,2]; D[1,2] <- gamma.mod(sqrt(dist%*%dist))
dist <- p.coords[,1]-p.coords[,3]; D[1,3] <- gamma.mod(sqrt(dist%*%dist))
dist <- p.coords[,2]-p.coords[,1]; D[2,1] <- gamma.mod(sqrt(dist%*%dist))
dist <- p.coords[,2]-p.coords[,2]; D[2,2] <- gamma.mod(sqrt(dist%*%dist))
dist <- p.coords[,2]-p.coords[,3]; D[2,3] <- gamma.mod(sqrt(dist%*%dist))
dist <- p.coords[,3]-p.coords[,1]; D[3,1] <- gamma.mod(sqrt(dist%*%dist))
dist <- p.coords[,3]-p.coords[,2]; D[3,2] <- gamma.mod(sqrt(dist%*%dist))
dist <- p.coords[,3]-p.coords[,3]; D[3,3] <- gamma.mod(sqrt(dist%*%dist))
D[4,] <- 1; D[,4] <- 1; D[4,4] <- 0
#
# Die rechte Seite des Gleichungssystems
#
rs <- 1:4
dist <- p.coords[,1]-y.coords; rs[1] <- gamma.mod(sqrt(dist%*%dist))
dist <- p.coords[,2]-y.coords; rs[2] <- gamma.mod(sqrt(dist%*%dist))
dist <- p.coords[,3]-y.coords; rs[3] <- gamma.mod(sqrt(dist%*%dist))
rs[4] <- 1
#
# Berechnung der Wichtungsfaktoren
#
D.inv <- solve(D)
lambda <- D.inv%*%rs
cat("\n\n  Ordinary Kriging: Lösung des LGS:   \n")
cat("--------------------------------------\n")
cat("\n Wichtungsfaktoren: \n")
print(as.matrix(lambda[1:3]))
cat("\n Lagrange multiplier: \n")
print(lambda[4])
cat("\n Summe der Wichtungsfaktoren: \n")
print(round(sum(lambda[1:3]),3)) 
#
# Berechnung des Schätzwertes und der Schätzvarianz:
#
T <- c(18,20,17) # wenn n.points vergrößert wird, müssen hier mehr Werte stehen
T.val <- as.real(T%*%lambda[1:3])
cat("\n Schätzwert: \n")
print(T.val)
#
T.var <- as.real(rs%*%lambda)
cat("\n Schätzvarianz: \n")
print(T.var)
#
cat("\n++++ Ende Ordinary Kriging  ++++\n")