crossprod function

http://www.endmemo.com/program/R/crossprod.php

R crossprod Function

crossprod() function returns matrix cross-product.

crossprod(x, y = NULL)
tcrossprod(x, y = NULL)

x: numeric matrix
y: numeric matrix, if y=NULL, y is the same as x
...

> x <- matrix(1:9,3,3)
> x

     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

> crossprod(x)

     [,1] [,2] [,3]
[1,]   14   32   50
[2,]   32   77  122
[3,]   50  122  194

> tcrossprod(x)

     [,1] [,2] [,3]
[1,]   66   78   90
[2,]   78   93  108
[3,]   90  108  126

R Tutorials

R Data Types

Loop, Condition Statements

Plotting and Graphics

String Manipulations

Math Functions

Matrix Manipulations

Read & Write Data

Statistical Analysis

All Functions List

:: Popular ::

»	R PCH Symbols
»	R Color Names
»	R Regular Expression
»	R tapply Function
»	R String Functions
»	R Plot Function
»	R Builtin Datasets List
	Python Tutorials
	HTML Tutorials
	JavaScript Tutorials
	Statistics
	News, Events Worldwide
	Unit Conversions
	Top Visited Websites Directory

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)