ningshuxia
2022-10-08 afedb8fd4e17a5a911deee3dae04a10a93e6a39a
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
using IStation.Numerics.LinearAlgebra;
 
namespace IStation.Numerics.Optimization
{
    public class NonlinearMinimizationResult
    {
        public IObjectiveModel ModelInfoAtMinimum { get; }
 
        /// <summary>
        /// Returns the best fit parameters.
        /// </summary>
        public Vector<double> MinimizingPoint => ModelInfoAtMinimum.Point;
 
        /// <summary>
        /// Returns the standard errors of the corresponding parameters
        /// </summary>
        public Vector<double> StandardErrors { get; private set; }
 
        /// <summary>
        /// Returns the y-values of the fitted model that correspond to the independent values.
        /// </summary>
        public Vector<double> MinimizedValues => ModelInfoAtMinimum.ModelValues;
 
        /// <summary>
        /// Returns the covariance matrix at minimizing point.
        /// </summary>
        public Matrix<double> Covariance { get; private set; }
 
        /// <summary>
        ///  Returns the correlation matrix at minimizing point.
        /// </summary>
        public Matrix<double> Correlation { get; private set; }
 
        public int Iterations { get; }
 
        public ExitCondition ReasonForExit { get; }
 
        public NonlinearMinimizationResult(IObjectiveModel modelInfo, int iterations, ExitCondition reasonForExit)
        {
            ModelInfoAtMinimum = modelInfo;
            Iterations = iterations;
            ReasonForExit = reasonForExit;
 
            EvaluateCovariance(modelInfo);
        }
 
        private void EvaluateCovariance(IObjectiveModel objective)
        {
            objective.EvaluateAt(objective.Point); // Hessian may be not yet updated.
 
            var Hessian = objective.Hessian;
            if (Hessian == null || objective.DegreeOfFreedom < 1)
            {
                Covariance = null;
                Correlation = null;
                StandardErrors = null;
                return;
            }
 
            Covariance = Hessian.PseudoInverse() * objective.Value / objective.DegreeOfFreedom;
 
            if (Covariance != null)
            {
                StandardErrors = Covariance.Diagonal().PointwiseSqrt();
 
                var correlation = Covariance.Clone();
                var d = correlation.Diagonal().PointwiseSqrt();
                var dd = d.OuterProduct(d);
                Correlation = correlation.PointwiseDivide(dd);
            }
            else
            {
                StandardErrors = null;
                Correlation = null;
            }
        }
    }
}