campo-sirio/gfm/linest.c
alex ba237a9d91 Patch level : no patch
Files correlati     :
Ricompilazione Demo : [ ]
Commento            :
Aggiunti i sorgenti per Greenleaf Math Library (gfm.dll)


git-svn-id: svn://10.65.10.50/trunk@10079 c028cbd2-c16b-5b4b-a496-9718f37d4682
2002-02-26 12:19:02 +00:00

146 lines
3.9 KiB
C
Executable File

/* DEC *LinearEstimate(pCorr, pLinA, pLinB, pSrcX, pSrcY, wSize);
*
* ARGUMENT
* DEC *pCorr, *pLinA, *pLinB;
* DEC **pSrcX, **pSrcY;
* int wSize;
*
* DESCRIPTION
* Determines the least squares line y = A + B * x and the
* linear correlation coefficient of a given set of (X, Y) values.
* The wSize (X, Y) values, where wSize > 1, are stored in arrays
* pSrcX and pSrcY. The correlation coefficient is put in pCorr.
* The value A is put in the global pGMStatA, along with pLinA,
* while B is put in the global pGMStatB and pLinB. Storing A and
* B in globals allow their later usage in PredictX and PredictY.
* The correlation coefficient is returned.
*
* SIDE EFFECTS
* None.
*
* RETURNS
* Returns pointer to the correlation coefficient if successful,
* otherwise returns GM_NULL.
*
* POSSIBLE ERROR CODES
*
* GM_NULLPOINTER
* GM_ARGVAL
* GM_OVERFLOW
* GM_UNDERFLOW
*
* AUTHOR
* Jared Levy
* Copyright (C) 1987-1990 Greenleaf Software Inc. All rights reserved.
*
* MODIFICATIONS
*
*/
#include <stdio.h>
#include "gm.h"
#include "gmsystem.h"
DEC *LinearEstimate(pCorr, pLinA, pLinB, pSrcX, pSrcY, wSize)
DEC *pCorr, *pLinA, *pLinB;
DEC **pSrcX, **pSrcY;
int wSize;
{
int i;
DEC dSumX, *pSumX=&dSumX, dSumY, *pSumY=&dSumY;
DEC dSumXX, *pSumXX=&dSumXX, dSumYY, *pSumYY=&dSumYY;
DEC dSumXY, *pSumXY=&dSumXY, dSize, *pSize=&dSize;
DEC dProd, *pProd=&dProd, dProd2, *pProd2=&dProd2;
DEC dVarX, *pVarX=&dVarX, dVarY, *pVarY=&dVarY;
DEC dNum, *pNum=&dNum, dDen, *pDen=&dDen;
_MacStart(GM_LINEST);
if (!pSrcX||!pSrcY||!pCorr||!pLinA||!pLinB) {
_MacErr(GM_NULLPOINTER);
_MacRet(GM_NULL);
}
if (wSize<=1) {
_MacErr(GM_ARGVAL);
_MacRet(GM_NULL);
}
for (i=0; i<wSize; i++) {
_MacInVarD(pSrcX[i]);
_MacInVarD(pSrcY[i]);
}
/* |pSrcX[i]| && |pSrcY[i]| must be < decMaxLinEst (1e6) to avoid overflow */
/* values less than 1e-6 may also lead to incorrect results */
for (i=0; i<wSize; i++)
if ((CompareDecimal(pSrcX[i],&decMaxLinEst)==1) ||
(CompareDecimal(pSrcY[i],&decMaxLinEst)==1) ||
(CompareDecimal(pSrcX[i],&decMinLinEst)==-1) ||
(CompareDecimal(pSrcY[i],&decMinLinEst)==-1)) {
_MacErr(GM_ARGVAL);
_MacRet(GM_NULL);
}
/* calculate the sum of X, Y, X^2, Y^2, and XY */
_MacDZero(pSumX);
_MacDZero(pSumY);
_MacDZero(pSumXX);
_MacDZero(pSumYY);
_MacDZero(pSumXY);
for (i=0;i<wSize;i++) {
(void) _AddDec80Bit(pSumX, pSumX, pSrcX[i]);
(void) _AddDec80Bit(pSumY, pSumY, pSrcY[i]);
(void) _MulDec80Bit(pProd, pSrcX[i], pSrcX[i]);
(void) _AddDec80Bit(pSumXX, pSumXX, pProd);
(void) _MulDec80Bit(pProd, pSrcY[i], pSrcY[i]);
(void) _AddDec80Bit(pSumYY, pSumYY, pProd);
(void) _MulDec80Bit(pProd, pSrcX[i], pSrcY[i]);
(void) _AddDec80Bit(pSumXY, pSumXY, pProd);
}
/* Compute n*sum(X^2) - sum(X)*sum(X) */
(void) ConvLongToDecimal(pSize, (long) wSize);
(void) _MulDec80Bit(pProd, pSumXX, pSize);
(void) _MulDec80Bit(pProd2, pSumX, pSumX);
(void) _SubDec80Bit(pVarX, pProd, pProd2);
if (_MacIsDecZ(pVarX)) {
_MacErr(GM_UNDERFLOW);
_MacRet(GM_NULL);
}
/* Compute A */
(void) _MulDec80Bit(pProd, pSumXX, pSumY);
(void) _MulDec80Bit(pProd2, pSumX, pSumXY);
(void) _SubDec80Bit(pNum, pProd, pProd2);
if (!DivideDecimal(pLinA,pNum,pVarX))
_MacRet(GM_NULL);
_MacDCopy(pGMStatA, pLinA);
/* Compute B */
(void) _MulDec80Bit(pProd, pSumXY, pSize);
(void) _MulDec80Bit(pProd2, pSumX, pSumY);
(void) _SubDec80Bit(pNum, pProd, pProd2);
if (!DivideDecimal(pLinB,pNum,pVarX))
_MacRet(GM_NULL);
_MacDCopy(pGMStatB, pLinB);
/* Compute correlation coefficient */
(void) _MulDec80Bit(pProd, pSumYY, pSize);
(void) _MulDec80Bit(pProd2, pSumY, pSumY);
(void) _SubDec80Bit(pVarY, pProd, pProd2);
(void) _MulDec80Bit(pProd, pVarX, pVarY);
_SqrtDec80Bit(pDen, pProd);
if (_MacIsDecZ(pDen)) {
_MacErr(GM_UNDERFLOW);
_MacRet(GM_NULL);
}
if (!DivideDecimal(pCorr, pNum, pDen))
_MacRet(GM_NULL);
_MacRet(pCorr);
}