martes, 21 de diciembre de 2010

Simulación de procesos estocásticos en R.

Teoría sobre generación de números pseudo-aleatorios
simulación de procesos estocásticos 5

Generación de números aleatorios, Generación de variables chi-cuadrado, normales, etc. por método de transformación inversa
simulación de procesos estocásticos 4

Procesos de Poisson, Movimiento Browniano, Sistemas de colas
simulación de procesos estocásticos 3

Simulación Monte Carlo
Simulacion_MonteCarlo2

Sistemas de colas
simulación de procesos estocásticos 2

Monte Carlo y cálculo de integrales
simulación de procesos estocásticos

Suavizado de series temporales con R.

Procesamiento preliminar de una serie temporal.

FUENTE: http://www.stat.pitt.edu/stoffer/tsa2/textRcode.htm#ch

##################################################
#Example 2.2: matriz de dispersión (scatterplot)
mort = ts(scan("/mydata/cmort.dat"),start=1970, frequency=52)
temp = ts(scan("/mydata/temp.dat"), start=1970, frequency=52)
part = ts(scan("/mydata/part.dat"), start=1970, frequency=52)
par(mfrow=c(3,1), mar=c(3,4,3,2))
plot(mort, main="Cardiovascular Mortality", xlab="", ylab="")
plot(temp, main="Temperature", xlab="", ylab="")
plot(part, main="Particulates", xlab="", ylab="")
x11() # open another graphic device
pairs(cbind(mort, temp, part)) # scatterplot matrix
temp = temp-mean(temp)
temp2 = temp^2
trend = time(mort) # time
fit = lm(mort~ trend + temp + temp2 + part, na.action=NULL)
summary(fit) # regression results
anova(lm(mort~cbind(trend, temp, temp2, part))) # pretty ANOVA table
num = length(mort) # sample size
AIC(fit)/num - log(2*pi) # AIC (corresponds to def 2.1)
AIC(fit, k=log(num))/num - log(2*pi) # BIC (corresponds to def 2.3)
(AICc = log(sum(resid(fit)^2)/num)+(num+5)/(num-5-2)) # AICc (as in def 2.2)

#Example 2.10: suavizado por media móvil
ma5 = filter(mort, sides=2, rep(1,5)/5)
ma53 = filter(mort, sides=2, rep(1,53)/53)
plot(mort, type="p", xlab="week", ylab="mortality")
lines(ma5)
lines(ma53)

#Example 2.11: suavizado por regresión polinómica y periódica
wk = time(mort) - mean(time(mort)) # week (centered) - wk is basically t/52 ...
wk2 = wk^2 # ... because frequency=52 for mort
wk3 = wk^3
c = cos(2*pi*wk)
s = sin(2*pi*wk)
reg1 = lm(mort~wk + wk2 + wk3, na.action=NULL)
reg2 = lm(mort~wk + wk2 + wk3 + c + s, na.action=NULL)
plot(mort, type="p", ylab="mortality")
lines(fitted(reg1))
lines(fitted(reg2))

# here's the original code where mort isn't a ts object
mort = scan("/mydata/cmort.dat")
t = 1:length(mort)
t2 = t^2
t3 = t^3
c = cos(2*pi*t/52)
s = sin(2*pi*t/52)
fit1 = lm(mort~t + t2 + t3)
fit2 = lm(mort~t + t2 + t3 + c + s)
plot(t, mort)
lines(fit1$fit)
lines(fit2$fit)

#Example 2.12: suavizado kernel
plot(mort, type="p", ylab="mortality")
lines(ksmooth(time(mort), mort, "normal", bandwidth=10/52)) # or try bandwidth=5/52
lines(ksmooth(time(mort), mort, "normal", bandwidth=2))
# - original code when mort wasn't a time series with frequency=52
plot(t, mort)
lines(ksmooth(t, mort, "normal", bandwidth=10))
lines(ksmooth(t, mort, "normal", bandwidth=104))
# Red indicates a change; Fig. 2.15 shows a bandwidth of 10.

#Example 2.13: Regresión ponderada localmente y vecino más cercano
par(mfrow=c(2,1))
plot(mort,type="p", ylab="mortality", main="nearest neighbor")
lines(supsmu(time(mort), mort, span=.5))
lines(supsmu(time(mort), mort, span=.01))
plot(mort, type="p", ylab="mortality", main="lowess")
lines(lowess(mort, f=.02))
lines(lowess(mort, f=2/3))

