Changeset 2568


Ignore:
Timestamp:
Sep 26, 2011, 3:38:30 AM (10 years ago)
Author:
Peter
Message:

refs #671

Location:
trunk/yat/statistics
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/yat/statistics/AveragerPairWeighted.cc

    r2567 r2568  
    3535
    3636  AveragerPairWeighted::AveragerPairWeighted(void)
    37     : wxy_(0)
     37    : wxy_centered_(0)
    3838  {
    3939  }
     
    5151    assert(!std::isnan(wy) && "wy is nan");
    5252    double w=wx*wy;
     53    xy_add(x, y, 0, w);
    5354    x_.add(x,w);
    5455    y_.add(y,w);
    55     wxy_ += w*x*y;
    5656  }
    5757
     
    7171  double AveragerPairWeighted::msd(void) const
    7272  {
    73     return (x_.sum_wxx()+y_.sum_wxx()-2*wxy_)/sum_w();
     73    return (x_.sum_wxx()+y_.sum_wxx()-2*sum_xy())/sum_w();
    7474  }
    7575
     
    8383  void AveragerPairWeighted::reset(void)
    8484  {
    85     x_.reset(); y_.reset(); wxy_=0;
     85    x_.reset(); y_.reset(); wxy_centered_=0;
    8686  }
    8787
     
    9595  double AveragerPairWeighted::sum_xy(void) const
    9696  {
    97     return wxy_;
     97    return sum_xy_centered() + x_.sum_wx()*y_.mean();
    9898  }
    9999
     
    101101  double AveragerPairWeighted::sum_xy_centered(void) const
    102102  {
    103     return sum_xy() - x_.sum_wx()*y_.mean();
     103    return wxy_centered_;
     104  }
     105
     106
     107  void AveragerPairWeighted::xy_add(double mx, double my, double wxy_centered,
     108                                    double w)
     109  {
     110    if (!w)
     111      return;
     112    /*
     113      For derivation of update formula refer to
     114      http://www.janinebennett.org/index_files/ParallelStatisticsAlgorithms.pdf
     115    */
     116    wxy_centered_ += wxy_centered;
     117    if (sum_w())
     118      wxy_centered_ +=
     119        w*sum_w() / (w+sum_w()) * (mx-x_.mean()) * (my-y_.mean());
    104120  }
    105121
  • trunk/yat/statistics/AveragerPairWeighted.h

    r2567 r2568  
    134134    AveragerWeighted x_; // weighted averager with w = w_x*w_y
    135135    AveragerWeighted y_; // weighted averager with w = w_x*w_y
    136     double wxy_;
     136    double wxy_centered_;
     137
     138    void xy_add(double mx, double my, double xy_centered, double w);
    137139  };
    138140
Note: See TracChangeset for help on using the changeset viewer.