OLS στη C++

Έστω ο πίνακας δεδομένων με δεδομένα x,y:

Να εκτιμήσετε με OLS τους συντελεστές:

$$ y = \beta_0 + \beta_1 x + u $$

#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <fstream>

#define MAX 100000

using namespace std;

double sum      (double *x, unsigned int n);
double sum2     (double *x, unsigned int n);
double average  (double *x, unsigned int n);
double variance (double *x, unsigned int n);
double stdev    (double *x, unsigned int n);
double sumxy    (double *x, double *y, unsigned int n);
double covar    (double *x, double *y, unsigned int n);
double correl   (double *x, double *y, unsigned int n);
double pearson  (double *x, double *y, unsigned int n);
void   ols      (double *x, double *y, unsigned int n);

int main()
{
    ifstream XY;
    XY.open("XYdata.txt");
    
    double x[MAX], y[MAX];
    double sumx;
    unsigned int i = 0, n;
    
    while (XY)
    {
        XY >> x[i];
        XY >> y[i];
        ++i;
    }
    n = --i;
 
    ols(x, y, n);

    return 0;
}

double sum(double *x, unsigned int n)
{
    double sum = 0.0;   
    for (int i=0; i < n; i++)
    {
        sum += x[i];
    }
    return sum; 
}


double sum2(double *x, unsigned int n)
{
    double sum2 = 0.0;   
    for (int i=0; i < n; i++)
    {
        sum2 += pow(x[i], 2.0);
    }
    return sum2; 
}


double sumxy(double *x, double *y, unsigned int n)
{
    double sumxy = 0.0;   
    for (int i=0; i < n; i++)
    {
        sumxy += x[i] * y[i];
    }
    return sumxy; 
}

double average(double *x, unsigned int n)
{
    double ave = sum(x,n) / n;   
    return ave; 
}


double variance(double *x, unsigned int n)
{
    double s  = sum(x, n);
    double s2 = sum2(x, n);   
    double avex = average(x,n);
    double var = 0.0;
    for (int i=0; i < n; i++)
    {
        var += pow( (x[i] - avex), 2.0 );
    }
    var /= (n-1);    
    return var; 
}


double stdev(double *x, unsigned int n)
{
    double var = variance(x,n);
    double std = sqrt(var);
    return std; 
}



double covar(double *x, double *y, unsigned int n)
{
    double avex = average(x,n);
    double avey = average(y,n);
    double covar = 0.0;
    for (int i=0; i < n; i++)
    {
        covar += (x[i] - avex) * (y[i] - avey);
    }
    covar /= n;
    return covar;   
}

double correl(double *x, double *y, unsigned int n)
{
    double cov = covar (x,y,n);
    double stdx = stdev(x,n);
    double stdy = stdev(y,n);
    double correl = cov / (stdx*stdy);
    return correl;
}


double pearson(double *x, double *y, unsigned int n)
{
    double sx = sum (x,n);
    double sy = sum (y,n);
    double s2x = sum2 (x,n);
    double s2y = sum2 (y,n);
    double sxy = sumxy (x,y,n);   
    double pearson = (n*sxy - sx*sy) / 
        sqrt( (n*s2x -pow(sx,2.0) ) * (n*s2y - pow(sy,2.0)) );
    return pearson;
}


void ols(double *x, double *y, unsigned int n)
{
    int i;
    double sy    = sum (y, n);
    double sxy   = sumxy (x, y, n);
    double sx2   = sum2 (x, n);
    double avex  = average (x, n);
    double avey  = average (y, n);
    double r     = correl (x, y, n);
    double b0, b1;       

    b1 = (sxy - n*avex*avey) / (sx2 - n*pow(avex,2.0));
    b0 = avey - b1*avex;

    printf (" correl    : %12.6fn" , r);
    printf (" intercept : %12.6fn" , b0);
    printf (" slope     : %12.6fn" , b1);
}

Μεταγλώττιση και δοκιμαστική εκτέλεση:

astavrak@pegasus:~$ g++ program.cpp
astavrak@pegasus:~$ ./a.out 
 correl    :     0.732554
 intercept :  -138.458784
 slope     :     1.185135

Συνδεθείτε για περισσότερες δυνατότητες αλληλεπίδρασης,
σχολιασμοί, εξωτερικοί σύνδεσμοι, βοήθεια, ψηφοφορίες, αρχεία, κτλ.

Εκπαιδευτικό υλικό από τον Αθανάσιο Σταυρακούδη σας παρέχετε κάτω από την άδεια Creative Commons Attribution-NonCommercial-ShareAlike 4.0 License.
Σας παρακαλώ να ενημερωθείτε για κάποιους επιπλέον περιορισμούς
http://stavrakoudis.econ.uoi.gr/stavrakoudis/?iid=401.

Athanassios Stavrakoudis, Sat 23 Nov 2024