function(d0,v0,C0,M0,y,group)
{J<-max(group)
 N<-numeric(J)
 ybar<-numeric(J)
 Sd<-0
 C1<-C0
 M1<-C0%*%M0
 for (j in 1:J)
     {N[j]<-sum(group==j)
      yj<-y[group==j]
      ybar[j]<-mean(yj)
      Sd<-Sd+sum((yj-ybar[j])^2)
      C1[j,j]<-C1[j,j]+N[j]
      }
 M1<-M1+N*ybar
 M1<-solve(C1,M1)
 R<-t(M0)%*%C0%*%M0+sum(N*ybar*ybar)-t(M1)%*%C1%*%M1
 N<-sum(N)
 Nvd<-Sd+R
 d1<-d0+N
 v1<-(d0*v0+Nvd)/d1
 list(d1=d1,v1=v1,C1=C1,M1=M1)
 }