Lab3

Respuestas Laboratorio 3

Exercise 1

library(spatstat)
K <- Kest(swedishpines)
plot(K, main = "K-function")

plot(K, . - pi*r^2 ~ r, legendpos = "bottomright")

#se calcula la función de correlación por pares
pcorf <- pcf(swedishpines)
plot(pcorf)

El gráfico de la función de correlación por pares y la función K muestran que el patrón es regular dado que la función k estimada está por debajo del modelo teórico y en la correlación por debajo de uno.

Exercise 2

help("rThomas")

## starting httpd help server ... done

replicate(3, plot(rThomas(10, 0.05, 8), main = ""))

## [[1]]
## Symbol map with no parameters
## 
## [[2]]
## Symbol map with no parameters
## 
## [[3]]
## Symbol map with no parameters

Se generaron diferentes realizaciones de un proceso de conglomerados, en promedio 10 conglomerados con

Se generaron 3 realizaciones de un proceso thomas con 10 cluster en promedio y cada uno con 8 puntos, debido a que este proceso genera los “offsprings” con una normal, la desviación estándar usada fue de 0.05.

#Clusters separados y al interior muy poca variabilidad
plot(rThomas(5, scale = 0.02, mu = 8), main = "") 

#Muy similar a un proceso poisson (aleatorio)
plot(rThomas(100, scale = 1.1, mu = 1), main = "")

Exercise 3

help("kppm")

#Se ajusta un modelo thomas, en kppm la opción por defecto es la estimación por mínimo contraste
fit <- kppm(redwood ~ 1, clusters="Thomas")
fit

## Stationary cluster point process model
## Fitted to point pattern dataset 'redwood'
## Fitted by minimum contrast
##  Summary statistic: K-function
## 
## Uniform intensity:   62
## 
## Cluster model: Thomas process
## Fitted cluster parameters:
##      kappa      scale 
## 23.5511449  0.0470461 
## Mean cluster size:  2.632568 points

plot(fit)

p <- parameters(fit)
rt2 <- rThomas(kappa = p$kappa, scale = p$scale, mu = p$mu)
plot(rt2, main = "")

plot(simulate(fit))

fit2 <- kppm(redwood ~ 1, clusters="Thomas", startpar=c(kappa=10, scale=0.1))

#variación del kappa y de la desviación
#kappa grande y desviación pequeña
kppm(redwood ~ 1, clusters="Thomas", startpar=c(kappa=100, scale=0.01))

## Stationary cluster point process model
## Fitted to point pattern dataset 'redwood'
## Fitted by minimum contrast
##  Summary statistic: K-function
## 
## Uniform intensity:   62
## 
## Cluster model: Thomas process
## Fitted cluster parameters:
##       kappa       scale 
## 23.54962713  0.04705395 
## Mean cluster size:  2.632738 points

#kappa pequeño y desviación grande
kppm(redwood ~ 1, clusters="Thomas", startpar=c(kappa=0.1, scale=10))

## Stationary cluster point process model
## Fitted to point pattern dataset 'redwood'
## Fitted by minimum contrast
##  Summary statistic: K-function
## 
## Uniform intensity:   62
## 
## Cluster model: Thomas process
## Fitted cluster parameters:
##       kappa       scale 
## 0.001824317 9.427722055 
## Mean cluster size:  33985.33 points

Se observa que para un kappa grande y una desviación pequeña los valores estimados son estables, sin embargo, para una kappa muy pequeño y una escala grande los resultados cambian con un promedio de puntos de 33985.

sim <- simulate(fit, nsim = 11)

## Generating 11 simulations... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,  11.
## Done.

sim[[12]] <- redwood
plot(sim, main = "", main.panel = "")

pcffit <- pcfmodel(fit)
plot(pcffit, xlim = c(0, 0.3)) #correlación por pares

plot(envelope(fit, Lest, nsim = 39, global = TRUE))

## Generating 78 simulated realisations of fitted cluster model (39 to estimate 
## the mean and 39 to calculate envelopes) ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77,  78.
## 
## Done.

El gráfico de los envelopes muestra que el modelo si parece plausible dado que se encuentra dentro de las bandas de confianza.

Exercise 4

mfit <- kppm(redwood ~ 1, clusters = "MatClust")
vcov(mfit)

##             (Intercept)
## (Intercept)  0.05304008

vcov(ppm(redwood ~ 1))

##             log(lambda)
## log(lambda)  0.01612903

Exercise 5

plot(Kest(swedishpines), abs(iso - theo) ~ r, main = "")

