[R] help with bivariate density plot question

Spencer Graves spencer.graves at pdf.com
Sun Jul 13 01:16:57 CEST 2003


I've never used KernSur, but the documentation says it, "returns two 
vectors and a matrix:

xvals vector of ordinates at which the density
	has been estimated in the x dimension
yvals vector of ordinates at which the density
	has been estimated in the y dimension
zden matrix of density for f(x,y) with
	dimensions xgridsize, ygridsize"

	  If we can ignore the probability outside the plot is negligible, then 
the following should produce what you want:

zd <- rev(sort(op$zden))
Surv.zd <- cumsum(zd)
Surv.zd <- (Surv.zd/Surv.zd[lenght(zd)])

plot(zd, Surv.zd)

Now pick the probability levels you want from Surv.zd and translate them 
into levels for zd and use those numbers for the "levels" argument in a 
call to

	 countour(op$xvals, op$yvals, op$zden, levels = ... )

hope this helps.
spencer graves

liping wrote:
 > Dear Spencer Graves:
 >
 > thank you for your help and actually I did something but I
do not know if what I did is fine. I know op$zden is the
estimated density so I  first sort op$zden( I doubt this
step but if i did not sort op$zden, the following step will
go wrong) and then I will sum op$zden to reach the limit i
set up(.68,0.84,.97... ) and output the last value from
op$zden, take it as the contour curve level.  thanks again
for your response.
 >
 > regards,
 >
 > liping
 >
 >
 >
 > xval<-op$xvals
 > yval<-op$yvals
 > z.den<-sort(op$zden)
 >
 > z.distr<-0
 > prob<-c(0.68, 0.84, 0.97, 0.99)
 >
 > for ( j in 1:4) {
 >       for (i in 1:length(op$zden)){
 >           z.distr<-z.den[i]+z.distr
 >           if (z.distr>prob[j]) {
 > 
c<-list(iteration=i,z.distr=z.distr,zden=z.den[i])
 >                         print(c)
 >                         break
 >
 >                              }
 >                  }
 > label<-c("within 68%", "within 84%", "within 97%", "within 99%")
 > contour(op$xvals, op$yvals, op$zden,
 > levels=z.den[i],col=j+5,labels=label[j],add=TRUE)
 > }
 >
 > legend(2,6,legend=c("WITHIN 68%", "WITHIN 84%","WITHIN 97%", "WITHIN 
99%"),
 >                lty=c(1,2,3,4),col=c(6,7,8,9))
########################################################################
Did you look at "?contour"?  The "contour" command has a "levels"
argument.  To translate op into percentage, you could produce an
empirical CDF of op$zden.

hope this helps.  spencer graves

liping wrote:
 > Dear R users:
 >
 > I have a dataset with two variables (>20000 observations, two samples 
from same subject) and I used "kernSur" from library(Genkern) to
 > get a estimated bivariate density and corresponding plots as follows:
 >
 > new.data.normal<-data.normal[!is.na(data.normal[,2]),]
 > x<-new.data.normal[,2]
 > y<-new.data.normal[,3]
 >
 > op <- KernSur(x,y, xgridsize=50, ygridsize=50, correlation=0.4968023,
 >               xbandwidth=1, ybandwidth=1)
 >
 > #3D density plot
 > persp(op$xvals, op$yvals, op$zden,
 >         theta=30,phi=10,expand=0.5,ltheta=120,
 >         xlab="TECH3661.A",ylab="TECH3661.B",zlab="Prob",col="pink",
 >       , main="3D DENSITY PLOT-TECH3661 ", sub=" TECH3661.A AND 
TECH3661.B",
 >         box = T, axes = TRUE,ticktype = "detailed", )
 >
 > #countour plot
 > image(op$xvals, op$yvals, op$zden, col=terrain.colors(100), 
axes=TRUE,xlab="TECH3661.A",ylab="TECH3661.B")
 > points(x,y,pch="*")
 >
 > Now after above step, how can I use 'contour' or other commands to 
draw ellipse curves over above plots indicating "including about 68% 
data", "including about 84% data", etc. similar to the (-std,std), 
(-2*std,2*std),(-3*std, 3*std) intervals for univariate variable.
 >
 > any suggestin will be appreciated.
 >
 > liping
 >
 >
 >
 >
 >
 >
 >
 >
 > 	[[alternative HTML version deleted]]
 >
 > ______________________________________________
 > R-help at stat.math.ethz.ch mailing list
 > https://www.stat.math.ethz.ch/mailman/listinfo/r-help

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help




More information about the R-help mailing list