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.e30; //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/(1p))=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,dim1)=expa/(1.+sum(expa)); 

64 
//p(dim)=1.sum(p(1,dim1)); //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 
Also available in: Unified diff