2022年4月27日星期三

MLE and REML in mixed models

 In mixed models, if use MLE, we have to assume that fixed effects are known, rather than being estimated from the data.  If we use REML (residual maximum likelihood) it will automatically adjust for the degree of freedom corresponding to estimated fixed effects. This is the key difference between MLE and REML in mixed models,   however, the difference will be very small when the sample size is big.

2021年3月1日星期一

character string as function argument r

 https://stackoverflow.com/questions/25989627/character-string-as-function-argument-r

2020年12月20日星期日

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.

2020年11月2日星期一

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