#Example 2.14: Suavizado por splines
plot(mort, type="p", ylab="mortality")
lines(smooth.spline(time(mort), mort))
lines(smooth.spline(time(mort), mort, spar=1))
# - original code when mort wasn't a time series with frequency=52
plot(t, mort)
lines(smooth.spline(t, mort, spar=.0000001))
lines(smooth.spline(t, mort, spar=.1))

#Example 2.15: Suavizado de una serie como función de otra
par(mfrow=c(2,1))
plot(temp, mort, main="lowess")
lines(lowess(temp,mort))
plot(temp, mort, main="smoothing splines")
lines(smooth.spline(temp,mort))

Estimación de modelos de series temporales con R.

Una vez que se ha identificado la estructura del modelo, calcularemos los valores de los parámetros involucrados en el mismo. El procedimiento de estimación más usado es el de máxima verosimilitud, aunque para modelos autorregresivos, la estimación por el método de los momentos es bastante simple.
La maximización de la función de verosimilitud es no lineal y por tanto se realiza numéricamente. Por ello, la convergencia al máximo será más rápida si se parte de un valor inicial de los parámetros próximo al valor de convergencia. A continuación se proponen 4 métodos para el cálculo de estos valores iniciales:
  1. Para el caso autoregresivo (AR)
  • Método de Yule-Walker
  • Algoritmo de Burg
2. Para el caso general (ARMA)
  • Algoritmo de las innovaciones
  • Algoritmo de Hannan-Rissanen
Estimamos los valores iniciales por métodos de máxima verosimilitud.
Una vez presentados los métodos para la estimación preliminar, pasamos a plantear el método de máxima verosimilitud. En este caso se supone que los residuos tienen distribución normal con media cero y varianza sigma cuadrado e independientes.

Referencias

lunes, 20 de diciembre de 2010

Validación de modelos de series temporales con R.

Aquí se concluyen las fases de modelización ARMA usando el mérodo de Box-Jenkins. Se utilizan herramientas para decidir si el modelo es o no adecuado para describir el comportamiento de los datos. En el caso de que no se considere adecuado, debemos volver a la fase de identificación y proporner un nuevo modelo.
Este análisis de basa en la evaluación de los residuos, y comparando su comportamiento con el del ruido blanco.
  • Estudio de incorrelación: los residuos deben estar incorrelados. Se utilizan test de tipo Pomanteau: test de Box-Pierce, Test de Ljung-Box, Test de McLEod y/o Test McLeod-Li. Otros test estadísticos generales son el Test de Durbin-Watson, Test de Wallis, etc..
  • Estudio de aleatoriedad: los residuos deben ser aleatorios. Evaluamos los residuos con test estadísticos, como el Test de puntos de cambios, el Test de rachas y/o el Test de rangos.
  • Estudio de normalidad: los residuos deben tener una distribución cercana a la normalidad. La ausencia de normalidad sugiere la necesidad de realizar transformaciones Box-Cox. Este análisis se realiza con gráficos Q-Qplot, y Test estadísticos como el Test DÁgostino, Test Shapiro-Wilk, etc..
También analizamos:
  • Sobreajuste: se trata de ajustar un modelo más complejo (con más parámetros) y estudiar la significación de los nuevos términos.
  • Correlación entre los estimadores: debe ser pequeña, sino las estimaciones serán inestables.
  • Poder predictivo: un método para elegir entre modelos es evaluar el menor error cuadrático medio de predicción. Otros estadísticos útiles son: AIC, BIC, etc..

Referencias

Identificación de modelos de series temporales con R.



En la fase de identificación aplicaremos distintos procedimientos para presentar modelos tentativos que puedan resumir el comportamiento de la serie temporal.

Incorporamos toda la información disponible sobre la serie:
  • estudios relacionados (la misma serie en otros períodos u otras series que se sospeche tengan comportamiento similar)
  • fenómeno subyacente
  • intervenciones externas
  • variables que modifiquen el comportamiento de nuestra serie
  • intervalo de observación (puede indicar la presencia de periodicidades)
Los pasos para el análisis preliminar de la serie son:

