|
|
Hi,
I have been trying to use the kde2d function in R to evaluate whether a pair of points lie within 95% contour band. I calculate whether the points lie in 2.5% contour and whether it lies in 97.5% contour and then take the difference. But I a getting negative values (variable 'within' below). All points in the inner contour should also be in the out contour which doesnt seem to be the case.
library(MASS)
library(car)
set.seed(1234)
require(mixtools)
p1<- c(0.363285618,0.136169868,0.073852521,0.603122717,0.120757527,0.110571766,0.254116432,
0.236207423,0.247417764,0.293465209,0.243177358,0.317391132,0.275942446,
0.158756938,0.077459794,0.175104215,0.221563864,0.293347044,0.290263433,0.025130803)
p2<- c(0.163566501,0.069177271,0.041183757,0.292957706,0.063548779,0.063027447,0.110613414,
0.117951705,0.113279684,0.143697269,0.101504557,0.14956093,0.132989584,
0.079083639,0.047644728,0.091237543,0.122406249,0.148372105,0.135008746,0.017191886)
#Set working directory
setwd("C://Tina/USB_Backup_042213/Paper II/MLN Automation/csvs_equal_20s")
bugs.output <- list()
for (k in 1:1){
Y <- read.csv(file=paste0("MVN",k,".csv"))
Y<-as.matrix(Y)
print(Y)
Y <- ifelse(Y==0,Y+.5,Y)
library("R2WinBUGS")
bugs.output[[1]] <- bugs(data = list(Y=as.matrix(Y), Nf=20, n=60, mn=c(-1.59, -2.44),
prec=matrix(c(.001,0,0,.001),nrow=2,ncol=2),
R=matrix(c(1,0,0,1),nrow=2,ncol=2)),
inits = list(list(gamma=c(0,0), T=matrix(c(0.01,0,0,0.01),nrow=2,ncol=2))),
param = c("p","p1","T","sigma2","rho"),
model = "M-LN_model_ell.txt",
n.iter = 12000, n.burnin = 5000, n.chains = 1, bugs.directory="C://Tina/USB_Backup_042213/winbugs14/WinBUGS14",
working.directory=NULL)
#save mcmc
#write.csv(bugs.output[[i]]$sims.list$ell,paste0("mcmc",i,".csv"))}
# Read CSV into R
#MyData <- read.csv(file="c:/Tina/USB_Backup_042213/Paper II/MLN Automation/csvs_equal_20s/mcmc1.csv", header=TRUE)
for (i in 1:20) # change for number of sites
{
x1<-bugs.output[[1]]$sims.list$p[,i,1]
y1<-bugs.output[[1]]$sims.list$p[,i,2]
## Points within polygons
library(MASS)
dens <- kde2d(x1, y1, n=100) # don't clip the contour
image(dens)
#image(dens, zlim=c(-7,15))
#filled.contour(dens,color.palette=colorRampPalette(c('white','blue','yellow','red','darkred')))
prob <- c(0.975, 0.025)
dx <- diff(dens$x[1:2])
dy <- diff(dens$y[1:2])
sz <- sort(dens$z)
c1 <- cumsum(sz) * dx * dy
levels <- sapply(prob, function(x) {
approx(c1, sz, xout = 1 - x)$y
})
#plot(p1,p2)
contour(dens, level=levels, labels=prob, add=T)
ls <- contourLines(dens, level=levels)
library(sp)
inner <- point.in.polygon(p1[i], p2[i], ls[[2]]$x, ls[[2]]$y)
out <- point.in.polygon(p1[i], p2[i], ls[[1]]$x, ls[[1]]$y)
within<-out-inner
#print(out)
#print(inner)
print(within)
}
}
Y1 Y2
[1,] 27 8
[2,] 8 4
[3,] 2 2
[4,] 36 24
[5,] 5 4
[6,] 10 1
[7,] 14 7
[8,] 12 8
[9,] 12 6
[10,] 22 11
[11,] 11 7
[12,] 24 7
[13,] 13 6
[14,] 6 3
[15,] 2 1
[16,] 10 1
[17,] 12 12
[18,] 11 8
[19,] 16 11
[20,] 2 0
within
[1] 1
[1] -1
[1] 1
[1] 1
[1] -1
[1] -1
[1] 0
[1] 1
[1] -1
[1] 0
[1] 1
[1] 1
[1] 1
[1] 1
[1] 1
[1] -1
[1] 1
[1] 1
[1] 1
[1] 0
|
Locked
