Skip to content

Commit 1306738

Browse files
diluculodiluculo
committed
Add ParameterStatistics class for computing parameter inference statistics
1 parent 9e0f585 commit 1306738

File tree

3 files changed

+706
-83
lines changed

3 files changed

+706
-83
lines changed

‎src/Numerics.Tests/StatisticsTests/ParameterStatisticsTests.cs‎

Lines changed: 240 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,240 @@
1+
// <copyright file="ParameterStatisticsTests.cs" company="Math.NET">
2+
// Math.NET Numerics, part of the Math.NET Project
3+
// https://numerics.mathdotnet.com
4+
// https://github.com/mathnet/mathnet-numerics
5+
//
6+
// Copyright (c) 2009-$CURRENT_YEAR$ Math.NET
7+
//
8+
// Permission is hereby granted, free of charge, to any person
9+
// obtaining a copy of this software and associated documentation
10+
// files (the "Software"), to deal in the Software without
11+
// restriction, including without limitation the rights to use,
12+
// copy, modify, merge, publish, distribute, sublicense, and/or sell
13+
// copies of the Software, and to permit persons to whom the
14+
// Software is furnished to do so, subject to the following
15+
// conditions:
16+
//
17+
// The above copyright notice and this permission notice shall be
18+
// included in all copies or substantial portions of the Software.
19+
//
20+
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
21+
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
22+
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
23+
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
24+
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
25+
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26+
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
27+
// OTHER DEALINGS IN THE SOFTWARE.
28+
// </copyright>
29+
30+
using MathNet.Numerics.LinearAlgebra;
31+
using MathNet.Numerics.Statistics;
32+
using NUnit.Framework;
33+
using System;
34+
using System.Linq;
35+
36+
namespace MathNet.Numerics.Tests.StatisticsTests
37+
{
38+
[TestFixture]
39+
public class ParameterStatisticsTests
40+
{
41+
#region Polynomial Regression Tests
42+
43+
[Test]
44+
public void PolynomialRegressionTest()
45+
{
46+
// https://github.com/mathnet/mathnet-numerics/discussions/801
47+
48+
// Y = B0 + B1*X + B2*X^2
49+
// Parameter Value Error t-value Pr(>|t|) LCL UCL CI half_width
50+
// --------------------------------------------------------------------------------------------
51+
// B0 -0.24 3.07019 -0.07817 0.94481 -13.44995 12.96995 13.20995
52+
// B1 3.46286 2.33969 1.48005 0.27700 -6.60401 13.52972 10.06686
53+
// B2 2.64286 0.38258 6.90799 0.02032 0.99675 4.28897 1.64611
54+
// --------------------------------------------------------------------------------------------
55+
//
56+
// Fit statistics
57+
// -----------------------------------------
58+
// Degree of freedom 2
59+
// Reduced Chi-Sqr 2.04914
60+
// Residual Sum of Sqaures 4.09829
61+
// R Value 0.99947
62+
// R-Square(COD) 0.99893
63+
// Adj. R-Square 0.99786
64+
// Root-MSE(SD) 1.43148
65+
// -----------------------------------------
66+
67+
double[] x = { 1, 2, 3, 4, 5 };
68+
double[] y = { 6.2, 16.9, 33, 57.5, 82.5 };
69+
var order = 2;
70+
71+
var Ns = x.Length;
72+
var k = order + 1; // number of parameters
73+
var dof = Ns - k; // degree of freedom
74+
75+
// Create the [Ns X k] design matrix
76+
// This matrix transforms the polynomial regression problem into a linear system
77+
// Each row represents one data point, and columns represent polynomial terms:
78+
// - First column: constant term (x^0 = 1)
79+
// - Second column: linear term (x^1)
80+
// - Third column: quadratic term (x^2)
81+
// The matrix looks like:
82+
// [ 1 x1 x1^2 ]
83+
// [ 1 x2 x2^2 ]
84+
// [ ... ]
85+
// [ 1 xN xN^2 ]
86+
var X = Matrix<double>.Build.Dense(Ns, k, (i, j) => Math.Pow(x[i], j));
87+
88+
// Create the Y vector
89+
var Y = Vector<double>.Build.DenseOfArray(y);
90+
91+
// Calculate best-fitted parameters using normal equations
92+
var XtX = X.TransposeThisAndMultiply(X);
93+
var XtXInv = XtX.Inverse();
94+
var Xty = X.TransposeThisAndMultiply(Y);
95+
var parameters = XtXInv.Multiply(Xty);
96+
97+
// Calculate the residuals
98+
var residuals = X.Multiply(parameters) - Y;
99+
100+
// Calculate residual variance (RSS/dof)
101+
var RSS = residuals.DotProduct(residuals);
102+
var residualVariance = RSS / dof;
103+
104+
var covariance = ParameterStatistics.CovarianceMatrixForLinearRegression(X, residualVariance);
105+
var standardErrors = ParameterStatistics.StandardErrors(covariance);
106+
var tStatistics = ParameterStatistics.TStatistics(parameters, standardErrors);
107+
var pValues = ParameterStatistics.PValues(tStatistics, dof);
108+
var confIntervals = ParameterStatistics.ConfidenceIntervalHalfWidths(standardErrors, dof, 0.95);
109+
110+
// Calculate total sum of squares for R-squared
111+
var yMean = Y.Average();
112+
var TSS = Y.Select(y_i => Math.Pow(y_i - yMean, 2)).Sum();
113+
var rSquared = 1.0 - RSS / TSS;
114+
var adjustedRSquared = 1 - (1 - rSquared) * (Ns - 1) / dof;
115+
var rootMSE = Math.Sqrt(residualVariance);
116+
117+
// Check parameters
118+
Assert.That(parameters[0], Is.EqualTo(-0.24).Within(0.001));
119+
Assert.That(parameters[1], Is.EqualTo(3.46286).Within(0.001));
120+
Assert.That(parameters[2], Is.EqualTo(2.64286).Within(0.001));
121+
122+
// Check standard errors
123+
Assert.That(standardErrors[0], Is.EqualTo(3.07019).Within(0.001));
124+
Assert.That(standardErrors[1], Is.EqualTo(2.33969).Within(0.001));
125+
Assert.That(standardErrors[2], Is.EqualTo(0.38258).Within(0.001));
126+
127+
// Check t-statistics
128+
Assert.That(tStatistics[0], Is.EqualTo(-0.07817).Within(0.001));
129+
Assert.That(tStatistics[1], Is.EqualTo(1.48005).Within(0.001));
130+
Assert.That(tStatistics[2], Is.EqualTo(6.90799).Within(0.001));
131+
132+
// Check p-values
133+
Assert.That(pValues[0], Is.EqualTo(0.94481).Within(0.001));
134+
Assert.That(pValues[1], Is.EqualTo(0.27700).Within(0.001));
135+
Assert.That(pValues[2], Is.EqualTo(0.02032).Within(0.001));
136+
137+
// Check confidence intervals
138+
Assert.That(confIntervals[0], Is.EqualTo(13.20995).Within(0.001));
139+
Assert.That(confIntervals[1], Is.EqualTo(10.06686).Within(0.001));
140+
Assert.That(confIntervals[2], Is.EqualTo(1.64611).Within(0.001));
141+
142+
// Check fit statistics
143+
Assert.That(dof, Is.EqualTo(2));
144+
Assert.That(residualVariance, Is.EqualTo(2.04914).Within(0.001));
145+
Assert.That(RSS, Is.EqualTo(4.09829).Within(0.001));
146+
Assert.That(Math.Sqrt(rSquared), Is.EqualTo(0.99947).Within(0.001)); // R value
147+
Assert.That(rSquared, Is.EqualTo(0.99893).Within(0.001));
148+
Assert.That(adjustedRSquared, Is.EqualTo(0.99786).Within(0.001));
149+
Assert.That(rootMSE, Is.EqualTo(1.43148).Within(0.001));
150+
}
151+
152+
#endregion
153+
154+
#region Matrix Utility Tests
155+
156+
[Test]
157+
public void CorrelationFromCovarianceTest()
158+
{
159+
var covariance = Matrix<double>.Build.DenseOfArray(new double[,] {
160+
{4.0, 1.2, -0.8},
161+
{1.2, 9.0, 0.6},
162+
{-0.8, 0.6, 16.0}
163+
});
164+
165+
var correlation = ParameterStatistics.CorrelationFromCovariance(covariance);
166+
167+
Assert.That(correlation.RowCount, Is.EqualTo(3));
168+
Assert.That(correlation.ColumnCount, Is.EqualTo(3));
169+
170+
// Diagonal elements should be 1
171+
for (var i = 0; i < correlation.RowCount; i++)
172+
{
173+
Assert.That(correlation[i, i], Is.EqualTo(1.0).Within(1e-10));
174+
}
175+
176+
// Off-diagonal elements should be between -1 and 1
177+
for (var i = 0; i < correlation.RowCount; i++)
178+
{
179+
for (var j = 0; j < correlation.ColumnCount; j++)
180+
{
181+
if (i != j)
182+
{
183+
Assert.That(correlation[i, j], Is.GreaterThanOrEqualTo(-1.0).And.LessThanOrEqualTo(1.0));
184+
}
185+
}
186+
}
187+
188+
// Check specific values (manually calculated)
189+
Assert.That(correlation[0, 1], Is.EqualTo(0.2).Within(1e-10));
190+
Assert.That(correlation[0, 2], Is.EqualTo(-0.1).Within(1e-10));
191+
Assert.That(correlation[1, 2], Is.EqualTo(0.05).Within(1e-10));
192+
}
193+
194+
#endregion
195+
196+
#region Special Cases Tests
197+
198+
[Test]
199+
public void DependenciesTest()
200+
{
201+
// Create a correlation matrix with high multicollinearity
202+
var correlation = Matrix<double>.Build.DenseOfArray(new double[,] {
203+
{1.0, 0.95, 0.3},
204+
{0.95, 1.0, 0.2},
205+
{0.3, 0.2, 1.0}
206+
});
207+
208+
var dependencies = ParameterStatistics.DependenciesFromCorrelation(correlation);
209+
210+
Assert.That(dependencies.Count, Is.EqualTo(3));
211+
212+
// First two parameters should have high dependency values
213+
Assert.That(dependencies[0], Is.GreaterThan(0.8));
214+
Assert.That(dependencies[1], Is.GreaterThan(0.8));
215+
216+
// Third parameter should have lower dependency
217+
Assert.That(dependencies[2], Is.LessThan(0.3));
218+
}
219+
220+
[Test]
221+
public void ConfidenceIntervalsTest()
222+
{
223+
var standardErrors = Vector<double>.Build.Dense(new double[] { 0.1, 0.2, 0.5 });
224+
var df = 10; // Degrees of freedom
225+
var confidenceLevel = 0.95; // 95% confidence
226+
227+
var halfWidths = ParameterStatistics.ConfidenceIntervalHalfWidths(standardErrors, df, confidenceLevel);
228+
229+
Assert.That(halfWidths.Count, Is.EqualTo(3));
230+
231+
// t-critical for df=10, 95% confidence (two-tailed) is approximately 2.228
232+
var expectedFactor = 2.228;
233+
Assert.That(halfWidths[0], Is.EqualTo(standardErrors[0] * expectedFactor).Within(0.1));
234+
Assert.That(halfWidths[1], Is.EqualTo(standardErrors[1] * expectedFactor).Within(0.1));
235+
Assert.That(halfWidths[2], Is.EqualTo(standardErrors[2] * expectedFactor).Within(0.1));
236+
}
237+
238+
#endregion
239+
}
240+
}