1. Graficar la serie: valores de la serie vs. tiempo. Puede poner de manifiesto:
  • autocorrelación (positiva -persistencia- o negativa -grandes fluctuaciones-)
  • estacionalidad (patrón que se repite aproximadamente transcurridos unos instantes de tiempo)
  • estacionariedad (puede existir no-estacionariedad en la media o en la varianza)
  • tendencias (caso particular de no-estacionariedad en la media)
  • necesidad de transformación de la serie (e.g. Box-Cox)
  • valores extremos (anómalos u outliers)
  • ciclos de período largo
  • intervenciones (cambios en el nivel de la serie por ingerencias externas a la misma)
Realizaremos un ejemplo en R dado por Shumway & Stoffer en su página web (descarcar desde allí el archivo de datos).
#Let's play with the Johnson & Johnson data set. Download jj.dat  to a directory called mydata (or wherever you choose ... the examples below and in the text assume the data are in that directory).
jj = scan("/mydata/jj.dat") # read the data
jj <- scan("/mydata/jj.dat") # read the data another way
scan("/mydata/jj.dat") -> jj # and another
#Now, let's make jj a time series object.
jj = ts(jj, start=1960, frequency=4) #Note that the data are quarterly earnings, hence the frequency=4 statement.

#Now try a plot of the data:
plot(jj, ylab="Earnings per Share", main="J & J")
#Try these and see what happens:
plot(jj, type="o", col="blue", lty="dashed")
plot(diff(log(jj)), main="logged and diffed")

Componentes de la variación: consideraremos los siguientes modelos de descomposición de series temporales.
  • aditivo: Yt=Tt+Ct+St+It. La varianza de ls serie original y la tendencia aumentan con el tiempo. Es decir, St e It son magnitudes que se agregan a Tt y Ct independientemente del valor que tenda.
  • multiplicativo: Yt=Tt*Ct*St+It. El efecto estacional tiene de aumentar al aumentar la tendencia. Es decir, St e It son proporciones de la tendencia Tt.
donde Tt es la tendencia, Ct es la ciclicidad (fluctuaciones en intervalos mayores al año, irregular o al azar), St es la estacionalidad (fluctuaciones periódicas de menso de un año) e It la irregularidad o aleatoriedad o error.

  1. Identificación de la estructura estacionaria: identificar el modelo a ajustar. Lo conseguimos con los pasos anteriores.
  2. Identificación de la estructura no-estacionaria: detectar qué transformaciones hay que aplicar a la serie para conseguir una meda y varianza constantes.
  • Sea xt=mu+yt, podemos eliminar la no-estacionaridad de la tendencia: eliminamos la tendencia al restarle a la serie su media estimada (detrending, yt=xt-mut) o diferenciando la serie (differencing, xt-xt-1). Una ventaja de diferenciar en lugar de eliminar la tendencia es que no estimamos parámetros en la diferenciación, pero como desventaja, la diferenciación no produce una estimación del proceso estocástico yt. Para aplicar diferencias definimos el operador "backshift": Bxt=xt-1, o para cualquierorden d B^dxt=xt-k o diff(d)=(1-B^d).
    Estimar la tendencia: por media móvil, suavizado/filtrado o descomposición. suavizado o alisado (smoothing): por media móvil, kernel de suavizado, regresión ponderado localmente y vecino más cercano, splines.
Ver aquí varios ejemplos de suavizado
  • no-estacionaridad de la varianza: I) transformación logarítimica de la serie: yt=ln(xt). Tiende a eliminar fluctuaciones largas que ocurren sobre porciones de la serie donde los valores subyacentes sin grandes. II) transformación Box-Cox o potencial de la serie: yt= 1) (xt^lambda-1)/lambda si lambda es distinto que cero y 2) ln(xt) si lambda es igual a cero. La transformación Box-Cox se utiliza para mejorar la aproximación a la normalidad o mejorar la linealidad en los valores predichos.
También el análisis de los componentes de variación, el análisis exploratorio de los datos (EDA), el el análisis de regresión y periodograma nos permiten describir la serie.
FILTRADO/SUAVIZADO
#How about filtering/smoothing the Johnson & Johnson series using a two-sided moving average? Let's try this:

