####### Additive Coefficient Estimation Using Polynomial Spline #######





######### need to input X--- the tuning variable matrix: n by d2 matrix 

#########               Tm---- linear variable matix:   n by d1 matrix

#########               Y---- predictor variable vector: length n vector 

#########               p----- order of polynomial spline: e.g. p=1 linear spline, p=3, cubic spline  

######### Output contains  (1) coefficients----- the estimated coefficients

#########                  (2) yhat

########               (3) the estimated functions

                       (4) the estimated constants



addspline<-function(X,Tm,Y,p)

{

n<-length(Y)                 ##### n is the sample size	

N<-ceiling(n^(1/(2*p+3)))	 ###### N is the number of interior knots

d1<-dim(Tm)[2] 

d2<-dim(X)[2]

knot<-matrix(0,d2,(N+1))          ###### knot stors the knot sequences of each direction of X    

for(i in 1:d2)

	knot[i,]<-min(X[,i])+(seq(0:N)-1)*((max(X[,i])-min(X[,i]))/(N+1))



xspline<-rep(1,n)       ##### xspline: additive spline matrix of X

if(p>1)

{

for(j in 1:(p-1))               

xspline<-cbind(xspline,X^j)

}



trun<-((1/2)*((kronecker(t(rep(1,(N+1))),X))-t(kronecker(c(knot),t(rep(1,n))))+abs(kronecker(t(rep(1,(N+1))),X)-t(kronecker(c(knot),t(rep(1,n)))))))^(p)



xspline<-cbind(xspline,trun)



designmatrix<-matrix(0,n,(1+(N+p)*d2)*d1)

for(j in 1:d1)		

designmatrix[,((1+(N+p)*d2)*(j-1)+1):((1+(N+p)*d2)*j)]<-xspline*Tm[,j]



a<-solve(t(designmatrix)%*%designmatrix)%*%t(designmatrix)%*%Y     ##### a: is the estimator of the coefficients.  (vector with length: (1+(N+p)*d2)*d1



yhat<-designmatrix%*%a                                             ##### yhat: gives the predicted value.  (vector with length n)



alphahat<-matrix(0,n,d1*d2)                                        #### alphahat contains the function estimats evaluated at the data points ( n by d1*d2 ) 

for(i in 1:d1)

{

	for(j in 1:d2)

	{	

	ind<-seq(((1+(N+p)*d2)*(i-1)+j+1),((1+(N+p)*d2)*i),d2)	

   alphahat[,(i-1)*d1+j]<-xspline[,seq((j+1),(1+(N+p)*d2),d2)]%*%a[ind]

   }

}



meanvector<-apply(alphahat,2,mean)

alphahat<-alphahat-matrix(rep(apply(alphahat,2,mean),n),n,byrow=T)





chat<-rep(0,d1)                                                  ###### chat contains the estimator of the functions



for(i in 1:d1)

	chat[i]<-a[1+((1+(N+p)*d2)*(i-1))]+sum(meanvector[(1+d1*(i-1)):(d1*i)])



return(list(a,yhat,alphahat,chat))

	

}