#### The Behrens-Fisher example # Log-likelihood function for one normal sample loglik<-function(mu,sig,y){ return(sum(log(dnorm(y,mu,sig)))) } # SD estimate with mu known (one normal sample) sigh<- function(mu,y){ return(sqrt(mean((y-mu)^2))) } # Marginal contour function on (mu1,mu2) for two normal samples (y1,y2) pl2 <- function(mu1,mu2,y1,y2){ n1<-length(y1) n2<-length(y2) L1<-loglik(mu1,sigh(mu1,y1),y1) L2<-loglik(mu2,sigh(mu2,y2),y2) L3<-loglik(mean(y1),sqrt((n1-1)/n1*var(y1)),y1) L4<-loglik(mean(y2),sqrt((n2-1)/n2*var(y2)),y2) return(exp(L1+L2-L3-L4)) } # Marginal contour function on (sig1,sig2) for two normal samples (y1,y2) pl2sig <- function(sig1,sig2,y1,y2){ n1<-length(y1) n2<-length(y2) L1<-loglik(mean(y1),sig1,y1) L2<-loglik(mean(y2),sig2,y2) L3<-loglik(mean(y1),sqrt((n1-1)/n1*var(y1)),y1) L4<-loglik(mean(y2),sqrt((n2-1)/n2*var(y2)),y2) return(exp(L1+L2-L3-L4)) } # plausibility of (mu1,mu1-delta) pl <- function(mu,delta,y1,y2){ return(pl2(mu,mu-delta,y1,y2)) } # The data y1<- c(6.5,6.8,7.1,7.3,10.2) y2<-c (5.8, 5.8, 5.9, 6.0, 6.0, 6.0, 6.3, 6.3, 6.4, 6.5, 6.5) n1<-length(y1) n2<-length(y2) # Contour plot of pl(mu1,mu2) xx<-seq(5.5,9.5,0.01) yy<-seq(5.9,6.4,0.01) nx=length(xx) ny=length(yy) z <- matrix(0,nrow=nx,ncol=ny) for(i in 1:nx){ for(j in 1:ny){ z[i,j]=pl2(xx[i],yy[j],y1,y2) } } contour(x=xx,y=yy,z,xlab=expression(mu[1]),ylab=expression(mu[2])) #,levels=c(-1.5,-1,0,2,10,20,50,100,500,1000, 2000, 3000,5000,7000)) points(mean(y1),mean(y2),pch=3) lines(c(min(c(xx,yy)),max(c(xx,yy))),c(min(c(xx,yy)),max(c(xx,yy))),lty="dashed") # Contour plot of pl(sig1,sig2) xx<-seq(0.2,5.5,0.01) yy<-seq(0.1,0.6,0.01) nx=length(xx) ny=length(yy) z <- matrix(0,nrow=nx,ncol=ny) for(i in 1:nx){ for(j in 1:ny){ z[i,j]=pl2sig(xx[i],yy[j],y1,y2) } } contour(x=xx,y=yy,z,xlab=expression(sigma[1]),ylab=expression(sigma[2]), levels=c(0.01,seq(0.1,0.9,0.1))) points(sqrt((n1-1)/n1*var(y1)),sqrt((n2-1)/n2*var(y2)),pch=3) lines(c(min(c(xx,yy)),max(c(xx,yy))),c(min(c(xx,yy)),max(c(xx,yy))),lty="dashed") # Plot of Pl(H_delta) as a function of delta i<-0 delta<-seq(-0.5,3.5,0.01) N<-length(delta) PL<-vector(length=N) for(d in delta){ i<-i+1 optim <- optimize(pl, c(mean(y1)-4, mean(y1)+4),maximum=TRUE,delta=d,y1=y1,y2=y2) PL[i] <- optim$objective } plot(delta,PL,type="l",xlab=expression(delta),ylab=expression(Pl(mu[1]-mu[2]==delta))) lines(c(0,0),c(0,1),lty="dashed") lines(range(delta),c(0.15,0.15),lty="dotdash") #The threshold corresponding to a 95% confidence level