fjj(t) = 1/8 jj(t-2) + ¼ jj(t-1) + ¼ jj(t) + ¼ jj(t+1) + 1/8 jj(t+2) #and we'll add a lowess fit for fun.
k = c(.5,1,1,1,.5) # k is the vector of weights
fjj = filter(jj, sides=2, k) # ?filter for help [but you knew that already]
plot(jj)
lines(fjj, col="red") # adds a line to the existing plot
lines(lowess(jj), col="blue", lty="dashed")

DIFERENCIACIÓN
#Let's difference the logged data and call it dljj. Then we'll play with dljj:
dljj = diff(log(jj)) # difference the logged data
plot(dljj) # plot it if you haven't already
shapiro.test(dljj) # test for normality
#Now a histogram and a Q-Q plot, one on top of the other:
par(mfrow=c(2,1)) # set up the graphics
hist(dljj, prob=TRUE, 12) # histogram
lines(density(dljj)) # smooth it - ?density for details
qqnorm(dljj) # normal Q-Q plot
qqline(dljj) # add a line

DESCOMPOSICIÓN
#let's try a structural decomposition of log(jj) = trend + season + error using lowess. Note, this example works only if jj is dimensionless (i.e., you didn't read it in using read.table ... thanks to Jon Moore of the University of Reading, U.K., for pointing this out.)

plot(dog <- stl(log(jj), "per"))
#You can do a little (very little) better using a local seasonal window, plot(dog <- stl(log(jj), s.win=4)), as opposed to the global one used by specifying "per". Type ?stl for details.

REGRESIÓN
#We're going to fit the regression log(jj)= beta*time + alfa1*Q1 + alfa2*Q2 + alfa3*Q3 + alfa4*Q4 + epsilon where Qi is an indicator of the quarter i = 1,2,3,4. Then we'll inspect the residuals.

Q = factor(rep(1:4,21)) # make (Q)uarter factors [that's repeat 1,2,3,4, 21 times]
trend = time(jj)-1970 # not necessary to "center" time, but the results look nicer
reg = lm(log(jj)~0+trend+Q, na.action=NULL) # run the regression without an intercept
#-- the na.action statement is to retain time series attributes
summary(reg)
model.matrix(reg) #You can view the model matrix (with the dummy variables) this way:
#Now check out what happened. Look at a plot of the observations and their fitted values:
plot(log(jj), type="o") # the data in black with little dots
lines(fitted(reg), col=2) # the fitted values in bloody red - or use lines(reg$fitted, col=2)
#and a plot of the residuals and the ACF of the residuals:
par(mfrow=c(2,1))
plot(resid(reg)) # residuals - reg$resid is same as resid(reg)
acf(resid(reg),20) # acf of the resids. Do those residuals look white? [Ignore the 0-lag correlation, it's always 1.]
#You have to be careful when you regress one time series on lagged components of another using lm(). There is a package called dynlm that makes it easy to fit lagged regressions, and I'll discuss that right after this example. If you use lm(), then what you have to do is "tie" the series together using ts.intersect. If you don't tie the series together, they won't be aligned properly.

2. Función de correlación (simple -FAC- y parcial -FACP-). Sugieren:
  • no-estacionariedad: si el FAC presenta persistencia o decaimiento lento (no exponencial)
  • los cortes/decaimientos en FAC y FACP sugiere la estructura del modelo que debe ajustarse a los datos.
  • matrices de dispersión (scatterplot): observamos las correlaciones a distintos lags.
Conceptos:
  • función media
  • función autocovarianza
  • función de autocorrelación (ACF o FAC)
  • función covarianza cruzada
  • función correlación cruzada (CCF o FCC)
Modelos: ruido blanco, MA(q), AR(p), ARMA()

#Let's check out the correlation structure of dljj using various techniques. First, we'll look at a grid of scatterplots of dljj(t-lag) vs dljj(t) for lag=1,2,...,9.
lag.plot(dljj, 9, do.lines=FALSE) # why the do.lines=FALSE? ... try leaving it out
#Notice the large positive correlation at lags 4 and 8 and the negative correlations at a few other lags:

