OriginC: Regression Analysis
Version Info
Minimum Origin Version Required: Origin 8 SR0
Linear Regression
#include <oc_nag.h>
//#include <nag\OC_nag_ex.h> //not need, only oc_nag.h required.
void nag_linear_regression_ex1()
{
// data input directly by vector
vector vx = {1,0,4,7.5,2.5,0,10,5};
vector vy = {20,15.5,28.3,45,24.5,10,99,31.2};
vector vw = {1,1,1,1,1,1,1,1};
Nag_SumSquare mean = Nag_AboutMean;
int n;
n = vx.GetSize();
double a, b, err_a, err_b, rsq, rss, df;
g02cac(mean, n, vx, vy, vw, &a, &b, &err_a, &err_b, &rsq, &rss,
&df, NAGERR_DEFAULT);
// output intepretation: a+b*x b is the slope, a is the intercept
if (mean == Nag_AboutMean)
{
printf("\nRegression constant a = %6.4f\n\n", a);
printf("Standard error of the regression constant a = %6.4f\n\n",err_a);
}
printf("Regression coefficient b = %6.4f\n\n", b);
printf("Standard error of the regression coefficient b = %6.4f\n\n",err_b);
printf("The regression coefficient of determination = %6.4f\n\n", rsq);
printf("The sum of squares of the residuals about the " "regression = %6.4f\n\n", rss);
printf("Number of degrees of freedom about the " "regression = %6.4f\n\n",df);
//Expected Result:
//Regression constant a = 7.5982
//Standard error of the regression constant a = 6.6858
//
//Regression coefficient b = 7.0905
//
//Standard error of the regression coefficient b = 1.3224
//
//The regression coefficient of determination = 0.8273
//
//The sum of squares of the residuals about the regression = 965.2454
//
//Number of degrees of freedom about the regression = 6.0000
}
#include <oc_nag.h>
//#include <nag\OC_nag_ex.h> //not need, only oc_nag.h required.
void nag_linear_regression_ex2()
{
// obtain regression data from an Origin Worksheet
Worksheet wks=Project.ActiveLayer();
if(!wks)
printf("no active worksheet found\n");
Dataset dsSrcX(wks,0); // assume that the x, y data and weight data are in the first three columns
Dataset dsSrcY(wks,1);
Dataset dsSrcW(wks,2);
int nSize;
vector vx(nSize);
vector vy(nSize);
vector vw(nSize);
vx=dsSrcX;
vy=dsSrcY;
vw=dsSrcW;
int n;
n = vx.GetSize();
Nag_SumSquare mean = Nag_AboutMean;
double a, b, err_a, err_b, rsq, rss, df;
g02cac(mean, n, vx, vy, vw, &a, &b, &err_a, &err_b, &rsq, &rss,
&df, NAGERR_DEFAULT);
// output intepretation: a+b*x b is the slope, a is the intercept
if (mean == Nag_AboutMean)
{
printf("\nRegression constant a = %6.4f\n\n", a);
printf("Standard error of the regression constant a = %6.4f\n\n",err_a);
}
printf("Regression coefficient b = %6.4f\n\n", b);
printf("Standard error of the regression coefficient b = %6.4f\n\n",err_b);
printf("The regression coefficient of determination = %6.4f\n\n", rsq);
printf("The sum of squares of the residuals about the " "regression = %6.4f\n\n", rss);
printf("Number of degrees of freedom about the " "regression = %6.4f\n\n",df);
}
Multiple Linear Regression
#include <oc_nag.h>
//#include <nag\OC_nag_ex.h> //not need, only oc_nag.h required.
void nag_regsn_mult_linear_ex()
{
double rss, tol;
int i, ip, rank, j, m, n;
double df;
Boolean svd;
char weight, meanc;
Nag_IncludeMean mean;
double b[20];
double cov[210], h[20], p[440];
double res[20], se[20];
double com_ar[495];
double wt[20];
matrix q;
q.SetSize(20, 21);
double *wtptr;
n = 12;
m = 4;
mean = Nag_MeanInclude;
wtptr =NULL;
// input data is only the first four columns, the design matrix
matrix x = {{1.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 0.0, 0.0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 1.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 0.0, 1.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 0.0, 0.0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 1.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 0.0, 0.0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 0.0, 1.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 0.0, 1.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 1.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
};
double y[20] = {33.63, 39.62, 38.18, 41.46, 38.02, 35.83,
35.99, 36.58, 42.92, 37.80, 40.43, 37.89}; // the y data
int sx[20] = {1, 1, 1, 1};
printf("The input data are as follows\n");
printf("n = 12, m = 4, tol = 0.000001e0\n");
printf("x\n");
for(i = 0; i < n; i++)
{
for(j = 0; j < m; j++)
printf("%2.1f, ",x[i][j]);
printf("\n");
}
printf("\ny\n");
for(i = 0; i < n; i ++)
{
printf("%5.3f, ", y[i]);
if((i + 1) % 7 == 0)
printf("\n");
}
ip = 0;
if (mean==Nag_MeanInclude)
ip += 1;
for (i=0; i<m; i++)
if (sx[i]>0) ip += 1;
tol = 0.00001e0;
nag_regsn_mult_linear(mean, n, x, 20, m, sx, ip, y,
wtptr, &rss, &df, b, se, cov, res, h, q, 21, &svd, &rank, p, tol, com_ar, NAGERR_DEFAULT);
if (svd)
{ printf("\nModel not of full rank, rank = %4ld\n", rank);
ASSERT( rank == 4 );
}
printf("Residual sum of squares = %12.4e\n", rss);
printf("Degrees of freedom = %3.1f\n\n", df);
printf("Variable Parameter Estimate Standard error\n\n"); // the standard error
for (j=0; j<ip; j++)
printf("%6ld%20.4e%20.4e\n", j+1, b[j], se[j]);
printf("\n");
printf("Obs Residuals h\n\n"); // the residuals
for (i=0; i<n; i++)
printf("%6ld%20.4e%20.4e\n", i+1, res[i], h[i]);
}
Nonlinear Regression
#include <OC_nag.h>
//Callback to calculate objective subfunction values, and (optionally) its Jacobian
void NAG_CALL objfun_callback(int m, int n, double p[], double f[],
double fjac[], int tdfjac,Nag_Comm *comm)
{
#define FJAC(I,J) fjac[(I)*tdfjac +(J)]
static vector x = {
8.0, 8.0, 10.0, 10.0, 10.0, 10.0, 12.0, 12.0, 12.0, 12.0, 14.0, 14.0,
14.0, 16.0, 16.0, 16.0, 18.0, 18.0, 20.0, 20.0, 20.0, 22.0, 22.0, 22.0,
24.0, 24.0, 24.0, 26.0, 26.0, 26.0, 28.0, 28.0, 30.0, 30.0, 30.0, 32.0,
32.0, 34.0, 36.0, 36.0, 38.0, 38.0, 40.0, 42.0};
double temp;
for(int i = 0; i < m; ++i)
{
//evaluate objective subfunction f(i+1) only if required
if(comm->needf == i+1 || comm->needf == 0)
{
temp = exp(-p[1] * (x[i] - 8.0));
f[i] = p[0] + (.49 - p[0]) * temp;
}
//evaluate the Jacobian if required
if(comm->flag == 2)
{
FJAC(i, 0) = 1.0 - temp;
FJAC(i, 1) = -(0.49 - p[0]) * (x[i] - 8.0) * temp;
}
}
}
//Callback to calculate nonlinear constraints and (optinally) its Jacobian
static void NAG_CALL confun_callback(int n, int ncnlin, int needc[], double p[],
double conf[], double cjac[], Nag_Comm *comm)
{
#define CJAC(I, J) cjac[(I)*n + (J)]
//First call to confun. Set all Jacobian elements to zero.
//Note that this will only work when options.obj_deriv = TRUE(default)
if(comm->first == true)
CJAC(0, 0) = CJAC(0, 1) = 0.0;
if(needc[0] > 0)
{
conf[0] = -.09 - p[0]*p[1] + 0.49*p[1];
if(comm->flag == 2)
{
CJAC(0, 0) = -p[1];
CJAC(0, 1) = -p[0] + .49;
}
}
}
void nag_opt_nlin_lsq_ex()
{
vector y = {0.49, 0.49, 0.48, 0.47, 0.48, 0.47, 0.46, 0.46, 0.45, 0.43, 0.45,
0.43, 0.43, 0.44, 0.43, 0.43, 0.46, 0.45, 0.42, 0.42, 0.43, 0.41,
0.41, 0.40, 0.42, 0.40, 0.40, 0.41, 0.40, 0.41, 0.41, 0.40, 0.40,
0.40, 0.38, 0.41, 0.40, 0.40, 0.41, 0.38, 0.40, 0.40, 0.39, 0.39};
int m = y.GetSize(), n = 2;
int tda = n, tdfjac = n;
int nclin = 1, ncnlin = 1;
vector p = {0.4, 0.0};//initial parameters
vector a = {1.0, 1.0};//linear constraints
vector bl = {0.4, -4.0, 1.0, 0.0}; //lower bounds
vector bu = {1.0e+25, 1.0e+25, 1.0e+25, 1.0e+25}; //upper bounds
vector f(m), fjac(m*n);
double objf;
Nag_E04_Opt options;
nag_opt_init(&options);
strcpy(options.outfile, "c:\\nag_opt_nlin_lsq_ex.txt");
NagError fail;
e04unc(m, n, nclin, ncnlin, a, tda, bl, bu, y, objfun_callback, confun_callback,
p, &objf, f, fjac, tdfjac, &options, NAGCOMM_NULL, &fail);
if(fail.code == NE_NOERROR)
{
printf("fitting paramters: %lf %lf \n", p[0], p[1]);
printf("value of objective function: %lf \n", objf);
}
else
{
printf("fitting failed, error code = %d", fail.code);
}
//Expected Result:
//fitting paramters: 0.419953 1.284845
//value of objective function: 0.014230
}
|