##FUMCTIONS Gamma_hat=function(xi,Z,n=Basics$n,G=Basics$G,B=Basics$B) { ### requires: n=nrow(X); G=solve(t(X)%*% X); B=G%*% t(X)%*% Y #use values if(class(xi)=="numeric") xi=t(xi) rX = c(1+ (n/(n+1)) * xi %*% G%*% t(xi) ) ; #equation 1987 2.9 Gamma_hat= ( S+ (n/((n+1)*rX)) * t ( Z - xi%*% B) %*% ( Z- xi%*% B ) )/ (n+1) Gamma_hat } ## Sundberg's inconsistency function and profile functions Rb=function(Z,xhat,B=Basics$B,Gammahat=Basics$Gammahat) { xhat= solve (B%*% solve(Gammahat)%*% t(B)) %*% B %*% solve(Gammahat)%*% as.matrix(Z);xhat if (class (Z)=="numeric") l=1 else l=nrow(Z) Rb= (Z - xhat%*% B) %*% solve(Gammahat/l) %*% t(Z- xhat%*% B) Rb } #this corresponds to equation 1.3 in Brown and Sundberg 1989 #large value of R attests to internal inconsistency of q components of Z in inference about p components of xhat #xhat above corresponds to eta_hat in 1.3 #also expressed in 3.1 in different notation: Gammahat^{-1} norm norm2=function(x,A) t(x) %*% A %*% as.matrix(x) #Brown profle likelihood given these facts brown.lik=function(x,Z,q=Basics$q,G=Basics$G,n=Basics$n,B=Basics$B,Sinv=Basics$Sinv,detS=Basics$detS) { a= (n+1)/2 *q* log(2*pi) ;a # 154.3817 b= (n+1)/2 * q;b G=c(G) rX = c(1+ (n/(n+1)) * x * G * x ) ;rX #equation 1987 2.9 k= (n/((n+1)*rX)) ;k # 0.9323255 A= sqrt(k) * ( Z - x * B ) ; A det_profile0= detS * c(1 + A %*% Sinv %*% t(A) )/ ( (n+1)^nrow(Sinv)) c= (n+1)/2 * log (det_profile0 );c # 27.83766 brown.lik=a+b-c # -441.6909 brown.lik } #alternative formulation which is being reconciled to brown.lik profile.lik=function(x,Z,X=Basics$X,Y=Basics$Y,n=Basics$n,G=Basics$G,B=Basics$B,l=Basics$l,Splus=Basics$S) { #requires n, G,Y,X #X=Basics$X;Y=Basics$Y;n=Basics$n;G=Basics$G;B=Basics$B;l=Basics$l;Splus=Basics$S xi=as.matrix(x) if (class (Z)=="numeric") l=1 else l=nrow(Z) sigmasq=c(1/l+1/n+xi%*%G%*%t(xi));sigmasq# 0.1570370 with l=9 # 1.045926 with l=1 #Y0=rbind(Y,Z) #defined after 1987 eq 2.2 #3X0=rbind(X,xi); #X0=scale(X0,center= xi/(n+1),scale=FALSE) #defined after following 1987 eq 2.2 #this re-centers X matrix #fm0=lm(Y0~X0) #S0=t(fm0$resid)%*% fm0$resid ;S0 #the RHS of 1987 2.5 #conf Jul 28, 2008 profile.lik=as.numeric( gg=((n+1)/2)*log (sigmasq/ ( sigmasq+( Z-xi%*%B)%*% solve(Splus) %*% t( Z-xi%*%B) ) ) ) ##profile.lik=list(gg=as.numeric(gg),sigmasq=sigmasq,test=Z-xi%*%B,B=B,Splus=Splus,xi=xi,Z=Z) profile.lik} # -6.264439 #this is benchmarked on Z=c(1.68,38.64)-Basics$ybar = -0.06777778 0.7037037 #Zb=t(c(1.68,38.64)-Basics$ybar) #profile.lik(x= -0.4,Zb) # -1.239506 cibrown=function(R,g=Profile$g) { test1=(1+g*(muhat^2+R))^2-4*g*R;test1 # 0.4028476 test2= 1+g*(muhat^2+R) # 1.41522 cibrown=(test2+c(-1,1)*sqrt( test1))/(2*g) # 3.902586 10.249614 names(cibrown)=c("L","U") cibrown} #cibrown(R=4);# 3.468871 11.531129 cstar=function(R) (m+f(muhat,R)) * exp(kgamma/(n+1)) -m # 1987 equation 4.3 cstar.brown=function(R,gamma=.95,cstar0) { #eq 4.4 kgamma=qchisq(gamma, df=p);kgamma # 3.841459 test1=muhat^2-(1-g*cstar0)*(R+muhat^2-cstar0) ;test1 # 1.614040 test2= muhat;test2#1 cstar.brown=(test2+c(-1,1)*sqrt( test1))/(1-g*cstar0) # 3.902586 10.249614 names(cstar.brown)=c("L","U") cstar.brown} brown.basics=function(X,Y,Z0,Xver) { #basics ybar=apply(Y,2,mean);sd0=apply(Y,2,sd) Y=scale(Y,center=ybar,scale=FALSE) xbar=apply(X,2,mean);X=as.matrix(scale(X,center=xbar,scale=FALSE)) n=nrow(Y);l=nrow(Z0);q=ncol(Y);p=ncol(X); c(n,p,q) # 27 2 6 n0=n-p-q #2.11 m=n0*(n+1)/n ;m# 24.88889 fm=lm(Y~X);round(coef(fm),3) alpha=as.matrix(fm$coef[1,]); alpha if (length(xbar)==1) B=t(fm$coef[2,]) else B=fm$coef[2:nrow(fm$coef),] ;round(B,3) E=fm$resid S= t(fm$resid) %*% fm$resid ; Sinv=solve(S) G=solve(t(X)%*% X);round(G,5) # 0.05556 ##Gammahat_inv=solve(S/n) #1987 different for l>1 S(n+l-1) detS=det (S) Gammahat=S/(n-p-q) # equation 2.11 Gammahat # Y1 Y4 #Y1 0.006107407 -0.008109722 #Y4 -0.008109722 1.823430864 brown.basics=list(X=X,Y=Y,xbar=xbar,ybar=ybar,n=n,l=l,m=m,p=p,q=q,n0=n0,fm=fm,alpha=alpha,B=B,S=S,Sinv=Sinv,detS=detS,G=G,Gammahat=Gammahat) brown.basics } #diagram here showing OLS FIT in calibration to Y-vars brown.profile=function(Z0,Basics,Xver){ #GLS FIT IN VERIFCIATION PERIOD ybar=Basics$ybar Z=as.matrix(scale(Z0,center=ybar,scale=FALSE));Z GammahatInv=solve(Basics$Gammahat) xhat= solve (Basics$B%*% GammahatInv%*% t(Basics$B)) %*% Basics$B %*% GammahatInv%*% t(Z);xhat # 0.1717605 H=Basics$B%*% GammahatInv %*% t(Basics$B);H # 4.585796 #simultanesous diagonalization: find algorithm U=1/sqrt(H);U # 0.4669739 #t(U) %*% H %*% U #1 D=(n/(n+1))* t(U)%*% Basics$G %*% U ;D # 0.01168203 #page 50 g=diag(D);g# 0.01168203 kgamma=qchisq(c(.95), df=Basics$p);kgamma#3.841459 #diagram here showing GLS FIT in prediction K=nrow(Z); stat=array(NA,dim=c(K,5));stat=data.frame(stat);names(stat)=c("X","xhat","mle","doublehat","R") stat$xhat=c(xhat);stat stat$X=Xver-Basics$xbar; xrange=3*range(stat$xhat) for (i in 1:K) { test=optimize(brown.lik,xrange,tol=0.001,Z=Z[i,],maximum=TRUE); stat$mle[i]=unlist(test)[1] #bingo stat$R[i]=Rb(Z[i,],xhat[i]) f=function(x,R=stat$R[i],muhat=stat$xhat[i] ) (R+ (x-muhat)^2)/(1+g*x^2) ##post 3.10 test=optimize(f,xrange,tol=0.001,R=stat$R[i]);unlist(test) stat$doublehat[i]=test[[1]]; } brown.profile=list(Z=Z,stat=stat,H=H,U=U,D=D,g=g,kgamma=kgamma,xrange=xrange) brown.profile } cstar=function(R,mu.double) (m+f(mu.double,R)) * exp(kgamma/(n+1)) -m # #eq 4.3 #cstar defined page 54 cibrown=function(R) { test1=(1+g*(muhat^2+R))^2-4*g*R;test1 # 0.4028476 test2= 1+g*(muhat^2+R) # 1.41522 cibrown=(test2+c(-1,1)*sqrt( test1))/(2*g) # 3.902586 10.249614 names(cibrown)=c("L","U") cibrown} #cibrown(R=4);# 3.468871 11.531129 cstar.brown=function(R,cstar0) { #eq 4.4 test1=muhat^2-(1-g*cstar0)*(R+muhat^2-cstar0) ;test1 # 1.614040 test2= muhat;test2#1 cstar.brown=(test2+c(-1,1)*sqrt( test1))/(1-g*cstar0) # 3.902586 10.249614 names(cstar.brown)=c("L","U") cstar.brown} #### ##2.14 profile likelihood #Basics=brown.basics(X,Y,Z0,Xver) #Profile=brown.profile(Z0,Basics,Xver)