# R/ldt-analytic.apply.R In spatialCovariance: Computation of Spatial Covariance Matrices for Data on Rectangles

```f.anal.ldt <- function(coords,info)
{
## Copied over from ldt-analytic.r on Dec 10th
## trying to incorporate apply into my code

## created May 12th
## analytic results for the ldt model
## see write up and notes in binder for more details.

a <- info\$rowwidth
b <- info\$colwidth
ax <- coords[1]
bx <- coords[2]
i <- coords[3]
j <- coords[4]
## 5
## 6
## 7 - which rows to evaluate - evalFactor
## 8 - lengths
## 9 - indicator for which rows to evaluate analytic results
evaluate <- coords[9]
## 10 - lower limit for numerical integration
## 11 - upper limit for numerical integration

ret <- 0
if(evaluate) {
if(i !=1 && j !=1)
{
## most common situation:  different row, different column
## simplified different row, column result
## this is bad due to 0/0 and log(0) see below
## also its old, note the double for loops
if(FALSE) {
as <- c(-ax,ax-a,ax,ax+a)
bs <- c(-bx,bx-b,bx,bx+b)
ret <- 100*a^2*b^2
for(ii in 1:4)
for(jj in 1:4)
{
ret <- ret + (-1)^(ii+jj+1)*(8*as[ii]*bs[jj]*(bs[jj]^2*atan(as[ii]/bs[jj])+as[ii]^2*atan(bs[jj]/as[ii]))-(as[ii]^4-6*as[ii]^2*bs[jj]^2+bs[jj]^4)*log(as[ii]^2+bs[jj]^2))
}
ret <- ret/(48*a^2*b^2)
}

## different version, takes care of 0/0, or log(0)
## doesn't have for loops
ret <- 100*a^2*b^2
as <- c(-ax,ax-a,ax,ax+a)
bs <- c(-bx,bx-b,bx,bx+b)
asbs <- cbind(as.vector(matrix(as,4,4)),as.vector(matrix(1:4,4,4)),as.vector(t(matrix(bs,4,4))),as.vector(t(matrix(1:4,4,4))))

ret <- 100*a^2*b^2
ret <- ret + sum(apply(asbs,1,sepRowsepCol.ldt.res))
ret <- ret/(48*a^2*b^2)
} else if(i == 1 && j !=1)
{
## same column, different row
ax <- bx
bx <- 0
a <- b
b <- info\$rowwidth
ret <- (50*a^2*b^2-8*b^3*(ax-a)*atan((ax-a)/b)+
16*ax*b^3*atan(ax/b)-8*b^3*(ax+a)*atan((ax+a)/b)+
16*ax^3*b*atan(b/ax)-8*b*(ax-a)^3*atan(b/(ax-a))-
8*b*(ax+a)^3*atan(b/(ax+a))+4*ax^4*log(ax)-
2*(ax+a)^4*log(ax+a)+
(((ax-a)^2-b^2)^2-4*b^2*(ax-a)^2)*log((ax-a)^2+b^2)-
2*((ax^2-b^2)^2-4*b^2*ax^2)*log(ax^2+b^2)+
(((ax+a)^2-b^2)^2-4*b^2*(ax+a)^2)*log((ax+a)^2+b^2))
if(ax>a) ret <- ret - 2*(ax-a)^4*log(ax-a)
ret <- ret/(24*a^2*b^2)
} else if(i!=1 && j==1)
{
## simplified same row result using mathematica
## problems when ax=a with log(ax-a) dealt with below
ret <- (50*a^2*b^2-8*b^3*(ax-a)*atan((ax-a)/b)+
16*ax*b^3*atan(ax/b)-8*b^3*(ax+a)*atan((ax+a)/b)+
16*ax^3*b*atan(b/ax)-8*b*(ax-a)^3*atan(b/(ax-a))-
8*b*(ax+a)^3*atan(b/(ax+a))+4*ax^4*log(ax)-
2*(ax+a)^4*log(ax+a)+
(((ax-a)^2-b^2)^2-4*b^2*(ax-a)^2)*log((ax-a)^2+b^2)-
2*((ax^2-b^2)^2-4*b^2*ax^2)*log(ax^2+b^2)+
(((ax+a)^2-b^2)^2-4*b^2*(ax+a)^2)*log((ax+a)^2+b^2))
if(ax>a) ret <- ret - 2*(ax-a)^4*log(ax-a)
ret <- ret/(24*a^2*b^2)
} else {
## diagonal result, the most rare
ret <- (25*a^2*b^2-8*a*b^3*atan(a/b)-8*a^3*b*atan(b/a)-2*a^4*log(a)-2*b^4*log(b)+(a^4-6*a^2*b^2+b^4)*log(a^2+b^2))/(12*a^2*b^2)
}
}
ret
}

sepRowsepCol.ldt.res <- function(abVal)
{
aVal <- abVal[1]
ii <- abVal[2]
bVal <- abVal[3]
jj <- abVal[4]

numSteps <- 2-((aVal==0) + (bVal==0))
ret <- 0
if(numSteps>=1) ret <- ret + (-1)^(ii+jj+2)*(aVal^4-6*aVal^2*bVal^2+bVal^4)*log(aVal^2+bVal^2)
if(numSteps>1) ret <- ret + (-1)^(ii+jj+1)*(8*aVal*bVal*(bVal^2*atan(aVal/bVal)+aVal^2*atan(bVal/aVal)))
ret
}
```

## Try the spatialCovariance package in your browser

Any scripts or data that you put into this service are public.

spatialCovariance documentation built on May 2, 2019, 9:43 a.m.