Revision 306 trunk/contrib/qfclib/qfc_est.cpp

qfc_est.cpp (revision 306)
7 7
  * \ingroup QFC
8 8
  *
9 9
  *  Following user defined functions are more useful for doing estimation model in admb,
10
  *  the easy way to use these functions is put this cpp file in the same directory with your tpl file
11
  *  and under the GLOBALS_SECTION , add a line on top as   #include "qfc_est.cpp" ,
10
  *  the easy way to use these functions is with your tpl file
11
  *  under the GLOBALS_SECTION , add a line on top as   #include "qfclib.h",
12 12
  *  there is a testfunction.tpl file within this folder show how to use these functions.
13 13
  *
14 14
  *  ================  brief list for usage  ============================\n
15 15
  *  1. want to constrain the vector of probability sum as 1 and each as (0,1),  
16 16
  *     see logitProp() and its inverse invLogitProp()\n
17
  *  2. overloaded posfun() for vector and matrix  \n
17
  *  2. overloaded posfun(), boundp() for vector and matrix  \n
18 18
  *  3. constrain parameter within some upper bound only, see mf_upper_bound2(), 
19 19
  *     constrain the sum as 1 use normalize_p() with penalty\n
20 20
  *  4. nll(neg. log. likelihood) for inverse gamma, see nllInverseGamma()  \n 
......
30 30
  *  Quantitative Fisheries Center(QFC), Michigan State University
31 31
  */
32 32

  
33
  #define __QFC_EST__
34
  #if !defined(__QFC_SIM__) //only include the following header file and constants once
35
    #include <math.h>
36
    #include <admodel.h>
37
    #include <df1b2fun.h>
38
    #include <adrndeff.h>
39
    // define constant variable
40
    const double EPS = 1.e-30;          //tiny number to avoid 0 in log
41
  #endif
42 33

  
34
  #include "qfclib.h"
43 35

  
44 36

  
45 37

  
38

  
46 39
  /** constrain probability vector as 1 and 
47 40
   * I forgot who create this first, may give credit to Punt., 
48 41
   * let logit(p)=log(p/(1-p))=a, so p=exp(a)/(1+exp(a)) ~[0,1]
......
61 54
    p(dim)=1./(1.+sum(expa));
62 55
    return p;
63 56
  }
57
  df1b2vector logitProp(const df1b2vector& a)
58
  {
59
    int dim; 
60
    dim=size_count(a)+1;
61
    df1b2vector p(1,dim);
62
    df1b2vector expa=mfexp(a);
63
    p(1,dim-1)=expa/(1.+sum(expa));
64
    //p(dim)=1.-sum(p(1,dim-1));  //original version was buggy for the last term
65
    p(dim)=1./(1.+sum(expa));
66
    return p;
67
  }
64 68

  
69

  
70

  
71

  
65 72
  /** reverse function for LogitProp 
66 73
   * \ingroup QFC
67 74
   * \param p : the prob. vector
......
78 85
    a=lp(1,dim)-lp(dim+1);//each subtract the last one
79 86
    return a;
80 87
  }
88
  df1b2vector invLogitProp(const df1b2vector& p)
89
  {
90
    int dim; 
91
    dim=size_count(p)-1;
92
    df1b2vector a(1,dim); 
93
    df1b2vector lp(1,dim+1);
94
  
95
    lp=log(p+EPS);//take log for each p
96
    a=lp(1,dim)-lp(dim+1);//each subtract the last one
97
    return a;
98
  }
81 99

  
82 100

  
101

  
102

  
103

  
83 104
  /** normailize p as sum(p)=1, return p and penalty if sum(p)!=1 
84 105
   * \ingroup QFC
85 106
   * \param p : the prob. vector
......
93 114
    fpen+=1000.*square(log(psum)); //penalty
94 115
    return p;
95 116
  }
117
  df1b2vector normalize_p(df1b2vector& p, df1b2variable fpen)
118
  {
119
    df1b2variable psum=sum(p);
120
    p/=psum; // Now the p will sum to 1
121
    fpen+=1000.*square(log(psum)); //penalty
122
    return p;
123
  }
96 124

  
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff