Changeset 1643 for trunk/yat


Ignore:
Timestamp:
Dec 13, 2008, 1:23:39 AM (12 years ago)
Author:
Jari Häkkinen
Message:

Addresses #466 and #425. Added cubic spline interpolation (not complete interface to GSL interpolation stuff).

Location:
trunk/yat/regression
Files:
1 edited
4 copied

Legend:

Unmodified
Added
Removed

  • trunk/yat/regression/CSpline.cc

    r1637 r1643  
    22
    33/*
    4   Copyright (C) 2004, 2005, 2006, 2007 Jari Häkkinen, Peter Johansson
    5   Copyright (C) 2008 Peter Johansson
     4  Copyright (C) 2008 Jari Häkkinen
    65
    76  This file is part of the yat library, http://dev.thep.lu.se/yat
    … …  
    2120*/
    2221
    23 #include "Linear.h"
    24 #include "yat/statistics/AveragerPair.h"
    25 #include "yat/utility/VectorBase.h"
     22#include "CSpline.h"
     23
     24#include <gsl/gsl_interp.h>
    2625
    2726namespace theplu {
    … …  
    2928namespace regression {
    3029
    31   Linear::Linear(void)
    32     : OneDimensional(), alpha_(0), alpha_var_(0), beta_(0), beta_var_(0)
     30
     31  CSpline::CSpline(const utility::VectorConstView& x,
     32                   const utility::VectorConstView& y)
     33    : GSLInterpolation(gsl_interp_cspline,x,y)
    3334  {
    3435  }
    3536
    36   Linear::~Linear(void)
     37
     38  CSpline::~CSpline(void)
    3739  {
    3840  }
    3941
    40   double Linear::alpha(void) const
    41   {
    42     return alpha_;
    43   }
    44 
    45   double Linear::alpha_var(void) const
    46   {
    47     return alpha_var_;
    48   }
    49 
    50   double Linear::beta(void) const
    51   {
    52     return beta_;
    53   }
    54 
    55   double Linear::beta_var(void) const
    56   {
    57     return beta_var_;
    58   }
    59 
    60   void Linear::fit(const utility::VectorBase& x, const utility::VectorBase& y)
    61   {
    62     ap_.reset();
    63     for (size_t i=0; i<x.size(); i++)
    64       ap_.add(x(i),y(i));
    65 
    66     alpha_ = ap_.y_averager().mean();
    67     beta_ = ap_.sum_xy_centered() / ap_.x_averager().sum_xx_centered();
    68 
    69     // calculating deviation between data and model
    70     chisq_ = (ap_.y_averager().sum_xx_centered() - ap_.sum_xy_centered()*
    71               ap_.sum_xy_centered()/ap_.x_averager().sum_xx_centered() );
    72     alpha_var_ = s2() / x.size();
    73     beta_var_ = s2() / ap_.x_averager().sum_xx_centered();
    74   }
    75 
    76   double Linear::predict(const double x) const
    77   {
    78     return alpha_ + beta_ * (x - ap_.x_averager().mean());
    79   }
    80 
    81   double Linear::s2(void) const
    82   {
    83     return chisq()/(ap_.n()-2);
    84   }
    85 
    86   double Linear::standard_error2(const double x) const
    87   {
    88     return alpha_var_+beta_var_*(x-ap_.x_averager().mean())*
    89       (x-ap_.x_averager().mean());
    90   }
    9142
    9243}}} // of namespaces regression, yat, and theplu

  • trunk/yat/regression/CSpline.h

    r1637 r1643  
    1 #ifndef _theplu_yat_regression_linear_
    2 #define _theplu_yat_regression_linear_
     1#ifndef _theplu_yat_regression_cspline_
     2#define _theplu_yat_regression_cspline_
    33
    44// $Id$
    55
    66/*
    7   Copyright (C) 2004, 2005, 2006, 2007 Jari Häkkinen, Peter Johansson
    8   Copyright (C) 2008 Peter Johansson
     7  Copyright (C) 2008 Jari Häkkinen
    98
    109  This file is part of the yat library, http://dev.thep.lu.se/yat
    … …  
    2423*/
    2524
    26 #include "OneDimensional.h"
    27 
    28 #include <cmath>
     25#include "GSLInterpolation.h"
    2926
    3027namespace theplu {
    3128namespace yat {
    32   namespace utility {
    33     class VectorBase;
    34   }
     29namespace utility {
     30  class VectorConstView;
     31}
    3532namespace regression {
    3633
    3734  /**
    38      @brief linear regression.   
    39      
    40      Data are modeled as \f$ y_i = \alpha + \beta (x_i-m_x) +
    41      \epsilon_i \f$.
     35     @brief Documentation please.
    4236  */
    43   class Linear : public OneDimensional
    44   {
    45  
    46   public:
    47     ///
    48     /// @brief The default constructor
    49     ///
    50     Linear(void);
     37  class CSpline : public GSLInterpolation
     38  {
    5139
    52     ///
    53     /// @brief The destructor
    54     ///
    55     virtual ~Linear(void);
    56          
     40  public:
    5741    /**
    58        The parameter \f$ \alpha \f$ is estimated as \f$
    59        \frac{1}{n}\sum y_i \f$
    60        
    61        @return the parameter \f$ \alpha \f$
     42       @brief The default constructor
    6243    */
    63     double alpha(void) const;
     44    CSpline(const utility::VectorConstView& x,
     45            const utility::VectorConstView& y);
    6446
    6547    /**
    66        The standard deviation is estimated as \f$ \sqrt{\frac{s^2}{n}}
    67        \f$ where \f$ s^2 = \frac{\sum \epsilon^2}{n-2} \f$
    68        
    69        @return standard deviation of parameter \f$ \alpha \f$
     48       @brief The destructor
    7049    */
    71     double alpha_var(void) const;
     50    ~CSpline(void);
    7251
     52  private:
    7353    /**
    74        The parameter \f$ \beta \f$ is estimated as \f$
    75        \frac{\textrm{Cov}(x,y)}{\textrm{Var}(x)} \f$
    76        
    77        @return the parameter \f$ \beta \f$
     54       Copy Constructor. (not implemented)
    7855    */
    79     double beta(void) const;
     56    CSpline(const CSpline&);
    8057
    81     /**
    82        The standard deviation is estimated as \f$ \frac{s^2}{\sum
    83        (x-m_x)^2} \f$ where \f$ s^2 = \frac{\sum \epsilon^2}{n-2} \f$
    84 
    85        @return standard deviation of parameter \f$ \beta \f$
    86     */
    87     double beta_var(void) const;
    88 
    89     /**
    90        Model is fitted by minimizing \f$ \sum{(y_i - \alpha - \beta
    91        (x-m_x))^2} \f$. By construction \f$ \alpha \f$ and \f$ \beta \f$
    92        are independent.
    93     */
    94     void fit(const utility::VectorBase& x, const utility::VectorBase& y) ;
    95    
    96     ///
    97     /// @return \f$ \alpha + \beta x \f$
    98     ///
    99     double predict(const double x) const;
    100 
    101     /**
    102        \f$ \frac{\sum \epsilon_i^2}{N-2} \f$
    103 
    104        @return variance of residuals
    105     */
    106     double s2(void) const;
    107 
    108     /**
    109        The error of the model is estimated as \f$
    110        \textrm{alpha\_err}^2+\textrm{beta\_err}^2*(x-m_x)*(x-m_x)\f$
    111    
    112        @return estimated error of model in @a x
    113     */
    114     double standard_error2(const double x) const;
    115 
    116 
    117   private:
    118     ///
    119     /// Copy Constructor. (not implemented)
    120     ///
    121     Linear(const Linear&);
    122 
    123     double alpha_;
    124     double alpha_var_;
    125     double beta_;
    126     double beta_var_;
    12758  };
    12859

  • trunk/yat/regression/GSLInterpolation.cc

    r1637 r1643  
    22
    33/*
    4   Copyright (C) 2004, 2005, 2006, 2007 Jari Häkkinen, Peter Johansson
    5   Copyright (C) 2008 Peter Johansson
     4  Copyright (C) 2008 Jari Häkkinen
    65
    76  This file is part of the yat library, http://dev.thep.lu.se/yat
    … …  
    2120*/
    2221
    23 #include "Linear.h"
    24 #include "yat/statistics/AveragerPair.h"
    25 #include "yat/utility/VectorBase.h"
     22#include "GSLInterpolation.h"
     23
     24#include "yat/utility/VectorConstView.h"
     25
     26#include <cassert>
    2627
    2728namespace theplu {
    … …  
    2930namespace regression {
    3031
    31   Linear::Linear(void)
    32     : OneDimensional(), alpha_(0), alpha_var_(0), beta_(0), beta_var_(0)
     32
     33  GSLInterpolation::GSLInterpolation(const gsl_interp_type* interpolator_type,
     34                                     const utility::VectorConstView& x,
     35                                     const utility::VectorConstView& y)
    3336  {
     37    assert(x.size()==y.size());
     38    size_t size=x.size();
     39    x_=new double[size];
     40    y_=new double[size];
     41    for (size_t i=0; i<size; ++i) {
     42      x_[i]=x(i);
     43      y_[i]=y(i);
     44    }
     45    interpolator_=gsl_interp_alloc (interpolator_type, size);
     46    gsl_interp_init(interpolator_, x_, y_, size);
     47    accelerator_ = gsl_interp_accel_alloc();
    3448  }
    3549
    36   Linear::~Linear(void)
     50
     51  GSLInterpolation::~GSLInterpolation(void)
    3752  {
     53    gsl_interp_free(interpolator_);
     54    gsl_interp_accel_free(accelerator_);
     55    delete x_;
     56    delete y_;
    3857  }
    3958
    40   double Linear::alpha(void) const
    41   {
    42     return alpha_;
     59
     60  double GSLInterpolation::evaluate(const double x) const
     61  {
     62    return gsl_interp_eval(interpolator_,x_,y_,x,accelerator_);
    4363  }
    4464
    45   double Linear::alpha_var(void) const
    46   {
    47     return alpha_var_;
    48   }
    49 
    50   double Linear::beta(void) const
    51   {
    52     return beta_;
    53   }
    54 
    55   double Linear::beta_var(void) const
    56   {
    57     return beta_var_;
    58   }
    59 
    60   void Linear::fit(const utility::VectorBase& x, const utility::VectorBase& y)
    61   {
    62     ap_.reset();
    63     for (size_t i=0; i<x.size(); i++)
    64       ap_.add(x(i),y(i));
    65 
    66     alpha_ = ap_.y_averager().mean();
    67     beta_ = ap_.sum_xy_centered() / ap_.x_averager().sum_xx_centered();
    68 
    69     // calculating deviation between data and model
    70     chisq_ = (ap_.y_averager().sum_xx_centered() - ap_.sum_xy_centered()*
    71               ap_.sum_xy_centered()/ap_.x_averager().sum_xx_centered() );
    72     alpha_var_ = s2() / x.size();
    73     beta_var_ = s2() / ap_.x_averager().sum_xx_centered();
    74   }
    75 
    76   double Linear::predict(const double x) const
    77   {
    78     return alpha_ + beta_ * (x - ap_.x_averager().mean());
    79   }
    80 
    81   double Linear::s2(void) const
    82   {
    83     return chisq()/(ap_.n()-2);
    84   }
    85 
    86   double Linear::standard_error2(const double x) const
    87   {
    88     return alpha_var_+beta_var_*(x-ap_.x_averager().mean())*
    89       (x-ap_.x_averager().mean());
    90   }
    9165
    9266}}} // of namespaces regression, yat, and theplu

  • trunk/yat/regression/GSLInterpolation.h

    r1637 r1643  
    1 #ifndef _theplu_yat_regression_linear_
    2 #define _theplu_yat_regression_linear_
     1#ifndef _theplu_yat_regression_gslinterpolation_
     2#define _theplu_yat_regression_gslinterpolation_
    33
    44// $Id$
    55
    66/*
    7   Copyright (C) 2004, 2005, 2006, 2007 Jari Häkkinen, Peter Johansson
    8   Copyright (C) 2008 Peter Johansson
     7  Copyright (C) 2008 Jari Häkkinen
    98
    109  This file is part of the yat library, http://dev.thep.lu.se/yat
    … …  
    2423*/
    2524
    26 #include "OneDimensional.h"
    27 
    28 #include <cmath>
     25#include <gsl/gsl_interp.h>
    2926
    3027namespace theplu {
    3128namespace yat {
    32   namespace utility {
    33     class VectorBase;
    34   }
     29namespace utility {
     30  class VectorConstView;
     31}
    3532namespace regression {
    3633
    3734  /**
    38      @brief linear regression.   
    39      
    40      Data are modeled as \f$ y_i = \alpha + \beta (x_i-m_x) +
    41      \epsilon_i \f$.
     35     @brief Documentation please.
    4236  */
    43   class Linear : public OneDimensional
    44   {
    45  
    46   public:
    47     ///
    48     /// @brief The default constructor
    49     ///
    50     Linear(void);
     37  class GSLInterpolation
     38  {
    5139
    52     ///
    53     /// @brief The destructor
    54     ///
    55     virtual ~Linear(void);
    56          
    57     /**
    58        The parameter \f$ \alpha \f$ is estimated as \f$
    59        \frac{1}{n}\sum y_i \f$
    60        
    61        @return the parameter \f$ \alpha \f$
    62     */
    63     double alpha(void) const;
     40  public:
    6441
    6542    /**
    66        The standard deviation is estimated as \f$ \sqrt{\frac{s^2}{n}}
    67        \f$ where \f$ s^2 = \frac{\sum \epsilon^2}{n-2} \f$
    68        
    69        @return standard deviation of parameter \f$ \alpha \f$
     43       @brief Documentation please.
    7044    */
    71     double alpha_var(void) const;
     45    double evaluate(const double x) const;
     46
     47  protected:
     48    /**
     49       @brief The default constructor
     50    */
     51    GSLInterpolation(const gsl_interp_type*, const utility::VectorConstView& x,
     52                     const utility::VectorConstView& y);
    7253
    7354    /**
    74        The parameter \f$ \beta \f$ is estimated as \f$
    75        \frac{\textrm{Cov}(x,y)}{\textrm{Var}(x)} \f$
    76        
    77        @return the parameter \f$ \beta \f$
     55       @brief The destructor
    7856    */
    79     double beta(void) const;
     57    virtual ~GSLInterpolation(void);
    8058
     59  private:
    8160    /**
    82        The standard deviation is estimated as \f$ \frac{s^2}{\sum
    83        (x-m_x)^2} \f$ where \f$ s^2 = \frac{\sum \epsilon^2}{n-2} \f$
     61       Copy Constructor. (not implemented)
     62    */
     63    GSLInterpolation(const GSLInterpolation&);
    8464
    85        @return standard deviation of parameter \f$ \beta \f$
    86     */
    87     double beta_var(void) const;
    88 
    89     /**
    90        Model is fitted by minimizing \f$ \sum{(y_i - \alpha - \beta
    91        (x-m_x))^2} \f$. By construction \f$ \alpha \f$ and \f$ \beta \f$
    92        are independent.
    93     */
    94     void fit(const utility::VectorBase& x, const utility::VectorBase& y) ;
    95    
    96     ///
    97     /// @return \f$ \alpha + \beta x \f$
    98     ///
    99     double predict(const double x) const;
    100 
    101     /**
    102        \f$ \frac{\sum \epsilon_i^2}{N-2} \f$
    103 
    104        @return variance of residuals
    105     */
    106     double s2(void) const;
    107 
    108     /**
    109        The error of the model is estimated as \f$
    110        \textrm{alpha\_err}^2+\textrm{beta\_err}^2*(x-m_x)*(x-m_x)\f$
    111    
    112        @return estimated error of model in @a x
    113     */
    114     double standard_error2(const double x) const;
    115 
    116 
    117   private:
    118     ///
    119     /// Copy Constructor. (not implemented)
    120     ///
    121     Linear(const Linear&);
    122 
    123     double alpha_;
    124     double alpha_var_;
    125     double beta_;
    126     double beta_var_;
     65    gsl_interp_accel* accelerator_;
     66    gsl_interp* interpolator_;
     67    double* x_;
     68    double* y_;
    12769  };
    12870

  • trunk/yat/regression/Makefile.am

    r1493 r1643  
    55# Copyright (C) 2006 Jari Häkkinen
    66# Copyright (C) 2007 Jari Häkkinen, Peter Johansson
     7# Copyright (C) 2008 Jari Häkkinen
    78#
    89# This file is part of the yat library, http://dev.thep.lu.se/yat
    … …  
    2223
    2324noinst_LTLIBRARIES = libregression.la
    24 libregression_la_SOURCES = KernelBox.cc KernelTriCube.cc Linear.cc  \
     25libregression_la_SOURCES = CSpline.cc GSLInterpolation.cc KernelBox.cc \
     26  KernelTriCube.cc Linear.cc  \
    2527  LinearWeighted.cc Local.cc MultiDimensional.cc      \
    2628  MultiDimensionalWeighted.cc Naive.cc NaiveWeighted.cc   \
    … …  
    3032include_regressiondir = $(includedir)/yat/regression
    3133
    32 include_regression_HEADERS = Kernel.h KernelBox.h KernelTriCube.h \
     34include_regression_HEADERS = CSpline.h GSLInterpolation.h Kernel.h \
     35  KernelBox.h KernelTriCube.h \
    3336  Linear.h LinearWeighted.h Local.h MultiDimensional.h    \
    3437  MultiDimensionalWeighted.h Naive.h NaiveWeighted.h    \
Note: See TracChangeset for help on using the changeset viewer.