Dear all, I created a bivariate normal distribution:
set.seed(138813) n<-100 x<-rnorm(n); y<-rnorm(n) and plotted a scatterplot of it: plot(x,y) Now I'd like to add the 2D-standard deviation. I found a thread regarding plotting arbitrary confidence boundaries from Pascal Hänggi http://www.mail-archive.com/r-help@r-project.org/msg24013.html which cites the even older thread http://tolstoy.newcastle.edu.au/R/help/03b/5384.html As I am unfortunately only a very poor R programmer, the code of Pascal Hänggi is a myth to me and I am not sure whether I was able to translate the recommendation of Brain Ripley in the later thread (which provides no code) into the the correct R code. Brain wrote: You need a 2D density estimate (e.g. kde2d in MASS) then compute the density values at the points and draw the contour of the density which includes 95% of the points (at a level computed from the sorted values via quantile()). [95% confidence interval was desired in thread instead of standard deviation...] So I tried this... den<-kde2d(x, y, n=n) #as I chose n to be the same as during creating the distributions x and y (see above), a z-value is assigned to every combination of x and y. # create a sorted vector of z-values (instead of the matrix stored inside the den object den.z <-sort(den$z) # set desired confidence border to draw and store it in variable confidence.border <- quantile(den.z, probs=0.6827, na.rm = TRUE) # draw a line representing confidence.border on the existing scatterplot par(new=TRUE) contour(den, levels=confidence.border, col = "red", add = TRUE) Unfortunately I doubt very much this is correct :( In fact I am sure this is wrong, because the border for probs=0.05 is drawn outside the values.... So please help and check. Pascal Hänggis code seems to work, but I don't understand the magic he does with pp <- array() for (i in 1:1000){ z.x <- max(which(den$x < x[i])) z.y <- max(which(den$y < y[i])) pp[i] <- den$z[z.x, z.y] } before doing the very same as I did above: confidencebound <- quantile(pp, 0.05, na.rm = TRUE) plot(x, y) contour(den, levels = confidencebound, col = "red", add = TRUE) My problems: 1.) setting probs=0.6827 is somehow a dirty trick which I can only use by simply knowing that this is the percentage of values inside +-1sd when a distribution is normal. Is there a way doing this with "native" sd function? sd(den.z) is not correct, as den.z is in contrast to x and y not normal any more. So ecdf(den.z)(sd(den.z)) results in a percentile of 0.5644 in this example instead of the desired 0.6827. 2.) I would like to have code that works with any desired confidence. Unfortunately setting probs to the desired confidence would probably be wrong (?!) as it relates to den.z instead of x and y, which are the underlying distributions I am interested in. To put it short I want the confidence of x/y and not of den.z. I am really completely stuck. Please help me out of this! Felix
______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