#Now let's take a look at the ACF and PACF of dljj:
par(mfrow=c(2,1)) # The power of accurate observation is commonly called cynicism
# by those who have not got it. - George Bernard Shaw
acf(dljj, 20) # ACF to lag 20 - no graph shown... keep reading. Ignore the 0-lag correlation, it's always 1
pacf(dljj, 20) # PACF to lag 20 - no graph shown... keep reading
#Note that the LAG axis is in terms of frequency, so 1,2,3,4,5 correspond to lags 4,8,12,16,20 because frequency=4 here. If you don't like this type of labeling, you can replace dljj in any of the above by ts(dljj, freq=1); e.g., acf(ts(dljj, freq=1), 20)


3. Modelos preliminares: a partir de los gráficos anteriores, analizamos los componentes de la variación y planteamos todos los modelos posibles que se puedan ajustar a los datos.


Referencias

viernes, 10 de diciembre de 2010

Análisis geoestadístico con R. Paquete geoR. parte II.

Más funciones para el análisis de datos geoestadístico según el libro"Model-based Geostatistics" de Peter J. Diggle & Paulo Justiniano Ribeiro Jr, 2007.
Iré agregando material teórico y práctico sobre el tema.
Book data set: http://www.leg.ufpr.br/doku.php/pessoais:paulojus:mbgbook:datasets

##############################################################
#### Diggle Model Based Geoestadistics. Paquete geoR en R ####
##############################################################
library(geoR)
data(elevation)

###2.8 perspectiva general de modelos geoestadísiticos:
#análisis exploratorio no-espacial
plot(elevation)

with(elevation, hist(data, main = "", xlab = "elevation"))
with(elevation, plot(coords[, 1], data, xlab = "W-E", ylab = "elevation data", pch = 20, cex = 0.7))
lines(lowess(elevation$data ~ elevation$coords[, 1]))
with(elevation, plot(coords[, 2], data, xlab = "S-N", ylab = "elevation data", pch = 20, cex = 0.7))
lines(with(elevation, lowess(data ~ coords[, 2])))
#modelo lineal con covariables. The values "1st" and "2nd" passed to the argument trend are aliases to indicate first-degree and second-degree polynomials on the coordinates.
points(elevation, cex.max = 2.5)
points(elevation, trend = "1st", pt.div = 2, abs = T, cex.max = 2.5)
points(elevation, trend = "2nd", pt.div = 2, abs = T, cex.max = 2.5)
#calculo de los variogramas empíricos para los datos y los residuales
plot(variog(elevation, uvec = seq(0, 5, by = 0.5)), type = "b")
res1.v <- variog(elevation, trend = "1st", uvec = seq(0, 5, by = 0.5))
plot(res1.v, type = "b")
res2.v <- variog(elevation, trend = "2nd", uvec = seq(0, 5, by = 0.5))
lines(res2.v, type = "b", lty = 2)
#variograma residual y envolventes de simulación bajo permutación aleatoria de los residuales.
set.seed(231)
mc1 <- variog.mc.env(elevation, obj = res1.v)
plot(res1.v, env = mc1, xlab = "u")
mc2 <- variog.mc.env(elevation, obj = res2.v)
plot(res2.v, env = mc2, xlab = "u")

#estimaciones por ML (maximum likelihood) del modelo Gaussiano sin y con el término de tendencia.
ml0 <- likfit(elevation, ini = c(3000, 2), cov.model = "matern", kappa = 1.5); ml0
ml1 <- likfit(elevation, trend = "1st", ini = c(1300, 2), cov.model = "matern", kappa = 1.5); ml1
#interpolación espacial con kriging
locs <- pred_grid(c(0, 6.3), c(0, 6.3), by = 0.1)
KC <- krige.control(type = "sk", obj.mod = ml0)
sk <- krige.conv(elevation, krige = KC, loc = locs)
KCt <- krige.control(type = "sk", obj.mod = ml1, trend.d = "1st", trend.l = "1st")
skt <- krige.conv(elevation, krige = KCt, loc = locs)
#selección de funciones gráficas para producir mapas usando argumentos opcionales para asegurar que los pares usan la misma escala gris
pred.lim <- range(c(sk$pred, skt$pred))
sd.lim <- range(sqrt(c(sk$kr, skt$kr)))
image(sk, col = gray(seq(1, 0, l = 51)), zlim = pred.lim)
contour(sk, add = T, nlev = 6)
points(elevation, add = TRUE, cex.max = 2)
image(skt, col = gray(seq(1, 0, l = 51)), zlim = pred.lim)
contour(skt, add = T, nlev = 6)
points(elevation, add = TRUE, cex.max = 2)
image(sk, value = sqrt(sk$krige.var), col = gray(seq(1, 0, l = 51)), zlim = sd.lim)
contour(sk, value = sqrt(sk$krige.var), levels = seq(10, 27, by = 2), add = T)
points(elevation$coords, pch = "+")
image(skt, value = sqrt(skt$krige.var), col = gray(seq(1, 0, l = 51)), zlim = sd.lim)
contour(skt, value = sqrt(skt$krige.var), levels = seq(10, 27, by = 2), add = T)
points(elevation$coords, pch = "+")
#dejar fuera el efecto de la covariable: remoción de tendencia
plot(elevation, low = TRUE, trend = ~coords[, 2], qt.col = 1)

