rm(list=ls())
n.points <- 4
#
# das Kovarianzmodell
#
covmod <- function(h){
  covmod <- exp(-abs(h)/5)
}
#
plot(1,1,type="n",xlim=c(0,10),ylim=c(0,10),xlab="x",ylab="y")
# Koordinaten des zu schätzenden Punktes:
y.coords <- c(4,5)
points(t(y.coords),cex=2.5,pch=20,col="red")
#
# Anklicken der Punkte des "Messgitters"
#
p.x <- 1:n.points
p.y <- 1:n.points
for (i in 1:n.points){
  p <- locator(1)
  p.x[i] <- p$x
  p.y[i] <- p$y
  points(p,cex=2.5,pch=20,col="darkgreen")
}
p.coords <- matrix(c(p.x,p.y),nrow=2,byrow=T)
#
# Aufbau der Krigematrix und der rechten Seite des LGS
#
D <- matrix(1:n.points^2,ncol=n.points)
rs <- 1:n.points
for(i in 1:n.points) {
  for(j in 1:n.points) {
    dist <- p.coords[,i]-p.coords[,j]
    D[i,j] <- covmod(sqrt(dist%*%dist))
  }
  dist <- y.coords-p.coords[,i]
  rs[i] <- covmod(sqrt(dist%*%dist))
}
#
# Berechnung der Wichtungsfaktoren
#
D.inv <- solve(D)
lamda <- D.inv%*%rs
cat("\nSimple Kriging\nGröße des Messgitters: ",n.points,"\n")
cat("\nWichtungsfaktoren: \n",lamda,"\n")
cat("\nSumme der Wichtungsfaktoren: \n",sum(lamda),"\n")
# Übrigens: bis hierher wurden noch keine Messwerte benötigt!!
# Einbeziehen von Messwerten:
T.mean = 19
T.long <- c(18,20,17,13,18,23,19,22,16,25) 
# wenn n.points vergrößert wird, müssen hier mehr Werte stehen
T <- T.long[1:n.points]
#
# Berechnung des Schätzwertes:
#
T.y <- (T-T.mean)%*%lamda + T.mean
cat("\ngeschätzter Wert: ",T.y,"\n")
#
# Berechnung der Schätzvarianz
#
est.var <- covmod(0) - as.vector(lamda)%*%as.vector(rs)
cat("\nSchätzvarianz: ",as.real(est.var),"\n")
