Gaussian elimination to solve AX=b

2020-11-03

library(matlib)

## Warning: package 'matlib' was built under R version 4.0.3

A<-matrix(c(2,3,-2,1,2,2,4,2,3),3,3)
b<-c(1,2,3)

gaussian_elimiation<-function(A,b){
  for(J in  1 : nrow(A)){
    C = 1 / A[J,J]
    for(K in 1:nrow(A)){
      A[J,K] = A[J,K] * C;
    }
    b[J] = b[J] * C
    
    for(I in 1 : nrow(A)){
      if(I != J){ C = A[I,J]
      for(K in 1:nrow(A)){
        A[I,K] = A[I,K] - A[J,K] * C;
      }
      b[I] = b[I] - b[J] * C;
      }
      
    }
  } 
  return (list("A"=A,"b"=b))
}

gaussian_elimiation(A,b)

## $A
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
## 
## $b
## [1] -0.19354839  1.25806452  0.03225806

# check with other methods
solve(A)%*%b

##             [,1]
## [1,] -0.19354839
## [2,]  1.25806452
## [3,]  0.03225806

gaussianElimination(A, b)

##      [,1] [,2] [,3]        [,4]
## [1,]    1    0    0 -0.19354839
## [2,]    0    1    0  1.25806452
## [3,]    0    0    1  0.03225806

Iterative re-weighted least squares(IRLS) for Poisson regression with identity link and canonical log link

 

enwul

2020-10-25

y<-c(2,3,6,7,8,9,10,12,15)
x<-matrix(c(1,1,1,1,1,1,1,1,1,-1, -1, 0, 0, 0, 0, 1, 1, 1), nrow=9,ncol=2)
b.init<-c(7,5)
x1<-c(-1, -1, 0, 0, 0, 0, 1, 1, 1)
data_4.3<-data.frame(y,x1)

poisson_IRLS<-function(x,y,b.init,tol=1e-8){
  change <- Inf
  b.old <- b.init
  while(change > tol) {
    eta <- x %*% b.old  # linear predictor
    w<-diag(1/as.vector(eta),nrow=nrow(x),ncol=nrow(x))
    z<-y
    b.new<-solve(t(x)%*%w%*%x)%*%t(x)%*%w%*%z
    change <- sqrt(sum((b.new - b.old)^2))
    b.old <- b.new
  }
  b.new
}

poisson_IRLS(x,y,b.init)

##          [,1]
## [1,] 7.451633
## [2,] 4.935300

m1<-glm(y ~ x1, family=poisson(link="identity"),data=data_4.3)
summary(m1)

## 
## Call:
## glm(formula = y ~ x1, family = poisson(link = "identity"), data = data_4.3)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7019  -0.3377  -0.1105   0.2958   0.7184  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   7.4516     0.8841   8.428  < 2e-16 ***
## x1            4.9353     1.0892   4.531 5.86e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 18.4206  on 8  degrees of freedom
## Residual deviance:  1.8947  on 7  degrees of freedom
## AIC: 40.008
## 
## Number of Fisher Scoring iterations: 3

inv_link<-function(x){exp(x)}

poisson_IRLS_log_link<-function(x,y,b.init,inv_link,tol=1e-8){
  change <- Inf
  b.old <- b.init
  while(change > tol){
    eta <- x %*% b.old  # linear predictor
    y.hat<-inv_link(eta)  #y.hat is mu in fact
    w<-diag(as.vector(y.hat),nrow=nrow(x),ncol=nrow(x))
    z<-solve(w)%*%(y-y.hat)
    b.new<-b.old+solve(t(x)%*%w%*%x)%*%t(x)%*%w%*%z
    change <- sqrt(sum((b.new - b.old)^2))
    b.old <- b.new
  }
  b.new
}

poisson_IRLS_log_link(x,y,rep(1,2),inv_link)

##           [,1]
## [1,] 1.8892720
## [2,] 0.6697856

m2<-glm(y ~ x1, family=poisson(link="identity"),data=data_4.3)
summary(m2)

## 
## Call:
## glm(formula = y ~ x1, family = poisson(link = "identity"), data = data_4.3)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7019  -0.3377  -0.1105   0.2958   0.7184  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   7.4516     0.8841   8.428  < 2e-16 ***
## x1            4.9353     1.0892   4.531 5.86e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 18.4206  on 8  degrees of freedom
## Residual deviance:  1.8947  on 7  degrees of freedom
## AIC: 40.008
## 
## Number of Fisher Scoring iterations: 3

SAS Macro functions – %eval() and %sysevalf()

The %EVAL function evaluates integer arithmetic or logical expressions. %EVAL operates by converting its argument from a character value to a numeric or logical expression. Then, it performs the evaluation. Finally, %EVAL converts the result back to a character value and returns that value.

Operands that contain a period character cause an error when they are part of an integer arithmetic expression. The following examples show correct and incorrect usage, respectively:

%let d=%eval(10+20);     /* Correct usage */

%let d=%eval(10.0+20.0); /* Incorrect usage */

When %EVAL encounters a value containing a period, it displays an error message about finding a character operand where a numeric operand is required.

view source

print?

%let a=1+2;

%let b=10*3;

%let c=5/3;

%let eval_a=%eval(&a);

%let eval_b=%eval(&b);

%let eval_c=%eval(&c);

%put &a is &eval_a;

%put &b is &eval_b;

%put &c is &eval_c; The %SYSEVALF function evaluates arithmetic and logical expressions using floating-point arithmetic and returns a value that is formatted using the BEST32. format. The result of the evaluation is always text. % SYSEVALF is the only macro function that can evaluate logical expressions that contain floating point or missing values. Specifying a conversion type can prevent problems when %SYSEVALF returns missing or floating point values to macro expressions or macro variables that are used in other macro expressions that require an integer value. Syntax: %SYSEVALF(expression<, conversion-type>) Taking an example… view source print? %let a=100; %let b=101.590; %let x=%sysevalf(&a+&b); %let y=%sysevalf(&a +&b, boolean); %let z=%sysevalf(&a +&b, ceil); %let v=%sysevalf(&a +&b, floor); %let w=%sysevalf(&a +&b, int); %put The result with SYSEVALF is: &x; %put  The BOOLEAN value is: &y; %put  The CEIL value is: &z; %put  The FLOOR value is: &v; %put  The INTEGER value is: &w; http://technico.qnownow.com/sas-macro-functions-eval-and-sysevalf-2/