mrwgmm

来自「LastWave」· 代码 · 共 99 行

#
#
#  Main estimation procedure
#
#
setproc mrwgmm {u {scales <1:20,21:2:40,41:5:100,150>} {&int iterMax 0} {&int tau 1} {lambda2Range 0.01:#30:0.04} {logTRange ln(100):#30:ln(1000)}} \
"{{{<mrw> <lags> [<iterMax>=0] [<tau>=1] [<lambda2Range>] [<lnTRange>]} {Perform GMM parameter estimation of the <mrw>. It is based on the increments of the \
mrw at scale <tau>. It uses correlaton function of the logarithm of these increments at different lags given by the signal <lags> (should not contain 0). \
The GMM function is evaluated for all the possible parameters lambda2 and ln(T) given by a specified range (default is  0.01:#30:0.04 for lambda2 \
and ln(100):#30:ln(1000) for ln(T)) only the parameters with the minimum value is stored. Then dichotomic iterations are performed (maximum number is <iterMax). \
It returns a listv of the form {{ln(sigma) lambda^2 T} covariance_matrix}}}}" {
  
#  sigma = [mrwinit u <scales> tau -s]
#  sigma = 1

  mrwinitgmm u <scales> tau :: >3

  # Starting value for sigma  
  s = 1
  
  # Signal values for lambda2
  lambda2 = lambda2Range 
  
  # Signal values for ln(T)
  T = exp(logTRange)
  
  #
  # Main Loop
  #
  M = 1e20
  last = 1  
  foreach l lambda2 {
    foreach t T {
    
      # Get first iteration estimation
      p = [gmm start {<ln(s),l,ln(t)>}]
      
      # Get quadratic form value for this estimation
      m = [gmm cut p[0] p[1] p[2]] 
      
      # If is minimum then save the estimation in variable P and the quadratic form value in M
      if (m < M && p[2]<ln(20000)) {
        P = [copy p]
        M = m
        if (last == 0) {printf "\n"}
        printf "%g %g %g --> %g\n" exp(p[0]) p[1] exp(p[2]) m 
        last = 1
      } else {
        last = 0
        printf "."
      }
    }
  }
  
  # Get the best estimation in variable p
  if (last == 0) {printf "\n"}
  p = P
  printf "\n"
  
  # Perform iteration if needed
  stop = 0.003
  if (iterMax != 0) {
    foreach iter 1:iterMax {
  
      printf "**** Iteration #%d\n" iter
      m = [gmm matrix P]
      p = [gmm start {P} 1 m]
      printf "%g %g %g\n" exp(p[0]) p[1] exp(p[2]) 
  
      if (abs((p[0]-P[0])/p[0]) < stop &&  abs((p[1]-P[1])/p[1]) < stop && abs((p[2]-P[2])/p[2]) < stop) {
        P = p
        break
      }
    
      P = p
    }
  }
  

  # The result is in P now
  # We get the varfiance of the estimators  
  var1 = [gmm therror P]
  
  # The the standard deviation
  var = 1.96*sqrt(var1[<0,1,2;0,1,2>])
  
  # We print everything
  printf "[%g;%g] [%g;%g] [%g;%g]\n" exp(P[0]-var[0]) exp(P[0]+var[0]) P[1]-var[1] P[1]+var[1] exp(P[2]-var[2]) exp(P[2]+var[2])
  
  # And then we return
  return {P var1}
}