Fabien Ors wrote:Nice to hear that !
Best regards
Hi Fabien,
Like I noted, the anim.fit works well in RGoestats 12.0.1 BUT the Gibbs sampling is not giving good results as the histogram of the truncated Gaussian NASC still contain a substantial proportion of zero values (attachement).
Could you please find out what the problem could be (I have also attached the code I am using)
- Code: Select all
##############################################################################################
areaPoly <- read.csv("lakeVictoriaRGeostatPolygon.csv", header=T)
areaPoly
dat <- read.csv("dagaa_NASC_by_EDSU3.csv", header=T)# Fish density data
densDb <- db.create(dat, flag.grid=F)#flag.grid=F indicates do not overlay the coastline
densDb = db.locate(densDb, "long", "x",1)#attribute long corresponds to the first (x1) coordinate
densDb = db.locate(densDb, "lat", "x",2)#attribute lat corresponds to the second (x2) coordinate
densDb = db.locate(densDb, "NASC","z")#attribute of interest NASC is desiginated locator "Z"
poly.data <- polygon.create(areaPoly)
projec.query()# the projection is turned off
plot(densDb, zmin=0.001, pch.low=3, cex.low=0.25, las=1,
pch=21, col=1, inches=5, title="dagaa_2015", asp=1)
densDb <- infl(densDb,
nodes=c(400,400),
origin=c(min(areaPoly$X),min(areaPoly$Y)),
extend=c(abs(max(areaPoly$X)-min(areaPoly$X)),
abs(max(areaPoly$Y))+abs(min(areaPoly$Y))),
polygon=poly.data,
plot=T,
asp=1)
projec.define(projection="mean", db=densDb ,flag.update=T)
plot(densDb, zmin=0.001, pch.low=3, cex.low=0.25, las=1,
pch=21, col=1, inches=5,
title="dagaa_2015 (Mean Projection)", asp=1)
pvals <- projec.operate(densDb@items$lon,
densDb@items$lat)
llag <- mean(nndist(cbind(pvals$x,pvals$y)))
nlag <- round((max(pairdist(cbind(pvals$x,pvals$y)))/2)/llag)
max(pairdist(cbind(pvals$x,pvals$y)))
hist(densDb@items$NASC,
main="",
xlab="", # Leave x-label blank
breaks=50,
col="#7f7f7f",
border="#7f7f7f")
mtext(expression(atop(s[a]~m^2~nmi^-2)),
side=1,
line=4.5)
mean(densDb@items$NASC)
var(densDb@items$NASC)
model.anam <- anam.fit(densDb,
type="emp", # Empirical anamorphosis with logarithmic dispersion
draw=T,
title="Model of anamorphosis",
ndisc=densDb$nech,
sigma2e=40)
db.data.trans <- anam.z2y(densDb, anam=model.anam)
db.data.trans <- db.rename(db.data.trans,
name="Gaussian.NASC",
newname="Yp")
ycut <- qnorm(sum(db.extract(db.data.trans,"NASC") == 0) /
length(db.extract(db.data.trans, name="NASC",
flag.compress=T)))
Y <- db.extract(db.data.trans, "Yp", flag.compress=T)
hist(Y, breaks = 50, main = "",
col="#b3b3b3", border="#b3b3b3",
xlab="Anamorphosis Transformed NASC")
nasc <- db.extract(db.data.trans, "NASC")
Y[nasc == 0] <- ycut
db.data.trans <- db.replace(db.data.trans, "Yp", Y)
vario.Yp <- vario.calc(db.data.trans, lag=llag, nlag=nlag)
n.H <- 50
vario.Y <- vario.trans.cut(vario.Yp, ycut, n.H)#
plot(vario.Y, xlab="Distance (n.mi)",col="red",lwd=2)
model.vario.Y <- model.auto(vario.Y, struc=melem.name(c(1,2)),col="red", draw=T)
plot(vario.Y ,npairdw=T, npairpt=F, inches=0.06, col="red", xlab="Distance (n.mi.)")
plot(vario.Yp, npairdw=T, npairpt=F, inches=0.06, col="#000000", add=T)
plot(model.vario.Y, add=T, col="blue", lwd=3)
legend("bottomright", cex=0.75, y.intersp=1.2,
legend=c("Gaussian Transformed Variogram",
"Truncated Variable Variogram",
"Model Variogram"),
col=c("red","black","blue"), lty=1, pch=c(19,19,NA), lwd=c(1,1,2))# pch controls
## a Gibbs Sampler
Ymax <- db.extract(db.data.trans, name="Yp", flag.compress=F)
Ymin <- db.extract(db.data.trans, name="Yp", flag.compress=F)
Ymin[Ymin <= ycut] <- -10
db.data.trans <- db.add(db.data.trans,Ymax)
db.data.trans <- db.locate(db.data.trans, db.data.trans$natt,"upper")
db.data.trans <- db.add(db.data.trans, Ymin)
db.data.trans <- db.locate(db.data.trans,db.data.trans$natt,"lower")
db.data.trans <- db.locate(db.data.trans, "Yp", "z")
db.data.trans <- gibbs(db = db.data.trans,
model = model.vario.Y,
seed = 1337, # Seed for reproducibility
nboot = 1000, # Number of iterations performed for bootstrap
niter = 1000, # Maximum iterations for calculating the
# conditional expectation
flag.norm = FALSE, # For normalizing the model
percent = 0, # Nugget effect added to model (0 = model unchanged)
toleps = 1, # Relative difference between the previous and the
# new value at a sample location below which the sample
# is considered as immobile
radix = "Gibbs", # For naming results variable
modify.target = TRUE) # Define locator in new variable?
db.data.trans <- db.rename(db.data.trans, "Gibbs.G1", "Y2")
# Truncated variable
histYp <- hist(db.data.trans[,"Yp"],plot=F,breaks=seq(-4,4,.1))
hist(db.data.trans[,"Yp"], proba=T, breaks=seq(-4,4,.1),
col="#b3b3b3", border="#b3b3b3",
main="", xlab="Lower-cut Gaussian Variable",
xlim=c(-4,4), ylim=c(0,ceiling(max(histYp$density))))
# Truncated variable after Gibbs sampling
hist(db.extract(db.data.trans,"Y2",flag.compress=T), proba=T, breaks=seq(-4,4,.1),
xlab="Gaussian Transformed Variable", main="",
col="#b3b3b3", border="#b3b3b3",
xlim=c(-4,4), ylim=c(0,ceiling(max(histYp$density))))