Changeset 730 for trunk/yat/regression/LinearWeighted.cc
- Timestamp:
- Jan 6, 2007, 12:02:21 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/yat/regression/LinearWeighted.cc
r729 r730 26 26 #include "yat/utility/vector.h" 27 27 28 #include <cassert> 29 28 30 namespace theplu { 29 31 namespace yat { … … 32 34 LinearWeighted::LinearWeighted(void) 33 35 : OneDimensionalWeighted(), alpha_(0), alpha_var_(0), beta_(0), 34 beta_var_(0) , m_x_(0), s2_(0)36 beta_var_(0) 35 37 { 36 38 } … … 47 49 double LinearWeighted::alpha_var(void) const 48 50 { 49 return sqrt(alpha_var_);51 return alpha_var_; 50 52 } 51 53 … … 57 59 double LinearWeighted::beta_var(void) const 58 60 { 59 return sqrt(beta_var_);61 return beta_var_; 60 62 } 61 63 … … 64 66 const utility::vector& w) 65 67 { 68 assert(x.size()==y.size()); 69 assert(x.size()==w.size()); 70 66 71 // AveragerPairWeighted requires 2 weights but works only on the 67 72 // product wx*wy, so we can send in w and a dummie to get what we … … 70 75 ap_.add_values(x,y,utility::vector(x.size(),1),w); 71 76 72 // estimating the noise level. see attached document for motivation73 // of the expression.74 s2_= (syy()-sxy()*sxy()/sxx())/(w.sum()-2*(w*w)/w.sum()) ;75 76 77 alpha_ = m_y(); 77 78 beta_ = sxy()/sxx(); 78 alpha_var_ = ap_.y_averager().standard_error() * 79 ap_.y_averager().standard_error(); 80 beta_var_ = s2_/sxx(); 81 m_x_=m_x(); 79 80 chisq_=0; 81 for (size_t i=0; i<x.size(); ++i){ 82 double res = predict(x(i))-y(i); 83 chisq_ += w(i)*res*res; 84 } 85 86 alpha_var_ = s2()/ap_.y_averager().sum_w(); 87 beta_var_ = s2()/sxx(); 82 88 } 83 89 … … 92 98 } 93 99 94 double LinearWeighted::predict ion_error2(const double x, const double w) const95 { 96 return sqrt(alpha_var_ + beta_var_*(x-m_x_)*(x-m_x_)+s2(w));100 double LinearWeighted::predict(const double x) const 101 { 102 return alpha_ + beta_ * (x-m_x()); 97 103 } 98 104 105 99 106 double LinearWeighted::s2(double w) const 100 107 { 101 return s2_/w;108 return chisq_/(w*(ap_.y_averager().n()-2)); 102 109 } 103 110 104 111 double LinearWeighted::standard_error2(const double x) const 105 112 { 106 return sqrt(alpha_var_ + beta_var_*(x-m_x_)*(x-m_x_));113 return alpha_var_ + beta_var_*(x-m_x())*(x-m_x()); 107 114 } 115 108 116 109 117 double LinearWeighted::sxx(void) const … … 112 120 } 113 121 122 114 123 double LinearWeighted::sxy(void) const 115 124 { 116 125 return ap_.sum_xy_centered(); 117 126 } 127 118 128 119 129 double LinearWeighted::syy(void) const
Note: See TracChangeset
for help on using the changeset viewer.