### 3.13 Modelos Gaussianos para datos geoestadísticos
#cálculo y gráfico de la función de correlación estándar. Se elige la familia de correlaciones con el argumento cov.model.
x <- seq(0, 1, l = 101)
plot(x, cov.spatial(x, cov.model = "mat", kappa = 0.5, cov.pars = c(1, 0.25)), type = "l", xlabel = "u", ylabel = expression(rho(u)), ylim = c(0, 1))
lines(x, cov.spatial(x, cov.model = "mat", kappa = 1.5, cov.pars = c(1, 0.16)), lty = 2)
lines(x, cov.spatial(x, cov.model = "mat", kappa = 2.5, cov.pars = c(1, 0.13)), lty = 3)

#generar simulaciones de procesos Gaussianos 2D. La función grf() especifica el modelo y las localizaciones para lso valores simulados requeridos.
set.seed(159) #asegura que las simulaciones se gneran con el mismo número aleatorio de semilla, y por tanto las realizaciones simuladas tienen diferencias solo por los valores diferentes de los parámetros del modelo.
image(grf(100^2, grid = "reg", cov.pars = c(1, 0.25)), col = gray(seq(1, 0, l = 51)), xlab = "", ylab = "")
set.seed(159)
image(grf(100^2, grid = "reg", cov.pars = c(1, 0.13), cov.model = "mat", kappa = 2.5), col = gray(seq(1, 0, l = 51)), xlab = "", ylab = "")
#simulaciones de modelos anisotrópicos
set.seed(421)
image(grf(201^2, grid = "reg", cov.pars = c(1, 0.25), aniso.pars = c(pi/3, 4)), col = gray(seq(1, 0, l = 51)), xlab = "", ylab = "")
set.seed(421)
image(grf(201^2, grid = "reg", cov.pars = c(1, 0.25), aniso.pars = c(3 * pi/4, 2)), col = gray(seq(1, 0, l = 51)), xlab = "", ylab = "")
#In this case, we have used a stationary Gaussian model with mean mu = 850, nugget variance tau^2 = 100, signal variance sigma^2 = 3500 and MatLern correlation function with k = 2.5 and phi= 0.8.
sim1 <- grf(100, cov.pars = c(1, 0.15), cov.model = "matern", kappa = 1.5)
points(sim1)
data(elevation)
sim2 <- grf(grid = elevation$coords, cov.pars = c(3500, 0.8), nugget = 100)
sim2$data <- sim2$data + 850
points(sim2)
### 4.6 Modelos lineales generalizados (GLM) para datos geoestadísticos
 #simulación con el GLM de un modelo de Poisson
#simulate a realisation of the Gaussian process at these locations with mu= 0.5, sigma^2 = 2 and MatLern correlation function with k = 1.5, phi = 0.2.
set.seed(371)
cp <- expand.grid(seq(0, 1, l = 10), seq(0, 1, l = 10))
s <- grf(grid = cp, cov.pars = c(2, 0.2), cov.model = "mat", kappa = 1.5)
image(s, col = gray(seq(1, 0.2, l = 21)))
lambda <- exp(0.5 + s$data)
y <- rpois(length(s$data), lambda = lambda)
text(cp[, 1], cp[, 2], y, cex = 1.5)
#extensión del modelo, incluyendo variación no-espacial extra de Poisson
lambda <- exp(s$data + tau * rnorm(length(s$data)))
#simulación con GLM del modelo Bernoulli
set.seed(34)
locs <- seq(0, 1, l = 401)
s <- grf(grid = cbind(locs, 1), cov.pars = c(5, 0.1), cov.model = "matern", kappa = 1.5)
p <- exp(s$data)/(1 + exp(s$data))
ind <- seq(1, 401, by = 8)
y <- rbinom(p[ind], size = 1, prob = p)
plot(locs[ind], y, xlab = "locations", ylab = "data")
lines(locs, p)

 #simulación con GLM del modelo Binomial. Con [Y (x)|S] ¡« Bin{n, p(x)} with n = 5 and p(x) = exp{mu + S(x)}/[1 + exp{mu + S(x)}],