discrepancia <- function(r) {
  return(abs(as.function(Kest(swedishpines))(r) - pi*r^2))
}
res <- optimise(discrepancia, interval = c(0.1, 20), maximum = TRUE)
print(res)

## $maximum
## [1] 9.84372
## 
## $objective
## [1] 150.6897

R <- res$maximum

el valor de r que maximiza la discrepancia es 9.84

ppm(swedishpines ~ 1, Strauss(R))

## Stationary Strauss process
## 
## First order term:  beta = 0.08310951
## 
## Interaction distance:    9.84372
## Fitted interaction parameter gamma:   0.2407279
## 
## Relevant coefficients:
## Interaction 
##   -1.424088 
## 
## For standard errors, type coef(summary(x))

fit <- ppm(swedishpines ~ 1, Strauss(R))
print(fit) ##debido a que el parámetro gamma, existe una relación fuerte negativa entre los puntos a una distancia R.

## Stationary Strauss process
## 
## First order term:  beta = 0.08310951
## 
## Interaction distance:    9.84372
## Fitted interaction parameter gamma:   0.2407279
## 
## Relevant coefficients:
## Interaction 
##   -1.424088 
## 
## For standard errors, type coef(summary(x))

plot(fitin(fit))

Exercise 6

plot(Kest(swedishpines))

rvals <- seq(5, 12, by = 0.1)
D <- data.frame(r = rvals)
fitp <- profilepl(D, Strauss, swedishpines ~ 1)

## (computing rbord)

## comparing 71 models...

## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,  71.

## fitting optimal model...

## done.

print(fitp)

## profile log pseudolikelihood
## for model:  ppm(swedishpines ~ 1,  interaction = Strauss)
## fitted with rbord = 12
## interaction: Strauss process
## irregular parameter: r in [5, 12]
## optimum value of irregular parameter:  r = 9.8

plot(fitp)

(Ropt <- reach(as.ppm(fitp)))

## [1] 9.8

La estimación de 9.8 es consistente con la estimación previa.

bestfit <- as.ppm(fitp)

Exercise 7

fit.mpl <- ppm(swedishpines ~ 1, Strauss(R), method = "mpl")
fit.ho  <- ppm(swedishpines ~ 1, Strauss(R), method = "ho")

## Simulating... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,  100.
## Done.

print(fit.ho)

## Stationary Strauss process
## 
## First order term:  beta = 0.1054425
## 
## Interaction distance:    9.84372
## Fitted interaction parameter gamma:   0.2297704
## 
## Relevant coefficients:
## Interaction 
##   -1.470675 
## 
## For standard errors, type coef(summary(x))

print(fit.mpl)

## Stationary Strauss process
## 
## First order term:  beta = 0.08310951
## 
## Interaction distance:    9.84372
## Fitted interaction parameter gamma:   0.2407279
## 
## Relevant coefficients:
## Interaction 
##   -1.424088 
## 
## For standard errors, type coef(summary(x))

Los ajustes por este método Huang-Ogata y maximo seudoverosimilitud son similares.

Exercise 8

fitp <- profilepl(D, Strauss, swedishpines ~ polynom(x,y,2))

## (computing rbord)

## comparing 71 models...

## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,  71.

## fitting optimal model...

## done.

fitp2 <- profilepl(D, Strauss, swedishpines ~ polynom(x,y,2))

## (computing rbord)

## comparing 71 models...

## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,  71.

## fitting optimal model...

## done.

print(fitp2)

## profile log pseudolikelihood
## for model:  ppm(swedishpines ~ polynom(x,  y,  2),  interaction = Strauss)
## fitted with rbord = 12
## interaction: Strauss process
## irregular parameter: r in [5, 12]
## optimum value of irregular parameter:  r = 9.8

(bestfit <- as.ppm(fitp2))

## Nonstationary Strauss process
## 
## Log trend:  ~x + y + I(x^2) + I(x * y) + I(y^2)
## 
## Fitted trend coefficients:
##   (Intercept)             x             y        I(x^2)      I(x * y) 
## -5.156605e+00  4.269469e-02  8.582459e-02 -5.324026e-06 -7.955428e-04 
##        I(y^2) 
## -5.506409e-04 
## 
## Interaction distance:    9.8
## Fitted interaction parameter gamma:   0.2390264
## 
## Relevant coefficients:
## Interaction 
##   -1.431181 
## 
## For standard errors, type coef(summary(x))

reach(bestfit)

## [1] 9.8