‎src/Numerics/Optimization/NonlinearMinimizationResult.cs‎

Lines changed: 22 additions & 83 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
using MathNet.Numerics.Distributions;
22
using MathNet.Numerics.LinearAlgebra;
3+
using MathNet.Numerics.Statistics;
34
using System;
45
using System.Linq;
56

@@ -138,6 +139,10 @@ void EvaluateCovariance(IObjectiveModel objective)
138139
Covariance = null;
139140
Correlation = null;
140141
StandardErrors = null;
142+
TStatistics = null;
143+
PValues = null;
144+
ConfidenceIntervalHalfWidths = null;
145+
Dependencies = null;
141146
return;
142147
}
143148

@@ -148,85 +153,19 @@ void EvaluateCovariance(IObjectiveModel objective)
148153

149154
if (Covariance != null)
150155
{
151-
StandardErrors = Covariance.Diagonal().PointwiseSqrt();
152-
153-
// Calculate t-statistics and p-values
154-
TStatistics = Vector<double>.Build.Dense(StandardErrors.Count);
155-
for (var i = 0; i < StandardErrors.Count; i++)
156-
{
157-
TStatistics[i] = StandardErrors[i] > double.Epsilon
158-
? objective.Point[i] / StandardErrors[i]
159-
: double.NaN;
160-
}
161-
162-
// Calculate p-values based on t-distribution with DegreeOfFreedom
163-
PValues = Vector<double>.Build.Dense(TStatistics.Count);
164-
var tDist = new StudentT(0, 1, objective.DegreeOfFreedom);
165-
166-
// Calculate confidence interval half-widths (for 95% confidence)
167-
ConfidenceIntervalHalfWidths = Vector<double>.Build.Dense(StandardErrors.Count);
168-
var tCritical = tDist.InverseCumulativeDistribution(0.975); // Two-tailed, 95% confidence
169-
170-
for (var i = 0; i < TStatistics.Count; i++)
171-
{
172-
// Two-tailed p-value from t-distribution
173-
var tStat = Math.Abs(TStatistics[i]);
174-
PValues[i] = 2 * (1 - tDist.CumulativeDistribution(tStat));
175-
176-
// Calculate confidence interval half-width
177-
ConfidenceIntervalHalfWidths[i] = tCritical * StandardErrors[i];
178-
}
179-
180-
var correlation = Covariance.Clone();
181-
var d = correlation.Diagonal().PointwiseSqrt();
182-
var dd = d.OuterProduct(d);
183-
Correlation = correlation.PointwiseDivide(dd);
184-
185-
// Calculate dependencies (measure of multicollinearity)
186-
Dependencies = Vector<double>.Build.Dense(Correlation.RowCount);
187-
for (var i = 0; i < Correlation.RowCount; i++)
188-
{
189-
// For each parameter, calculate 1 - 1/VIF where VIF is the variance inflation factor
190-
// VIF is calculated as 1/(1-R²) where R² is the coefficient of determination when
191-
// parameter i is regressed against all other parameters
192-
193-
// Extract the correlation coefficients related to parameter i (excluding self-correlation)
194-
var correlations = Vector<double>.Build.Dense(Correlation.RowCount - 1);
195-
var index = 0;
196-
for (var j = 0; j < Correlation.RowCount; j++)
197-
{
198-
if (j != i)
199-
{
200-
correlations[index++] = Correlation[i, j];
201-
}
202-
}
203-
204-
// Calculate the square of the multiple correlation coefficient
205-
// In the simple case, this is the maximum of the squared individual correlations
206-
var maxSquaredCorrelation = 0.0;
207-
for (var j = 0; j < correlations.Count; j++)
208-
{
209-
var squaredCorr = correlations[j] * correlations[j];
210-
if (squaredCorr > maxSquaredCorrelation)
211-
{
212-
maxSquaredCorrelation = squaredCorr;
213-
}
214-
}
215-
216-
// Calculate dependency = 1 - 1/VIF
217-
var maxSquaredCorrelationCapped = Math.Min(maxSquaredCorrelation, 0.9999);
218-
var vif = 1.0 / (1.0 - maxSquaredCorrelationCapped);
219-
Dependencies[i] = 1.0 - 1.0 / vif;
220-
}
221-
}
222-
else
223-
{
224-
StandardErrors = null;
225-
TStatistics = null;
226-
PValues = null;
227-
ConfidenceIntervalHalfWidths = null;
228-
Correlation = null;
229-
Dependencies = null;
156+
// Use ParameterStatistics class to compute all statistics at once
157+
var stats = ParameterStatistics.ComputeStatistics(
158+
objective.Point,
159+
Covariance,
160+
objective.DegreeOfFreedom
161+
);
162+
163+
StandardErrors = stats.StandardErrors;
164+
TStatistics = stats.TStatistics;
165+
PValues = stats.PValues;
166+
ConfidenceIntervalHalfWidths = stats.ConfidenceIntervalHalfWidths;
167+
Correlation = stats.Correlation;
168+
Dependencies = stats.Dependencies;
230169
}
231170
}
232171

@@ -256,12 +195,12 @@ void EvaluateGoodnessOfFit(IObjectiveModel objective)
256195
{
257196
return;
258197
}
259-
198+
260199
var n = hasObservations ? objective.ObservedY.Count : objective.Residuals.Count;
261200
var dof = objective.DegreeOfFreedom;
262201

263-
// Calculate sum of squared residuals using vector operations
264-
var ssRes = objective.Residuals.DotProduct(objective.Residuals);
202+
// Calculate sum of squared residuals
203+
var ssRes = 2.0 * objective.Value;
265204

266205
// Guard against zero or negative SSR
267206
if (ssRes <= 0)
@@ -348,7 +287,7 @@ void EvaluateGoodnessOfFit(IObjectiveModel objective)
348287
}
349288

350289
// Only calculate correlation coefficient if we have model values
351-
if (hasModelValues)
290+
if (hasModelValues && hasObservations)
352291
{
353292
// Calculate correlation coefficient between observed and predicted
354293
CorrelationCoefficient = GoodnessOfFit.R(objective.ModelValues, objective.ObservedY);

0 commit comments

Comments
 (0)