set.seed(23)
s <- grf(60, cov.pars = c(5, 0.25))
p <- exp(2 + s$data)/(1 + exp(2 + s$data))
y <- rbinom(length(p), size = 5, prob = p)
points(s)
text(s$coords, label = y, pos = 3, offset = 0.3)





### Estimación de parámetros clásica: estimaciónd e tendencia, variogramas, estimación de la estructura de covarianza (OLS, WLS, etc), ML, estimación paramétrica por modelos lineales generalizados, MCML, etc.
#calculo de variograma con diferentes opciones del tamaño bin
data(elevation)
plot(variog(elevation, option = "cloud"), xlab = "u", ylab = "V(u)")
plot(variog(elevation, uvec = seq(0, 8, by = 0.5)), xlab = "u", ylab = "V(u)")
plot(variog(elevation, trend = "2nd", max.dist = 6.5), xlab = "u", ylab = "V(u)")


    data(s100) #datos simulados con un modelo con media=cero, varianza=uno, nugget=0, y función de correlación exponencial con phi=.3.
v1 <- variog(s100); plot(v1)
v2 <- variog(s100, uvec = seq(0, 1, by = 0.1)); plot(v2)
v3 <- variog(s100, max.dist = 1); plot(v3)
round(v1$u, dig = 2) #The midpoints of the bins obtained from the commands listed above
round(v2$u, dig = 2)
round(v3$u, dig = 2)

    data(ca20)
plot(ca20,trend= area+altitude, low=T) #sugiere tendencia cuadratica
plot(variog(ca20, max.dist = 510)) #construimos distintos variogramas: observamos que incorporando sub-regoses al modelo este se comporta aproximadamente estacionario en los residuales
plot(variog(ca20, trend = ~area, max.dist = 510))
plot(variog(ca20, trend = ~area + altitude, max.dist = 510))
t.all <- trend.spatial(ca20, trend = ~area + altitude, add = "2nd")
plot(variog(ca20, trend = ~t.all, max.dist = 510))


set.seed(83) #variogramas de las simulaciones del proceso Gaussiano, cada gráfico indica el variograma verdadero y la curva suavizada, así como los variogramas empíricos de las tres simulaciones del proceso en 100 localizaciones.
sim1 <- grf(100, cov.pars = c(1, 0.05), cov.model = "mat", kap = 1.5, nsim = 3) #Parameter values are k 1.5, sigma^2= 1, tau^2 = 0 in both panels. In the left-hand panel, phi= 0.05, in the right-hand panel phi= 0.2.
plot(variog(sim1, max.dist = 1), type = "l", lty = 1:3, col = 1)
lines.variomodel(seq(0, 1, l = 100), cov.model = "mat", kap = 1.5, cov.pars = c(1, 0.05), nug = 0)
set.seed(83)
sim2 <- grf(100, cov.pars = c(1, 0.2), cov.model = "mat", kap = 1.5, nsim = 3)
plot(variog(sim2, max.dist = 1), type = "l", lty = 1:3, col = 1)
lines.variomodel(seq(0, 1, l = 100), cov.model = "mat", kap = 1.5, cov.pars = c(1, 0.2), nug = 0)






#estimación de parámetros en modelos gaussianos con particular foco en los parámetros que definen la estructura de covarianza del modelo.
s100.v <- variog(s100, max.dist = 1) #estimación mediante variograma
plot(s100.v); lines.variomodel(seq(0, 1, l = 100), cov.pars = c(0.9, 0.2), cov.model = "mat", kap = 1.5, nug = 0.2)






