Scaling.cpp

/*          
 *             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;
    }
*/

---