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.
______________________________________________
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.
