Migration from SVN

eo/contrib/mathsym/regression/Scaling.cpp

@@ -0,0 +1,417 @@
+/*	    
+ *             Copyright (C) 2005 Maarten Keijzer
+ *
+ *          This program is free software; you can redistribute it and/or modify
+ *          it under the terms of version 2 of the GNU General Public License as 
+ *          published by the Free Software Foundation. 
+ *
+ *          This program is distributed in the hope that it will be useful,
+ *          but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *          MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *          GNU General Public License for more details.
+ *
+ *          You should have received a copy of the GNU General Public License
+ *          along with this program; if not, write to the Free Software
+ *          Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#include "Scaling.h"
+#include "TargetInfo.h"
+
+using namespace std;
+
+Scaling slope(const std::valarray<double>& x, const TargetInfo& targets) {
+    
+    double xx = 0.0;
+    double xy = 0.0;
+    
+    const valarray<double>& y = targets.targets();
+    
+    for (unsigned i = 0; i < x.size(); ++i) {
+	xx += x[i] * x[i];
+	xy += x[i] * y[i];
+    }
+    
+    if (xx < 1e-7) return Scaling(new LinearScaling(0.0,0.0));
+    
+    double b = xy / xx;
+   
+    return Scaling(new LinearScaling(0.0, b));
+    
+}
+
+// Still needs proper testing with non-trivial lambda
+Scaling regularized_least_squares(const std::valarray<double>& inputs, const TargetInfo& targets, double lambda) {
+    
+    double n = inputs.size();
+    
+    valarray<double> x = inputs;
+
+    double a,b,d;
+    a=b=d=0;
+
+    for (unsigned i = 0; i < n; ++i) {
+	a += 1 + lambda;
+	b += x[i];
+	d += x[i] * x[i] + lambda;
+    }
+    
+    //invert
+    
+    double ad_bc = a*d - b * b;
+    // if ad_bc equals zero there's a problem
+   
+    if (ad_bc < 1e-17) return Scaling(new LinearScaling);
+    
+    double ai = d/ad_bc;
+    double bi = -b/ad_bc;
+    double di = a/ad_bc;
+    double ci = bi;
+    
+    // Now multiply this inverted covariance matrix (C^-1) with x' * t
+    
+    std::valarray<double> ones = x;
+    
+    // calculate C^-1 * x' )
+    for (unsigned i = 0; i < n; ++i) 
+    {
+	ones[i] = (ai + bi * x[i]);
+	x[i]    = (ci + di * x[i]);
+    }
+
+    // results are in [ones, x], now multiply with y
+
+    a = 0.0; // intercept
+    b = 0.0; // slope
+    
+    const valarray<double>& t = targets.targets();
+    
+    for (unsigned i = 0; i < n; ++i)
+    {
+	a += ones[i] * t[i];
+	b += x[i] * t[i];
+    }
+    
+    return Scaling(new LinearScaling(a,b));
+}
+
+Scaling ols(const std::valarray<double>& y, const std::valarray<double>& t) {
+    double n = y.size();
+    
+    double y_mean = y.sum() / n;
+    double t_mean = t.sum() / n;
+   
+    std::valarray<double> y_var = (y - y_mean);
+    std::valarray<double> t_var = (t - t_mean);
+    std::valarray<double> cov = t_var * y_var;
+    
+    y_var *= y_var;
+    t_var *= t_var;
+    
+    double sumvar = y_var.sum();
+    
+    if (sumvar == 0. || sumvar/n < 1e-7 || sumvar/n > 1e+7) // breakout when numerical problems are likely
+	return Scaling(new LinearScaling(t_mean,0.));
+
+    
+    double b = cov.sum() / sumvar;
+    double a = t_mean - b * y_mean;
+    
+    Scaling s = Scaling(new LinearScaling(a,b));
+
+    return s;
+}
+
+Scaling ols(const std::valarray<double>& y, const TargetInfo& targets) {
+    double n = y.size();
+    
+    double y_mean = y.sum() / n;
+    
+    std::valarray<double> y_var = (y - y_mean);
+    std::valarray<double> cov = targets.tcov_part() * y_var;
+    
+    y_var *= y_var;
+
+    double sumvar = y_var.sum();
+    
+    if (sumvar == 0. || sumvar/n < 1e-7 || sumvar/n > 1e+7) // breakout when numerical problems are likely
+	return Scaling(new LinearScaling(targets.tmean(),0.));
+
+    
+    double b = cov.sum() / sumvar;
+    double a = targets.tmean() - b * y_mean;
+    
+    if (!finite(b)) {
+	
+	cout << a << ' ' << b << endl;
+	cout << sumvar << endl;
+	cout << y_mean << endl;
+	cout << cov.sum() << endl;
+	exit(1);
+    }
+	
+    Scaling s = Scaling(new LinearScaling(a,b));
+
+    return s;
+}
+
+
+Scaling wls(const std::valarray<double>& inputs, const TargetInfo& targets) {
+    
+    std::valarray<double> x = inputs;
+    const std::valarray<double>& w = targets.weights();
+    
+    unsigned n = x.size();
+    // First calculate x'*W (as W is a diagonal matrix it's simply elementwise multiplication
+    std::valarray<double> wx = targets.weights() * x;
+    
+    // Now x'*W is contained in [w,wx], calculate x' * W * x (the covariance)
+    double a,b,d;
+    a=b=d=0.0;
+    
+    for (unsigned i = 0; i < n; ++i)
+    {
+	a += w[i];
+	b += wx[i];
+	d += x[i] * wx[i];
+    }
+ 
+    //invert
+    
+    double ad_bc = a*d - b * b;
+    // if ad_bc equals zero there's a problem
+   
+    if (ad_bc < 1e-17) return Scaling(new LinearScaling);
+    
+    double ai = d/ad_bc;
+    double bi = -b/ad_bc;
+    double di = a/ad_bc;
+    double ci = bi;
+    
+    // Now multiply this inverted covariance matrix (C^-1) with x' * W * y
+    
+    // create alias to reuse the wx we do not need anymore
+    std::valarray<double>& ones = wx;
+    
+    // calculate C^-1 * x' * W (using the fact that W is diagonal)
+    for (unsigned i = 0; i < n; ++i) 
+    {
+	ones[i] = w[i]*(ai + bi * x[i]);
+	x[i]    = w[i]*(ci + di * x[i]);
+    }
+
+    // results are in [ones, x], now multiply with y
+
+    a = 0.0; // intercept
+    b = 0.0; // slope
+    
+    const valarray<double>& t = targets.targets();
+    
+    for (unsigned i = 0; i < n; ++i)
+    {
+	a += ones[i] * t[i];
+	b += x[i] * t[i];
+    }
+    
+    return Scaling(new LinearScaling(a,b));
+}
+
+
+//Scaling med(const std::valarray<double>& inputs, const TargetInfo& targets);
+
+double mse(const std::valarray<double>& y, const TargetInfo& t) {
+
+    valarray<double> residuals = t.targets()-y;
+    residuals *= residuals;
+    double sz = residuals.size();
+    if (t.has_weights()) {
+	residuals *= t.weights();
+	sz = 1.0;
+    }
+	
+    return residuals.sum() / sz;
+}
+
+double rms(const std::valarray<double>& y, const TargetInfo& t) {
+    return sqrt(mse(y,t));
+}
+    
+double mae(const std::valarray<double>& y, const TargetInfo& t) {
+    valarray<double> residuals = abs(t.targets()-y);
+    if (t.has_weights()) residuals *= t.weights();
+    return residuals.sum() / residuals.size();
+}
+
+
+/*
+    double standard_error(const std::valarray<double>& y, const std::pair<double,double>& scaling) {
+	double a = scaling.first;
+	double b = scaling.second;
+	double n = y.size();
+	double se = sqrt( pow(a+b*y-current_set->targets,2.0).sum() / (n-2));
+	
+	double mean_y = y.sum() / n;
+	double sxx = pow( y - mean_y, 2.0).sum();
+
+	return se / sqrt(sxx);
+    }
+  
+    double scaled_mse(const std::valarray<double>& y){
+	std::pair<double,double> scaling;
+	return scaled_mse(y,scaling);
+    }
+    
+    double scaled_mse(const std::valarray<double>& y, std::pair<double, double>& scaling)
+    {
+	scaling = scale(y);
+	
+	double a = scaling.first;
+	double b = scaling.second;
+	
+	std::valarray<double> tmp = current_set->targets - a - b * y;
+	tmp *= tmp;
+	
+	if (weights.size())
+	    return (weights * tmp).sum();
+
+	return tmp.sum() / tmp.size();
+    }
+   
+    double robust_mse(const std::valarray<double>& ny, std::pair<double, double>& scaling) {
+	
+	double smse = scaled_mse(ny,scaling);
+
+	std::valarray<double> y = ny;
+	// find maximum covariance case 
+	double n = y.size();
+
+	int largest = 0;
+	
+	{
+	    double y_mean = y.sum() / n;
+	
+	    std::valarray<double> y_var = (y - y_mean);
+	    std::valarray<double> cov = tcov * y_var;
+	
+	    std::valarray<bool> maxcov = cov == cov.max();
+	
+	    for (unsigned i = 0; i < maxcov.size(); ++i) {
+		if (maxcov[i]) {
+		    largest = i;
+		    break;
+		}
+	    }
+	}
+	
+	double y_mean = (y.sum() - y[largest]) / (n-1);
+	y[largest] = y_mean; // dissappears from covariance calculation
+	
+	std::valarray<double> y_var = (y - y_mean);
+	std::valarray<double> cov = tcov * y_var;
+	y_var *= y_var;
+
+	double sumvar = y_var.sum();
+	
+	if (sumvar == 0. || sumvar/n < 1e-7 || sumvar/n > 1e+7) // breakout when numerical problems are likely
+	    return worst_performance();
+	
+	double b = cov.sum() / sumvar;
+	double a = tmean - b * y_mean;
+	
+	std::valarray<double> tmp = current_set->targets - a - b * y;
+	tmp[largest] = 0.0;
+	tmp *= tmp;
+	
+	double smse2 = tmp.sum() / (tmp.size()-1);
+	
+	static std::ofstream os("smse.txt");
+	os << smse << ' ' << smse2 << '\n';
+	
+	if (smse2 > smse) {
+	    return worst_performance();
+	    //std::cerr << "overfit? " << smse << ' ' << smse2 << '\n';
+	}
+	
+	scaling.first = a;
+	scaling.second = b;
+	
+	return smse2;
+    }
+    
+    class Sorter {
+	const std::valarray<double>& scores;
+	public:
+	    Sorter(const std::valarray<double>& _scores) : scores(_scores) {}
+	    
+	    bool operator()(unsigned i, unsigned j) const {
+		return scores[i] < scores[j];
+	    }
+    };
+    
+    double coc(const std::valarray<double>& y) {
+	std::vector<unsigned> indices(y.size());
+	for (unsigned i = 0; i < y.size(); ++i) indices[i] = i;
+	std::sort(indices.begin(), indices.end(), Sorter(y));
+	
+	const std::valarray<double>& targets = current_set->targets;
+	
+	double neg = 1.0 - targets[indices[0]];
+	double pos = targets[indices[0]];
+	
+	double cumpos = 0;
+	double cumneg = 0;
+	double sum=0;
+	
+	double last_score = y[indices[0]];
+	
+	for(unsigned i = 1; i < targets.size(); ++i) {
+	        
+	    if (fabs(y[indices[i]] - last_score) < 1e-9) { // we call it tied
+		pos += targets[indices[i]];
+		neg += 1.0 - targets[indices[i]];
+		
+		if (i < targets.size()-1)
+		    continue;
+	    }
+	    sum += pos * cumneg + (pos * neg) * 0.5;
+	    cumneg += neg;
+	    cumpos += pos;
+	    pos = targets[indices[i]];
+	    neg = 1.0 - targets[indices[i]];
+	    last_score = y[indices[i]];
+	}
+	
+	return sum / (cumneg * cumpos);
+    }
+   
+    // iterative re-weighted least squares (for parameters.classification)
+    double irls(const std::valarray<double>& scores, std::pair<double,double>& scaling) {
+	const std::valarray<double>& t = current_set->targets;
+	
+	std::valarray<double> e(scores.size());
+	std::valarray<double> u(scores.size()); 
+	std::valarray<double> w(scores.size());
+	std::valarray<double> z(scores.size());
+	
+	parameters.use_irls = false; parameters.classification=false;
+	scaling = scale(scores);
+	parameters.use_irls=true;parameters.classification=true;
+	
+	if (scaling.second == 0.0) return worst_performance(); 
+	
+	for (unsigned i = 0; i < 10; ++i) {
+	    e = exp(scaling.first + scaling.second*scores);
+	    u = e / (e + exp(-(scaling.first + scaling.second * scores)));
+	    w = u*(1.-u);
+	    z = (t-u)/w;
+	    scaling = wls(scores, u, w);
+	    //double ll = (log(u)*t + (1.-log(u))*(1.-t)).sum();
+	    //std::cout << "Scale " << i << ' ' << scaling.first << " " << scaling.second << " LL " << 2*ll << std::endl;
+	}
+
+	// log-likelihood
+	u = exp(scaling.first + scaling.second*scores) / (1 + exp(scaling.first + scaling.second*scores));
+	double ll = (log(u)*t + (1.-log(u))*(1.-t)).sum();
+	return 2*ll;
+    }
+*/