[Libreoffice-commits] .: 5 commits - sc/source
Kohei Yoshida
kohei at kemper.freedesktop.org
Fri Nov 5 10:16:39 PDT 2010
sc/source/core/inc/interpre.hxx | 15
sc/source/core/tool/interpr4.cxx | 2
sc/source/core/tool/interpr5.cxx | 1720 ++++++++++++++++++++++++---------------
3 files changed, 1096 insertions(+), 641 deletions(-)
New commits:
commit 44931e585d56e0d31f811a03f3d41d204b1acf49
Merge: 43a8963... 96f07b9...
Author: Kohei Yoshida <kyoshida at novell.com>
Date: Fri Nov 5 13:15:31 2010 -0400
Merge branch 'feature/calc-function-linest-logest'
commit 96f07b9ca435d5641214c8d65cc407495d4fb3be
Author: Kohei Yoshida <kyoshida at novell.com>
Date: Thu Nov 4 23:12:15 2010 -0400
Removed unnecessary return statement.
diff --git a/sc/source/core/tool/interpr5.cxx b/sc/source/core/tool/interpr5.cxx
index 757dfce..759614a 100644
--- a/sc/source/core/tool/interpr5.cxx
+++ b/sc/source/core/tool/interpr5.cxx
@@ -2874,7 +2874,6 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
PushMatrix(pResMat);
}
}
- return;
}
void ScInterpreter::ScTrend()
commit d77b1bd4c25818e587420d79f755f341ca9b8ee1
Author: Kohei Yoshida <kyoshida at novell.com>
Date: Thu Nov 4 18:51:24 2010 -0400
Fixed wrong indentations.
For some reason the indentations were all out of wrack.
diff --git a/sc/source/core/tool/interpr5.cxx b/sc/source/core/tool/interpr5.cxx
index 93763ec..757dfce 100644
--- a/sc/source/core/tool/interpr5.cxx
+++ b/sc/source/core/tool/interpr5.cxx
@@ -1856,29 +1856,29 @@ namespace {
// Multiply n x m Mat A with m x l Mat B to n x l Mat R
void lcl_MFastMult(ScMatrixRef pA, ScMatrixRef pB, ScMatrixRef pR,
- SCSIZE n, SCSIZE m, SCSIZE l)
- {
+ SCSIZE n, SCSIZE m, SCSIZE l)
+{
double sum;
for (SCSIZE row = 0; row < n; row++)
- {
+ {
for (SCSIZE col = 0; col < l; col++)
{ // result element(col, row) =sum[ (row of A) * (column of B)]
- sum = 0.0;
+ sum = 0.0;
for (SCSIZE k = 0; k < m; k++)
sum += pA->GetDouble(k,row) * pB->GetDouble(col,k);
pR->PutDouble(sum, col, row);
- }
- }
}
+ }
+}
// <A;B> over all elements; uses the matrices as vectors of length M
double lcl_GetSumProduct(ScMatrixRef pMatA, ScMatrixRef pMatB, SCSIZE nM)
- {
+{
double fSum = 0.0;
for (SCSIZE i=0; i<nM; i++)
fSum += pMatA->GetDouble(i) * pMatB->GetDouble(i);
return fSum;
- }
+}
// Special version for use within QR decomposition.
// Euclidean norm of column index C starting in row index R;
@@ -1889,7 +1889,7 @@ double lcl_GetColumnEuclideanNorm(ScMatrixRef pMatA, SCSIZE nC, SCSIZE nR, SCSIZ
for (SCSIZE row=nR; row<nN; row++)
fNorm += (pMatA->GetDouble(nC,row)) * (pMatA->GetDouble(nC,row));
return sqrt(fNorm);
- }
+}
// Euclidean norm of row index R starting in column index C;
// matrix A has count N columns.
@@ -1899,7 +1899,7 @@ double lcl_TGetColumnEuclideanNorm(ScMatrixRef pMatA, SCSIZE nR, SCSIZE nC, SCSI
for (SCSIZE col=nC; col<nN; col++)
fNorm += (pMatA->GetDouble(col,nR)) * (pMatA->GetDouble(col,nR));
return sqrt(fNorm);
- }
+}
// Special version for use within QR decomposition.
// Maximum norm of column index C starting in row index R;
@@ -1911,7 +1911,7 @@ double lcl_GetColumnMaximumNorm(ScMatrixRef pMatA, SCSIZE nC, SCSIZE nR, SCSIZE
if (fNorm < fabs(pMatA->GetDouble(nC,row)))
fNorm = fabs(pMatA->GetDouble(nC,row));
return fNorm;
- }
+}
// Maximum norm of row index R starting in col index C;
// matrix A has count N columns.
@@ -1922,14 +1922,14 @@ double lcl_TGetColumnMaximumNorm(ScMatrixRef pMatA, SCSIZE nR, SCSIZE nC, SCSIZE
if (fNorm < fabs(pMatA->GetDouble(col,nR)))
fNorm = fabs(pMatA->GetDouble(col,nR));
return fNorm;
- }
+}
// Special version for use within QR decomposition.
// <A(Ca);B(Cb)> starting in row index R;
// Ca and Cb are indices of columns, matrices A and B have count N rows.
double lcl_GetColumnSumProduct(ScMatrixRef pMatA, SCSIZE nCa,
- ScMatrixRef pMatB, SCSIZE nCb, SCSIZE nR, SCSIZE nN)
- {
+ ScMatrixRef pMatB, SCSIZE nCb, SCSIZE nR, SCSIZE nN)
+{
double fResult = 0.0;
for (SCSIZE row=nR; row<nN; row++)
fResult += pMatA->GetDouble(nCa,row) * pMatB->GetDouble(nCb,row);
@@ -1939,8 +1939,8 @@ double lcl_GetColumnSumProduct(ScMatrixRef pMatA, SCSIZE nCa,
// <A(Ra);B(Rb)> starting in column index C;
// Ra and Rb are indices of rows, matrices A and B have count N columns.
double lcl_TGetColumnSumProduct(ScMatrixRef pMatA, SCSIZE nRa,
- ScMatrixRef pMatB, SCSIZE nRb, SCSIZE nC, SCSIZE nN)
- {
+ ScMatrixRef pMatB, SCSIZE nRb, SCSIZE nC, SCSIZE nN)
+{
double fResult = 0.0;
for (SCSIZE col=nC; col<nN; col++)
fResult += pMatA->GetDouble(col,nRa) * pMatB->GetDouble(col,nRb);
@@ -1948,13 +1948,13 @@ double lcl_TGetColumnSumProduct(ScMatrixRef pMatA, SCSIZE nRa,
}
double lcl_GetSign(double fValue)
- {
+{
if (fValue < 0.0)
return -1.0;
else if (fValue > 0.0)
- return 1.0;
- else
- return 0.0;
+ return 1.0;
+ else
+ return 0.0;
}
/* Calculates a QR decomposition with Householder reflection.
@@ -1971,8 +1971,8 @@ double lcl_GetSign(double fValue)
* errors singularity is often not detected.
*/
bool lcl_CalculateQRdecomposition(ScMatrixRef pMatA,
- ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN)
- {
+ ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN)
+{
double fScale ;
double fEuclid ;
double fFactor ;
@@ -1984,10 +1984,10 @@ bool lcl_CalculateQRdecomposition(ScMatrixRef pMatA,
// calculate vector u of the householder transformation
fScale = lcl_GetColumnMaximumNorm(pMatA, col, col, nN);
if (fScale == 0.0)
- {
+ {
// A is singular
return false;
- }
+ }
for (SCSIZE row = col; row <nN; row++)
pMatA->PutDouble( pMatA->GetDouble(col,row)/fScale, col, row);
@@ -2002,16 +2002,15 @@ bool lcl_CalculateQRdecomposition(ScMatrixRef pMatA,
{
fSum =lcl_GetColumnSumProduct(pMatA, col, pMatA, c, col, nN);
for (SCSIZE row = col; row <nN; row++)
- pMatA->PutDouble( pMatA->GetDouble(c,row)
- - fSum * fFactor * pMatA->GetDouble(col,row), c, row);
- }
- }
- return true;
+ pMatA->PutDouble( pMatA->GetDouble(c,row) - fSum * fFactor * pMatA->GetDouble(col,row), c, row);
}
+ }
+ return true;
+}
// same with transposed matrix A, N is count of columns, K count of rows
bool lcl_TCalculateQRdecomposition(ScMatrixRef pMatA,
- ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN)
+ ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN)
{
double fScale ;
double fEuclid ;
@@ -2024,10 +2023,10 @@ bool lcl_TCalculateQRdecomposition(ScMatrixRef pMatA,
// calculate vector u of the householder transformation
fScale = lcl_TGetColumnMaximumNorm(pMatA, row, row, nN);
if (fScale == 0.0)
- {
+ {
// A is singular
return false;
- }
+ }
for (SCSIZE col = row; col <nN; col++)
pMatA->PutDouble( pMatA->GetDouble(col,row)/fScale, col, row);
@@ -2042,12 +2041,12 @@ bool lcl_TCalculateQRdecomposition(ScMatrixRef pMatA,
{
fSum =lcl_TGetColumnSumProduct(pMatA, row, pMatA, r, row, nN);
for (SCSIZE col = row; col <nN; col++)
- pMatA->PutDouble( pMatA->GetDouble(col,r)
- - fSum * fFactor * pMatA->GetDouble(col,row), col, r);
+ pMatA->PutDouble(
+ pMatA->GetDouble(col,r) - fSum * fFactor * pMatA->GetDouble(col,row), col, r);
}
}
return true;
- }
+}
/* Applies a Householder transformation to a column vector Y with is given as
@@ -2057,15 +2056,15 @@ bool lcl_TCalculateQRdecomposition(ScMatrixRef pMatA,
* lcl_CaluclateQRdecomposition.
*/
void lcl_ApplyHouseholderTransformation(ScMatrixRef pMatA, SCSIZE nC,
- ScMatrixRef pMatY, SCSIZE nN)
- {
+ ScMatrixRef pMatY, SCSIZE nN)
+{
// ScMatrix matrices are zero based, index access (column,row)
double fDenominator = lcl_GetColumnSumProduct(pMatA, nC, pMatA, nC, nC, nN);
double fNumerator = lcl_GetColumnSumProduct(pMatA, nC, pMatY, 0, nC, nN);
double fFactor = 2.0 * (fNumerator/fDenominator);
for (SCSIZE row = nC; row < nN; row++)
pMatY->PutDouble(
- pMatY->GetDouble(row) - fFactor * pMatA->GetDouble(nC,row), row);
+ pMatY->GetDouble(row) - fFactor * pMatA->GetDouble(nC,row), row);
}
// Same with transposed matrices A and Y.
@@ -2115,8 +2114,8 @@ void lcl_SolveWithUpperRightTriangle(ScMatrixRef pMatA,
* zero elements, no check is done.
*/
void lcl_SolveWithLowerLeftTriangle(ScMatrixRef pMatA,
- ::std::vector< double>& pVecR, ScMatrixRef pMatT,
- SCSIZE nK, bool bIsTransposed)
+ ::std::vector< double>& pVecR, ScMatrixRef pMatT,
+ SCSIZE nK, bool bIsTransposed)
{
// ScMatrix matrices are zero based, index access (column,row)
double fSum;
@@ -2127,12 +2126,12 @@ void lcl_SolveWithLowerLeftTriangle(ScMatrixRef pMatA,
{
if (bIsTransposed)
fSum -= pMatA->GetDouble(col,row) * pMatT->GetDouble(col);
- else
+ else
fSum -= pMatA->GetDouble(row,col) * pMatT->GetDouble(col);
- }
- pMatT->PutDouble( fSum / pVecR[row] , row);
- }
}
+ pMatT->PutDouble( fSum / pVecR[row] , row);
+ }
+}
/* Calculates Z = R * B
* R is given in matrix A and vector VecR as obtained from the QR
@@ -2141,18 +2140,18 @@ void lcl_SolveWithLowerLeftTriangle(ScMatrixRef pMatA,
* not used.
*/
void lcl_ApplyUpperRightTriangle(ScMatrixRef pMatA,
- ::std::vector< double>& pVecR, ScMatrixRef pMatB,
- ScMatrixRef pMatZ, SCSIZE nK, bool bIsTransposed)
+ ::std::vector< double>& pVecR, ScMatrixRef pMatB,
+ ScMatrixRef pMatZ, SCSIZE nK, bool bIsTransposed)
{
// ScMatrix matrices are zero based, index access (column,row)
double fSum;
for (SCSIZE row = 0; row < nK; row++)
- {
+ {
fSum = pVecR[row] * pMatB->GetDouble(row);
for (SCSIZE col = row+1; col < nK; col++)
if (bIsTransposed)
fSum += pMatA->GetDouble(row,col) * pMatB->GetDouble(col);
- else
+ else
fSum += pMatA->GetDouble(col,row) * pMatB->GetDouble(col);
pMatZ->PutDouble( fSum, row);
}
@@ -2161,33 +2160,33 @@ void lcl_ApplyUpperRightTriangle(ScMatrixRef pMatA,
double lcl_GetMeanOverAll(ScMatrixRef pMat, SCSIZE nN)
- {
+{
double fSum = 0.0;
for (SCSIZE i=0 ; i<nN; i++)
fSum += pMat->GetDouble(i);
return fSum/static_cast<double>(nN);
- }
+}
// Calculates means of the columns of matrix X. X is a RxC matrix;
// ResMat is a 1xC matrix (=row).
void lcl_CalculateColumnMeans(ScMatrixRef pX, ScMatrixRef pResMat,
- SCSIZE nC, SCSIZE nR)
- {
+ SCSIZE nC, SCSIZE nR)
+{
double fSum = 0.0;
for (SCSIZE i=0; i < nC; i++)
- {
+ {
fSum =0.0;
for (SCSIZE k=0; k < nR; k++)
fSum += pX->GetDouble(i,k); // GetDouble(Column,Row)
pResMat ->PutDouble( fSum/static_cast<double>(nR),i);
- }
+ }
}
// Calculates means of the rows of matrix X. X is a RxC matrix;
// ResMat is a Rx1 matrix (=column).
void lcl_CalculateRowMeans(ScMatrixRef pX, ScMatrixRef pResMat,
- SCSIZE nC, SCSIZE nR)
- {
+ SCSIZE nC, SCSIZE nR)
+{
double fSum = 0.0;
for (SCSIZE k=0; k < nR; k++)
{
@@ -2195,32 +2194,32 @@ void lcl_CalculateRowMeans(ScMatrixRef pX, ScMatrixRef pResMat,
for (SCSIZE i=0; i < nC; i++)
fSum += pX->GetDouble(i,k); // GetDouble(Column,Row)
pResMat ->PutDouble( fSum/static_cast<double>(nC),k);
- }
- }
+ }
+}
void lcl_CalculateColumnsDelta(ScMatrixRef pMat, ScMatrixRef pColumnMeans,
- SCSIZE nC, SCSIZE nR)
+ SCSIZE nC, SCSIZE nR)
{
for (SCSIZE i = 0; i < nC; i++)
for (SCSIZE k = 0; k < nR; k++)
pMat->PutDouble( ::rtl::math::approxSub
- (pMat->GetDouble(i,k) , pColumnMeans->GetDouble(i) ) , i, k);
- }
+ (pMat->GetDouble(i,k) , pColumnMeans->GetDouble(i) ) , i, k);
+}
void lcl_CalculateRowsDelta(ScMatrixRef pMat, ScMatrixRef pRowMeans,
- SCSIZE nC, SCSIZE nR)
+ SCSIZE nC, SCSIZE nR)
{
for (SCSIZE k = 0; k < nR; k++)
for (SCSIZE i = 0; i < nC; i++)
pMat->PutDouble( ::rtl::math::approxSub
- ( pMat->GetDouble(i,k) , pRowMeans->GetDouble(k) ) , i, k);
+ ( pMat->GetDouble(i,k) , pRowMeans->GetDouble(k) ) , i, k);
}
// Case1 = simple regression
// MatX = X - MeanX, MatY = Y - MeanY, y - haty = (y - MeanY) - (haty - MeanY)
// = (y-MeanY)-((slope*x+a)-(slope*MeanX+a)) = (y-MeanY)-slope*(x-MeanX)
double lcl_GetSSresid(ScMatrixRef pMatX, ScMatrixRef pMatY, double fSlope,
- SCSIZE nN)
+ SCSIZE nN)
{
double fSum = 0.0;
double fTemp = 0.0;
@@ -2341,6 +2340,7 @@ bool ScInterpreter::CheckMatrix(bool _bLOG, BYTE& nCase, SCSIZE& nCX,
}
return true;
}
+
// -----------------------------------------------------------------------------
// LINEST
@@ -2360,7 +2360,7 @@ void ScInterpreter::ScRKP()
void ScInterpreter::CalulateRGPRKP(bool _bRKP)
{
BYTE nParamCount = GetByte();
- if ( !MustHaveParamCount( nParamCount, 1, 4 ) )
+ if (!MustHaveParamCount( nParamCount, 1, 4 ))
return;
bool bConstant, bStats;
@@ -2382,7 +2382,7 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
// return;
}
else
- bConstant = GetBool();
+ bConstant = GetBool();
}
else
bConstant = true;
@@ -2397,7 +2397,7 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
}
else
{
- pMatX = GetMatrix();
+ pMatX = GetMatrix();
}
}
else
@@ -2410,46 +2410,46 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
PushIllegalParameter();
return;
}
-
+
// 1 = simple; 2 = multiple with Y as column; 3 = multiple with Y as row
BYTE nCase;
SCSIZE nCX, nCY; // number of columns
SCSIZE nRX, nRY; //number of rows
SCSIZE K = 0, N = 0; // K=number of variables X, N=number of data samples
- if ( !CheckMatrix(_bRKP,nCase,nCX,nCY,nRX,nRY,K,N,pMatX,pMatY) )
- {
+ if (!CheckMatrix(_bRKP,nCase,nCX,nCY,nRX,nRY,K,N,pMatX,pMatY))
+ {
PushIllegalParameter();
- return;
- }
+ return;
+ }
// Enough data samples?
- if ( (bConstant && (N<K+1)) || (!bConstant && (N<K)) || (N<1) || (K<1) )
- {
+ if ((bConstant && (N<K+1)) || (!bConstant && (N<K)) || (N<1) || (K<1))
+ {
PushIllegalParameter();
return;
- }
+ }
ScMatrixRef pResMat;
if (bStats)
pResMat = GetNewMat(K+1,5);
- else
+ else
pResMat = GetNewMat(K+1,1);
if (!pResMat)
- {
+ {
PushError(errCodeOverflow);
return;
- }
+ }
// Fill unused cells in pResMat; order (column,row)
- if (bStats)
- {
+ if (bStats)
+ {
for (SCSIZE i=2; i<K+1; i++)
- {
+ {
pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 2 );
pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 3 );
pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 4 );
- }
- }
+ }
+ }
// Uses sum(x-MeanX)^2 and not [sum x^2]-N * MeanX^2 in case bConstant.
// Clone constant matrices, so that Mat = Mat - Mean is possible.
@@ -2468,30 +2468,30 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
// DeltaY is possible here; DeltaX depends on nCase, so later
fMeanY = lcl_GetMeanOverAll(pMatY, N);
for (SCSIZE i=0; i<N; i++)
- {
+ {
pMatY->PutDouble( ::rtl::math::approxSub(pMatY->GetDouble(i),fMeanY), i );
}
}
if (nCase==1)
- {
+ {
// calculate simple regression
double fMeanX = 0.0;
if (bConstant)
{ // Mat = Mat - Mean
fMeanX = lcl_GetMeanOverAll(pMatX, N);
for (SCSIZE i=0; i<N; i++)
- {
+ {
pMatX->PutDouble( ::rtl::math::approxSub(pMatX->GetDouble(i),fMeanX), i );
- }
- }
+ }
+ }
double fSumXY = lcl_GetSumProduct(pMatX,pMatY,N);
double fSumX2 = lcl_GetSumProduct(pMatX,pMatX,N);
if (fSumX2==0.0)
{
PushNoValue(); // all x-values are identical
return;
- }
+ }
double fSlope = fSumXY / fSumX2;
double fIntercept = 0.0;
if (bConstant)
@@ -2522,11 +2522,11 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
else
pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 1, 1);
pResMat->PutDouble(1.0, 0, 2); // R^2
- }
- else
- {
+ }
+ else
+ {
double fFstatistic = (fSSreg / static_cast<double>(K))
- / (fSSresid / fDegreesFreedom);
+ / (fSSresid / fDegreesFreedom);
pResMat->PutDouble(fFstatistic, 0, 3);
// standard error of estimate
@@ -2537,24 +2537,24 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
pResMat->PutDouble(fSigmaSlope, 0, 1);
if (bConstant)
- {
+ {
double fSigmaIntercept = fRMSE
- * sqrt(fMeanX*fMeanX/fSumX2 + 1.0/static_cast<double>(N));
+ * sqrt(fMeanX*fMeanX/fSumX2 + 1.0/static_cast<double>(N));
pResMat->PutDouble(fSigmaIntercept, 1, 1);
- }
+ }
else
{
pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 1, 1);
- }
-
+ }
+
double fR2 = fSSreg / (fSSreg + fSSresid);
pResMat->PutDouble(fR2, 0, 2);
}
+ }
+ PushMatrix(pResMat);
}
- PushMatrix(pResMat);
-}
else // calculate multiple regression;
-{
+ {
// Uses a QR decomposition X = QR. The solution B = (X'X)^(-1) * X' * Y
// becomes B = R^(-1) * Q' * Y
if (nCase ==2) // Y is column
@@ -2569,12 +2569,12 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
pMatZ = pMatY; // Y can be overwritten
ScMatrixRef pSlopes = GetNewMat(1,K); // from b1 to bK
if (!pMeans || !pMatZ || !pSlopes)
-{
+ {
PushError(errCodeOverflow);
return;
}
- if (bConstant)
- {
+ if (bConstant)
+ {
lcl_CalculateColumnMeans(pMatX, pMeans, K, N);
lcl_CalculateColumnsDelta(pMatX, pMeans, K, N);
}
@@ -2589,21 +2589,21 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
for (SCSIZE row=0; row < K && !bIsSingular; row++)
bIsSingular = bIsSingular || aVecR[row]==0.0;
if (bIsSingular)
- {
+ {
PushNoValue();
return;
- }
+ }
// Z = Q' Y;
for (SCSIZE col = 0; col < K; col++)
- {
+ {
lcl_ApplyHouseholderTransformation(pMatX, col, pMatZ, N);
- }
+ }
// B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z
// result Z should have zeros for index>=K; if not, ignore values
for (SCSIZE col = 0; col < K ; col++)
- {
+ {
pSlopes->PutDouble( pMatZ->GetDouble(col), col);
-}
+ }
lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, false);
double fIntercept = 0.0;
if (bConstant)
@@ -2612,11 +2612,11 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, K, 0 );
for (SCSIZE i = 0; i < K; i++)
pResMat->PutDouble(_bRKP ? exp(pSlopes->GetDouble(i))
- : pSlopes->GetDouble(i) , K-1-i, 0);
+ : pSlopes->GetDouble(i) , K-1-i, 0);
if (bStats)
- {
+ {
double fSSreg = 0.0;
double fSSresid = 0.0;
// re-use memory of Z;
@@ -2660,9 +2660,9 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
pResMat->PutDouble(1.0, 0, 2); // R^2
}
else
- {
+ {
double fFstatistic = (fSSreg / static_cast<double>(K))
- / (fSSresid / fDegreesFreedom);
+ / (fSSresid / fDegreesFreedom);
pResMat->PutDouble(fFstatistic, 0, 3);
// standard error of estimate = root mean SSE
@@ -2692,28 +2692,28 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
pResMat->PutDouble(fSigmaSlope, K-1-col, 1);
// (R' R) ^(-1) is symmetric
if (bConstant)
- {
+ {
fPart = lcl_GetSumProduct(pMeans, pMatZ, K);
fSigmaIntercept += fPart * pMeans->GetDouble(col);
- }
- }
+ }
+ }
if (bConstant)
{
fSigmaIntercept = fRMSE
- * sqrt(fSigmaIntercept + 1.0 / static_cast<double>(N));
+ * sqrt(fSigmaIntercept + 1.0 / static_cast<double>(N));
pResMat->PutDouble(fSigmaIntercept, K, 1);
- }
- else
- {
+ }
+ else
+ {
pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1);
- }
+ }
double fR2 = fSSreg / (fSSreg + fSSresid);
pResMat->PutDouble(fR2, 0, 2);
+ }
}
+ PushMatrix(pResMat);
}
- PushMatrix(pResMat);
- }
else // nCase == 3, Y is row, all matrices are transposed
{
::std::vector< double> aVecR(N); // for QR decomposition
@@ -2748,10 +2748,10 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
for (SCSIZE row=0; row < K && !bIsSingular; row++)
bIsSingular = bIsSingular || aVecR[row]==0.0;
if (bIsSingular)
- {
+ {
PushNoValue();
return;
- }
+ }
// Z = Q' Y
for (SCSIZE row = 0; row < K; row++)
{
@@ -2762,7 +2762,7 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
for (SCSIZE col = 0; col < K ; col++)
{
pSlopes->PutDouble( pMatZ->GetDouble(col), col);
- }
+ }
lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, true);
double fIntercept = 0.0;
if (bConstant)
@@ -2771,7 +2771,7 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, K, 0 );
for (SCSIZE i = 0; i < K; i++)
pResMat->PutDouble(_bRKP ? exp(pSlopes->GetDouble(i))
- : pSlopes->GetDouble(i) , K-1-i, 0);
+ : pSlopes->GetDouble(i) , K-1-i, 0);
if (bStats)
@@ -2817,11 +2817,11 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
// R^2 = SSreg / (SSreg + SSresid) = 1.0
pResMat->PutDouble(1.0, 0, 2); // R^2
- }
- else
- {
+ }
+ else
+ {
double fFstatistic = (fSSreg / static_cast<double>(K))
- / (fSSresid / fDegreesFreedom);
+ / (fSSresid / fDegreesFreedom);
pResMat->PutDouble(fFstatistic, 0, 3);
// standard error of estimate = root mean SSE
@@ -2854,24 +2854,24 @@ void ScInterpreter::CalulateRGPRKP(bool _bRKP)
{
fPart = lcl_GetSumProduct(pMeans, pMatZ, K);
fSigmaIntercept += fPart * pMeans->GetDouble(row);
- }
-}
+ }
+ }
if (bConstant)
- {
+ {
fSigmaIntercept = fRMSE
- * sqrt(fSigmaIntercept + 1.0 / static_cast<double>(N));
+ * sqrt(fSigmaIntercept + 1.0 / static_cast<double>(N));
pResMat->PutDouble(fSigmaIntercept, K, 1);
- }
+ }
else
- {
+ {
pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1);
- }
+ }
double fR2 = fSSreg / (fSSreg + fSSresid);
pResMat->PutDouble(fR2, 0, 2);
+ }
}
- }
- PushMatrix(pResMat);
+ PushMatrix(pResMat);
}
}
return;
@@ -2892,7 +2892,7 @@ void ScInterpreter::ScGrowth()
void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
{
BYTE nParamCount = GetByte();
- if ( !MustHaveParamCount( nParamCount, 1, 4 ) )
+ if (!MustHaveParamCount( nParamCount, 1, 4 ))
return;
// optional forth parameter
@@ -2913,7 +2913,7 @@ void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
pMatNewX = NULL;
}
else
- pMatNewX = GetMatrix();
+ pMatNewX = GetMatrix();
}
else
pMatNewX = NULL;
@@ -2930,7 +2930,7 @@ void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
}
else
{
- pMatX = GetMatrix();
+ pMatX = GetMatrix();
}
}
else
@@ -2950,19 +2950,19 @@ void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
SCSIZE nCX, nCY; // number of columns
SCSIZE nRX, nRY; //number of rows
SCSIZE K = 0, N = 0; // K=number of variables X, N=number of data samples
- if ( !CheckMatrix(_bGrowth,nCase,nCX,nCY,nRX,nRY,K,N,pMatX,pMatY) )
+ if (!CheckMatrix(_bGrowth,nCase,nCX,nCY,nRX,nRY,K,N,pMatX,pMatY))
{
PushIllegalParameter();
return;
}
-
+
// Enough data samples?
- if ( (bConstant && (N<K+1)) || (!bConstant && (N<K)) || (N<1) || (K<1) )
+ if ((bConstant && (N<K+1)) || (!bConstant && (N<K)) || (N<1) || (K<1))
{
PushIllegalParameter();
return;
}
-
+
// Set default pMatNewX if necessary
SCSIZE nCXN, nRXN;
SCSIZE nCountXN;
@@ -2982,7 +2982,7 @@ void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
return;
}
nCountXN = nCXN * nRXN;
- for ( SCSIZE i = 0; i < nCountXN; i++ )
+ for (SCSIZE i = 0; i < nCountXN; i++)
if (!pMatNewX->IsValue(i))
{
PushIllegalArgument();
@@ -2991,7 +2991,7 @@ void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
}
ScMatrixRef pResMat; // size depends on nCase
if (nCase == 1)
- pResMat = GetNewMat(nCXN,nRXN);
+ pResMat = GetNewMat(nCXN,nRXN);
else
{
if (nCase==2)
@@ -3000,15 +3000,15 @@ void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
pResMat = GetNewMat(nCXN,1);
}
if (!pResMat)
- {
+ {
PushError(errCodeOverflow);
return;
- }
+ }
// Uses sum(x-MeanX)^2 and not [sum x^2]-N * MeanX^2 in case bConstant.
// Clone constant matrices, so that Mat = Mat - Mean is possible.
double fMeanY = 0.0;
if (bConstant)
- {
+ {
ScMatrixRef pCopyX = pMatX->CloneIfConst();
ScMatrixRef pCopyY = pMatY->CloneIfConst();
if (!pCopyX || !pCopyY)
@@ -3030,7 +3030,7 @@ void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
{
// calculate simple regression
double fMeanX = 0.0;
- if (bConstant)
+ if (bConstant)
{ // Mat = Mat - Mean
fMeanX = lcl_GetMeanOverAll(pMatX, N);
for (SCSIZE i=0; i<N; i++)
@@ -3056,9 +3056,9 @@ void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
fHelp = pMatNewX->GetDouble(i)*fSlope + fIntercept;
pResMat->PutDouble(_bGrowth ? exp(fHelp) : fHelp, i);
}
- }
- else
- {
+ }
+ else
+ {
for (SCSIZE i = 0; i < nCountXN; i++)
{
fHelp = pMatNewX->GetDouble(i)*fSlope;
@@ -3067,7 +3067,7 @@ void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
}
}
else // calculate multiple regression;
- {
+ {
if (nCase ==2) // Y is column
{
::std::vector< double> aVecR(N); // for QR decomposition
@@ -3102,13 +3102,13 @@ void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
// Z := Q' Y; Y is overwritten with result Z
for (SCSIZE col = 0; col < K; col++)
{
- lcl_ApplyHouseholderTransformation(pMatX, col, pMatY, N);
+ lcl_ApplyHouseholderTransformation(pMatX, col, pMatY, N);
}
// B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z
// result Z should have zeros for index>=K; if not, ignore values
for (SCSIZE col = 0; col < K ; col++)
{
- pSlopes->PutDouble( pMatY->GetDouble(col), col);
+ pSlopes->PutDouble( pMatY->GetDouble(col), col);
}
lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, false);
@@ -3124,9 +3124,9 @@ void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
{
for (SCSIZE i = 0; i < nRXN; i++)
pResMat->PutDouble(exp(pResMat->GetDouble(i)), i);
+ }
}
- }
- else
+ else
{ // nCase == 3, Y is row, all matrices are transposed
::std::vector< double> aVecR(N); // for QR decomposition
@@ -3134,12 +3134,12 @@ void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
ScMatrixRef pMeans = GetNewMat(1, K); // mean of each row
ScMatrixRef pSlopes = GetNewMat(K,1); // row from b1 to bK
if (!pMeans || !pSlopes)
- {
+ {
PushError(errCodeOverflow);
- return;
+ return;
}
if (bConstant)
- {
+ {
lcl_CalculateRowMeans(pMatX, pMeans, N, K);
lcl_CalculateRowsDelta(pMatX, pMeans, N, K);
}
@@ -3161,13 +3161,13 @@ void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
// Z := Q' Y; Y is overwritten with result Z
for (SCSIZE row = 0; row < K; row++)
{
- lcl_TApplyHouseholderTransformation(pMatX, row, pMatY, N);
- }
+ lcl_TApplyHouseholderTransformation(pMatX, row, pMatY, N);
+ }
// B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z
// result Z should have zeros for index>=K; if not, ignore values
for (SCSIZE col = 0; col < K ; col++)
- {
- pSlopes->PutDouble( pMatY->GetDouble(col), col);
+ {
+ pSlopes->PutDouble( pMatY->GetDouble(col), col);
}
lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, true);
commit 666817d535aa5b56733fc6aad0e344fa68bdef66
Author: Kohei Yoshida <kyoshida at novell.com>
Date: Thu Nov 4 18:01:45 2010 -0400
Put local functions inside anonymous namespace.
diff --git a/sc/source/core/tool/interpr5.cxx b/sc/source/core/tool/interpr5.cxx
index f3dcd5c..93763ec 100644
--- a/sc/source/core/tool/interpr5.cxx
+++ b/sc/source/core/tool/interpr5.cxx
@@ -1843,6 +1843,9 @@ void ScInterpreter::ScFrequency()
pResMat->PutDouble(static_cast<double>(nDataSize-i), j);
PushMatrix(pResMat);
}
+
+namespace {
+
// -----------------------------------------------------------------------------
// Helper methods for LINEST/LOGEST and TREND/GROWTH
// All matrices must already exist and have the needed size, no control tests
@@ -2229,6 +2232,8 @@ double lcl_GetSSresid(ScMatrixRef pMatX, ScMatrixRef pMatY, double fSlope,
return fSum;
}
+}
+
// Fill default values in matrix X, transform Y to log(Y) in case LOGEST|GROWTH,
// determine sizes of matrices X and Y, determine kind of regression, clone
// Y in case LOGEST|GROWTH, if constant.
commit 94b0dade627b89935530fe6f23cc2f859a510f4c
Author: Regina Henschel <rb.henschel at t-online.de>
Date: Thu Nov 4 17:45:33 2010 -0400
Patch for LINEST|LOGEST and TREND|GROWTH.
diff --git a/sc/source/core/inc/interpre.hxx b/sc/source/core/inc/interpre.hxx
index f356dbe..07f95ba 100644
--- a/sc/source/core/inc/interpre.hxx
+++ b/sc/source/core/inc/interpre.hxx
@@ -664,7 +664,6 @@ void ScLCM();
void ScMatValue();
void MEMat(ScMatrix* mM, SCSIZE n);
-void MFastMult(ScMatrix* pA, ScMatrix* pB, ScMatrix* pR, SCSIZE n, SCSIZE m, SCSIZE l);
void ScMatDet();
void ScMatInv();
void ScMatMult();
@@ -677,13 +676,6 @@ void ScSumX2MY2();
void ScSumX2DY2();
void ScSumXMY2();
void ScGrowth();
-// multiple Regression: Varianzen der Koeffizienten
-BOOL RGetVariances( ScMatrix* pV, ScMatrix* pX, SCSIZE nC, SCSIZE nR,
- BOOL bSwapColRow, BOOL bZeroConstant );
-void Calculate(ScMatrixRef& pResMat,ScMatrixRef& pE,ScMatrixRef& pQ,ScMatrixRef& pV,ScMatrixRef& pMatX,BOOL bConstant,SCSIZE N,SCSIZE M,BYTE nCase);
-ScMatrixRef Calculate2(const BOOL bConstant,const SCSIZE M ,const SCSIZE N,ScMatrixRef& pMatX,ScMatrixRef& pMatY,BYTE nCase);
-bool Calculate3(const SCSIZE M ,ScMatrixRef& pQ);
-bool Calculate4(BOOL _bExp,ScMatrixRef& pResMat,ScMatrixRef& pQ,BOOL bConstant,SCSIZE N,SCSIZE M);
bool CalculateSkew(double& fSum,double& fCount,double& vSum,std::vector<double>& values);
void CalculateSlopeIntercept(BOOL bSlope);
void CalculateSmallLarge(BOOL bSmall);
@@ -695,12 +687,11 @@ bool CalculateTest( BOOL _bTemplin
void CalculateLookup(BOOL HLookup);
bool FillEntry(ScQueryEntry& rEntry);
void CalculateAddSub(BOOL _bSub);
-void CalculateTrendGrowth(BOOL _bGrowth);
-void CalulateRGPRKP(BOOL _bRKP);
+void CalculateTrendGrowth(bool _bGrowth);
+void CalulateRGPRKP(bool _bRKP);
void CalculateSumX2MY2SumX2DY2(BOOL _bSumX2DY2);
void CalculateMatrixValue(const ScMatrix* pMat,SCSIZE nC,SCSIZE nR);
-bool CheckMatrix(BOOL _bLOG,BOOL _bTrendGrowth,BYTE& nCase,SCSIZE& nCX,SCSIZE& nCY,SCSIZE& nRX,SCSIZE& nRY,SCSIZE& M,SCSIZE& N,ScMatrixRef& pMatX,ScMatrixRef& pMatY);
-
+bool CheckMatrix(bool _bLOG,BYTE& nCase,SCSIZE& nCX,SCSIZE& nCY,SCSIZE& nRX,SCSIZE& nRY,SCSIZE& M,SCSIZE& N,ScMatrixRef& pMatX,ScMatrixRef& pMatY);
void ScRGP();
void ScRKP();
void ScForecast();
diff --git a/sc/source/core/tool/interpr4.cxx b/sc/source/core/tool/interpr4.cxx
index a1be445..e0ecf26 100644
--- a/sc/source/core/tool/interpr4.cxx
+++ b/sc/source/core/tool/interpr4.cxx
@@ -1623,7 +1623,7 @@ bool ScInterpreter::ConvertMatrixParameters()
for ( USHORT i=1; i <= nParams && i <= sp; ++i )
{
FormulaToken* p = pStack[ sp - i ];
- if ( p->GetOpCode() != ocPush )
+ if ( p->GetOpCode() != ocPush && p->GetOpCode() != ocMissing)
{
DBG_ERRORFILE( "ConvertMatrixParameters: not a push");
}
diff --git a/sc/source/core/tool/interpr5.cxx b/sc/source/core/tool/interpr5.cxx
index ff72817..f3dcd5c 100644
--- a/sc/source/core/tool/interpr5.cxx
+++ b/sc/source/core/tool/interpr5.cxx
@@ -660,25 +660,6 @@ void ScInterpreter::MEMat(ScMatrix* mM, SCSIZE n)
mM->PutDouble(1.0, i, i);
}
-void ScInterpreter::MFastMult(ScMatrix* pA, ScMatrix* pB, ScMatrix* pR,
- SCSIZE n, SCSIZE m, SCSIZE l)
- // Multipliziert n x m Mat a mit m x l Mat b nach Mat r
-{
- RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::MFastMult" );
- double sum;
- for (SCSIZE i = 0; i < n; i++)
- {
- for (SCSIZE j = 0; j < l; j++)
- {
- sum = 0.0;
- for (SCSIZE k = 0; k < m; k++)
- sum += pA->GetDouble(i,k)*pB->GetDouble(k,j);
- pR->PutDouble(sum, i, j);
- }
- }
-}
-
-
/* Matrix LUP decomposition according to the pseudocode of "Introduction to
* Algorithms" by Cormen, Leiserson, Rivest, Stein.
*
@@ -783,6 +764,13 @@ static int lcl_LUP_decompose( ScMatrix* mA, const SCSIZE n,
fprintf( stderr, "%5u ", (unsigned)P[j]);
fprintf( stderr, "\n%s\n", "");
#endif
+
+ bool bSingular=false;
+ for (SCSIZE i=0; i<n && !bSingular; i++)
+ bSingular = bSingular || ((mA->GetDouble(i,i))==0.0);
+ if (bSingular)
+ nSign = 0;
+
return nSign;
}
@@ -950,7 +938,7 @@ void ScInterpreter::ScMatInv()
if (xR)
{
ScMatrix* pR = xR;
- MFastMult( pMat, pY, pR, nR, nR, nR);
+ lcl_MFastMult( pMat, pY, pR, nR, nR, nR);
#if OSL_DEBUG_LEVEL > 1
fprintf( stderr, "\n%s\n", "ScMatInv(): mult-identity");
#endif
@@ -1855,213 +1843,400 @@ void ScInterpreter::ScFrequency()
pResMat->PutDouble(static_cast<double>(nDataSize-i), j);
PushMatrix(pResMat);
}
+// -----------------------------------------------------------------------------
+// Helper methods for LINEST/LOGEST and TREND/GROWTH
+// All matrices must already exist and have the needed size, no control tests
+// done. Those methodes, which names start with lcl_T, are adapted to case 3,
+// where Y (=observed values) is given as row.
+// Remember, ScMatrix matrices are zero based, index access (column,row).
+// -----------------------------------------------------------------------------
-BOOL ScInterpreter::RGetVariances( ScMatrix* pV, ScMatrix* pX,
- SCSIZE nC, SCSIZE nR, BOOL bSwapColRow, BOOL bZeroConstant )
-{ // multiple Regression: Varianzen der Koeffizienten
- // bSwapColRow==TRUE : Koeffizienten in Zeilen statt Spalten angeordnet
- SCSIZE i, j, k;
- double sum;
- ScMatrixRef pC = GetNewMat(nC, nC);
- if ( !pC )
- return FALSE;
- // X transformiert mit X multipziert, X'X Matrix
- if ( !bZeroConstant )
- { // in der X-Designmatrix existiert ein gedachtes X0j==1
- if ( bSwapColRow )
- {
- for ( i=0; i<nC; i++ )
+// Multiply n x m Mat A with m x l Mat B to n x l Mat R
+void lcl_MFastMult(ScMatrixRef pA, ScMatrixRef pB, ScMatrixRef pR,
+ SCSIZE n, SCSIZE m, SCSIZE l)
{
- for ( j=0; j<nC; j++ )
+ double sum;
+ for (SCSIZE row = 0; row < n; row++)
{
+ for (SCSIZE col = 0; col < l; col++)
+ { // result element(col, row) =sum[ (row of A) * (column of B)]
sum = 0.0;
- for ( k=0; k<nR; k++ )
- {
- sum += (j==0 ? 1 : pX->GetDouble(k,j-1))
- * (i==0 ? 1 : pX->GetDouble(k,i-1));
- }
- pC->PutDouble(sum, i, j);
+ for (SCSIZE k = 0; k < m; k++)
+ sum += pA->GetDouble(k,row) * pB->GetDouble(col,k);
+ pR->PutDouble(sum, col, row);
}
}
}
- else
- {
- for ( i=0; i<nC; i++ )
- {
- for ( j=0; j<nC; j++ )
- {
- sum = 0.0;
- for ( k=0; k<nR; k++ )
+
+// <A;B> over all elements; uses the matrices as vectors of length M
+double lcl_GetSumProduct(ScMatrixRef pMatA, ScMatrixRef pMatB, SCSIZE nM)
{
- sum += (j==0 ? 1 : pX->GetDouble(j-1,k))
- * (i==0 ? 1 : pX->GetDouble(i-1,k));
+ double fSum = 0.0;
+ for (SCSIZE i=0; i<nM; i++)
+ fSum += pMatA->GetDouble(i) * pMatB->GetDouble(i);
+ return fSum;
}
- pC->PutDouble(sum, i, j);
+
+// Special version for use within QR decomposition.
+// Euclidean norm of column index C starting in row index R;
+// matrix A has count N rows.
+double lcl_GetColumnEuclideanNorm(ScMatrixRef pMatA, SCSIZE nC, SCSIZE nR, SCSIZE nN)
+{
+ double fNorm = 0.0;
+ for (SCSIZE row=nR; row<nN; row++)
+ fNorm += (pMatA->GetDouble(nC,row)) * (pMatA->GetDouble(nC,row));
+ return sqrt(fNorm);
}
+
+// Euclidean norm of row index R starting in column index C;
+// matrix A has count N columns.
+double lcl_TGetColumnEuclideanNorm(ScMatrixRef pMatA, SCSIZE nR, SCSIZE nC, SCSIZE nN)
+{
+ double fNorm = 0.0;
+ for (SCSIZE col=nC; col<nN; col++)
+ fNorm += (pMatA->GetDouble(col,nR)) * (pMatA->GetDouble(col,nR));
+ return sqrt(fNorm);
}
+
+// Special version for use within QR decomposition.
+// Maximum norm of column index C starting in row index R;
+// matrix A has count N rows.
+double lcl_GetColumnMaximumNorm(ScMatrixRef pMatA, SCSIZE nC, SCSIZE nR, SCSIZE nN)
+{
+ double fNorm = 0.0;
+ for (SCSIZE row=nR; row<nN; row++)
+ if (fNorm < fabs(pMatA->GetDouble(nC,row)))
+ fNorm = fabs(pMatA->GetDouble(nC,row));
+ return fNorm;
}
- }
- else
- {
- if ( bSwapColRow )
+
+// Maximum norm of row index R starting in col index C;
+// matrix A has count N columns.
+double lcl_TGetColumnMaximumNorm(ScMatrixRef pMatA, SCSIZE nR, SCSIZE nC, SCSIZE nN)
+{
+ double fNorm = 0.0;
+ for (SCSIZE col=nC; col<nN; col++)
+ if (fNorm < fabs(pMatA->GetDouble(col,nR)))
+ fNorm = fabs(pMatA->GetDouble(col,nR));
+ return fNorm;
+ }
+
+// Special version for use within QR decomposition.
+// <A(Ca);B(Cb)> starting in row index R;
+// Ca and Cb are indices of columns, matrices A and B have count N rows.
+double lcl_GetColumnSumProduct(ScMatrixRef pMatA, SCSIZE nCa,
+ ScMatrixRef pMatB, SCSIZE nCb, SCSIZE nR, SCSIZE nN)
+ {
+ double fResult = 0.0;
+ for (SCSIZE row=nR; row<nN; row++)
+ fResult += pMatA->GetDouble(nCa,row) * pMatB->GetDouble(nCb,row);
+ return fResult;
+}
+
+// <A(Ra);B(Rb)> starting in column index C;
+// Ra and Rb are indices of rows, matrices A and B have count N columns.
+double lcl_TGetColumnSumProduct(ScMatrixRef pMatA, SCSIZE nRa,
+ ScMatrixRef pMatB, SCSIZE nRb, SCSIZE nC, SCSIZE nN)
{
- for ( i=0; i<nC; i++ )
+ double fResult = 0.0;
+ for (SCSIZE col=nC; col<nN; col++)
+ fResult += pMatA->GetDouble(col,nRa) * pMatB->GetDouble(col,nRb);
+ return fResult;
+}
+
+double lcl_GetSign(double fValue)
{
- for ( j=0; j<nC; j++ )
+ if (fValue < 0.0)
+ return -1.0;
+ else if (fValue > 0.0)
+ return 1.0;
+ else
+ return 0.0;
+}
+
+/* Calculates a QR decomposition with Householder reflection.
+ * For each NxK matrix A exists a decomposition A=Q*R with an orthogonal
+ * NxN matrix Q and a NxK matrix R.
+ * Q=H1*H2*...*Hk with Householder matrices H. Such a householder matrix can
+ * be build from a vector u by H=I-(2/u'u)*(u u'). This vectors u are returned
+ * in the columns of matrix A, overwriting the old content.
+ * The matrix R has a quadric upper part KxK with values in the upper right
+ * triangle and zeros in all other elements. Here the diagonal elements of R
+ * are stored in the vector R and the other upper right elements in the upper
+ * right of the matrix A.
+ * The function returns false, if calculation breaks. But because of round-off
+ * errors singularity is often not detected.
+ */
+bool lcl_CalculateQRdecomposition(ScMatrixRef pMatA,
+ ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN)
{
- sum = 0.0;
- for ( k=0; k<nR; k++ )
+ double fScale ;
+ double fEuclid ;
+ double fFactor ;
+ double fSignum ;
+ double fSum ;
+ // ScMatrix matrices are zero based, index access (column,row)
+ for (SCSIZE col = 0; col <nK; col++)
+ {
+ // calculate vector u of the householder transformation
+ fScale = lcl_GetColumnMaximumNorm(pMatA, col, col, nN);
+ if (fScale == 0.0)
{
- sum += pX->GetDouble(k,j) * pX->GetDouble(k,i);
+ // A is singular
+ return false;
}
- pC->PutDouble(sum, i, j);
+ for (SCSIZE row = col; row <nN; row++)
+ pMatA->PutDouble( pMatA->GetDouble(col,row)/fScale, col, row);
+
+ fEuclid = lcl_GetColumnEuclideanNorm(pMatA, col, col, nN);
+ fFactor = 1.0/fEuclid/(fEuclid + fabs(pMatA->GetDouble(col,col)));
+ fSignum = lcl_GetSign(pMatA->GetDouble(col,col));
+ pMatA->PutDouble( pMatA->GetDouble(col,col) + fSignum*fEuclid, col,col);
+ pVecR[col] = -fSignum * fScale * fEuclid;
+
+ // apply Householder transformation to A
+ for (SCSIZE c=col+1; c<nK; c++)
+ {
+ fSum =lcl_GetColumnSumProduct(pMatA, col, pMatA, c, col, nN);
+ for (SCSIZE row = col; row <nN; row++)
+ pMatA->PutDouble( pMatA->GetDouble(c,row)
+ - fSum * fFactor * pMatA->GetDouble(col,row), c, row);
}
}
+ return true;
}
- else
- {
- for ( i=0; i<nC; i++ )
- {
- for ( j=0; j<nC; j++ )
- {
- sum = 0.0;
- for ( k=0; k<nR; k++ )
+
+// same with transposed matrix A, N is count of columns, K count of rows
+bool lcl_TCalculateQRdecomposition(ScMatrixRef pMatA,
+ ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN)
+{
+ double fScale ;
+ double fEuclid ;
+ double fFactor ;
+ double fSignum ;
+ double fSum ;
+ // ScMatrix matrices are zero based, index access (column,row)
+ for (SCSIZE row = 0; row <nK; row++)
+ {
+ // calculate vector u of the householder transformation
+ fScale = lcl_TGetColumnMaximumNorm(pMatA, row, row, nN);
+ if (fScale == 0.0)
{
- sum += pX->GetDouble(j,k) * pX->GetDouble(i,k);
- }
- pC->PutDouble(sum, i, j);
- }
+ // A is singular
+ return false;
}
+ for (SCSIZE col = row; col <nN; col++)
+ pMatA->PutDouble( pMatA->GetDouble(col,row)/fScale, col, row);
+
+ fEuclid = lcl_TGetColumnEuclideanNorm(pMatA, row, row, nN);
+ fFactor = 1.0/fEuclid/(fEuclid + fabs(pMatA->GetDouble(row,row)));
+ fSignum = lcl_GetSign(pMatA->GetDouble(row,row));
+ pMatA->PutDouble( pMatA->GetDouble(row,row) + fSignum*fEuclid, row,row);
+ pVecR[row] = -fSignum * fScale * fEuclid;
+
+ // apply Householder transformation to A
+ for (SCSIZE r=row+1; r<nK; r++)
+ {
+ fSum =lcl_TGetColumnSumProduct(pMatA, row, pMatA, r, row, nN);
+ for (SCSIZE col = row; col <nN; col++)
+ pMatA->PutDouble( pMatA->GetDouble(col,r)
+ - fSum * fFactor * pMatA->GetDouble(col,row), col, r);
}
}
- // X'X Inverse
- BOOL bOk = TRUE;
- USHORT nErr = nGlobalError;
- PushMatrix(pC);
- BYTE nTmp = cPar;
- cPar = 1;
- ScMatInv();
- cPar = nTmp;
- if ( nGlobalError )
- {
- nGlobalError = nErr;
- bOk = FALSE;
+ return true;
}
- else
- {
- // #i61216# ScMatInv no longer modifies the original matrix, so just calling Pop() doesn't work
- pC = PopMatrix();
- if ( pC.Is() )
- {
- // Varianzen auf der Diagonalen, andere sind Kovarianzen
- for (i = 0; i < nC; i++)
- pV->PutDouble(pC->GetDouble(i, i), i);
- }
+
+
+/* Applies a Householder transformation to a column vector Y with is given as
+ * Nx1 Matrix. The Vektor u, from which the Householder transformation is build,
+ * is the column part in matrix A, with column index C, starting with row
+ * index C. A is the result of the QR decomposition as obtained from
+ * lcl_CaluclateQRdecomposition.
+ */
+void lcl_ApplyHouseholderTransformation(ScMatrixRef pMatA, SCSIZE nC,
+ ScMatrixRef pMatY, SCSIZE nN)
+ {
+ // ScMatrix matrices are zero based, index access (column,row)
+ double fDenominator = lcl_GetColumnSumProduct(pMatA, nC, pMatA, nC, nC, nN);
+ double fNumerator = lcl_GetColumnSumProduct(pMatA, nC, pMatY, 0, nC, nN);
+ double fFactor = 2.0 * (fNumerator/fDenominator);
+ for (SCSIZE row = nC; row < nN; row++)
+ pMatY->PutDouble(
+ pMatY->GetDouble(row) - fFactor * pMatA->GetDouble(nC,row), row);
+}
+
+// Same with transposed matrices A and Y.
+void lcl_TApplyHouseholderTransformation(ScMatrixRef pMatA, SCSIZE nR,
+ ScMatrixRef pMatY, SCSIZE nN)
+{
+ // ScMatrix matrices are zero based, index access (column,row)
+ double fDenominator = lcl_TGetColumnSumProduct(pMatA, nR, pMatA, nR, nR, nN);
+ double fNumerator = lcl_TGetColumnSumProduct(pMatA, nR, pMatY, 0, nR, nN);
+ double fFactor = 2.0 * (fNumerator/fDenominator);
+ for (SCSIZE col = nR; col < nN; col++)
+ pMatY->PutDouble(
+ pMatY->GetDouble(col) - fFactor * pMatA->GetDouble(col,nR), col);
+}
+
+/* Solve for X in R*X=S using back substitution. The solution X overwrites S.
+ * Uses R from the result of the QR decomposition of a NxK matrix A.
+ * S is a column vector given as matrix, with at least elements on index
+ * 0 to K-1; elements on index>=K are ignored. Vector R must not have zero
+ * elements, no check is done.
+ */
+void lcl_SolveWithUpperRightTriangle(ScMatrixRef pMatA,
+ ::std::vector< double>& pVecR, ScMatrixRef pMatS,
+ SCSIZE nK, bool bIsTransposed)
+{
+ // ScMatrix matrices are zero based, index access (column,row)
+ double fSum;
+ SCSIZE row;
+ // SCSIZE is never negative, therefore test with rowp1=row+1
+ for (SCSIZE rowp1 = nK; rowp1>0; rowp1--)
+ {
+ row = rowp1-1;
+ fSum = pMatS->GetDouble(row);
+ for (SCSIZE col = rowp1; col<nK ; col++)
+ if (bIsTransposed)
+ fSum -= pMatA->GetDouble(row,col) * pMatS->GetDouble(col);
+ else
+ fSum -= pMatA->GetDouble(col,row) * pMatS->GetDouble(col);
+ pMatS->PutDouble( fSum / pVecR[row] , row);
}
- return bOk;
}
-// -----------------------------------------------------------------------------
-void ScInterpreter::Calculate(ScMatrixRef& pResMat,ScMatrixRef& pE,ScMatrixRef& pQ,ScMatrixRef& pV,ScMatrixRef& pMatX,BOOL bConstant,SCSIZE N,SCSIZE M,BYTE nCase)
+
+/* Solve for X in R' * X= T using forward substitution. The solution X
+ * overwrites T. Uses R from the result of the QR decomposition of a NxK
+ * matrix A. T is a column vectors given as matrix, with at least elements on
+ * index 0 to K-1; elements on index>=K are ignored. Vector R must not have
+ * zero elements, no check is done.
+ */
+void lcl_SolveWithLowerLeftTriangle(ScMatrixRef pMatA,
+ ::std::vector< double>& pVecR, ScMatrixRef pMatT,
+ SCSIZE nK, bool bIsTransposed)
{
- RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::RGetVariances" );
- // pE[0] := Sigma i=1...n (Yi)
- // pE[k] := Sigma i=1...n (Xki*Yi)
- // pE[M+1] := Sigma i=1...n (Yi**2)
- // pQ[0,M+1]:= B
- // pQ[k,M+1]:= Mk
- double fSQR, fSQT, fSQE;
- fSQT = pE->GetDouble(M+1)
- - pE->GetDouble(0) * pE->GetDouble(0) / (double)N;
- fSQR = pE->GetDouble(M+1);
- SCSIZE i, j;
- for (i = 0; i < M+1; i++)
- fSQR -= pQ->GetDouble(i, M+1) * pE->GetDouble(i);
- fSQE = fSQT-fSQR;
- // r2 (Bestimmtheitsmass, 0...1)
- if (fSQT == 0.0)
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), 0, 2);
- else
- pResMat->PutDouble (fSQE/fSQT, 0, 2);
- // ssReg (Regressions-Quadratsumme)
- pResMat->PutDouble(fSQE, 0, 4);
- // ssResid (Residual-Quadratsumme, Summe der Abweichungsquadrate)
- pResMat->PutDouble(fSQR, 1, 4);
- for (i = 2; i < 5; i++)
- for (j = 2; j < M+1; j++)
- pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), j, i);
- if (bConstant)
+ // ScMatrix matrices are zero based, index access (column,row)
+ double fSum;
+ for (SCSIZE row = 0; row < nK; row++)
{
- if (N-M-1 == 0)
+ fSum = pMatT -> GetDouble(row);
+ for (SCSIZE col=0; col < row; col++)
{
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), 1, 2);
- for (i = 0; i < M+1; i++)
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i, 1);
- }
+ if (bIsTransposed)
+ fSum -= pMatA->GetDouble(col,row) * pMatT->GetDouble(col);
else
- {
- double fSE2 = fSQR/(N-M-1);
- // sey (Standardfehler des Schaetzwertes y)
- pResMat->PutDouble(sqrt(fSE2), 1, 2);
- // sen...se1 (Standardfehler der Koeffizienten mn...m1)
- // seb (Standardfehler der Konstanten b)
- if ( RGetVariances( pV, pMatX, M+1, N, nCase != 2, FALSE ) )
- {
- for (i = 0; i < M+1; i++)
- pResMat->PutDouble( sqrt(fSE2 * pV->GetDouble(i)), M-i, 1 );
+ fSum -= pMatA->GetDouble(row,col) * pMatT->GetDouble(col);
}
- else
- {
- for (i = 0; i < M+1; i++)
- pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 1);
+ pMatT->PutDouble( fSum / pVecR[row] , row);
}
}
- // F (F-Statistik)
- if (fSQR == 0.0)
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), 0, 3);
+
+/* Calculates Z = R * B
+ * R is given in matrix A and vector VecR as obtained from the QR
+ * decompostion in lcl_CalculateQRdecomposition. B and Z are column vectors
+ * given as matrix with at least index 0 to K-1; elements on index>=K are
+ * not used.
+ */
+void lcl_ApplyUpperRightTriangle(ScMatrixRef pMatA,
+ ::std::vector< double>& pVecR, ScMatrixRef pMatB,
+ ScMatrixRef pMatZ, SCSIZE nK, bool bIsTransposed)
+{
+ // ScMatrix matrices are zero based, index access (column,row)
+ double fSum;
+ for (SCSIZE row = 0; row < nK; row++)
+ {
+ fSum = pVecR[row] * pMatB->GetDouble(row);
+ for (SCSIZE col = row+1; col < nK; col++)
+ if (bIsTransposed)
+ fSum += pMatA->GetDouble(row,col) * pMatB->GetDouble(col);
else
- pResMat->PutDouble(((double)(N-M-1))*fSQE/fSQR/((double)M),0, 3);
- // df (Freiheitsgrad)
- pResMat->PutDouble(((double)(N-M-1)), 1, 3);
+ fSum += pMatA->GetDouble(col,row) * pMatB->GetDouble(col);
+ pMatZ->PutDouble( fSum, row);
}
- else
- {
- if (N-M == 0)
+}
+
+
+
+double lcl_GetMeanOverAll(ScMatrixRef pMat, SCSIZE nN)
{
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), 1, 2);
- for (i = 0; i < M+1; i++)
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i, 1);
+ double fSum = 0.0;
+ for (SCSIZE i=0 ; i<nN; i++)
+ fSum += pMat->GetDouble(i);
+ return fSum/static_cast<double>(nN);
}
- else
+
+// Calculates means of the columns of matrix X. X is a RxC matrix;
+// ResMat is a 1xC matrix (=row).
+void lcl_CalculateColumnMeans(ScMatrixRef pX, ScMatrixRef pResMat,
+ SCSIZE nC, SCSIZE nR)
{
- double fSE2 = fSQR/(N-M);
- pResMat->PutDouble(sqrt(fSE2), 1, 2);
- if ( RGetVariances( pV, pMatX, M, N, nCase != 2, TRUE ) )
+ double fSum = 0.0;
+ for (SCSIZE i=0; i < nC; i++)
{
- for (i = 0; i < M; i++)
- pResMat->PutDouble( sqrt(fSE2 * pV->GetDouble(i)), M-i-1, 1 );
- pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), M, 1);
+ fSum =0.0;
+ for (SCSIZE k=0; k < nR; k++)
+ fSum += pX->GetDouble(i,k); // GetDouble(Column,Row)
+ pResMat ->PutDouble( fSum/static_cast<double>(nR),i);
}
- else
+}
+
+// Calculates means of the rows of matrix X. X is a RxC matrix;
+// ResMat is a Rx1 matrix (=column).
+void lcl_CalculateRowMeans(ScMatrixRef pX, ScMatrixRef pResMat,
+ SCSIZE nC, SCSIZE nR)
{
- for (i = 0; i < M+1; i++)
- pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 1);
+ double fSum = 0.0;
+ for (SCSIZE k=0; k < nR; k++)
+ {
+ fSum =0.0;
+ for (SCSIZE i=0; i < nC; i++)
+ fSum += pX->GetDouble(i,k); // GetDouble(Column,Row)
+ pResMat ->PutDouble( fSum/static_cast<double>(nC),k);
}
}
- if (fSQR == 0.0)
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), 0, 3);
- else
- pResMat->PutDouble(((double)(N-M))*fSQE/fSQR/((double)M),0, 3);
- pResMat->PutDouble(((double)(N-M)), 1, 3);
+
+void lcl_CalculateColumnsDelta(ScMatrixRef pMat, ScMatrixRef pColumnMeans,
+ SCSIZE nC, SCSIZE nR)
+{
+ for (SCSIZE i = 0; i < nC; i++)
+ for (SCSIZE k = 0; k < nR; k++)
+ pMat->PutDouble( ::rtl::math::approxSub
+ (pMat->GetDouble(i,k) , pColumnMeans->GetDouble(i) ) , i, k);
}
+
+void lcl_CalculateRowsDelta(ScMatrixRef pMat, ScMatrixRef pRowMeans,
+ SCSIZE nC, SCSIZE nR)
+{
+ for (SCSIZE k = 0; k < nR; k++)
+ for (SCSIZE i = 0; i < nC; i++)
+ pMat->PutDouble( ::rtl::math::approxSub
+ ( pMat->GetDouble(i,k) , pRowMeans->GetDouble(k) ) , i, k);
}
-void ScInterpreter::ScRGP()
+// Case1 = simple regression
+// MatX = X - MeanX, MatY = Y - MeanY, y - haty = (y - MeanY) - (haty - MeanY)
+// = (y-MeanY)-((slope*x+a)-(slope*MeanX+a)) = (y-MeanY)-slope*(x-MeanX)
+double lcl_GetSSresid(ScMatrixRef pMatX, ScMatrixRef pMatY, double fSlope,
+ SCSIZE nN)
{
- RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScRGP" );
- CalulateRGPRKP(FALSE);
+ double fSum = 0.0;
+ double fTemp = 0.0;
+ for (SCSIZE i=0; i<nN; i++)
+ {
+ fTemp = pMatY->GetDouble(i) - fSlope * pMatX->GetDouble(i);
+ fSum += fTemp * fTemp;
+ }
+ return fSum;
}
-bool ScInterpreter::CheckMatrix(BOOL _bLOG,BOOL _bTrendGrowth,BYTE& nCase,SCSIZE& nCX,SCSIZE& nCY,SCSIZE& nRX,SCSIZE& nRY,SCSIZE& M,SCSIZE& N,ScMatrixRef& pMatX,ScMatrixRef& pMatY)
+
+// Fill default values in matrix X, transform Y to log(Y) in case LOGEST|GROWTH,
+// determine sizes of matrices X and Y, determine kind of regression, clone
+// Y in case LOGEST|GROWTH, if constant.
+bool ScInterpreter::CheckMatrix(bool _bLOG, BYTE& nCase, SCSIZE& nCX,
+ SCSIZE& nCY, SCSIZE& nRX, SCSIZE& nRY, SCSIZE& M,
+ SCSIZE& N, ScMatrixRef& pMatX, ScMatrixRef& pMatY)
{
+
nCX = 0;
nCY = 0;
nRX = 0;
@@ -2107,7 +2282,11 @@ bool ScInterpreter::CheckMatrix(BOOL _bLOG,BOOL _bTrendGrowth,BYTE& nCase,SCSIZE
return false;
}
if (nCX == nCY && nRX == nRY)
- nCase = 1; // einfache Regression
+ {
+ nCase = 1; // simple regression
+ M = 1;
+ N = nCountY;
+ }
else if (nCY != 1 && nRY != 1)
{
PushIllegalArgument();
@@ -2122,7 +2301,7 @@ bool ScInterpreter::CheckMatrix(BOOL _bLOG,BOOL _bTrendGrowth,BYTE& nCase,SCSIZE
}
else
{
- nCase = 2; // zeilenweise
+ nCase = 2; // Y is column
N = nRY;
M = nCX;
}
@@ -2134,7 +2313,7 @@ bool ScInterpreter::CheckMatrix(BOOL _bLOG,BOOL _bTrendGrowth,BYTE& nCase,SCSIZE
}
else
{
- nCase = 3; // spaltenweise
+ nCase = 3; // Y is row
N = nCY;
M = nRX;
}
@@ -2142,460 +2321,644 @@ bool ScInterpreter::CheckMatrix(BOOL _bLOG,BOOL _bTrendGrowth,BYTE& nCase,SCSIZE
else
{
pMatX = GetNewMat(nCY, nRY);
- if ( _bTrendGrowth )
- {
nCX = nCY;
nRX = nRY;
- }
if (!pMatX)
{
PushIllegalArgument();
return false;
}
for ( SCSIZE i = 1; i <= nCountY; i++ )
- pMatX->PutDouble((double)i, i-1);
+ pMatX->PutDouble(static_cast<double>(i), i-1);
nCase = 1;
+ N = nCountY;
+ M = 1;
}
return true;
}
-void ScInterpreter::CalulateRGPRKP(BOOL _bRKP)
+// -----------------------------------------------------------------------------
+
+// LINEST
+void ScInterpreter::ScRGP()
+{
+ RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScRGP" );
+ CalulateRGPRKP(false);
+}
+
+// LOGEST
+void ScInterpreter::ScRKP()
+{
+ RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScRKP" );
+ CalulateRGPRKP(true);
+}
+
+void ScInterpreter::CalulateRGPRKP(bool _bRKP)
{
- RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CheckMatrix" );
BYTE nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 1, 4 ) )
return;
- BOOL bConstant, bStats;
+ bool bConstant, bStats;
+
+ // optional forth parameter
if (nParamCount == 4)
bStats = GetBool();
else
- bStats = FALSE;
+ bStats = false;
+
+ // The third parameter may not be missing in ODF, if the forth parameter
+ // is present. But Excel allows it with default true, we too.
if (nParamCount >= 3)
+ {
+ if (IsMissing())
+ {
+ Pop();
+ bConstant = true;
+// PushIllegalParameter(); if ODF behavior is desired
+// return;
+ }
+ else
bConstant = GetBool();
+ }
else
- bConstant = TRUE;
+ bConstant = true;
+
ScMatrixRef pMatX;
- ScMatrixRef pMatY;
if (nParamCount >= 2)
+ {
+ if (IsMissing())
+ { //In ODF1.2 empty second parameter (which is two ;; ) is allowed
+ Pop();
+ pMatX = NULL;
+ }
+ else
+ {
pMatX = GetMatrix();
+ }
+ }
else
pMatX = NULL;
+
+ ScMatrixRef pMatY;
pMatY = GetMatrix();
if (!pMatY)
{
PushIllegalParameter();
return;
- } // if (!pMatY)
- BYTE nCase; // 1 = normal, 2,3 = mehrfach
- SCSIZE nCX, nCY;
- SCSIZE nRX, nRY;
- SCSIZE M = 0, N = 0;
- if ( !CheckMatrix(_bRKP,FALSE,nCase,nCX,nCY,nRX,nRY,M,N,pMatX,pMatY) )
- return;
+ }
- ScMatrixRef pResMat;
- if (nCase == 1)
- {
- if (!bStats)
- pResMat = GetNewMat(2,1);
- else
- pResMat = GetNewMat(2,5);
- if (!pResMat)
+ // 1 = simple; 2 = multiple with Y as column; 3 = multiple with Y as row
+ BYTE nCase;
+
+ SCSIZE nCX, nCY; // number of columns
+ SCSIZE nRX, nRY; //number of rows
+ SCSIZE K = 0, N = 0; // K=number of variables X, N=number of data samples
+ if ( !CheckMatrix(_bRKP,nCase,nCX,nCY,nRX,nRY,K,N,pMatX,pMatY) )
{
- PushIllegalArgument();
+ PushIllegalParameter();
return;
}
- double fCount = 0.0;
- double fSumX = 0.0;
- double fSumSqrX = 0.0;
- double fSumY = 0.0;
- double fSumSqrY = 0.0;
- double fSumXY = 0.0;
- double fValX, fValY;
- for (SCSIZE i = 0; i < nCY; i++)
- for (SCSIZE j = 0; j < nRY; j++)
- {
- fValX = pMatX->GetDouble(i,j);
- fValY = pMatY->GetDouble(i,j);
- fSumX += fValX;
- fSumSqrX += fValX * fValX;
- fSumY += fValY;
- fSumSqrY += fValY * fValY;
- fSumXY += fValX*fValY;
- fCount++;
- }
- if (fCount < 1.0)
- PushNoValue();
- else
- {
- double f1 = fCount*fSumXY-fSumX*fSumY;
- double fX = fCount*fSumSqrX-fSumX*fSumX;
- double b, m;
- if (bConstant)
+
+ // Enough data samples?
+ if ( (bConstant && (N<K+1)) || (!bConstant && (N<K)) || (N<1) || (K<1) )
{
- b = fSumY/fCount - f1/fX*fSumX/fCount;
- m = f1/fX;
+ PushIllegalParameter();
+ return;
}
+
+ ScMatrixRef pResMat;
+ if (bStats)
+ pResMat = GetNewMat(K+1,5);
else
+ pResMat = GetNewMat(K+1,1);
+ if (!pResMat)
{
- b = 0.0;
- m = fSumXY/fSumSqrX;
+ PushError(errCodeOverflow);
+ return;
}
- pResMat->PutDouble(_bRKP ? exp(m) : m, 0, 0);
- pResMat->PutDouble(_bRKP ? exp(b) : b, 1, 0);
+ // Fill unused cells in pResMat; order (column,row)
if (bStats)
{
- double fY = fCount*fSumSqrY-fSumY*fSumY;
- double fSyx = fSumSqrY-b*fSumY-m*fSumXY;
- double fR2 = f1*f1/(fX*fY);
- pResMat->PutDouble (fR2, 0, 2);
- if (fCount < 3.0)
- {
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), 0, 1 );
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), 1, 1 );
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), 1, 2 );
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), 0, 3 );
- }
- else
+ for (SCSIZE i=2; i<K+1; i++)
{
- pResMat->PutDouble(sqrt(fSyx*fCount/(fX*(fCount-2.0))), 0, 1);
- pResMat->PutDouble(sqrt(fSyx*fSumSqrX/fX/(fCount-2.0)), 1, 1);
- pResMat->PutDouble(
- sqrt((fCount*fSumSqrY - fSumY*fSumY - f1*f1/fX)/
- (fCount*(fCount-2.0))), 1, 2);
- if (fR2 == 1.0)
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), 0, 3 );
- else
- pResMat->PutDouble(fR2*(fCount-2.0)/(1.0-fR2), 0, 3);
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 2 );
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 3 );
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 4 );
}
- pResMat->PutDouble(((double)(nCY*nRY))-2.0, 1, 3);
- pResMat->PutDouble(fY/fCount-fSyx, 0, 4);
- pResMat->PutDouble(fSyx, 1, 4);
}
- }
- } // if (nCase == 1)
- if ( nCase != 1 )
+
+ // Uses sum(x-MeanX)^2 and not [sum x^2]-N * MeanX^2 in case bConstant.
+ // Clone constant matrices, so that Mat = Mat - Mean is possible.
+ double fMeanY = 0.0;
+ if (bConstant)
{
- SCSIZE i, j, k;
- if (!bStats)
- pResMat = GetNewMat(M+1,1);
- else
- pResMat = GetNewMat(M+1,5);
- if (!pResMat)
+ ScMatrixRef pNewX = pMatX->CloneIfConst();
+ ScMatrixRef pNewY = pMatY->CloneIfConst();
+ if (!pNewX || !pNewY)
{
- PushIllegalArgument();
+ PushError(errCodeOverflow);
return;
}
- ScMatrixRef pQ = GetNewMat(M+1, M+2);
- ScMatrixRef pE = GetNewMat(M+2, 1);
- ScMatrixRef pV = GetNewMat(M+1, 1);
- pE->PutDouble(0.0, M+1);
- pQ->FillDouble(0.0, 0, 0, M, M+1);
- if (nCase == 2)
- {
- for (k = 0; k < N; k++)
- {
- double Yk = pMatY->GetDouble(k);
- pE->PutDouble( pE->GetDouble(M+1)+Yk*Yk, M+1 );
- double sumYk = pQ->GetDouble(0, M+1) + Yk;
- pQ->PutDouble( sumYk, 0, M+1 );
- pE->PutDouble( sumYk, 0 );
- for (i = 0; i < M; i++)
+ pMatX = pNewX;
+ pMatY = pNewY;
+ // DeltaY is possible here; DeltaX depends on nCase, so later
+ fMeanY = lcl_GetMeanOverAll(pMatY, N);
+ for (SCSIZE i=0; i<N; i++)
+ {
+ pMatY->PutDouble( ::rtl::math::approxSub(pMatY->GetDouble(i),fMeanY), i );
+ }
+ }
+
+ if (nCase==1)
{
- double Xik = pMatX->GetDouble(i,k);
- double sumXik = pQ->GetDouble(0, i+1) + Xik;
- pQ->PutDouble( sumXik, 0, i+1);
- pQ->PutDouble( sumXik, i+1, 0);
- double sumXikYk = pQ->GetDouble(i+1, M+1) + Xik * Yk;
- pQ->PutDouble( sumXikYk, i+1, M+1);
- pE->PutDouble( sumXikYk, i+1);
- for (j = i; j < M; j++)
+ // calculate simple regression
+ double fMeanX = 0.0;
+ if (bConstant)
+ { // Mat = Mat - Mean
+ fMeanX = lcl_GetMeanOverAll(pMatX, N);
+ for (SCSIZE i=0; i<N; i++)
{
- const double fVal = pMatX->GetDouble(j,k);
- double sumXikXjk = pQ->GetDouble(j+1, i+1) +
- Xik * fVal;
- pQ->PutDouble( sumXikXjk, j+1, i+1);
- pQ->PutDouble( sumXikXjk, i+1, j+1);
+ pMatX->PutDouble( ::rtl::math::approxSub(pMatX->GetDouble(i),fMeanX), i );
}
}
+ double fSumXY = lcl_GetSumProduct(pMatX,pMatY,N);
+ double fSumX2 = lcl_GetSumProduct(pMatX,pMatX,N);
+ if (fSumX2==0.0)
+ {
+ PushNoValue(); // all x-values are identical
+ return;
}
+ double fSlope = fSumXY / fSumX2;
+ double fIntercept = 0.0;
+ if (bConstant)
+ fIntercept = fMeanY - fSlope * fMeanX;
+ pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, 1, 0); //order (column,row)
+ pResMat->PutDouble(_bRKP ? exp(fSlope) : fSlope, 0, 0);
+
+ if (bStats)
+ {
+ double fSSreg = fSlope * fSlope * fSumX2;
+ pResMat->PutDouble(fSSreg, 0, 4);
+
+ double fDegreesFreedom =static_cast<double>( (bConstant) ? N-2 : N-1 );
+ pResMat->PutDouble(fDegreesFreedom, 1, 3);
+
+ double fSSresid = lcl_GetSSresid(pMatX,pMatY,fSlope,N);
+ pResMat->PutDouble(fSSresid, 1, 4);
+
+ if (fDegreesFreedom == 0.0 || fSSresid == 0.0 || fSSreg == 0.0)
+ { // exact fit; test SSreg too, because SSresid might be
+ // unequal zero due to round of errors
+ pResMat->PutDouble(0.0, 1, 4); // SSresid
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 0, 3); // F
+ pResMat->PutDouble(0.0, 1, 2); // RMSE
+ pResMat->PutDouble(0.0, 0, 1); // SigmaSlope
+ if (bConstant)
+ pResMat->PutDouble(0.0, 1, 1); //SigmaIntercept
+ else
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 1, 1);
+ pResMat->PutDouble(1.0, 0, 2); // R^2
}
else
{
- for (k = 0; k < N; k++)
- {
- double Yk = pMatY->GetDouble(k);
- pE->PutDouble( pE->GetDouble(M+1)+Yk*Yk, M+1 );
- double sumYk = pQ->GetDouble(0, M+1) + Yk;
- pQ->PutDouble( sumYk, 0, M+1 );
- pE->PutDouble( sumYk, 0 );
- for (i = 0; i < M; i++)
- {
- double Xki = pMatX->GetDouble(k,i);
- double sumXki = pQ->GetDouble(0, i+1) + Xki;
- pQ->PutDouble( sumXki, 0, i+1);
- pQ->PutDouble( sumXki, i+1, 0);
- double sumXkiYk = pQ->GetDouble(i+1, M+1) + Xki * Yk;
- pQ->PutDouble( sumXkiYk, i+1, M+1);
- pE->PutDouble( sumXkiYk, i+1);
- for (j = i; j < M; j++)
+ double fFstatistic = (fSSreg / static_cast<double>(K))
+ / (fSSresid / fDegreesFreedom);
+ pResMat->PutDouble(fFstatistic, 0, 3);
+
+ // standard error of estimate
+ double fRMSE = sqrt(fSSresid / fDegreesFreedom);
+ pResMat->PutDouble(fRMSE, 1, 2);
+
+ double fSigmaSlope = fRMSE / sqrt(fSumX2);
+ pResMat->PutDouble(fSigmaSlope, 0, 1);
+
+ if (bConstant)
{
- const double fVal = pMatX->GetDouble(k,j);
- double sumXkiXkj = pQ->GetDouble(j+1, i+1) +
- Xki * fVal;
- pQ->PutDouble( sumXkiXkj, j+1, i+1);
- pQ->PutDouble( sumXkiXkj, i+1, j+1);
- }
- }
+ double fSigmaIntercept = fRMSE
+ * sqrt(fMeanX*fMeanX/fSumX2 + 1.0/static_cast<double>(N));
+ pResMat->PutDouble(fSigmaIntercept, 1, 1);
}
+ else
+ {
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 1, 1);
}
- if ( !Calculate4(_bRKP,pResMat,pQ,bConstant,N,M) )
- return;
- if (bStats)
- Calculate(pResMat,pE,pQ,pV,pMatX,bConstant,N,M,nCase);
+ double fR2 = fSSreg / (fSSreg + fSSresid);
+ pResMat->PutDouble(fR2, 0, 2);
+ }
}
PushMatrix(pResMat);
}
-
-void ScInterpreter::ScRKP()
+ else // calculate multiple regression;
{
- RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScRKP" );
- CalulateRGPRKP(TRUE);
-}
-// -----------------------------------------------------------------------------
-bool ScInterpreter::Calculate4(BOOL _bExp,ScMatrixRef& pResMat,ScMatrixRef& pQ,BOOL bConstant,SCSIZE N,SCSIZE M)
+ // Uses a QR decomposition X = QR. The solution B = (X'X)^(-1) * X' * Y
+ // becomes B = R^(-1) * Q' * Y
+ if (nCase ==2) // Y is column
+ {
+ ::std::vector< double> aVecR(N); // for QR decomposition
+ // Enough memory for needed matrices?
+ ScMatrixRef pMeans = GetNewMat(K, 1); // mean of each column
+ ScMatrixRef pMatZ; // for Q' * Y , inter alia
+ if (bStats)
+ pMatZ = pMatY->Clone(); // Y is used in statistic, keep it
+ else
+ pMatZ = pMatY; // Y can be overwritten
+ ScMatrixRef pSlopes = GetNewMat(1,K); // from b1 to bK
+ if (!pMeans || !pMatZ || !pSlopes)
{
- RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::Calculate4" );
- pQ->PutDouble((double)N, 0, 0);
+ PushError(errCodeOverflow);
+ return;
+ }
if (bConstant)
{
- SCSIZE S, L;
- for (S = 0; S < M+1; S++)
- {
- SCSIZE i = S;
- while (i < M+1 && pQ->GetDouble(i, S) == 0.0)
- i++;
- if (i >= M+1)
- {
- PushNoValue();
- return false;
+ lcl_CalculateColumnMeans(pMatX, pMeans, K, N);
+ lcl_CalculateColumnsDelta(pMatX, pMeans, K, N);
}
- double fVal;
- for (L = 0; L < M+2; L++)
+ if (!lcl_CalculateQRdecomposition(pMatX, aVecR, K, N))
{
- fVal = pQ->GetDouble(S, L);
- pQ->PutDouble(pQ->GetDouble(i, L), S, L);
- pQ->PutDouble(fVal, i, L);
+ PushNoValue();
+ return;
}
- fVal = 1.0/pQ->GetDouble(S, S);
- for (L = 0; L < M+2; L++)
- pQ->PutDouble(pQ->GetDouble(S, L)*fVal, S, L);
- for (i = 0; i < M+1; i++)
- {
- if (i != S)
+ // Later on we will divide by elements of aVecR, so make sure
+ // that they aren't zero.
+ bool bIsSingular=false;
+ for (SCSIZE row=0; row < K && !bIsSingular; row++)
+ bIsSingular = bIsSingular || aVecR[row]==0.0;
+ if (bIsSingular)
{
- fVal = -pQ->GetDouble(i, S);
- for (L = 0; L < M+2; L++)
- pQ->PutDouble(
- pQ->GetDouble(i,L)+fVal*pQ->GetDouble(S,L),i,L);
- }
- }
- }
+ PushNoValue();
+ return;
}
- else
+ // Z = Q' Y;
+ for (SCSIZE col = 0; col < K; col++)
{
- if ( !Calculate3(M,pQ) )
- return false;
-
+ lcl_ApplyHouseholderTransformation(pMatX, col, pMatZ, N);
}
- for (SCSIZE i = 0; i < M+1; i++)
+ // B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z
+ // result Z should have zeros for index>=K; if not, ignore values
+ for (SCSIZE col = 0; col < K ; col++)
{
- const double d = pQ->GetDouble(M-i,M+1);
- pResMat->PutDouble(_bExp ? exp(d) : d, i, 0);
- } // for (SCSIZE i = 0; i < M+1; i++)
- return true;
+ pSlopes->PutDouble( pMatZ->GetDouble(col), col);
}
+ lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, false);
+ double fIntercept = 0.0;
+ if (bConstant)
+ fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K);
+ // Fill first line in result matrix
+ pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, K, 0 );
+ for (SCSIZE i = 0; i < K; i++)
+ pResMat->PutDouble(_bRKP ? exp(pSlopes->GetDouble(i))
+ : pSlopes->GetDouble(i) , K-1-i, 0);
-ScMatrixRef ScInterpreter::Calculate2(const BOOL bConstant,const SCSIZE M ,const SCSIZE N,ScMatrixRef& pMatX,ScMatrixRef& pMatY,BYTE nCase)
-{
- RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::Calculate2" );
- SCSIZE i, j, k;
- ScMatrixRef pQ = GetNewMat(M+1, M+2);
- ScMatrixRef pE = GetNewMat(M+2, 1);
- pE->PutDouble(0.0, M+1);
- pQ->FillDouble(0.0, 0, 0, M, M+1);
- if (nCase == 2)
- {
- for (k = 0; k < N; k++)
- {
- pE->PutDouble(
- pE->GetDouble(M+1)+pMatY->GetDouble(k)*pMatY->GetDouble(k), M+1);
- pQ->PutDouble(pQ->GetDouble(0, M+1) + pMatY->GetDouble(k), 0, M+1);
- pE->PutDouble(pQ->GetDouble(0, M+1), 0);
- for (i = 0; i < M; i++)
- {
- pQ->PutDouble(pQ->GetDouble(0, i+1)+pMatX->GetDouble(i,k), 0, i+1);
- pQ->PutDouble(pQ->GetDouble(0, i+1), i+1, 0);
- pQ->PutDouble(pQ->GetDouble(i+1, M+1) +
- pMatX->GetDouble(i,k)*pMatY->GetDouble(k), i+1, M+1);
- pE->PutDouble(pQ->GetDouble(i+1, M+1), i+1);
- for (j = i; j < M; j++)
+
+ if (bStats)
+ {
+ double fSSreg = 0.0;
+ double fSSresid = 0.0;
+ // re-use memory of Z;
+ pMatZ->FillDouble(0.0, 0, 0, 0, N-1);
+ // Z = R * Slopes
+ lcl_ApplyUpperRightTriangle(pMatX, aVecR, pSlopes, pMatZ, K, false);
+ // Z = Q * Z, that is Q * R * Slopes = X * Slopes
+ for (SCSIZE colp1 = K; colp1 > 0; colp1--)
{
- pQ->PutDouble(pQ->GetDouble(j+1, i+1) +
- pMatX->GetDouble(i,k)*pMatX->GetDouble(j,k), j+1, i+1);
- pQ->PutDouble(pQ->GetDouble(j+1, i+1), i+1, j+1);
+ lcl_ApplyHouseholderTransformation(pMatX, colp1-1, pMatZ,N);
+ }
+ fSSreg =lcl_GetSumProduct(pMatZ, pMatZ, N);
+ // re-use Y for residuals, Y = Y-Z
+ for (SCSIZE row = 0; row < N; row++)
+ pMatY->PutDouble(pMatY->GetDouble(row) - pMatZ->GetDouble(row), row);
+ fSSresid = lcl_GetSumProduct(pMatY, pMatY, N);
+ pResMat->PutDouble(fSSreg, 0, 4);
+ pResMat->PutDouble(fSSresid, 1, 4);
+
+ double fDegreesFreedom =static_cast<double>( (bConstant) ? N-K-1 : N-K );
+ pResMat->PutDouble(fDegreesFreedom, 1, 3);
+
+ if (fDegreesFreedom == 0.0 || fSSresid == 0.0 || fSSreg == 0.0)
+ { // exact fit; incl. observed values Y are identical
+ pResMat->PutDouble(0.0, 1, 4); // SSresid
+ // F = (SSreg/K) / (SSresid/df) = #DIV/0!
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 0, 3); // F
+ // RMSE = sqrt(SSresid / df) = sqrt(0 / df) = 0
+ pResMat->PutDouble(0.0, 1, 2); // RMSE
+ // SigmaSlope[i] = RMSE * sqrt(matrix[i,i]) = 0 * sqrt(...) = 0
+ for (SCSIZE i=0; i<K; i++)
+ pResMat->PutDouble(0.0, K-1-i, 1);
+
+ // SigmaIntercept = RMSE * sqrt(...) = 0
+ if (bConstant)
+ pResMat->PutDouble(0.0, K, 1); //SigmaIntercept
+ else
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1);
+
+ // R^2 = SSreg / (SSreg + SSresid) = 1.0
+ pResMat->PutDouble(1.0, 0, 2); // R^2
}
+ else
+ {
+ double fFstatistic = (fSSreg / static_cast<double>(K))
+ / (fSSresid / fDegreesFreedom);
+ pResMat->PutDouble(fFstatistic, 0, 3);
+
+ // standard error of estimate = root mean SSE
+ double fRMSE = sqrt(fSSresid / fDegreesFreedom);
+ pResMat->PutDouble(fRMSE, 1, 2);
+
+ // standard error of slopes
+ // = RMSE * sqrt(diagonal element of (R' R)^(-1) )
+ // standard error of intercept
+ // = RMSE * sqrt( Xmean * (R' R)^(-1) * Xmean' + 1/N)
+ // (R' R)^(-1) = R^(-1) * (R')^(-1). Do not calculate it as
+ // a whole matrix, but iterate over unit vectors.
+ double fSigmaSlope = 0.0;
+ double fSigmaIntercept = 0.0;
+ double fPart; // for Xmean * single column of (R' R)^(-1)
+ for (SCSIZE col = 0; col < K; col++)
+ {
+ //re-use memory of MatZ
+ pMatZ->FillDouble(0.0,0,0,0,K-1); // Z = unit vector e
+ pMatZ->PutDouble(1.0, col);
+ //Solve R' * Z = e
+ lcl_SolveWithLowerLeftTriangle(pMatX, aVecR, pMatZ, K, false);
+ // Solve R * Znew = Zold
+ lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pMatZ, K, false);
+ // now Z is column col in (R' R)^(-1)
+ fSigmaSlope = fRMSE * sqrt(pMatZ->GetDouble(col));
+ pResMat->PutDouble(fSigmaSlope, K-1-col, 1);
+ // (R' R) ^(-1) is symmetric
+ if (bConstant)
+ {
+ fPart = lcl_GetSumProduct(pMeans, pMatZ, K);
+ fSigmaIntercept += fPart * pMeans->GetDouble(col);
}
}
+ if (bConstant)
+ {
+ fSigmaIntercept = fRMSE
+ * sqrt(fSigmaIntercept + 1.0 / static_cast<double>(N));
+ pResMat->PutDouble(fSigmaIntercept, K, 1);
}
else
{
- for (k = 0; k < N; k++)
- {
- pE->PutDouble(
- pE->GetDouble(M+1)+pMatY->GetDouble(k)*pMatY->GetDouble(k), M+1);
- pQ->PutDouble(pQ->GetDouble(0, M+1) + pMatY->GetDouble(k), 0, M+1);
- pE->PutDouble(pQ->GetDouble(0, M+1), 0);
- for (i = 0; i < M; i++)
- {
- pQ->PutDouble(pQ->GetDouble(0, i+1)+pMatX->GetDouble(k,i), 0, i+1);
- pQ->PutDouble(pQ->GetDouble(0, i+1), i+1, 0);
- pQ->PutDouble(pQ->GetDouble(i+1, M+1) +
- pMatX->GetDouble(k,i)*pMatY->GetDouble(k), i+1, M+1);
- pE->PutDouble(pQ->GetDouble(i+1, M+1), i+1);
- for (j = i; j < M; j++)
- {
- pQ->PutDouble(pQ->GetDouble(j+1, i+1) +
- pMatX->GetDouble(k, i)*pMatX->GetDouble(k, j), j+1, i+1);
- pQ->PutDouble(pQ->GetDouble(j+1, i+1), i+1, j+1);
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1);
}
+
+ double fR2 = fSSreg / (fSSreg + fSSresid);
+ pResMat->PutDouble(fR2, 0, 2);
}
}
+ PushMatrix(pResMat);
}
- pQ->PutDouble((double)N, 0, 0);
- if (bConstant)
- {
- SCSIZE S, L;
- for (S = 0; S < M+1; S++)
+ else // nCase == 3, Y is row, all matrices are transposed
{
- i = S;
- while (i < M+1 && pQ->GetDouble(i, S) == 0.0)
- i++;
- if (i >= M+1)
+ ::std::vector< double> aVecR(N); // for QR decomposition
+ // Enough memory for needed matrices?
+ ScMatrixRef pMeans = GetNewMat(1, K); // mean of each row
+ ScMatrixRef pMatZ; // for Q' * Y , inter alia
+ if (bStats)
+ pMatZ = pMatY->Clone(); // Y is used in statistic, keep it
+ else
+ pMatZ = pMatY; // Y can be overwritten
+ ScMatrixRef pSlopes = GetNewMat(K,1); // from b1 to bK
+ if (!pMeans || !pMatZ || !pSlopes)
{
- PushNoValue();
- return ScMatrixRef();
+ PushError(errCodeOverflow);
+ return;
}
- double fVal;
- for (L = 0; L < M+2; L++)
+ if (bConstant)
{
- fVal = pQ->GetDouble(S, L);
- pQ->PutDouble(pQ->GetDouble(i, L), S, L);
- pQ->PutDouble(fVal, i, L);
+ lcl_CalculateRowMeans(pMatX, pMeans, N, K);
+ lcl_CalculateRowsDelta(pMatX, pMeans, N, K);
}
- fVal = 1.0/pQ->GetDouble(S, S);
- for (L = 0; L < M+2; L++)
- pQ->PutDouble(pQ->GetDouble(S, L)*fVal, S, L);
- for (i = 0; i < M+1; i++)
+
+ if (!lcl_TCalculateQRdecomposition(pMatX, aVecR, K, N))
{
- if (i != S)
+ PushNoValue();
+ return;
+ }
+
+ // Later on we will divide by elements of aVecR, so make sure
+ // that they aren't zero.
+ bool bIsSingular=false;
+ for (SCSIZE row=0; row < K && !bIsSingular; row++)
+ bIsSingular = bIsSingular || aVecR[row]==0.0;
+ if (bIsSingular)
{
- fVal = -pQ->GetDouble(i, S);
- for (L = 0; L < M+2; L++)
- pQ->PutDouble(
- pQ->GetDouble(i,L)+fVal*pQ->GetDouble(S,L),i,L);
+ PushNoValue();
+ return;
}
+ // Z = Q' Y
+ for (SCSIZE row = 0; row < K; row++)
+ {
+ lcl_TApplyHouseholderTransformation(pMatX, row, pMatZ, N);
}
+ // B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z
+ // result Z should have zeros for index>=K; if not, ignore values
+ for (SCSIZE col = 0; col < K ; col++)
+ {
+ pSlopes->PutDouble( pMatZ->GetDouble(col), col);
}
+ lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, true);
+ double fIntercept = 0.0;
+ if (bConstant)
+ fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K);
+ // Fill first line in result matrix
+ pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, K, 0 );
+ for (SCSIZE i = 0; i < K; i++)
+ pResMat->PutDouble(_bRKP ? exp(pSlopes->GetDouble(i))
+ : pSlopes->GetDouble(i) , K-1-i, 0);
+
+
+ if (bStats)
+ {
+ double fSSreg = 0.0;
+ double fSSresid = 0.0;
+ // re-use memory of Z;
+ pMatZ->FillDouble(0.0, 0, 0, N-1, 0);
+ // Z = R * Slopes
+ lcl_ApplyUpperRightTriangle(pMatX, aVecR, pSlopes, pMatZ, K, true);
+ // Z = Q * Z, that is Q * R * Slopes = X * Slopes
+ for (SCSIZE rowp1 = K; rowp1 > 0; rowp1--)
+ {
+ lcl_TApplyHouseholderTransformation(pMatX, rowp1-1, pMatZ,N);
+ }
+ fSSreg =lcl_GetSumProduct(pMatZ, pMatZ, N);
+ // re-use Y for residuals, Y = Y-Z
+ for (SCSIZE col = 0; col < N; col++)
+ pMatY->PutDouble(pMatY->GetDouble(col) - pMatZ->GetDouble(col), col);
+ fSSresid = lcl_GetSumProduct(pMatY, pMatY, N);
+ pResMat->PutDouble(fSSreg, 0, 4);
+ pResMat->PutDouble(fSSresid, 1, 4);
+
+ double fDegreesFreedom =static_cast<double>( (bConstant) ? N-K-1 : N-K );
+ pResMat->PutDouble(fDegreesFreedom, 1, 3);
+
+ if (fDegreesFreedom == 0.0 || fSSresid == 0.0 || fSSreg == 0.0)
+ { // exact fit; incl. case observed values Y are identical
+ pResMat->PutDouble(0.0, 1, 4); // SSresid
+ // F = (SSreg/K) / (SSresid/df) = #DIV/0!
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 0, 3); // F
+ // RMSE = sqrt(SSresid / df) = sqrt(0 / df) = 0
+ pResMat->PutDouble(0.0, 1, 2); // RMSE
+ // SigmaSlope[i] = RMSE * sqrt(matrix[i,i]) = 0 * sqrt(...) = 0
+ for (SCSIZE i=0; i<K; i++)
+ pResMat->PutDouble(0.0, K-1-i, 1);
+
+ // SigmaIntercept = RMSE * sqrt(...) = 0
+ if (bConstant)
+ pResMat->PutDouble(0.0, K, 1); //SigmaIntercept
+ else
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1);
+
+ // R^2 = SSreg / (SSreg + SSresid) = 1.0
+ pResMat->PutDouble(1.0, 0, 2); // R^2
}
else
{
- if ( !Calculate3(M,pQ) )
- return ScMatrixRef();
+ double fFstatistic = (fSSreg / static_cast<double>(K))
+ / (fSSresid / fDegreesFreedom);
+ pResMat->PutDouble(fFstatistic, 0, 3);
+
+ // standard error of estimate = root mean SSE
+ double fRMSE = sqrt(fSSresid / fDegreesFreedom);
+ pResMat->PutDouble(fRMSE, 1, 2);
+
+ // standard error of slopes
+ // = RMSE * sqrt(diagonal element of (R' R)^(-1) )
+ // standard error of intercept
+ // = RMSE * sqrt( Xmean * (R' R)^(-1) * Xmean' + 1/N)
+ // (R' R)^(-1) = R^(-1) * (R')^(-1). Do not calculate it as
+ // a whole matrix, but iterate over unit vectors.
+ // (R' R) ^(-1) is symmetric
+ double fSigmaSlope = 0.0;
+ double fSigmaIntercept = 0.0;
+ double fPart; // for Xmean * single col of (R' R)^(-1)
+ for (SCSIZE row = 0; row < K; row++)
+ {
+ //re-use memory of MatZ
+ pMatZ->FillDouble(0.0,0,0,K-1,0); // Z = unit vector e
+ pMatZ->PutDouble(1.0, row);
+ //Solve R' * Z = e
+ lcl_SolveWithLowerLeftTriangle(pMatX, aVecR, pMatZ, K, true);
+ // Solve R * Znew = Zold
+ lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pMatZ, K, true);
+ // now Z is column col in (R' R)^(-1)
+ fSigmaSlope = fRMSE * sqrt(pMatZ->GetDouble(row));
+ pResMat->PutDouble(fSigmaSlope, K-1-row, 1);
+ if (bConstant)
+ {
+ fPart = lcl_GetSumProduct(pMeans, pMatZ, K);
+ fSigmaIntercept += fPart * pMeans->GetDouble(row);
}
- return pQ;
}
-bool ScInterpreter::Calculate3(const SCSIZE M ,ScMatrixRef& pQ)
-{
- RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::Calculate3" );
- SCSIZE S, L;
- for (S = 1; S < M+1; S++)
- {
- SCSIZE i = S;
- while (i < M+1 && pQ->GetDouble(i, S) == 0.0)
- i++;
- if (i >= M+1)
+ if (bConstant)
{
- PushNoValue();
- return ScMatrixRef();
+ fSigmaIntercept = fRMSE
+ * sqrt(fSigmaIntercept + 1.0 / static_cast<double>(N));
+ pResMat->PutDouble(fSigmaIntercept, K, 1);
}
- double fVal;
- for (L = 1; L < M+2; L++)
+ else
{
- fVal = pQ->GetDouble(S, L);
- pQ->PutDouble(pQ->GetDouble(i, L), S, L);
- pQ->PutDouble(fVal, i, L);
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1);
}
- fVal = 1.0/pQ->GetDouble(S, S);
- for (L = 1; L < M+2; L++)
- pQ->PutDouble(pQ->GetDouble(S, L)*fVal, S, L);
- for (i = 1; i < M+1; i++)
- {
- if (i != S)
- {
- fVal = -pQ->GetDouble(i, S);
- for (L = 1; L < M+2; L++)
- pQ->PutDouble(
- pQ->GetDouble(i,L)+fVal*pQ->GetDouble(S,L),i,L);
+
+ double fR2 = fSSreg / (fSSreg + fSSresid);
+ pResMat->PutDouble(fR2, 0, 2);
}
}
- pQ->PutDouble(0.0, 0, M+1);
- } // for (S = 1; S < M+1; S++)
- return true;
+ PushMatrix(pResMat);
+ }
+ }
+ return;
}
void ScInterpreter::ScTrend()
{
RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScTrend" );
- CalculateTrendGrowth(FALSE);
+ CalculateTrendGrowth(false);
+}
+
+void ScInterpreter::ScGrowth()
+{
+ RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScGrowth" );
+ CalculateTrendGrowth(true);
}
-void ScInterpreter::CalculateTrendGrowth(BOOL _bGrowth)
+
+void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
{
- RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CalculateTrendGrowth" );
BYTE nParamCount = GetByte();
if ( !MustHaveParamCount( nParamCount, 1, 4 ) )
return;
- BOOL bConstant;
+
+ // optional forth parameter
+ bool bConstant;
if (nParamCount == 4)
bConstant = GetBool();
else
- bConstant = TRUE;
- ScMatrixRef pMatX;
- ScMatrixRef pMatY;
+ bConstant = true;
+
+ // The third parameter may be missing in ODF, although the forth parameter
+ // is present. Default values depend on data not yet read.
ScMatrixRef pMatNewX;
if (nParamCount >= 3)
+ {
+ if (IsMissing())
+ {
+ Pop();
+ pMatNewX = NULL;
+ }
+ else
pMatNewX = GetMatrix();
+ }
else
pMatNewX = NULL;
+
+ //In ODF1.2 empty second parameter (which is two ;; ) is allowed
+ //Defaults will be set in CheckMatrix
+ ScMatrixRef pMatX;
if (nParamCount >= 2)
+ {
+ if (IsMissing())
+ {
+ Pop();
+ pMatX = NULL;
+ }
+ else
+ {
pMatX = GetMatrix();
+ }
+ }
else
pMatX = NULL;
+
+ ScMatrixRef pMatY;
pMatY = GetMatrix();
if (!pMatY)
{
PushIllegalParameter();
return;
- } // if (!pMatY)
+ }
+
+ // 1 = simple; 2 = multiple with Y as column; 3 = multiple with Y as row
+ BYTE nCase;
- BYTE nCase; // 1 = normal, 2,3 = mehrfach
- SCSIZE nCX, nCY;
- SCSIZE nRX, nRY;
- SCSIZE M = 0, N = 0;
- if ( !CheckMatrix(_bGrowth,TRUE,nCase,nCX,nCY,nRX,nRY,M,N,pMatX,pMatY) )
+ SCSIZE nCX, nCY; // number of columns
+ SCSIZE nRX, nRY; //number of rows
+ SCSIZE K = 0, N = 0; // K=number of variables X, N=number of data samples
+ if ( !CheckMatrix(_bGrowth,nCase,nCX,nCY,nRX,nRY,K,N,pMatX,pMatY) )
+ {
+ PushIllegalParameter();
return;
+ }
+ // Enough data samples?
+ if ( (bConstant && (N<K+1)) || (!bConstant && (N<K)) || (N<1) || (K<1) )
+ {
+ PushIllegalParameter();
+ return;
+ }
+ // Set default pMatNewX if necessary
... etc. - the rest is truncated
More information about the Libreoffice-commits
mailing list