Changeset 235

Ignore:
Timestamp:
Feb 21, 2005, 3:53:48 PM (18 years ago)
Message:

Major modifications

Location:
trunk/src
Files:
8 edited

Unmodified
Removed
• trunk/src/Regression.cc

 r221 Regression::Regression() : x_(0.0), y_(0.0), y_err_(0.0) { }
• trunk/src/Regression.h

 r221 #define _theplu_statistics_regression_ // C++ tools include ///////////////////// #include "vector.h" // Standard C++ includes //////////////////////// /// virtual void predict(const double x, double& y, double& y_err, const double w=1) const=0; const double w=1) =0; /// /// @return prediction value and parameters /// virtual std::ostream& print(std::ostream&) const=0; /// /// @return header for print() /// virtual std::ostream& print_header(std::ostream& s) const=0; protected: double x_; double y_; double y_err_; };
• trunk/src/RegressionLinear.cc

 r221 } std::ostream& RegressionLinear::print(std::ostream& s) const { s << x_ << "\t" << y_ << "\t" << y_err_ << "\t" << alpha_ << "\t" << alpha_err() << "\t" << beta_ << "\t" << beta_err() << "\t" << s2_ << "\t" << r2_; return s; } std::ostream& RegressionLinear::print_header(std::ostream& s) const { s << "# column 1: x\n" << "# column 2: y\n" << "# column 3: y_err\n" << "# column 4: alpha\n" << "# column 5: alpha_err\n" << "# column 6: beta\n" << "# column 7: beta_err\n" << "# column 8: s2 (var(y|x))\n" << "# column 9: r2 (coefficient of determination)"; return s; } }} // of namespace statistics and namespace theplu

 r221 /// @return standard deviation of parameter \f$\alpha \f$ /// inline double alpha_var(void) const { return sqrt(alpha_var_); } inline double alpha_err(void) const { return sqrt(alpha_var_); } /// /// @return the parameter \f$\beta /// @return the parameter \f$ \beta \f$/// inline double beta(void) const { return beta_; } /// /// @return standard devaition of parameter \f$ \beta \f$/// @return standard deviation of parameter \f$ \beta \f$/// inline double beta_var(void) const { return sqrt(beta_var_); } inline double beta_err(void) const { return sqrt(beta_var_); } /// /// coefficients \f$ (\alpha, \beta)\f$of the model \f$ y = /// \alpha + \beta (x-m_x) \f$from vectors \a x and \a y, by /// minimizing \f$ \sum{(y_i - \alpha - \beta (x-m_x))^2} \f$. /// minimizing \f$ \sum{(y_i - \alpha - \beta (x-m_x))^2} \f$. By /// construction \f$ \alpha \f$and \f$ \beta \f$are independent. /// void fit(const gslapi::vector& x, const gslapi::vector& y) ; /// coefficients \f$ (\alpha, \beta)\f$of the model \f$ y = /// \alpha + \beta (x-m_x) \f$from vectors \a x and \a y, by /// minimizing \f$ \sum{w_i(y_i - \alpha - \beta (x-m_x))^2} \f$. /// minimizing \f$ \sum{w_i(y_i - \alpha - \beta (x-m_x))^2} \f$, /// where \f$ m_x \f$is the weighted average. By construction \f$ /// \alpha \f$and \f$ \beta \f$are independent. /// void fit(const gslapi::vector& x, const gslapi::vector& y, /// inline void predict(const double x, double& y, double& y_err, const double w=1.0) const const double w=1.0) { y = alpha_ + beta_ * (x-m_x_); y_err = sqrt( alpha_var_+beta_var_*(x-m_x_)*(x-m_x_)+s2_/w ); } /// /// @return prediction value and parameters /// std::ostream& print(std::ostream&) const; /// /// @return header for print() /// std::ostream& print_header(std::ostream&) const; /// • trunk/src/RegressionLocal.cc  r221 namespace statistics { RegressionLocal::RegressionLocal(const gslapi::vector& x, const gslapi::vector& y, Regression& r, RegressionKernel& k, const size_t nof_predictions) : kernel_(&k), regression_(&r) RegressionLocal::RegressionLocal(Regression& r, RegressionKernel& k) : kernel_(&k), regressor_(&r) { estimated_x_ = gslapi::vector(nof_predictions); estimated_y_ = estimated_x_; estimated_y_err_ = estimated_x_; //sorting data with respect to x // Peter, could this be done without all the copying std::vector > data; for (size_t i=0; i tmp(x(i),y(i)); data_.push_back(tmp); } sort(data_.begin(), data_.end()); for (size_t i=0; i((i+0.5)/nof_predictions*data_.size())].first; } void RegressionLocal::fit(const double f) void RegressionLocal::fit(const double f, const u_int step_size) { sort(data_.begin(), data_.end()); // number of points on each side size_t nof_data_points = static_cast(f/2*data_.size()); size_t index_min = 0; size_t index_max = 0; for (size_t i=0; i(f*data_.size()); double width = 0; int min_index = 0; for (size_t i=0; iestimated_x_(i)){ index_min = 0; index_max = 2*nof_data_points-1; if (data_[i].first < (data_[nof_points-1].first+data_[nof_points-1].first)/2 ){ min_index=0; width = data_[nof_points-1].first-data_[i].first; } // extreme right case else if (data_[data_.size()-estimated_x_.size()-1].first< estimated_x_(i)){ index_min = data_.size()-2*nof_data_points; index_max = data_.size()-1; else if (data_[i].first > (data_[data_.size()-1].first+ data_[data_.size()-nof_points].first)/2 ){ min_index = data_.size() - nof_points; width = data_[i].first - data_[data_.size()-nof_points].first; } // normal case - window in the middle else { for (size_t j=nof_data_points; jestimated_x_(i)) { index_min = j-nof_data_points; index_max = j+nof_data_points-1; break; } // stepping forward until x is in the middle while(data_[min_index+nof_points-1].first-data_[i].first < data_[i].first-data_[min_index].first) min_index++; width = data_[min_index+nof_points-1].first-data_[i].first; // Peter, should check if x gets closer to centrum if we step // back one step } // copying data // Peter, this is a stupid solution which could be // improved by using a moving view (in the vector class). gslapi::vector x(index_max-index_min+1); gslapi::vector y(index_max-index_min+1); gslapi::vector w(index_max-index_min+1); // Peter, too much copying. Move the copying outside loop and // use a vector view (when available in gslapi::vector class). gslapi::vector x(nof_points); gslapi::vector y(nof_points); gslapi::vector w(nof_points); for (size_t j=0; j10) std::cout << x(j) << "\t" << y(j) << std::endl; y(j)=data_[min_index+j].second; } // calculating weights for (size_t j=0; jweight( 2*(x(j)-data_[index_min].first)/ (data_[index_max].first- data_[index_min].first)-1); w(j) = kernel_->weight( (x(j)-(data_[min_index].first+ data_[min_index+nof_points-1].first) /2) /width ); // fitting the regressor locally std::cout << "x: " << estimated_x_(i) << std::endl; std::cout << "x.begin: " << x(0) << std::endl; std::cout << "x.end: " << x(x.size()-1) << std::endl; regression_->fit(x,y,w); std::cout << "estimating value\n"; regression_->predict(estimated_x_(i), estimated_y_(i), estimated_y_err_(i)); regressor_->fit(x,y,w); regressor_->predict(x_[i], y_[i], y_err_[i]); } } std::ostream& operator<< (std::ostream& s, const RegressionLocal& r) std::ostream& RegressionLocal::print(std::ostream& s) const { gslapi::vector x=r.estimated_x(); gslapi::vector y(r.estimated_y()); gslapi::vector y_err(r.estimated_y_err()); for (size_t i=0; iprint(s); s << std::endl; } return s; } }} // of namespace cpptools and namespace theplu • trunk/src/RegressionLocal.h  r221 /// /// Constructor loading the object with data, type of regressor, /// type of kernel and in how many points to predict. /// Constructor taking type of \a regressor, /// type of \a kernel. /// RegressionLocal(const gslapi::vector& x, const gslapi::vector& y, Regression& r, RegressionKernel& k, const size_t); RegressionLocal(Regression& regressor, RegressionKernel& kernel); /// virtual ~RegressionLocal(void) {}; /// /// adding a data point /// inline void add(const double x, const double y) { data_.push_back(std::make_pair(x,y)); } /// /// Function returning the points where to predict /// inline gslapi::vector estimated_x(void) const { return estimated_x_; } inline const std::vector& x(void) const { return x_; } /// /// Function returning predicted values /// inline gslapi::vector estimated_y(void) const { return estimated_y_; } inline const std::vector& y(void) const { return y_; } /// /// Function returning error of predictions /// inline gslapi::vector estimated_y_err(void) const {return estimated_y_err_;} inline const std::vector& y_err(void) const { return y_err_; } /// /// Function performing the fit, using a \a fraction of the data /// point and regression method defined in the constructor. The /// algorithm uses equally many points around the point to /// predict. If this is not possible (because the point is too far /// left/right) the points to the extreme left/right is used. /// Performs the fit in data defined by add using a kernel and a /// regression method defined in the constructor. The algorithm /// selects boundaries for the kernel such that \a fraction of the /// data points are used and the point where the fit is done is in /// the middle. Starting with the smallest x, the function jumps /// \a step_size point in each iteration to do the next fit /// void fit(const double fraction); void fit(const double fraction, const u_int step_size=1); /// /// @return prediction values and parameters /// std::ostream& print(std::ostream&) const; /// /// @return header for print() /// inline std::ostream& print_header(std::ostream& s) const { return regressor_->print_header(s); } private: gslapi::vector data_y_; RegressionKernel* kernel_; Regression* regression_; gslapi::vector estimated_x_; gslapi::vector estimated_y_; gslapi::vector estimated_y_err_; Regression* regressor_; std::vector x_; std::vector y_; std::vector y_err_; }; /// /// The output operator for Local Regression Class. /// std::ostream& operator<< (std::ostream& s, const RegressionLocal&); }} // of namespace statistics and namespace theplu • trunk/src/RegressionNaive.cc  r204 //$Id$// Header #include "RegressionNaive.h" // C++_tools #include "Averager.h" #include "Regression.h" #include "WeightedAverager.h" #include RegressionNaive::RegressionNaive(void) : Regression() : Regression(), m_(0.0), m_err_(0.0) { } void RegressionNaive::fit(const gslapi::vector& x, const gslapi::vector& y) { Averager a; for (size_t i=0; i • trunk/src/RegressionNaive.h  r221 //////////////////////// //#include #include #include namespace theplu { /// /// This function computes the best-fit for the naive model \f$ /// y = m \f$from vectors \a x and \a y, by minimizing \f$ /// \sum{(y_i-m)^2} \f$. /// This function computes the best-fit for the naive model \f$ y /// = m \f$from vectors \a x and \a y, by minimizing \f$ /// \sum{(y_i-m)^2} \f\$. This function is the same as using the /// weighted version with unity weights. /// inline void fit(const gslapi::vector& x, const gslapi::vector& y) { Averager a; for (size_t i=0; i
Note: See TracChangeset for help on using the changeset viewer.