#	data piling example

dp = function(n=c(10,10),m1=rep(0,50),m2=rep(1,50),seed=976) {

    set.seed(seed); p = length(m1); n1 = n[1]; n2 = n[2]

    x1 = matrix( rnorm(n1*p,m1), ncol=p, byrow=T )
    x2 = matrix( rnorm(n2*p,m2), ncol=p, byrow=T )

    W = (n1-1)*var(x1) + (n2-1)*var(x2)

    ev = eigen(W)
    kk = seq(p)[ abs(ev$val)<1.e-10 ]
    A  = ev$vec[,kk]
    dx = x1[1,] - x2[1,]

    ev1 = eigen( t(A) %*% outer(dx,dx) %*% A )

    alpha = ev1$vec[,1]
    beta  = c( A %*% alpha )

    c1 = x1 %*% beta
    c2 = x2 %*% beta

	 ###############

    del = rep(1/sqrt(p),p)
    print( acos( sum(del*beta) ) * 180/pi )

plot(  rbind(x1[,1:2],x2[,1:2]), col=c(rep(1,10),rep(2,10)), pch=16,
       main="1st 2 variables of p=50  n1=n2=10" )
                 text(0.25,-0.8,"MEAN GROUP 1 = 0 0 0 ...",cex=1.5)
  if(m2[1]==0) { text(0.25,-1.3,"MEAN GROUP 2 = 0 0 0 ...",cex=1.5) } else {
                 text(0.25,-1.3,"MEAN GROUP 2 = 1 1 1 ...",cex=1.5) }

browser()

pairs( rbind(x1[,1:4],x2[,1:4]), col=c(rep(1,10),rep(2,10)), pch=16,
       main="1st 4 variables of 50" )

xp = c( rbind(x1,x2) %*% del )
plot( density( xp ),main="Projection onto 1 vector" )
      rug(xp[1:10]); rug(xp[11:20],col=2)

c1p = rep(0,32); c2p = c1p
for(k in 19:50) { v = ev$vec[,k]; c1p[k-18]=mean(x1 %*% v); c2p[k-18]=mean(x2 %*% v) }
plot(range(c1p,c2p),18+c(1,32),type="n",main="Project onto Eigenvectors")
points( c(c1p,c2p), 18+c(1:32,1:32), pch=16, col=2 )
for(k in 1:32) lines( c(c1p[k],c2p[k]), 18+rep(k,2),lty=2 )


plot(19:50,alpha,type="b",pch=16,main="Optimal Weights on 32 Eigenvectors")
  abline(h=0,lty=2)


plot(range(c1p,c2p,c1,c2),18+c(1,32),type="n",main="OPTIMAL PROJECTION")
points( c(c1p,c2p), 18+c(1:32,1:32), pch=16, col=2 )
for(k in 1:32) lines( c(c1p[k],c2p[k]), 18+rep(k,2),lty=2 )
points(c(c1[1],c2[1]),18+rep(16.5,2),pch=16,cex=1.5,col=4)
lines( c(c1[1],c2[1]),18+rep(16.5,2),lwd=2,col=4)



 browser()


}


dp()
dp(m2=rep(0,50))