Schur complement

https://en.wikipedia.org/wiki/Schur_complement

Output variance -covariance matrix from Genmod

proc genmod data=outmi;
   model Oxygen= RunTime RunPulse/covb;
   ods output ParameterEstimates=gmparms
              ParmInfo=gmpinfo
              CovB=gmcovb;
run;

Mean and Variance of epxonential family

https://math.stackexchange.com/a/1634096/166263

Piecewise linear basis, tent function, B-spline in GAM

require(gamair)
data("engine")
engine
attach(engine)
plot(size,wear,xlab='Engine capacity',ylab='Wear index')


tf<-function(x,xj,j){
  dj<-xj*0;
  dj[j]<-1
  approx(xj,dj,x)$y # approx return two subject, only return y here

}
tf.X<-function(x,xj){
  ##tent function basis matrix given data x
  ## and knot sequnce xj
  nk<-length(xj);
  n<-length(x)
  X<-matrix(NA,n,nk)
  for(j in 1:nk){
    X[,j]<-tf(x,xj,j)
  }
  X
}
sj<-seq(min(size),max(size),length=6)
sj
X<-tf.X(size,sj)
X

b<-lm(wear~X-1)
b
s<-seq(min(size),max(size),length=200)
Xp<-tf.X(s,sj)
plot(size,wear)
lines(s,Xp%*%coef(b))

Xp
library(splines)
B=bs(size,knots=sj,Boundary.knots=c(1.420,2.980),degree=1)
B

b2<-lm(wear~B-1)
b2

s2<-seq(min(size),max(size),length=200)
s2
Xp2<-bs(s2,knots=sj,Boundary.knots=c(1.420,2.980),degre=1)
Xp2[,1:6]

plot(size,wear)
lines(s2,Xp2[,1:6]%*%coef(b))

B=bs(s,knots=sj,Boundary.knots=c(1.420,2.980),degree =1)

matplot(s,B,type="l")

Triangular function and B spline

https://en.wikipedia.org/wiki/Triangular_function

https://stats.stackexchange.com/questions/163393/splines-basis-functions-clarification

A function for Gauss-Seidel Iterative methods

Gauss-Seidel Iterative method is essentially the back fitting method for GAM model.

Such as to solve equations

 5x1- 2x2+3x3=-1
-3x1+9x2+  x3= 2
 2x1 - x2 - 7x3= 3


library(MASS) # package needed for ginv() function
Gauss_Seidel<-function(A,b,x0){
 
B<-t(A)%*%A
C<-t(A)%*%b

iter<-100          # number of Gauss-Seidel iterations to run,100 iterations
L<-lower.tri(B)*B# lower triang. B,lower.tri() Returns a matrix of logicals,therefore need to *B
U<-upper.tri(B)*B  # upper triang. B
D<-diag(diag(B))   # diag of A
n<-1
while(n<=iter){
  print(x0)
 
  # Gauss-Seidel formula
  x1<-(ginv((L+D)))%*%((-U%*%x0)+C)
  x0<-x1
  n=n+1
 
}

return(x0)

}

A<-matrix(c(5,-3,2,-2,9,-1,3,1,-7),ncol=3,nrow=3)
b<-matrix(c(-1,2,3),ncol=1,nrow=3)
x0<-t(t(c(0,0,0)))

Gauss_Seidel(A,b,x0)

After 10 iterations:


> Gauss_Seidel(A,b,x0)
     [,1]
[1,]    0
[2,]    0
[3,]    0
           [,1]
[1,] -0.1315789
[2,]  0.1380049
[3,] -0.4007323
            [,1]
[1,] -0.01103352
[2,]  0.23926763
[3,] -0.41380921
            [,1]
[1,]  0.09220576
[2,]  0.28760601
[3,] -0.41850252
           [,1]
[1,]  0.1415692
[2,]  0.3105375
[3,] -0.4207159
           [,1]
[1,]  0.1649876
[2,]  0.3214148
[3,] -0.4217656
           [,1]
[1,]  0.1760960
[2,]  0.3265744
[3,] -0.4222636
           [,1]
[1,]  0.1813651
[2,]  0.3290218
[3,] -0.4224998
           [,1]
[1,]  0.1838645
[2,]  0.3301827
[3,] -0.4226118
           [,1]
[1,]  0.1850501
[2,]  0.3307334
[3,] -0.4226650
           [,1]
[1,]  0.1856124
[2,]  0.3309946
[3,] -0.4226902

#Check with solve function.
solve(A,b)

Convergence Analysis of Iterative Methods

https://www.maa.org/press/periodicals/loci/joma/iterative-methods-for-solving-iaxi-ibi-convergence-analysis-of-iterative-methods