Integral, Measure and Martingale

2019年11月4日星期一

Output variance -covariance matrix from Genmod

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

2019年10月26日星期六

Mean and Variance of epxonential family

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

2019年9月16日星期一

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")

2019年9月15日星期日

Triangular function and B spline

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

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

2019年9月11日星期三

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)