[R] Ellipse that Contains 95% of the Observed Data

[Tom La Bone]

Concisely, here is what I am trying to do:

#I take a random sample of 300 measurements. After I have the measurements 
#I post stratify them to 80 type A measurements and 220 type B measurements. 
#These measurements tend to be lognormally distributed so I fit them to 
#determine the geometric mean and geometric standard deviation of each
stratum. 
#The question is: are the geometric mean and geometric standard deviation of 
#the type A measurements the same as the geometric mean and geometric 
#standard deviation of the type B measurements?


library(MASS)
library(car)
setwd("C:/Documents and Settings/Tom/workspace/Work")
source("bagplot.r") 
#http://addictedtor.free.fr/graphiques/RGraphGallery.php?graph=112

#Here are the data. Unknown to me, they are indeed drawn from the same
#lognormal distribution
X <- cbind(1:300,rlnorm(300,log(30),log(3)))
X.A <- X[1:80,2]
X.B <- X[81:300,2]

#These are the data I see. Fit the type A and type B measurements
fit.XA <- fitdistr(X.A,"lognormal")
fit.XB <- fitdistr(X.B,"lognormal")

#Fit 2000 random samples of size 80 and 220 and calculate the difference in
#the GM and GSD for each
x <- numeric(2000)
y <- numeric(2000)
for (i in 1:2000) {
	k <- sample(X[,1],80,replace=FALSE)
	x.a <- X[k,2]
	x.b <- X[setdiff(1:300,k),2]
	fit.xa <- fitdistr(x.a,"lognormal")
	fit.xb <- fitdistr(x.b,"lognormal")
	x[i] <- coef(fit.xa)[1] - coef(fit.xb)[1]
	y[i] <- coef(fit.xa)[2] - coef(fit.xb)[2]
}

#Create the bagplot and superimpose the 95% joint normal confidence ellipse. 
#Does the difference in GM and GSD actually observed for the type A and type
B
#measurements look like the result of the random draw?
bagplot(x,y,show.whiskers=FALSE,approx.limit=1000)
data.ellipse(x,y,plot.points=FALSE,levels=c(0.95),col="black",center.cex=0)
box()
points(coef(fit.XA)[1]-coef(fit.XB)[1],coef(fit.XA)[2]-coef(fit.XB)[2],
		cex=1.5,col="black",pch=19)

	
-- 
View this message in context: http://n4.nabble.com/Ellipse-that-Contains-95-of-the-Observed-Data-tp1694538p1695134.html
Sent from the R help mailing list archive at Nabble.com.