wls <- variofit(s100.v, ini = c(0.9, 0.2), cov.model = "mat", kap = 1.5, nug = 0.2); wls #estimación mediante WLS
ml <- likfit(s100, cov.model = "mat", kap = 1.5, ini = c(0.9, 0.2), nug = 0.2); ml #estimación por ML
reml <- likfit(s100, cov.model = "mat", kap = 1.5, ini = c(0.9, 0.2), nug = 0.2, met = "reml"); reml #estimación por ML restringido
#ajuste de 5 modelos alternaticos a los datos de suelo y sus covariables.
data(ca20)
m1 <- likfit(ca20, ini = c(100, 200), nug = 50)
m2 <- likfit(ca20, trend = ~area, ini = c(60, 100), nug = 40)
m3 <- likfit(ca20, trend = ~area + altitude, ini = c(60, 100), nug = 40)
m4 <- likfit(ca20, trend = ~area + coords, ini = c(60, 100), nug = 40)
m5 <- likfit(ca20, trend = ~area + altitude + coords, ini = c(60, 100), nug = 40)

### Predicción espacial
#predicción por "conventional kriging"
data(SIC)
ml <- likfit(sic.all, ini = c(100, 40), nug = 10, lambda = 0.5, kappa = 1) #estimjación de los parámetros del modelo
gr <- pred_grid(sic.borders, by = 7.5) #definimos una grilla para lso puntos de predicción
KC <- krige.control(obj.model = ml) #pasa los parámetros del modelo mediante "ordinary krigning"
OC <- output.control(n.pred = 1000, simul = TRUE, thres = 250) #incluye opciones para generar y guardar 1000 simulaciones de las distribuciones condicionales de S dado Y y definir un valor umbral de 250 que será usado para calcular la probabilidad de excedencia en cada localización de la predicción.
set.seed(2419) #permite la reproducción de los resultados de la simulación
pred <- krige.conv(sic.all, loc = gr, borders = sic.borders,krige = KC, out = OC)
image(pred, col = gray(seq(1, 0.1, l = 21)), zlim = predlim, x.leg = c(0, 350), y.leg = c(-60, -30))
image(pred, loc = gr, val = sqrt(pred$krige.var), zlim = selim, col = gray(seq(1, 0.1, l = 21)), x.leg = c(0, 350), y.leg = c(-60, -30))
dim(pred$simulations)
A200 <- apply(pred$simul, 2, function(y) sum(y > 200)/length(y))
hist(A200, main = "")
image(pred, val = 1 - pred$prob, col = gray(seq(0.9, 0.1, l = 41)), x.leg = c(0, 350), y.leg = c(-60, -30))



#el mismo análisis para los datos de calcio en el suelo
data(ca20)
fit <- likfit(ca20, ini = c(100, 60), trend = ~area)
gr <- pred_grid(ca20$borders, by = 10)
gr0 <- polygrid(gr, borders = ca20$border, bound = T)
ind.reg <- numeric(nrow(gr0))
ind.reg[.geoR_inout(gr0, ca20$reg1)] <- 1
ind.reg[.geoR_inout(gr0, ca20$reg2)] <- 2
ind.reg[.geoR_inout(gr0, ca20$reg3)] <- 3
ind.reg <- as.factor(ind.reg)
KC <- krige.control(trend.d = ~area, trend.l = ~ind.reg, obj.model = fit)
ca20pred <- krige.conv(ca20, loc = gr, krige = KC)
par(mar = c(2.8, 2.5, 0.5, 0.5), mgp = c(1.8, 0.7, 0), mfrow = c(1, 2))
image(ca20pred, loc = gr, col = gray(seq(1, 0, l = 21)), x.leg = c(4930, 5350), y.leg = c(4790, 4840))
polygon(ca20$reg1)
polygon(ca20$reg2)
polygon(ca20$reg3)
image(ca20pred, loc = gr, val = sqrt(ca20pred$krige.var), col = gray(seq(1, 0, l = 21)), x.leg = c(4930, 5350), y.leg = c(4790, 4840))
polygon(ca20$reg1)
polygon(ca20$reg2)
polygon(ca20$reg3)
coords <- cbind(c(0.2, 0.25, 0.6, 0.7), c(0.1, 0.8, 0.9, 0.3))
KC <- krige.control(ty = "ok", cov.model = "mat", kap = 1.5, nug = 0.1, cov.pars = c(1, 0.1))
krweights(coords, c(0.5, 0.5), KC)