Semivariograma

datos=read.table("clipboard")
datos

##   V1 V2 V3 V4 V5 V6 V7 V8 V9
## 1 44 NA 40 42 40 39 37 36 NA
## 2 42 NA 43 42 39 39 41 40 38
## 3 37 37 37 35 38 37 37 33 34
## 4 35 38 NA 35 37 36 36 35 NA
## 5 36 35 36 35 34 33 32 29 28
## 6 38 37 35 NA 30 NA 29 30 32

datos[ ,1]

## [1] 44 42 37 35 36 38

diff(datos[,1],lag=1)^2

## [1]  4 25  4  1  4

apply(datos,1,diff)^2

##    [,1] [,2] [,3] [,4] [,5] [,6]
## V2   NA   NA    0    9    1    1
## V3   NA   NA    0   NA    1    4
## V4    4    1    4   NA    1   NA
## V5    4    9    9    4    1   NA
## V6    1    0    1    1    1   NA
## V7    4    4    0    0    1   NA
## V8    1    1   16    1    9    1
## V9   NA    4    1   NA    1    4

apply(datos,2,diff)^2

##      V1 V2 V3 V4 V5 V6 V7 V8 V9
## [1,]  4 NA  9  0  1  0 16 16 NA
## [2,] 25 NA 36 49  1  4 16 49 16
## [3,]  4  1 NA  0  1  1  1  4 NA
## [4,]  1  9 NA  0  9  9 16 36 NA
## [5,]  4  4  1 NA 16 NA  9  1 16

c(apply(datos,1,diff)^2,apply(datos,2,diff)^2)

##  [1] NA NA  4  4  1  4  1 NA NA NA  1  9  0  4  1  4  0  0  4  9  1  0 16  1  9
## [26] NA NA  4  1  0  1 NA  1  1  1  1  1  1  9  1  1  4 NA NA NA NA  1  4  4 25
## [51]  4  1  4 NA NA  1  9  4  9 36 NA NA  1  0 49  0  0 NA  1  1  1  9 16  0  4
## [76]  1  9 NA 16 16  1 16  9 16 49  4 36  1 NA 16 NA NA 16

diferencias=c(apply(datos,1,diff, lag=1)^2,apply(datos,2,diff, lag=1)^2)
diferencias=na.omit(diferencias)             
sum(diferencias)/(2*length(diferencias))

## [1] 3.402778

sum(diferencias)

## [1] 490

gamma=sum(diferencias)/(2*length(diferencias))
gamma

## [1] 3.402778

##funcion semivariograma
calc_semivariograma=function(lag){
diferencias=c(apply(datos,1,diff,lag=lag)^2,apply(datos,2,diff, lag=lag)^2)
diferencias=na.omit(diferencias)             
sum(diferencias)/(2*length(diferencias))
sum(diferencias)
gamma=sum(diferencias)/(2*length(diferencias))
return(gamma)
}
calc_semivariograma(lag=1)

## [1] 3.402778

calc_semivariograma(lag=2)

## [1] 6.258333

h=1:7
sapply(h,calc_semivariograma)

## [1]  3.402778  6.258333 10.687500 14.222222 12.795455 13.035714 15.555556

gamma=sapply(h,calc_semivariograma)
plot(h*100,gamma,pch=16, main="Semivariograma")