SHARE
TWEET

quadratic loess in C++

a guest May 5th, 2010 663 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <iostream>
  2. #include <iterator>
  3. #include <vector>
  4. #include <cmath>
  5. #include <gsl_multifit.h>
  6. #include "NumericUtils.h"
  7.  
  8. double tricube(double u, double t)
  9. {
  10.     // In our case, u represents a distance and hence should be strictly positive.
  11.     if (u < 0) {
  12.         std::cout << "In tricube(u, t), u < 0" << std::endl;
  13.         std::abort();
  14.     }
  15.  
  16.     // 0 <= u < t
  17.     if ( (closeAtTolerance(u, 0.0) || 0.0 < u) && (u < t) ) {
  18.         // (1 - (u/t)^3)^3
  19.         return pow( ( 1.0 - pow(u/t, 3.0)), 3.0 );
  20.     }
  21.     // u >= t
  22.     else {
  23.         return 0.0;
  24.     }
  25. }
  26.  
  27. double getDistance(double x, double xi)
  28. {
  29.     return fabs(x - xi);
  30. }
  31.  
  32. int main (int argc, char* argv[])
  33. {
  34.     // Loess algorithm parameters.
  35.     double alpha = 0.33;
  36.  
  37.     // Create test dataset.
  38.     const unsigned n = 21;
  39.  
  40.     double xValues[n] = { 0.5578196, 2.0217271, 2.5773252, 3.4140288, 4.3014084, 4.7448394, 5.1073781,
  41.                          6.5411662, 6.7216176, 7.2600583, 8.1335874, 9.1224379, 11.9296663, 12.3797674,
  42.                          13.2728619, 14.2767453, 15.3731026, 15.6476637, 18.5605355, 18.5866354, 18.7572812 };
  43.  
  44.     double yValues[n] = { 18.63654, 103.49646, 150.35391, 190.51031, 208.70115, 213.71135, 228.49353,
  45.                          233.55387, 234.55054, 223.89225, 227.68339, 223.91982, 168.01999, 164.95750,
  46.                          152.61107, 160.78742, 168.55567, 152.42658, 221.70702, 222.69040, 243.18828 };
  47.  
  48.     // Loess neighborhood parameter.
  49.     // XXX: Careful, alpha is assume to be <= 1 here.
  50.     const unsigned q = floor(n * alpha);
  51.  
  52.     // Working arrays. Yes we need both distances and sortedDistances.
  53.     double distances[n];
  54.     double sortedDistances[n];
  55.     double ySmooth[n];
  56.  
  57.     gsl_matrix *X = gsl_matrix_alloc(n, 3);
  58.     gsl_matrix *cov = gsl_matrix_alloc(3, 3);
  59.     gsl_vector *weights = gsl_vector_alloc(n);
  60.     gsl_vector *c = gsl_vector_alloc(3);
  61.     gsl_vector *x = gsl_vector_alloc(3);
  62.  
  63.     double y, yErr, chisq;
  64.    
  65.     // Setup the model matrix X for a quadratic fit.
  66.     for (unsigned i = 0; i < n; ++i) {
  67.         gsl_matrix_set(X, i, 0, 1.0);
  68.         gsl_matrix_set(X, i, 1, xValues[i]);
  69.         gsl_matrix_set(X, i, 2, xValues[i] * xValues[i]);
  70.     }
  71.  
  72.     for (unsigned i = 0; i < n; ++i) {
  73.  
  74.         // Compute distances.
  75.         for (unsigned j = 0; j < n; ++j) {
  76.             distances[j] = getDistance(xValues[i], xValues[j]);
  77.             sortedDistances[j] = distances[j];
  78.         }
  79.  
  80.         // Sort distances in order from smallest to largest.
  81.         std::sort(sortedDistances, sortedDistances + n);
  82.  
  83.         // Compute weights.
  84.         for (unsigned j = 0; j < n; ++j) {
  85.             gsl_vector_set(weights, j, tricube(distances[j], sortedDistances[q]));
  86.         }
  87.  
  88.         gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(n, 3);
  89.         gsl_multifit_wlinear(X, weights, &gsl_vector_view_array(yValues, n).vector, c, cov, &chisq, work);
  90.         gsl_multifit_linear_free(work);
  91.  
  92.         gsl_vector_set(x, 0, 1.0);
  93.         gsl_vector_set(x, 1, xValues[i]);
  94.         gsl_vector_set(x, 2, xValues[i] * xValues[i]);
  95.         gsl_multifit_linear_est (x, c, cov, &y, &yErr);
  96.        
  97.         ySmooth[i] = y;
  98.     }
  99.  
  100.     std::cout.setf(std::ios::fixed, std::ios::floatfield);    
  101.     for (unsigned i = 0; i < n; ++i) {
  102.         std::cout << xValues[i] << " "
  103.                   << yValues[i] << " "
  104.                   << ySmooth[i]
  105.                   << std::endl;
  106.     }
  107.  
  108.     gsl_matrix_free(X);
  109.     gsl_vector_free(weights);
  110.     gsl_vector_free(c);
  111.     gsl_matrix_free(cov);
  112.  
  113.     return 0;
  114. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Top