empirical p-value function

examples/exRand.cpp

@@ -32,9 +32,10 @@ int main(void)
     for (double s = 1.; s < 5.; ++s)
     {
         auto ci = h.confidenceInterval(s);
-        
+
         cout << static_cast<int>(s) << " sigma(s) interval= [";
         cout << ci.first << ", " << ci.second << "]" << endl;
+        cout << "P(X > " << s << ") = " << h.pValue(s) + 1. - h.pValue(-s) << endl;
     }
     p << PlotHistogram(h);
     p << PlotFunction(compile("return exp(-x_0^2/2)/sqrt(2*pi);", 1), -5., 5.);

lib/LatAnalyze/Core/Plot.cpp

@@ -511,6 +511,27 @@ PlotImpulses::PlotImpulses(const DVec &x, const DVec &y)
     setCommand("'" + tmpFileName + "' u 1:2 w impulses");
 }
 
+// PlotSteps constructor ////////////////////////////////////////////////////
+PlotSteps::PlotSteps(const DVec &x, const DVec &y)
+{
+    if (x.rows() != y.rows())
+    {
+        LATAN_ERROR(Size, "x and y vector does not have the same size");
+    }
+
+    DMat   d(x.rows(), 2);
+    string tmpFileName;
+
+    for (Index i = 0; i < x.rows(); ++i)
+    {
+        d(i, 0) = x(i);
+        d(i, 1) = y(i);
+    }
+    tmpFileName = dumpToTmpFile(d);
+    pushTmpFile(tmpFileName);
+    setCommand("'" + tmpFileName + "' u 1:2 w steps");
+}
+
 // PlotGrid constructor ////////////////////////////////////////////////////////
 PlotGrid::PlotGrid(const DVec &x, const DVec &y, const DMat &value)
 {

lib/LatAnalyze/Core/Plot.hpp

@@ -191,6 +191,15 @@ public:
     virtual ~PlotHistogram(void) = default;
 };
 
+class PlotSteps: public PlotObject
+{
+public:
+    // constructor
+    PlotSteps(const DVec &x, const DVec &y);
+    // destructor
+    virtual ~PlotSteps(void) = default;
+};
+
 class PlotImpulses: public PlotObject
 {
 public:

lib/LatAnalyze/Numerical/DWT.cpp

@@ -38,6 +38,10 @@ void filterConvolution(MatType &out, const MatType &data,
 {
     Index n = data.rows(), nf = n*filter.size();
 
+    if (&out == &data)
+    {
+        LATAN_ERROR(Argument, "filter convolution does not support in-place operation");
+    }
     out.resizeLike(data);
     out.fill(0.);
     for (unsigned int i = 0; i < filter.size(); ++i)

lib/LatAnalyze/Statistics/Histogram.cpp

@@ -162,6 +162,20 @@ double Histogram::operator()(const double x) const
     return (*this)[static_cast<Index>(i)];
 }
 
+// p-value P(x > x0) ///////////////////////////////////////////////////////////
+double Histogram::pValue(const double x0) const
+{
+    Index n = data_.size();
+    double count = 0.;
+
+    FOR_STAT_ARRAY(data_, s)
+    {
+        count += (data_[s] > x0) ? 1. : 0.;
+    }
+    
+    return count/n;
+} 
+
 // percentiles & confidence interval ///////////////////////////////////////////
 double Histogram::percentile(const double p) const
 {

lib/LatAnalyze/Statistics/Histogram.hpp

@@ -54,6 +54,8 @@ public:
     double                    getX(const Index i) const;
     double                    operator[](const Index i) const;
     double                    operator()(const double x) const;
+    // p-value P(x > x0)
+    double                    pValue(const double x0) const;
     // percentiles & confidence interval
     double                    percentile(const double p) const;
     double                    median(void) const;