Selected_real_kind

 

Selected_real_kind

Function Type

Selected_real_kind is a transformational function.

Function Purpose

Selected_real_kind gives a programer the ability to choose how many decimal digits of accuracy and the exponent range that he or she needs for a particular purpose. This is an attempt by Fortran 90 to deal with the problem of having different default configurations for the number of decimal digits of precision and exponent range on different computer systems. By giving the programmer the ability to explicitly choose the configuration that is needed to fulfill a task, the accuracy of results should be the same from system to system. There is something that one should be aware of before using this function. There are limitations to what a programmer can request. A programmer can ask for more digits of precision or a greater exponent range then is possible to have on their system.

Declaring variable types with Selected_real_kind

This process is very similar to declaring variable types using selected_int_kind. First thing that needs to be done is to declare a parameter that contains a number that can be used to declare the variables with the desired configuration (KIND). This is done by the following method.

		Integer, parameter :: r14 = selected_real_kind(14,30)

Right away, take note that this intrinsic function and selected_int_kind are the only two intrinsic functions that can be used inside a parameter statement. Next, unlike the selected_int_kind intrinsic function which only has one argument passed to it, this function requires that two arguments be passed to it. The first number after the opening parenthesis is the number of decimal digits of precision that the programmer wants. The next number is then the exponent range that is needed. That means in this example the programmer is looking for fourteen decimal digits of precision and an exponent range from 1.0e-30 up to 1.0e30. Finally, the parameter r14 contains the value needed to declare any real variables with at least fourteen decimal digits of precision and an exponent range of thirty.

Now that r14 is defined, real variables of the required type can be declared. This is simply done in one of two slightly different ways. First, it can be done by typing

		Real(kind=r14)Re,rho,f,roughness

or by typing

		real(r14)Re,rho,f,roughness

There is one cautionary note that any programmer should be aware of when using this function. It is possible to ask the computer for more decimal digits of precision or a greater power of ten than it represent in any configuration. When one or both of these limits are reached, the compiler will load a specific value into the parameter that is being declared. For instance, if a greater degree of precision is asked for then can be represented the parameter will have the value of negative one stored in it. If the exponent range is exceeded, then the value of negative two is stored inside the parameter. Finally, if both the precision and range are exceeded the compiler return the value negative three. To prevent this from screwing up a major program, it is a good idea to write a small test program to test to see if the required set up is possible on your machine. This could potentially save lots of frustration.

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/