simple.krig = function(x,z,x0,p=1,d=1) {
  if(length(x) != length(z)) stop("input vectors of different length")
  n=length(x)
  sigma=matrix(1,ncol=n,nrow=n)
  for(i in 2:n) {
    for(j in 1:(i-1)) {
      sigma[i,j] = exp( - (abs( (x[i] - x[j])/d ))^p )
      sigma[j,i] = sigma[i,j]
    }
  }
  sigma.inv=solve(sigma)
  m=length(x0)
  pred=numeric(m)
  se=numeric(m)
  for(i in 1:m) {
    sigma0 = exp( - (abs( (x0[i] - x)/d ))^p )
    lambda = sigma.inv%*%sigma0
    pred[i] = t(lambda)%*%z
    se[i] = 1 - t(lambda)%*%sigma0
  }
  return(list(pred=pred,se=se))
}

x = c(1:4)
z = c(-.5,.5,-.5,.5)
x0 = seq(0,5,by=.1)

plot.krig = function(x,z,x0,p=1,d=1) {
  temp.krig = simple.krig(x,z,x0,p,d)
  par(mfcol=c(2,1))
  plot(x0,temp.krig$pred,type="l")
  plot(x0,temp.krig$se,type="l")
  par(mfcol=c(1,1))
  return()
}

plot.krig(x,z,x0,1,1)
plot.krig(x,z,x0,1,0.5)
plot.krig(x,z,x0,1,0.25)
plot.krig(x,z,x0,1,0.1)
plot.krig(x,z,x0,1,0.05)
plot.krig(x,z,x0,1,1)
plot.krig(x,z,x0,1,4)
plot.krig(x,z,x0,1,20)

plot.krig(x,z,x0,1.5,1)


library(fields)
data(ozone)
?ozone

ozone
plot(ozone$x)
plot(ozone$x[,1],ozone$y)
plot(ozone$x[,2],ozone$y)

library(akima)
contour(interp(ozone$x[,1],ozone$x[,2],ozone$y))

fit = Krig(ozone$x,ozone$y, theta=20)
summary(fit)
plot(fit)
surface(fit)
surface(fit,type="C")

par(mfrow=c(1,2))
surface(fit,type="C")
contour(interp(ozone$x[,1],ozone$x[,2],ozone$y))

fit = Krig(ozone$x,ozone$y, theta=10)
surface(fit,type="C")
contour(interp(ozone$x[,1],ozone$x[,2],ozone$y))

par(mfrow=c(1,1))