当前位置： 文江博客话题详情

快速& C 中的有效最小二乘拟合算法？

发布于 2024-10-18 20:38:42 字数 157 浏览 5 评论 0原文

我正在尝试对两个数据数组实现线性最小二乘拟合：时间与幅度。到目前为止，我知道的唯一技术是测试 (y = m*x+b) 中所有可能的 m 和 b 点，然后找出最适合我的数据的组合，以便其误差最小。然而，我认为迭代这么多组合有时是没有用的，因为它测试了一切。有什么我不知道的技术可以加快这个过程吗？谢谢。

原文

I am trying to implement a linear least squares fit onto 2 arrays of data: time vs amplitude. The only technique I know so far is to test all of the possible m and b points in (y = m*x+b) and then find out which combination fits my data best so that it has the least error. However, I think iterating so many combinations is sometimes useless because it tests out everything. Are there any techniques to speed up the process that I don't know about? Thanks.

如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论， 你可以免费 注册 一个本站的账号。

评论（8）

生死何惧 2024-10-25 20:38:42

试试这个代码。它适合您的 (x,y) 数据。

linreg 的参数为

linreg(int n, REAL x[], REAL y[], REAL* b, REAL* m, REAL* r)

n = number of data points
x,y  = arrays of data
*b = output intercept
*m  = output slope
*r = output correlation coefficient (can be NULL if you don't want it)

成功时返回值为 0，失败时返回值为 !=0。

以下是代码

#include "linreg.h"
#include <stdlib.h>
#include <math.h>                           /* math functions */

//#define REAL float
#define REAL double


inline static REAL sqr(REAL x) {
    return x*x;
}


int linreg(int n, const REAL x[], const REAL y[], REAL* m, REAL* b, REAL* r){
    REAL   sumx = 0.0;                      /* sum of x     */
    REAL   sumx2 = 0.0;                     /* sum of x**2  */
    REAL   sumxy = 0.0;                     /* sum of x * y */
    REAL   sumy = 0.0;                      /* sum of y     */
    REAL   sumy2 = 0.0;                     /* sum of y**2  */

    for (int i=0;i<n;i++){ 
        sumx  += x[i];       
        sumx2 += sqr(x[i]);  
        sumxy += x[i] * y[i];
        sumy  += y[i];      
        sumy2 += sqr(y[i]); 
    } 

    REAL denom = (n * sumx2 - sqr(sumx));
    if (denom == 0) {
        // singular matrix. can't solve the problem.
        *m = 0;
        *b = 0;
        if (r) *r = 0;
            return 1;
    }

    *m = (n * sumxy  -  sumx * sumy) / denom;
    *b = (sumy * sumx2  -  sumx * sumxy) / denom;
    if (r!=NULL) {
        *r = (sumxy - sumx * sumy / n) /    /* compute correlation coeff */
              sqrt((sumx2 - sqr(sumx)/n) *
              (sumy2 - sqr(sumy)/n));
    }

    return 0; 
}

示例

您可以在线运行此示例。

int main()
{
    int n = 6;
    REAL x[6]= {1, 2, 4,  5,  10, 20};
    REAL y[6]= {4, 6, 12, 15, 34, 68};

    REAL m,b,r;
    linreg(n,x,y,&m,&b,&r);
    printf("m=%g b=%g r=%g\n",m,b,r);
    return 0;
}

这是输出

m=3.43651 b=-0.888889 r=0.999192

这是 Excel 绘图和线性拟合（用于验证）。

所有值与上面的 C 代码完全一致（注意 C 代码返回 r，而 Excel 返回 R**2）。

Try this code. It fits y = mx + b to your (x,y) data.

The arguments to linreg are

linreg(int n, REAL x[], REAL y[], REAL* b, REAL* m, REAL* r)

n = number of data points
x,y  = arrays of data
*b = output intercept
*m  = output slope
*r = output correlation coefficient (can be NULL if you don't want it)

The return value is 0 on success, !=0 on failure.

Here's the code

#include "linreg.h"
#include <stdlib.h>
#include <math.h>                           /* math functions */

//#define REAL float
#define REAL double


inline static REAL sqr(REAL x) {
    return x*x;
}


int linreg(int n, const REAL x[], const REAL y[], REAL* m, REAL* b, REAL* r){
    REAL   sumx = 0.0;                      /* sum of x     */
    REAL   sumx2 = 0.0;                     /* sum of x**2  */
    REAL   sumxy = 0.0;                     /* sum of x * y */
    REAL   sumy = 0.0;                      /* sum of y     */
    REAL   sumy2 = 0.0;                     /* sum of y**2  */

    for (int i=0;i<n;i++){ 
        sumx  += x[i];       
        sumx2 += sqr(x[i]);  
        sumxy += x[i] * y[i];
        sumy  += y[i];      
        sumy2 += sqr(y[i]); 
    } 

    REAL denom = (n * sumx2 - sqr(sumx));
    if (denom == 0) {
        // singular matrix. can't solve the problem.
        *m = 0;
        *b = 0;
        if (r) *r = 0;
            return 1;
    }

    *m = (n * sumxy  -  sumx * sumy) / denom;
    *b = (sumy * sumx2  -  sumx * sumxy) / denom;
    if (r!=NULL) {
        *r = (sumxy - sumx * sumy / n) /    /* compute correlation coeff */
              sqrt((sumx2 - sqr(sumx)/n) *
              (sumy2 - sqr(sumy)/n));
    }

    return 0; 
}

Example

You can run this example online.

int main()
{
    int n = 6;
    REAL x[6]= {1, 2, 4,  5,  10, 20};
    REAL y[6]= {4, 6, 12, 15, 34, 68};

    REAL m,b,r;
    linreg(n,x,y,&m,&b,&r);
    printf("m=%g b=%g r=%g\n",m,b,r);
    return 0;
}

Here is the output

m=3.43651 b=-0.888889 r=0.999192

Here is the Excel plot and linear fit (for verification).

All values agree exactly with the C code above (note C code returns r while Excel returns R**2).

Excel version of fit

回复 原文
云胡 2024-10-25 20:38:42

有有效的最小二乘拟合算法；有关详细信息，请参阅维基百科。还有一些库可以为您实现算法，可能比简单的实现更有效； GNU Scientific Library 就是一个例子，但还有其他一些具有更宽松许可的例子。

There are efficient algorithms for least-squares fitting; see Wikipedia for details. There are also libraries that implement the algorithms for you, likely more efficiently than a naive implementation would do; the GNU Scientific Library is one example, but there are others under more lenient licenses as well.

回复 原文
殤城〤 2024-10-25 20:38:42

这是我的 C/C++ 函数版本，它执行简单的线性回归。计算遵循有关简单线性回归的维基百科文章。这是在 github 上作为单头公共域 (MIT) 库发布的：simple_linear_regression。该库（.h 文件）经过测试，可以在 Linux 和 Windows 上运行，也可以使用 -Wall -Werror 和 clang/gcc 支持的所有 -std 版本从 C 和 C++ 运行。

#define SIMPLE_LINEAR_REGRESSION_ERROR_INPUT_VALUE -2
#define SIMPLE_LINEAR_REGRESSION_ERROR_NUMERIC     -3

int simple_linear_regression(const double * x, const double * y, const int n, double * slope_out, double * intercept_out, double * r2_out) {
    double sum_x = 0.0;
    double sum_xx = 0.0;
    double sum_xy = 0.0;
    double sum_y = 0.0;
    double sum_yy = 0.0;
    double n_real = (double)(n);
    int i = 0;
    double slope = 0.0;
    double denominator = 0.0;

    if (x == NULL || y == NULL || n < 2) {
        return SIMPLE_LINEAR_REGRESSION_ERROR_INPUT_VALUE;
    }

    for (i = 0; i < n; ++i) {
        sum_x += x[i];
        sum_xx += x[i] * x[i];
        sum_xy += x[i] * y[i];
        sum_y += y[i];
        sum_yy += y[i] * y[i];
    }

    denominator = n_real * sum_xx - sum_x * sum_x;
    if (denominator == 0.0) {
        return SIMPLE_LINEAR_REGRESSION_ERROR_NUMERIC;
    }
    slope = (n_real * sum_xy - sum_x * sum_y) / denominator;

    if (slope_out != NULL) {
        *slope_out = slope;
    }

    if (intercept_out != NULL) {
        *intercept_out = (sum_y  - slope * sum_x) / n_real;
    }

    if (r2_out != NULL) {
        denominator = ((n_real * sum_xx) - (sum_x * sum_x)) * ((n_real * sum_yy) - (sum_y * sum_y));
        if (denominator == 0.0) {
            return SIMPLE_LINEAR_REGRESSION_ERROR_NUMERIC;
        }
        *r2_out = ((n_real * sum_xy) - (sum_x * sum_y)) * ((n_real * sum_xy) - (sum_x * sum_y)) / denominator;
    }

    return 0;
}

使用示例：

#define SIMPLE_LINEAR_REGRESSION_IMPLEMENTATION
#include "simple_linear_regression.h"

#include <stdio.h>

/* Some data that we want to find the slope, intercept and r2 for */
static const double x[] = { 1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63, 1.65, 1.68, 1.70, 1.73, 1.75, 1.78, 1.80, 1.83 };
static const double y[] = { 52.21, 53.12, 54.48, 55.84, 57.20, 58.57, 59.93, 61.29, 63.11, 64.47, 66.28, 68.10, 69.92, 72.19, 74.46 };

int main() {
    double slope = 0.0;
    double intercept = 0.0;
    double r2 = 0.0;
    int res = 0;

    res = simple_linear_regression(x, y, sizeof(x) / sizeof(x[0]), &slope, &intercept, &r2);
    if (res < 0) {
        printf("Error: %s\n", simple_linear_regression_error_string(res));
        return res;
    }

    printf("slope: %f\n", slope);
    printf("intercept: %f\n", intercept);
    printf("r2: %f\n", r2);

    return 0;
}

Here is my version of a C/C++ function that does simple linear regression. The calculations follow the wikipedia article on simple linear regression. This is published as a single-header public-domain (MIT) library on github: simple_linear_regression. The library (.h file) is tested to work on Linux and Windows, and from C and C++ using -Wall -Werror and all -std versions supported by clang/gcc.

#define SIMPLE_LINEAR_REGRESSION_ERROR_INPUT_VALUE -2
#define SIMPLE_LINEAR_REGRESSION_ERROR_NUMERIC     -3

int simple_linear_regression(const double * x, const double * y, const int n, double * slope_out, double * intercept_out, double * r2_out) {
    double sum_x = 0.0;
    double sum_xx = 0.0;
    double sum_xy = 0.0;
    double sum_y = 0.0;
    double sum_yy = 0.0;
    double n_real = (double)(n);
    int i = 0;
    double slope = 0.0;
    double denominator = 0.0;

    if (x == NULL || y == NULL || n < 2) {
        return SIMPLE_LINEAR_REGRESSION_ERROR_INPUT_VALUE;
    }

    for (i = 0; i < n; ++i) {
        sum_x += x[i];
        sum_xx += x[i] * x[i];
        sum_xy += x[i] * y[i];
        sum_y += y[i];
        sum_yy += y[i] * y[i];
    }

    denominator = n_real * sum_xx - sum_x * sum_x;
    if (denominator == 0.0) {
        return SIMPLE_LINEAR_REGRESSION_ERROR_NUMERIC;
    }
    slope = (n_real * sum_xy - sum_x * sum_y) / denominator;

    if (slope_out != NULL) {
        *slope_out = slope;
    }

    if (intercept_out != NULL) {
        *intercept_out = (sum_y  - slope * sum_x) / n_real;
    }

    if (r2_out != NULL) {
        denominator = ((n_real * sum_xx) - (sum_x * sum_x)) * ((n_real * sum_yy) - (sum_y * sum_y));
        if (denominator == 0.0) {
            return SIMPLE_LINEAR_REGRESSION_ERROR_NUMERIC;
        }
        *r2_out = ((n_real * sum_xy) - (sum_x * sum_y)) * ((n_real * sum_xy) - (sum_x * sum_y)) / denominator;
    }

    return 0;
}

Usage example:

#define SIMPLE_LINEAR_REGRESSION_IMPLEMENTATION
#include "simple_linear_regression.h"

#include <stdio.h>

/* Some data that we want to find the slope, intercept and r2 for */
static const double x[] = { 1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63, 1.65, 1.68, 1.70, 1.73, 1.75, 1.78, 1.80, 1.83 };
static const double y[] = { 52.21, 53.12, 54.48, 55.84, 57.20, 58.57, 59.93, 61.29, 63.11, 64.47, 66.28, 68.10, 69.92, 72.19, 74.46 };

int main() {
    double slope = 0.0;
    double intercept = 0.0;
    double r2 = 0.0;
    int res = 0;

    res = simple_linear_regression(x, y, sizeof(x) / sizeof(x[0]), &slope, &intercept, &r2);
    if (res < 0) {
        printf("Error: %s\n", simple_linear_regression_error_string(res));
        return res;
    }

    printf("slope: %f\n", slope);
    printf("intercept: %f\n", intercept);
    printf("r2: %f\n", r2);

    return 0;
}
回复 原文
深海夜未眠 2024-10-25 20:38:42

来自数值食谱：科学计算的艺术中的(15.2)将数据拟合为直线：

线性回归：

考虑将一组 N 个数据点 (xi, yi) 拟合到直线模型的问题：

在此处输入图像描述

假设不确定性：sigmai 与每个 yi 相关并且 xi（因变量的值）是准确已知的。为了衡量模型与数据的吻合程度，我们使用卡方函数，在本例中为：

在此处输入图像描述

最小化上述方程以确定 a 和 b。这是通过求上述方程对 a 和 b 的导数，将它们等于零并求解 a 和 b 来完成的。然后我们估计 a 和 b 估计中的可能不确定性，因为显然数据中的测量误差必然会在这些参数的确定中引入一些不确定性。此外，我们必须估计数据与
模型。如果没有这个估计，我们没有丝毫迹象表明模型中的参数 a 和 b 有任何意义。

下面的结构执行上述计算：

struct Fitab {
// Object for fitting a straight line y = a + b*x to a set of 
// points (xi, yi), with or without available
// errors sigma i . Call one of the two constructors to calculate the fit.
// The answers are then available as the variables:
// a, b, siga, sigb, chi2, and either q or sigdat.
int ndata;
double a, b, siga, sigb, chi2, q, sigdat; // Answers.
vector<double> &x, &y, &sig;
// Constructor. 
Fitab(vector<double> &xx, vector<double> &yy, vector<double> &ssig)
    : ndata(xx.size()), x(xx), y(yy), sig(ssig), chi2(0.), q(1.), sigdat(0.) 
  {
    // Given a set of data points x[0..ndata-1], y[0..ndata-1] 
    // with individual standard deviations sig[0..ndata-1], 
    // sets a,b and their respective probable uncertainties
    // siga and sigb, the chi-square: chi2, and the goodness-of-fit
    //  probability: q  
    Gamma gam;
    int i;
    double ss=0., sx=0., sy=0., st2=0., t, wt, sxoss; b=0.0; 

    for (i=0;i < ndata; i++) { // Accumulate sums ...
        wt = 1.0 / SQR(sig[i]); //...with weights
        ss += wt;
        sx += x[i]*wt;
        sy += y[i]*wt;
    }
    sxoss = sx/ss;

    for (i=0; i < ndata; i++) {
        t = (x[i]-sxoss) / sig[i];
        st2 += t*t;
        b += t*y[i]/sig[i];
    }
    b /= st2; // Solve for a, b, sigma-a, and simga-b.
    a = (sy-sx*b) / ss;
    siga = sqrt((1.0+sx*sx/(ss*st2))/ss);
    sigb = sqrt(1.0/st2); // Calculate chi2.
    for (i=0;i<ndata;i++) chi2 += SQR((y[i]-a-b*x[i])/sig[i]);

    if (ndata>2) q=gam.gammq(0.5*(ndata-2),0.5*chi2); // goodness of fit
  }
// Constructor. 
Fitab(vector<double> &xx, vector<double> &yy)
    : ndata(xx.size()), x(xx), y(yy), sig(xx), chi2(0.), q(1.), sigdat(0.) 
  {
    // As above, but without known errors (sig is not used). 
    // The uncertainties siga and sigb are estimated by assuming
    // equal errors for all points, and that a straight line is
    // a good fit. q is returned as 1.0, the normalization of chi2
    // is to unit standard deviation on all points, and sigdat
    // is set to the estimated error of each point.
    int i;
    double ss,sx=0.,sy=0.,st2=0.,t,sxoss;
    b=0.0; // Accumulate sums ...
    for (i=0; i < ndata; i++) {
        sx += x[i]; // ...without weights.
        sy += y[i];
    }
    ss = ndata;
    sxoss = sx/ss;
    for (i=0;i < ndata; i++) {
        t = x[i]-sxoss;
        st2 += t*t;
        b += t*y[i];
    }
    b /= st2;  // Solve for a, b, sigma-a, and sigma-b.
    a = (sy-sx*b)/ss;
    siga=sqrt((1.0+sx*sx/(ss*st2))/ss);
    sigb=sqrt(1.0/st2); // Calculate chi2.
    for (i=0;i<ndata;i++) chi2 += SQR(y[i]-a-b*x[i]);

    if (ndata > 2) sigdat=sqrt(chi2/(ndata-2)); 
    // For unweighted data evaluate typical
    // sig using chi2, and adjust
    // the standard deviations.
    siga *= sigdat;
    sigb *= sigdat;
  }
};

其中 struct Gamma:

struct Gamma : Gauleg18 {  
// Object for incomplete gamma function. 
// Gauleg18 provides coefficients for Gauss-Legendre quadrature.
static const Int ASWITCH=100; When to switch to quadrature method.
static const double EPS;   // See end of struct for initializations.
static const double FPMIN; 
double gln;
double gammp(const double a, const double x) {
    // Returns the incomplete gamma function P(a,x)
    if (x < 0.0 || a <= 0.0) throw("bad args in gammp");
    if (x == 0.0) return 0.0;
    else if ((Int)a >= ASWITCH) return gammpapprox(a,x,1); // Quadrature.
    else if (x < a+1.0) return gser(a,x); // Use the series representation.
    else return 1.0-gcf(a,x); // Use the continued fraction representation.
}

double gammq(const double a, const double x) {
    // Returns the incomplete gamma function Q(a,x) = 1 - P(a,x)
    if (x < 0.0 || a <= 0.0) throw("bad args in gammq");
    if (x == 0.0) return 1.0;
    else if ((Int)a >= ASWITCH) return gammpapprox(a,x,0); // Quadrature.
    else if (x < a+1.0) return 1.0-gser(a,x); // Use the series representation.
    else return gcf(a,x); // Use the continued fraction representation.
}
double gser(const Doub a, const Doub x) {
    // Returns the incomplete gamma function P(a,x) evaluated by its series representation.
    // Also sets ln (gamma) as gln. User should not call directly.
    double sum,del,ap;
    gln=gammln(a);
    ap=a;
    del=sum=1.0/a;
    for (;;) {
        ++ap;
        del *= x/ap;
        sum += del;
        if (fabs(del) < fabs(sum)*EPS) {
            return sum*exp(-x+a*log(x)-gln);
        }
    }
}
double gcf(const Doub a, const Doub x) {
    // Returns the incomplete gamma function Q(a, x) evaluated 
    // by its continued fraction representation. 
    // Also sets ln (gamma) as gln. User should not call directly.
    int i;
    double an,b,c,d,del,h;
    gln=gammln(a);
    b=x+1.0-a; // Set up for evaluating continued fraction
               // by modified Lentz’s method with with b0 = 0.
    c=1.0/FPMIN;
    d=1.0/b;
    h=d;
    for (i=1;;i++) { 
        // Iterate to convergence.
        an = -i*(i-a);
        b += 2.0;
        d=an*d+b;

        if (fabs(d) < FPMIN) d=FPMIN;
        c=b+an/c;
        if (fabs(c) < FPMIN) c=FPMIN;
        d=1.0/d;
        del=d*c;
        h *= del;
        if (fabs(del-1.0) <= EPS) break;
    }
    return exp(-x+a*log(x)-gln)*h; Put factors in front.
}
double gammpapprox(double a, double x, int psig) {
    // Incomplete gamma by quadrature. Returns P(a,x) or Q(a, x), 
    // when psig is 1 or 0, respectively. User should not call directly.
    int j;
    double xu,t,sum,ans;
    double a1 = a-1.0, lna1 = log(a1), sqrta1 = sqrt(a1);
    gln = gammln(a);
    // Set how far to integrate into the tail:
    if (x > a1) xu = MAX(a1 + 11.5*sqrta1, x + 6.0*sqrta1);
    else xu = MAX(0.,MIN(a1 - 7.5*sqrta1, x - 5.0*sqrta1));
    sum = 0;
    for (j=0;j<ngau;j++) { // Gauss-Legendre.
        t = x + (xu-x)*y[j];
        sum += w[j]*exp(-(t-a1)+a1*(log(t)-lna1));
    }
    ans = sum*(xu-x)*exp(a1*(lna1-1.)-gln);
    return (psig?(ans>0.0? 1.0-ans:-ans):(ans>=0.0? ans:1.0+ans));
}
double invgammp(Doub p, Doub a);
// Inverse function on x of P(a,x) .
};
const Doub Gamma::EPS = numeric_limits<Doub>::epsilon();
const Doub Gamma::FPMIN = numeric_limits<Doub>::min()/EPS

和 stuct Gauleg18:

struct Gauleg18 {
// Abscissas and weights for Gauss-Legendre quadrature.
   static const Int ngau = 18;
   static const Doub y[18];
   static const Doub w[18];
};

const Doub Gauleg18::y[18] = {0.0021695375159141994,
   0.011413521097787704,0.027972308950302116,0.051727015600492421,
   0.082502225484340941, 0.12007019910960293,0.16415283300752470,
   0.21442376986779355, 0.27051082840644336, 0.33199876341447887,
   0.39843234186401943, 0.46931971407375483, 0.54413605556657973,
   0.62232745288031077, 0.70331500465597174, 0.78649910768313447,
   0.87126389619061517, 0.95698180152629142};

const Doub Gauleg18::w[18] = {0.0055657196642445571,
   0.012915947284065419,0.020181515297735382,0.027298621498568734,
   0.034213810770299537,0.040875750923643261,0.047235083490265582,
   0.053244713977759692,0.058860144245324798,0.064039797355015485
   0.068745323835736408,0.072941885005653087,0.076598410645870640,
   0.079687828912071670,0.082187266704339706,0.084078218979661945,
   0.085346685739338721,0.085983275670394821};

以及最后的 fuinction Gamma::invgamp():

double Gamma::invgammp(double p, double a) { 
    // Returns x such that P(a,x) =  p for an argument p between 0 and 1.
    int j;
    double x,err,t,u,pp,lna1,afac,a1=a-1;
    const double EPS=1.e-8; // Accuracy is the square of EPS.
    gln=gammln(a);
    if (a <= 0.) throw("a must be pos in invgammap");
    if (p >= 1.) return MAX(100.,a + 100.*sqrt(a));
    if (p <= 0.) return 0.0;
    if (a > 1.) {   
        lna1=log(a1);
        afac = exp(a1*(lna1-1.)-gln);
        pp = (p < 0.5)? p : 1. - p;
        t = sqrt(-2.*log(pp));
        x = (2.30753+t*0.27061)/(1.+t*(0.99229+t*0.04481)) - t;
        if (p < 0.5) x = -x;
        x = MAX(1.e-3,a*pow(1.-1./(9.*a)-x/(3.*sqrt(a)),3));
   } else {  
        t = 1.0 - a*(0.253+a*0.12); and (6.2.9).
        if (p < t) x = pow(p/t,1./a);
        else x = 1.-log(1.-(p-t)/(1.-t));
   }
   for (j=0;j<12;j++) {
        if (x <= 0.0) return 0.0; // x too small to compute accurately.
        err = gammp(a,x) - p;
        if (a > 1.) t = afac*exp(-(x-a1)+a1*(log(x)-lna1));
        else t = exp(-x+a1*log(x)-gln);
        u = err/t;
        // Halley’s method.
        x -= (t = u/(1.-0.5*MIN(1.,u*((a-1.)/x - 1)))); 
        // Halve old value if x tries to go negative.
        if (x <= 0.) x = 0.5*(x + t); 
        if (fabs(t) < EPS*x ) break;
    }
    return x;
}

From Numerical Recipes: The Art of Scientific Computing in (15.2) Fitting Data to a Straight Line:

Linear Regression:

Consider the problem of fitting a set of N data points (xi, yi) to a straight-line model:

enter image description here

Assume that the uncertainty: sigmai associated with each yi and that the xi’s (values of the dependent variable) are known exactly. To measure how well the model agrees with the data, we use the chi-square function, which in this case is:

enter image description here

The above equation is minimized to determine a and b. This is done by finding the derivative of the above equation with respect to a and b, equate them to zero and solve for a and b. Then we estimate the probable uncertainties in the estimates of a and b, since obviously the measurement errors in the data must introduce some uncertainty in the determination of those parameters. Additionally, we must estimate the goodness-of-fit of the data to the
model. Absent this estimate, we have not the slightest indication that the parameters a and b in the model have any meaning at all.

The below struct performs the mentioned calculations:

struct Fitab {
// Object for fitting a straight line y = a + b*x to a set of 
// points (xi, yi), with or without available
// errors sigma i . Call one of the two constructors to calculate the fit.
// The answers are then available as the variables:
// a, b, siga, sigb, chi2, and either q or sigdat.
int ndata;
double a, b, siga, sigb, chi2, q, sigdat; // Answers.
vector<double> &x, &y, &sig;
// Constructor. 
Fitab(vector<double> &xx, vector<double> &yy, vector<double> &ssig)
    : ndata(xx.size()), x(xx), y(yy), sig(ssig), chi2(0.), q(1.), sigdat(0.) 
  {
    // Given a set of data points x[0..ndata-1], y[0..ndata-1] 
    // with individual standard deviations sig[0..ndata-1], 
    // sets a,b and their respective probable uncertainties
    // siga and sigb, the chi-square: chi2, and the goodness-of-fit
    //  probability: q  
    Gamma gam;
    int i;
    double ss=0., sx=0., sy=0., st2=0., t, wt, sxoss; b=0.0; 

    for (i=0;i < ndata; i++) { // Accumulate sums ...
        wt = 1.0 / SQR(sig[i]); //...with weights
        ss += wt;
        sx += x[i]*wt;
        sy += y[i]*wt;
    }
    sxoss = sx/ss;

    for (i=0; i < ndata; i++) {
        t = (x[i]-sxoss) / sig[i];
        st2 += t*t;
        b += t*y[i]/sig[i];
    }
    b /= st2; // Solve for a, b, sigma-a, and simga-b.
    a = (sy-sx*b) / ss;
    siga = sqrt((1.0+sx*sx/(ss*st2))/ss);
    sigb = sqrt(1.0/st2); // Calculate chi2.
    for (i=0;i<ndata;i++) chi2 += SQR((y[i]-a-b*x[i])/sig[i]);

    if (ndata>2) q=gam.gammq(0.5*(ndata-2),0.5*chi2); // goodness of fit
  }
// Constructor. 
Fitab(vector<double> &xx, vector<double> &yy)
    : ndata(xx.size()), x(xx), y(yy), sig(xx), chi2(0.), q(1.), sigdat(0.) 
  {
    // As above, but without known errors (sig is not used). 
    // The uncertainties siga and sigb are estimated by assuming
    // equal errors for all points, and that a straight line is
    // a good fit. q is returned as 1.0, the normalization of chi2
    // is to unit standard deviation on all points, and sigdat
    // is set to the estimated error of each point.
    int i;
    double ss,sx=0.,sy=0.,st2=0.,t,sxoss;
    b=0.0; // Accumulate sums ...
    for (i=0; i < ndata; i++) {
        sx += x[i]; // ...without weights.
        sy += y[i];
    }
    ss = ndata;
    sxoss = sx/ss;
    for (i=0;i < ndata; i++) {
        t = x[i]-sxoss;
        st2 += t*t;
        b += t*y[i];
    }
    b /= st2;  // Solve for a, b, sigma-a, and sigma-b.
    a = (sy-sx*b)/ss;
    siga=sqrt((1.0+sx*sx/(ss*st2))/ss);
    sigb=sqrt(1.0/st2); // Calculate chi2.
    for (i=0;i<ndata;i++) chi2 += SQR(y[i]-a-b*x[i]);

    if (ndata > 2) sigdat=sqrt(chi2/(ndata-2)); 
    // For unweighted data evaluate typical
    // sig using chi2, and adjust
    // the standard deviations.
    siga *= sigdat;
    sigb *= sigdat;
  }
};

where struct Gamma:

struct Gamma : Gauleg18 {  
// Object for incomplete gamma function. 
// Gauleg18 provides coefficients for Gauss-Legendre quadrature.
static const Int ASWITCH=100; When to switch to quadrature method.
static const double EPS;   // See end of struct for initializations.
static const double FPMIN; 
double gln;
double gammp(const double a, const double x) {
    // Returns the incomplete gamma function P(a,x)
    if (x < 0.0 || a <= 0.0) throw("bad args in gammp");
    if (x == 0.0) return 0.0;
    else if ((Int)a >= ASWITCH) return gammpapprox(a,x,1); // Quadrature.
    else if (x < a+1.0) return gser(a,x); // Use the series representation.
    else return 1.0-gcf(a,x); // Use the continued fraction representation.
}

double gammq(const double a, const double x) {
    // Returns the incomplete gamma function Q(a,x) = 1 - P(a,x)
    if (x < 0.0 || a <= 0.0) throw("bad args in gammq");
    if (x == 0.0) return 1.0;
    else if ((Int)a >= ASWITCH) return gammpapprox(a,x,0); // Quadrature.
    else if (x < a+1.0) return 1.0-gser(a,x); // Use the series representation.
    else return gcf(a,x); // Use the continued fraction representation.
}
double gser(const Doub a, const Doub x) {
    // Returns the incomplete gamma function P(a,x) evaluated by its series representation.
    // Also sets ln (gamma) as gln. User should not call directly.
    double sum,del,ap;
    gln=gammln(a);
    ap=a;
    del=sum=1.0/a;
    for (;;) {
        ++ap;
        del *= x/ap;
        sum += del;
        if (fabs(del) < fabs(sum)*EPS) {
            return sum*exp(-x+a*log(x)-gln);
        }
    }
}
double gcf(const Doub a, const Doub x) {
    // Returns the incomplete gamma function Q(a, x) evaluated 
    // by its continued fraction representation. 
    // Also sets ln (gamma) as gln. User should not call directly.
    int i;
    double an,b,c,d,del,h;
    gln=gammln(a);
    b=x+1.0-a; // Set up for evaluating continued fraction
               // by modified Lentz’s method with with b0 = 0.
    c=1.0/FPMIN;
    d=1.0/b;
    h=d;
    for (i=1;;i++) { 
        // Iterate to convergence.
        an = -i*(i-a);
        b += 2.0;
        d=an*d+b;

        if (fabs(d) < FPMIN) d=FPMIN;
        c=b+an/c;
        if (fabs(c) < FPMIN) c=FPMIN;
        d=1.0/d;
        del=d*c;
        h *= del;
        if (fabs(del-1.0) <= EPS) break;
    }
    return exp(-x+a*log(x)-gln)*h; Put factors in front.
}
double gammpapprox(double a, double x, int psig) {
    // Incomplete gamma by quadrature. Returns P(a,x) or Q(a, x), 
    // when psig is 1 or 0, respectively. User should not call directly.
    int j;
    double xu,t,sum,ans;
    double a1 = a-1.0, lna1 = log(a1), sqrta1 = sqrt(a1);
    gln = gammln(a);
    // Set how far to integrate into the tail:
    if (x > a1) xu = MAX(a1 + 11.5*sqrta1, x + 6.0*sqrta1);
    else xu = MAX(0.,MIN(a1 - 7.5*sqrta1, x - 5.0*sqrta1));
    sum = 0;
    for (j=0;j<ngau;j++) { // Gauss-Legendre.
        t = x + (xu-x)*y[j];
        sum += w[j]*exp(-(t-a1)+a1*(log(t)-lna1));
    }
    ans = sum*(xu-x)*exp(a1*(lna1-1.)-gln);
    return (psig?(ans>0.0? 1.0-ans:-ans):(ans>=0.0? ans:1.0+ans));
}
double invgammp(Doub p, Doub a);
// Inverse function on x of P(a,x) .
};
const Doub Gamma::EPS = numeric_limits<Doub>::epsilon();
const Doub Gamma::FPMIN = numeric_limits<Doub>::min()/EPS

and stuct Gauleg18:

struct Gauleg18 {
// Abscissas and weights for Gauss-Legendre quadrature.
   static const Int ngau = 18;
   static const Doub y[18];
   static const Doub w[18];
};

const Doub Gauleg18::y[18] = {0.0021695375159141994,
   0.011413521097787704,0.027972308950302116,0.051727015600492421,
   0.082502225484340941, 0.12007019910960293,0.16415283300752470,
   0.21442376986779355, 0.27051082840644336, 0.33199876341447887,
   0.39843234186401943, 0.46931971407375483, 0.54413605556657973,
   0.62232745288031077, 0.70331500465597174, 0.78649910768313447,
   0.87126389619061517, 0.95698180152629142};

const Doub Gauleg18::w[18] = {0.0055657196642445571,
   0.012915947284065419,0.020181515297735382,0.027298621498568734,
   0.034213810770299537,0.040875750923643261,0.047235083490265582,
   0.053244713977759692,0.058860144245324798,0.064039797355015485
   0.068745323835736408,0.072941885005653087,0.076598410645870640,
   0.079687828912071670,0.082187266704339706,0.084078218979661945,
   0.085346685739338721,0.085983275670394821};

and, finally fuinction Gamma::invgamp():

double Gamma::invgammp(double p, double a) { 
    // Returns x such that P(a,x) =  p for an argument p between 0 and 1.
    int j;
    double x,err,t,u,pp,lna1,afac,a1=a-1;
    const double EPS=1.e-8; // Accuracy is the square of EPS.
    gln=gammln(a);
    if (a <= 0.) throw("a must be pos in invgammap");
    if (p >= 1.) return MAX(100.,a + 100.*sqrt(a));
    if (p <= 0.) return 0.0;
    if (a > 1.) {   
        lna1=log(a1);
        afac = exp(a1*(lna1-1.)-gln);
        pp = (p < 0.5)? p : 1. - p;
        t = sqrt(-2.*log(pp));
        x = (2.30753+t*0.27061)/(1.+t*(0.99229+t*0.04481)) - t;
        if (p < 0.5) x = -x;
        x = MAX(1.e-3,a*pow(1.-1./(9.*a)-x/(3.*sqrt(a)),3));
   } else {  
        t = 1.0 - a*(0.253+a*0.12); and (6.2.9).
        if (p < t) x = pow(p/t,1./a);
        else x = 1.-log(1.-(p-t)/(1.-t));
   }
   for (j=0;j<12;j++) {
        if (x <= 0.0) return 0.0; // x too small to compute accurately.
        err = gammp(a,x) - p;
        if (a > 1.) t = afac*exp(-(x-a1)+a1*(log(x)-lna1));
        else t = exp(-x+a1*log(x)-gln);
        u = err/t;
        // Halley’s method.
        x -= (t = u/(1.-0.5*MIN(1.,u*((a-1.)/x - 1)))); 
        // Halve old value if x tries to go negative.
        if (x <= 0.) x = 0.5*(x + t); 
        if (fabs(t) < EPS*x ) break;
    }
    return x;
}
回复 原文
天邊彩虹 2024-10-25 20:38:42

上面的原始示例对我来说在斜率和偏移方面效果很好，但我在校正系数方面遇到了困难。也许我的括号的作用与假定的优先级不同？不管怎样，在其他网页的帮助下，我终于得到了与 Excel 中的线性趋势线相匹配的值。我想我会使用 Mark Lakata 的变量名来分享我的代码。希望这有帮助。

double slope = ((n * sumxy) - (sumx * sumy )) / denom;
double intercept = ((sumy * sumx2) - (sumx * sumxy)) / denom;
double term1 = ((n * sumxy) - (sumx * sumy));
double term2 = ((n * sumx2) - (sumx * sumx));
double term3 = ((n * sumy2) - (sumy * sumy));
double term23 = (term2 * term3);
double r2 = 1.0;
if (fabs(term23) > MIN_DOUBLE)  // Define MIN_DOUBLE somewhere as 1e-9 or similar
    r2 = (term1 * term1) / term23;

The original example above worked well for me with slope and offset but I had a hard time with the corr coef. Maybe I don't have my parenthesis working the same as the assumed precedence? Anyway, with some help of other web pages I finally got values that match the linear trend-line in Excel. Thought I would share my code using Mark Lakata's variable names. Hope this helps.

double slope = ((n * sumxy) - (sumx * sumy )) / denom;
double intercept = ((sumy * sumx2) - (sumx * sumxy)) / denom;
double term1 = ((n * sumxy) - (sumx * sumy));
double term2 = ((n * sumx2) - (sumx * sumx));
double term3 = ((n * sumy2) - (sumy * sumy));
double term23 = (term2 * term3);
double r2 = 1.0;
if (fabs(term23) > MIN_DOUBLE)  // Define MIN_DOUBLE somewhere as 1e-9 or similar
    r2 = (term1 * term1) / term23;
回复 原文
转身泪倾城 2024-10-25 20:38:42

作为一项作业，我必须使用 RMSE 损失函数用 C 语言编写一个简单的线性回归。该程序是动态的，您可以输入自己的值并选择自己的损失函数，目前仅限于均方根误差。但首先是我使用的算法：

在此处输入图像描述
输入图片此处描述
输入图片此处描述

现在代码...你需要gnuplot来显示图表，sudo apt install gnuplot

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/types.h>

#define BUFFSIZE 64
#define MAXSIZE 100

static double vector_x[MAXSIZE] = {0};
static double vector_y[MAXSIZE] = {0};
static double vector_predict[MAXSIZE] = {0};

static double max_x;
static double max_y;
static double mean_x;
static double mean_y;
static double teta_0_intercept;
static double teta_1_grad;
static double RMSE;
static double r_square;
static double prediction;

static char intercept[BUFFSIZE];
static char grad[BUFFSIZE];
static char xrange[BUFFSIZE];
static char yrange[BUFFSIZE];
static char lossname_RMSE[BUFFSIZE] = "Simple Linear Regression using RMSE'";

static char cmd_gnu_0[BUFFSIZE] = "set title '";
static char cmd_gnu_1[BUFFSIZE] = "intercept = ";
static char cmd_gnu_2[BUFFSIZE] = "grad = ";
static char cmd_gnu_3[BUFFSIZE] = "set xrange [0:";
static char cmd_gnu_4[BUFFSIZE] = "set yrange [0:";
static char cmd_gnu_5[BUFFSIZE] = "f(x) = (grad * x) + intercept";
static char cmd_gnu_6[BUFFSIZE] = "plot f(x), 'data.temp' with points pointtype 7";

static char const *commands_gnuplot[] = {
    cmd_gnu_0,
    cmd_gnu_1,
    cmd_gnu_2,
    cmd_gnu_3,
    cmd_gnu_4,
    cmd_gnu_5,
    cmd_gnu_6,
};

static size_t size;

static void user_input()
{
    printf("Enter x,y vector size, MAX = 100\n");
    scanf("%lu", &size);
    if (size > MAXSIZE) {
        printf("Wrong input size is too big\n");
        user_input();
    }
    printf("vector's size is %lu\n", size);

    size_t i;
    for (i = 0; i < size; i++) {
        printf("Enter vector_x[%ld] values\n", i);
        scanf("%lf", &vector_x[i]);
    }

    for (i = 0; i < size; i++) {
        printf("Enter vector_y[%ld] values\n", i);
        scanf("%lf", &vector_y[i]);
    }
}

static void display_vector()
{
    size_t i;
    for (i = 0; i < size; i++){
        printf("vector_x[%lu] = %lf\t", i, vector_x[i]);
        printf("vector_y[%lu] = %lf\n", i, vector_y[i]);
    }
}

static void concatenate(char p[], char q[]) {
   int c;
   int d;
   c = 0;

   while (p[c] != '\0') {
      c++;
   }
   d = 0;

   while (q[d] != '\0') {
      p[c] = q[d];
      d++;
      c++;
   }
   p[c] = '\0';
}

static void compute_mean_x_y()
{
    size_t i;
    double tmp_x = 0.0;
    double tmp_y = 0.0;
    for (i = 0; i < size; i++) {
        tmp_x += vector_x[i];
        tmp_y += vector_y[i];
    }

    mean_x = tmp_x / size;
    mean_y = tmp_y / size;

    printf("mean_x = %lf\n", mean_x);
    printf("mean_y = %lf\n", mean_y);
}

static void compute_teta_1_grad()
{
    double numerator = 0.0;
    double denominator = 0.0;
    double tmp1 = 0.0;
    double tmp2 = 0.0;
    size_t i;

    for (i = 0; i < size; i++) {
        numerator += (vector_x[i] - mean_x) * (vector_y[i] - mean_y);
    }

    for (i = 0; i < size; i++) {
        tmp1 = vector_x[i] - mean_x;
        tmp2 = tmp1 * tmp1;
        denominator += tmp2;
    }

    teta_1_grad = numerator / denominator;
    printf("teta_1_grad = %lf\n", teta_1_grad);
}

static void compute_teta_0_intercept()
{
    teta_0_intercept = mean_y - (teta_1_grad * mean_x);
    printf("teta_0_intercept = %lf\n", teta_0_intercept);
}

static void compute_prediction()
{
    size_t i;
    for (i = 0; i < size; i++) {
        vector_predict[i] = teta_0_intercept + (teta_1_grad * vector_x[i]);
        printf("y^[%ld] = %lf\n", i, vector_predict[i]);
    }
    printf("\n");
}

static void compute_RMSE()
{
    compute_prediction();
    double error = 0;
    size_t i;
    for (i = 0; i < size; i++) {
        error = (vector_predict[i] - vector_y[i]) * (vector_predict[i] - vector_y[i]);
        printf("error y^[%ld] =  %lf\n", i, error);
        RMSE += error;
    }
    /* mean */
    RMSE = RMSE / size;
    /* square root mean */
    RMSE = sqrt(RMSE);
    printf("\nRMSE = %lf\n", RMSE);
}

static void compute_loss_function()
{
    int input = 0;
    printf("Which loss function do you want to use?\n");
    printf(" 1 - RMSE\n");
    scanf("%d", &input);
    switch(input) {
        case 1:
            concatenate(cmd_gnu_0, lossname_RMSE);
            compute_RMSE();
            printf("\n");
            break;
        default:
            printf("Wrong input try again\n");
            compute_loss_function(size);
    }
}

static void compute_r_square(size_t size)
{
    double num_err = 0.0;
    double den_err = 0.0;
    size_t i;

    for (i = 0; i < size; i++) {
        num_err += (vector_y[i] - vector_predict[i]) * (vector_y[i] - vector_predict[i]);
        den_err += (vector_y[i] - mean_y) * (vector_y[i] - mean_y);
    }
    r_square = 1 - (num_err/den_err);
    printf("R_square = %lf\n", r_square);
}

static void compute_predict_for_x()
{
    double x = 0.0;
    printf("Please enter x value\n");
    scanf("%lf", &x);
    prediction = teta_0_intercept + (teta_1_grad * x);
    printf("y^ if x = %lf -> %lf\n",x, prediction);
}

static void compute_max_x_y()
{
    size_t i;
    double tmp1= 0.0;
    double tmp2= 0.0;

    for (i = 0; i < size; i++) {
        if (vector_x[i] > tmp1) {
            tmp1 = vector_x[i];
            max_x = vector_x[i];

        }
        if (vector_y[i] > tmp2) {
            tmp2 = vector_y[i];
            max_y = vector_y[i];
        }
    }
    printf("vector_x max value %lf\n", max_x);
    printf("vector_y max value %lf\n", max_y);
}

static void display_model_line()
{
    sprintf(intercept, "%0.7lf", teta_0_intercept);
    sprintf(grad, "%0.7lf", teta_1_grad);
    sprintf(xrange, "%0.7lf", max_x + 1);
    sprintf(yrange, "%0.7lf", max_y + 1);

    concatenate(cmd_gnu_1, intercept);
    concatenate(cmd_gnu_2, grad);
    concatenate(cmd_gnu_3, xrange);
    concatenate(cmd_gnu_3, "]");
    concatenate(cmd_gnu_4, yrange);
    concatenate(cmd_gnu_4, "]");

    printf("grad = %s\n", grad);
    printf("intercept = %s\n", intercept);
    printf("xrange = %s\n", xrange);
    printf("yrange = %s\n", yrange);

    printf("cmd_gnu_0: %s\n", cmd_gnu_0);
    printf("cmd_gnu_1: %s\n", cmd_gnu_1);
    printf("cmd_gnu_2: %s\n", cmd_gnu_2);
    printf("cmd_gnu_3: %s\n", cmd_gnu_3);
    printf("cmd_gnu_4: %s\n", cmd_gnu_4);
    printf("cmd_gnu_5: %s\n", cmd_gnu_5);
    printf("cmd_gnu_6: %s\n", cmd_gnu_6);

    /* print plot */
    FILE *gnuplot_pipe = (FILE*)popen("gnuplot -persistent", "w");
    FILE *temp = (FILE*)fopen("data.temp", "w");

    /* create data.temp */
    size_t i;
    for (i = 0; i < size; i++)
    {
        fprintf(temp, "%f %f \n", vector_x[i], vector_y[i]);
    }

    /* display gnuplot */
    for (i = 0; i < 7; i++)
    {
        fprintf(gnuplot_pipe, "%s \n", commands_gnuplot[i]);
    }
}

int main(void)
{
    printf("===========================================\n");
    printf("INPUT DATA\n");
    printf("===========================================\n");
    user_input();
    display_vector();
    printf("\n");

    printf("===========================================\n");
    printf("COMPUTE MEAN X:Y, TETA_1 TETA_0\n");
    printf("===========================================\n");
    compute_mean_x_y();
    compute_max_x_y();
    compute_teta_1_grad();
    compute_teta_0_intercept();
    printf("\n");

    printf("===========================================\n");
    printf("COMPUTE LOSS FUNCTION\n");
    printf("===========================================\n");
    compute_loss_function();

    printf("===========================================\n");
    printf("COMPUTE R_square\n");
    printf("===========================================\n");
    compute_r_square(size);
    printf("\n");

    printf("===========================================\n");
    printf("COMPUTE y^ according to x\n");
    printf("===========================================\n");
    compute_predict_for_x();
    printf("\n");

    printf("===========================================\n");
    printf("DISPLAY LINEAR REGRESSION\n");
    printf("===========================================\n");
    display_model_line();
    printf("\n");

    return 0;
}

as an assignment I had to code in C a simple linear regression using RMSE loss function. The program is dynamic and you can enter your own values and choose your own loss function which is for now limited to Root Mean Square Error. But first here are the algorithms I used:

enter image description here
enter image description here
enter image description here
enter image description here

now the code... you need gnuplot to display the chart, sudo apt install gnuplot

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/types.h>

#define BUFFSIZE 64
#define MAXSIZE 100

static double vector_x[MAXSIZE] = {0};
static double vector_y[MAXSIZE] = {0};
static double vector_predict[MAXSIZE] = {0};

static double max_x;
static double max_y;
static double mean_x;
static double mean_y;
static double teta_0_intercept;
static double teta_1_grad;
static double RMSE;
static double r_square;
static double prediction;

static char intercept[BUFFSIZE];
static char grad[BUFFSIZE];
static char xrange[BUFFSIZE];
static char yrange[BUFFSIZE];
static char lossname_RMSE[BUFFSIZE] = "Simple Linear Regression using RMSE'";

static char cmd_gnu_0[BUFFSIZE] = "set title '";
static char cmd_gnu_1[BUFFSIZE] = "intercept = ";
static char cmd_gnu_2[BUFFSIZE] = "grad = ";
static char cmd_gnu_3[BUFFSIZE] = "set xrange [0:";
static char cmd_gnu_4[BUFFSIZE] = "set yrange [0:";
static char cmd_gnu_5[BUFFSIZE] = "f(x) = (grad * x) + intercept";
static char cmd_gnu_6[BUFFSIZE] = "plot f(x), 'data.temp' with points pointtype 7";

static char const *commands_gnuplot[] = {
    cmd_gnu_0,
    cmd_gnu_1,
    cmd_gnu_2,
    cmd_gnu_3,
    cmd_gnu_4,
    cmd_gnu_5,
    cmd_gnu_6,
};

static size_t size;

static void user_input()
{
    printf("Enter x,y vector size, MAX = 100\n");
    scanf("%lu", &size);
    if (size > MAXSIZE) {
        printf("Wrong input size is too big\n");
        user_input();
    }
    printf("vector's size is %lu\n", size);

    size_t i;
    for (i = 0; i < size; i++) {
        printf("Enter vector_x[%ld] values\n", i);
        scanf("%lf", &vector_x[i]);
    }

    for (i = 0; i < size; i++) {
        printf("Enter vector_y[%ld] values\n", i);
        scanf("%lf", &vector_y[i]);
    }
}

static void display_vector()
{
    size_t i;
    for (i = 0; i < size; i++){
        printf("vector_x[%lu] = %lf\t", i, vector_x[i]);
        printf("vector_y[%lu] = %lf\n", i, vector_y[i]);
    }
}

static void concatenate(char p[], char q[]) {
   int c;
   int d;
   c = 0;

   while (p[c] != '\0') {
      c++;
   }
   d = 0;

   while (q[d] != '\0') {
      p[c] = q[d];
      d++;
      c++;
   }
   p[c] = '\0';
}

static void compute_mean_x_y()
{
    size_t i;
    double tmp_x = 0.0;
    double tmp_y = 0.0;
    for (i = 0; i < size; i++) {
        tmp_x += vector_x[i];
        tmp_y += vector_y[i];
    }

    mean_x = tmp_x / size;
    mean_y = tmp_y / size;

    printf("mean_x = %lf\n", mean_x);
    printf("mean_y = %lf\n", mean_y);
}

static void compute_teta_1_grad()
{
    double numerator = 0.0;
    double denominator = 0.0;
    double tmp1 = 0.0;
    double tmp2 = 0.0;
    size_t i;

    for (i = 0; i < size; i++) {
        numerator += (vector_x[i] - mean_x) * (vector_y[i] - mean_y);
    }

    for (i = 0; i < size; i++) {
        tmp1 = vector_x[i] - mean_x;
        tmp2 = tmp1 * tmp1;
        denominator += tmp2;
    }

    teta_1_grad = numerator / denominator;
    printf("teta_1_grad = %lf\n", teta_1_grad);
}

static void compute_teta_0_intercept()
{
    teta_0_intercept = mean_y - (teta_1_grad * mean_x);
    printf("teta_0_intercept = %lf\n", teta_0_intercept);
}

static void compute_prediction()
{
    size_t i;
    for (i = 0; i < size; i++) {
        vector_predict[i] = teta_0_intercept + (teta_1_grad * vector_x[i]);
        printf("y^[%ld] = %lf\n", i, vector_predict[i]);
    }
    printf("\n");
}

static void compute_RMSE()
{
    compute_prediction();
    double error = 0;
    size_t i;
    for (i = 0; i < size; i++) {
        error = (vector_predict[i] - vector_y[i]) * (vector_predict[i] - vector_y[i]);
        printf("error y^[%ld] =  %lf\n", i, error);
        RMSE += error;
    }
    /* mean */
    RMSE = RMSE / size;
    /* square root mean */
    RMSE = sqrt(RMSE);
    printf("\nRMSE = %lf\n", RMSE);
}

static void compute_loss_function()
{
    int input = 0;
    printf("Which loss function do you want to use?\n");
    printf(" 1 - RMSE\n");
    scanf("%d", &input);
    switch(input) {
        case 1:
            concatenate(cmd_gnu_0, lossname_RMSE);
            compute_RMSE();
            printf("\n");
            break;
        default:
            printf("Wrong input try again\n");
            compute_loss_function(size);
    }
}

static void compute_r_square(size_t size)
{
    double num_err = 0.0;
    double den_err = 0.0;
    size_t i;

    for (i = 0; i < size; i++) {
        num_err += (vector_y[i] - vector_predict[i]) * (vector_y[i] - vector_predict[i]);
        den_err += (vector_y[i] - mean_y) * (vector_y[i] - mean_y);
    }
    r_square = 1 - (num_err/den_err);
    printf("R_square = %lf\n", r_square);
}

static void compute_predict_for_x()
{
    double x = 0.0;
    printf("Please enter x value\n");
    scanf("%lf", &x);
    prediction = teta_0_intercept + (teta_1_grad * x);
    printf("y^ if x = %lf -> %lf\n",x, prediction);
}

static void compute_max_x_y()
{
    size_t i;
    double tmp1= 0.0;
    double tmp2= 0.0;

    for (i = 0; i < size; i++) {
        if (vector_x[i] > tmp1) {
            tmp1 = vector_x[i];
            max_x = vector_x[i];

        }
        if (vector_y[i] > tmp2) {
            tmp2 = vector_y[i];
            max_y = vector_y[i];
        }
    }
    printf("vector_x max value %lf\n", max_x);
    printf("vector_y max value %lf\n", max_y);
}

static void display_model_line()
{
    sprintf(intercept, "%0.7lf", teta_0_intercept);
    sprintf(grad, "%0.7lf", teta_1_grad);
    sprintf(xrange, "%0.7lf", max_x + 1);
    sprintf(yrange, "%0.7lf", max_y + 1);

    concatenate(cmd_gnu_1, intercept);
    concatenate(cmd_gnu_2, grad);
    concatenate(cmd_gnu_3, xrange);
    concatenate(cmd_gnu_3, "]");
    concatenate(cmd_gnu_4, yrange);
    concatenate(cmd_gnu_4, "]");

    printf("grad = %s\n", grad);
    printf("intercept = %s\n", intercept);
    printf("xrange = %s\n", xrange);
    printf("yrange = %s\n", yrange);

    printf("cmd_gnu_0: %s\n", cmd_gnu_0);
    printf("cmd_gnu_1: %s\n", cmd_gnu_1);
    printf("cmd_gnu_2: %s\n", cmd_gnu_2);
    printf("cmd_gnu_3: %s\n", cmd_gnu_3);
    printf("cmd_gnu_4: %s\n", cmd_gnu_4);
    printf("cmd_gnu_5: %s\n", cmd_gnu_5);
    printf("cmd_gnu_6: %s\n", cmd_gnu_6);

    /* print plot */
    FILE *gnuplot_pipe = (FILE*)popen("gnuplot -persistent", "w");
    FILE *temp = (FILE*)fopen("data.temp", "w");

    /* create data.temp */
    size_t i;
    for (i = 0; i < size; i++)
    {
        fprintf(temp, "%f %f \n", vector_x[i], vector_y[i]);
    }

    /* display gnuplot */
    for (i = 0; i < 7; i++)
    {
        fprintf(gnuplot_pipe, "%s \n", commands_gnuplot[i]);
    }
}

int main(void)
{
    printf("===========================================\n");
    printf("INPUT DATA\n");
    printf("===========================================\n");
    user_input();
    display_vector();
    printf("\n");

    printf("===========================================\n");
    printf("COMPUTE MEAN X:Y, TETA_1 TETA_0\n");
    printf("===========================================\n");
    compute_mean_x_y();
    compute_max_x_y();
    compute_teta_1_grad();
    compute_teta_0_intercept();
    printf("\n");

    printf("===========================================\n");
    printf("COMPUTE LOSS FUNCTION\n");
    printf("===========================================\n");
    compute_loss_function();

    printf("===========================================\n");
    printf("COMPUTE R_square\n");
    printf("===========================================\n");
    compute_r_square(size);
    printf("\n");

    printf("===========================================\n");
    printf("COMPUTE y^ according to x\n");
    printf("===========================================\n");
    compute_predict_for_x();
    printf("\n");

    printf("===========================================\n");
    printf("DISPLAY LINEAR REGRESSION\n");
    printf("===========================================\n");
    display_model_line();
    printf("\n");

    return 0;
}
回复 原文
若有似无的小暗淡 2024-10-25 20:38:42

请参阅本文的第 1 节。本节将二维线性回归表示为矩阵乘法练习。只要您的数据表现良好，这种技术就应该允许您开发快速的最小二乘拟合。

根据数据的大小，可能值得将矩阵乘法从代数角度简化为简单的方程组，从而避免编写 matmult() 函数。 （预先警告，对于超过 4 或 5 个数据点，这是完全不切实际的！）

Look at Section 1 of this paper. This section expresses a 2D linear regression as a matrix multiplication exercise. As long as your data is well-behaved, this technique should permit you to develop a quick least squares fit.

Depending on the size of your data, it might be worthwhile to algebraically reduce the matrix multiplication to simple set of equations, thereby avoiding the need to write a matmult() function. (Be forewarned, this is completely impractical for more than 4 or 5 data points!)

回复 原文
隱形的亼 2024-10-25 20:38:42

据我所知，求解最小二乘法最快、最有效的方法是从参数向量中减去（梯度）/（二阶梯度）。 （二阶梯度 = 即 Hessian 矩阵的对角线。）

直觉如下：

假设您想要优化单个参数的最小二乘法。这相当于找到抛物线的顶点。然后，对于任意随机初始参数 x0，损失函数的顶点位于 x0 - f(1) / f (2)。这是因为将 - f(1) / f(2) 添加到 x 总是会将导数 f(1) 归零。

旁注：在 Tensorflow 中实现这一点，解决方案出现在 w0 - f(1) / f(2) / (权重数），但我不确定这是由于 Tensorflow 还是其他原因造成的。

The fastest, most efficient way to solve least squares, as far as I am aware, is to subtract (the gradient)/(the 2nd order gradient) from your parameter vector. (2nd order gradient = i.e. the diagonal of the Hessian.)

Here is the intuition:

Let's say you want to optimize least squares over a single parameter. This is equivalent to finding the vertex of a parabola. Then, for any random initial parameter, x0, the vertex of the loss function is located at x0 - f(1) / f(2). That's because adding - f(1) / f(2) to x will always zero out the derivative, f(1).

Side note: Implementing this in Tensorflow, the solution appeared at w0 - f(1) / f(2) / (number of weights), but I'm not sure if that's due to Tensorflow or if it's due to something else..

回复 原文
~没有更多了~

关于作者

梦行七里

暂无简介

0 文章
0 评论
24 人气
发私信

相关话题

更多

热门标签

操作系统 程序设计 IT运维 Linux系统管理 JavaScript 服务器应用 solaris C/C++ PHP Shell BSD Vue.js aix Oracle Python HTML 系统管理 HTML5 CSS 前端
更多

推荐作者

苦中寻乐

文章 0 评论 0

lueluelue

文章 0 评论 0

嗼ふ静

文章 0 评论 0

王权女流氓

文章 0 评论 0

与花如笺

文章 0 评论 0

残酷

文章 0 评论 0

更多

友情链接

文江博客
    我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站，即表示您同意使用 Cookies 和您的相关数据。
    原文