1 | #region License Information |
---|
2 | /* HeuristicLab |
---|
3 | * Copyright (C) 2002-2019 Heuristic and Evolutionary Algorithms Laboratory (HEAL) |
---|
4 | * |
---|
5 | * This file is part of HeuristicLab. |
---|
6 | * |
---|
7 | * HeuristicLab is free software: you can redistribute it and/or modify |
---|
8 | * it under the terms of the GNU General Public License as published by |
---|
9 | * the Free Software Foundation, either version 3 of the License, or |
---|
10 | * (at your option) any later version. |
---|
11 | * |
---|
12 | * HeuristicLab is distributed in the hope that it will be useful, |
---|
13 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
14 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
15 | * GNU General Public License for more details. |
---|
16 | * |
---|
17 | * You should have received a copy of the GNU General Public License |
---|
18 | * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>. |
---|
19 | */ |
---|
20 | #endregion |
---|
21 | |
---|
22 | using System; |
---|
23 | using System.Collections.Generic; |
---|
24 | using System.Linq; |
---|
25 | |
---|
26 | namespace HeuristicLab.Common { |
---|
27 | public class EmpiricalCumulativeDistributionFunction { |
---|
28 | private static readonly AbscissaComparer abscissaComparer = new AbscissaComparer(); |
---|
29 | private static readonly OrdinateComparer ordinateComparer = new OrdinateComparer(); |
---|
30 | |
---|
31 | private List<Point2D<double>> ecdf; |
---|
32 | public IEnumerable<Point2D<double>> SupportingPoints { |
---|
33 | get { return ecdf; } |
---|
34 | } |
---|
35 | |
---|
36 | public EmpiricalCumulativeDistributionFunction(IList<double> sample) { |
---|
37 | ecdf = new List<Point2D<double>>(); |
---|
38 | |
---|
39 | var len = sample.Count; |
---|
40 | var cumulative = 0; |
---|
41 | var localcumulative = 0; |
---|
42 | var prev = double.NaN; |
---|
43 | foreach (var p in sample.OrderBy(x => x)) { |
---|
44 | if (double.IsNaN(p) || double.IsInfinity(p)) continue; |
---|
45 | if (!double.IsNaN(prev) && prev < p) { |
---|
46 | cumulative += localcumulative; |
---|
47 | localcumulative = 0; |
---|
48 | ecdf.Add(Point2D<double>.Create(prev, cumulative / (double)len)); |
---|
49 | } |
---|
50 | prev = p; |
---|
51 | localcumulative++; |
---|
52 | } |
---|
53 | if (!double.IsNaN(prev)) { |
---|
54 | cumulative += localcumulative; |
---|
55 | ecdf.Add(Point2D<double>.Create(prev, cumulative / (double)len)); |
---|
56 | } |
---|
57 | } |
---|
58 | public EmpiricalCumulativeDistributionFunction(IEnumerable<Point2D<double>> ecdf) { |
---|
59 | this.ecdf = new List<Point2D<double>>(); |
---|
60 | var prev = Point2D<double>.Empty; |
---|
61 | foreach (var point in ecdf) { |
---|
62 | if (point.Y < 0 || point.Y > 1 || double.IsNaN(point.X) || double.IsInfinity(point.X) |
---|
63 | || point.IsEmpty || (!prev.IsEmpty && (point.X <= prev.X || point.Y <= prev.Y))) |
---|
64 | throw new ArgumentException("Invalid supporting points of a cumulative distribution function. Must be strictly monotonically increasing in both X and Y with X in R and Y in [0;1].", "ecdf"); |
---|
65 | |
---|
66 | this.ecdf.Add(point); |
---|
67 | prev = point; |
---|
68 | } |
---|
69 | } |
---|
70 | |
---|
71 | public double Evaluate(double x) { |
---|
72 | if (ecdf.Count == 0) return double.NaN; |
---|
73 | if (x < ecdf[0].X) return 0; |
---|
74 | var last = ecdf[ecdf.Count - 1]; |
---|
75 | if (x >= last.X) return last.Y; |
---|
76 | |
---|
77 | var index = ecdf.BinarySearch(Point2D<double>.Create(x, 0), abscissaComparer); |
---|
78 | if (index >= 0) return ecdf[index].Y; |
---|
79 | return ecdf[~index - 1].Y; |
---|
80 | } |
---|
81 | |
---|
82 | public double InterpolateLinear(double x) { |
---|
83 | if (ecdf.Count == 0) return double.NaN; |
---|
84 | if (x < ecdf[0].X) return 0; |
---|
85 | var last = ecdf[ecdf.Count - 1]; |
---|
86 | if (x >= last.X) return last.Y; |
---|
87 | |
---|
88 | var index = ecdf.BinarySearch(Point2D<double>.Create(x, 0), abscissaComparer); |
---|
89 | if (index >= 0) return ecdf[index].Y; |
---|
90 | var prev = ecdf[~index - 1]; |
---|
91 | var next = ecdf[~index]; |
---|
92 | return prev.Y + (next.Y - prev.Y) * ((x - prev.X) / (next.X - prev.X)); |
---|
93 | } |
---|
94 | |
---|
95 | public double InterpolateNearest(double x) { |
---|
96 | if (ecdf.Count == 0) return double.NaN; |
---|
97 | if (x < ecdf[0].X) return 0; |
---|
98 | var last = ecdf[ecdf.Count - 1]; |
---|
99 | if (x >= last.X) return last.Y; |
---|
100 | |
---|
101 | var index = ecdf.BinarySearch(Point2D<double>.Create(x, 0), abscissaComparer); |
---|
102 | if (index >= 0) return ecdf[index].Y; |
---|
103 | var prev = ecdf[~index - 1]; |
---|
104 | var next = ecdf[~index]; |
---|
105 | if (x - prev.X < next.X - x) return prev.Y; |
---|
106 | return next.Y; |
---|
107 | } |
---|
108 | |
---|
109 | public double Inverse(double y) { |
---|
110 | if (ecdf.Count == 0) return double.NaN; |
---|
111 | if (y < 0 || y > 1) throw new ArgumentException("parameter must be in interval [0;1]", "y"); |
---|
112 | if (ecdf[ecdf.Count - 1].Y < y) return double.PositiveInfinity; |
---|
113 | var index = ecdf.BinarySearch(Point2D<double>.Create(0, y), ordinateComparer); |
---|
114 | if (index >= 0) return ecdf[index].X; |
---|
115 | return ecdf[Math.Max(~index - 1, 0)].X; |
---|
116 | } |
---|
117 | |
---|
118 | private class AbscissaComparer : Comparer<Point2D<double>> { |
---|
119 | public override int Compare(Point2D<double> x, Point2D<double> y) { |
---|
120 | return x.X.CompareTo(y.X); |
---|
121 | } |
---|
122 | } |
---|
123 | |
---|
124 | private class OrdinateComparer : Comparer<Point2D<double>> { |
---|
125 | public override int Compare(Point2D<double> x, Point2D<double> y) { |
---|
126 | return x.Y.CompareTo(y.Y); |
---|
127 | } |
---|
128 | } |
---|
129 | } |
---|
130 | } |
---|