1 /*
2 * Copyright (C) 2012 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.numerical;
17
18 import com.irurueta.numerical.optimization.BracketedSingleOptimizer;
19 import com.irurueta.numerical.optimization.BrentSingleOptimizer;
20
21 /**
22 * Class to estimate the most likely value from a series of samples assumed to
23 * be normally distributed.
24 * This implementation will first use an internal HistogramMaximumLikelihood
25 * to estimate the most likely value using a histogram, and then it will
26 * refine the solution by using a BrentSingleOptimizer in order to find the
27 * maximum of the probability distribution function assuming that such function
28 * is computed by aggregating small Gaussians (of size gaussianSigma) centered
29 * at the location of each sample.
30 */
31 public class AccurateMaximumLikelihoodEstimator extends MaximumLikelihoodEstimator {
32
33 /**
34 * Boolean indicating if an initial solution should be obtained first by
35 * using the Histogram method. It is suggested to always enable this option.
36 */
37 public static final boolean DEFAULT_USE_HISTOGRAM_INITIAL_SOLUTION = true;
38
39 /**
40 * Value to be considered as the machine precision.
41 */
42 public static final double EPS = 1e-9;
43
44 /**
45 * Boolean that indicates that an initial coarse solution will be computed
46 * first by using an internal HistogramMaximumLikelihoodEstimator in order
47 * to initialize the internal BrentSingleOptimizer to obtain a more accurate
48 * solution.
49 */
50 private boolean useHistogramInitialSolution;
51
52 /**
53 * Internal maximum likelihood estimator based on the Histogram method.
54 */
55 private HistogramMaximumLikelihoodEstimator internalEstimator;
56
57 /**
58 * Internal optimizer to find the true maximum of the probability
59 * distribution function. Because a BrentSingleOptimizer is only guaranteed
60 * to obtain local minima/maxima, it is preferred to start the optimizer
61 * near the true solution to be found, for that reason it is suggested to
62 * always use the Histogram initial solution as coarse approximation to
63 * start the optimizer and get a more accurate solution.
64 */
65 private BrentSingleOptimizer optimizer;
66
67 /**
68 * Constructor.
69 *
70 * @param gaussianSigma Gaussian sigma to be used on each sample.
71 * @param useHistogramInitialSolution Boolean indicating whether an internal
72 * HistogramMaximumLikelihoodEstimator will be used to obtain a coarse
73 * initial solution to initialize the BrentSingleOptimizer. It is suggested
74 * to set this value always to true.
75 * @throws IllegalArgumentException Raised if provided Gaussian sigma is
76 * negative or zero.
77 */
78 public AccurateMaximumLikelihoodEstimator(final double gaussianSigma, final boolean useHistogramInitialSolution) {
79 super(gaussianSigma);
80 this.useHistogramInitialSolution = useHistogramInitialSolution;
81 internalEstimator = null;
82 optimizer = null;
83 }
84
85 /**
86 * Empty constructor.
87 */
88 public AccurateMaximumLikelihoodEstimator() {
89 super();
90 this.useHistogramInitialSolution = DEFAULT_USE_HISTOGRAM_INITIAL_SOLUTION;
91 internalEstimator = null;
92 optimizer = null;
93 }
94
95 /**
96 * Constructor
97 *
98 * @param inputData Array containing input data where most likely value must
99 * be estimated from.
100 * @param gaussianSigma Gaussian sigma to be used on each sample.
101 * @param useHistogramInitialSolution Boolean indicating whether an internal
102 * HistogramMaximumLikelihoodEstimator will be used to obtain a coarse
103 * initial solution to initialize the BrentSingleOptimizer. It is suggested
104 * to set this value always to true.
105 * @throws IllegalArgumentException Raised if provided Gaussian sigma is
106 * negative or zero.
107 */
108 public AccurateMaximumLikelihoodEstimator(
109 final double[] inputData, final double gaussianSigma, final boolean useHistogramInitialSolution) {
110 super(inputData, gaussianSigma);
111 this.useHistogramInitialSolution = useHistogramInitialSolution;
112 internalEstimator = null;
113 optimizer = null;
114 }
115
116 /**
117 * Constructor.
118 *
119 * @param minValue Minimum value assumed to be contained within input data
120 * array.
121 * @param maxValue Maximum value assumed to be contained within input data
122 * array.
123 * @param inputData Array containing input data where most likely value must
124 * be estimated from.
125 * @param gaussianSigma Gaussian sigma to be used on each sample.
126 * @param useHistogramInitialSolution Boolean indicating whether an internal
127 * HistogramMaximumLikelihoodEstimator will be used to obtain a coarse
128 * initial solution to initialize the BrentSingleOptimizer. It is suggested
129 * to set this value always to true.
130 * @throws IllegalArgumentException Raised if provided Gaussian sigma is
131 * negative or zero, or if minValue < maxValue.
132 */
133 public AccurateMaximumLikelihoodEstimator(
134 final double minValue, final double maxValue, final double[] inputData, final double gaussianSigma,
135 final boolean useHistogramInitialSolution) {
136 super(minValue, maxValue, inputData, gaussianSigma);
137 this.useHistogramInitialSolution = useHistogramInitialSolution;
138 internalEstimator = null;
139 optimizer = null;
140 }
141
142 /**
143 * Returns method to be used for maximum likelihood estimation, which for
144 * this class is MaximumLikelihoodEstimatorMethod.
145 * ACCURATE_MAXIMUM_LIKELIHOOD_ESTIMATOR.
146 *
147 * @return Method for maximum likelihood estimation.
148 */
149 @Override
150 public MaximumLikelihoodEstimatorMethod getMethod() {
151 return MaximumLikelihoodEstimatorMethod.ACCURATE_MAXIMUM_LIKELIHOOD_ESTIMATOR;
152 }
153
154 /**
155 * Returns boolean that indicates that an initial coarse solution will be
156 * computed first by using an internal HistogramMaximumLikelihoodEstimator
157 * in order to initialize the internal BrentSingleOptimizer to obtain a more
158 * accurate solution.
159 *
160 * @return True if an initial coarse solution if found by the Histogram
161 * method, false otherwise.
162 */
163 public boolean isHistogramInitialSolutionUsed() {
164 return useHistogramInitialSolution;
165 }
166
167 /**
168 * Sets boolean that indicates that an initial coarse solution will be
169 * computed first by using an internal HistogramMaximumLikelihoodEstimator
170 * in order to initialize the internal BrentSingleOptimizer to obtain a more
171 * accurate solution.
172 *
173 * @param used True if an initial coarse solution will be found by the
174 * Histogram method, false otherwise.
175 * @throws LockedException Exception raised if this instance is locked.
176 * This method can only be executed when computations finish and this
177 * instance becomes unlocked.
178 */
179 public void setHistogramInitialSolutionUsed(final boolean used) throws LockedException {
180 if (isLocked()) {
181 throw new LockedException();
182 }
183 useHistogramInitialSolution = used;
184 }
185
186 /**
187 * Starts the estimation of the most likely value contained within provided
188 * input data array.
189 *
190 * @return The most likely value.
191 * @throws LockedException Exception raised if this instance is locked.
192 * This method can only be executed when computations finish and this
193 * instance becomes unlocked.
194 * @throws NotReadyException Exception raised if this instance is not yet
195 * ready.
196 * @see #isReady()
197 */
198 @Override
199 public double estimate() throws LockedException, NotReadyException {
200 if (isLocked()) {
201 throw new LockedException();
202 }
203 if (!isReady()) {
204 throw new NotReadyException();
205 }
206
207 locked = true;
208
209 final double minEvalPoint;
210 final double middleEvalPoint;
211 final double maxEvalPoint;
212
213 if (useHistogramInitialSolution) {
214 if (internalEstimator == null) {
215 internalEstimator = new HistogramMaximumLikelihoodEstimator();
216 }
217 internalEstimator.setInputData(inputData);
218 internalEstimator.setGaussianSigma(gaussianSigma);
219 internalEstimator.computeMinMaxValues();
220
221 middleEvalPoint = internalEstimator.estimate();
222
223 var localMinValue = 0.0;
224 var localMaxValue = 0.0;
225
226 try {
227 localMinValue = internalEstimator.getMinValue();
228 localMaxValue = internalEstimator.getMaxValue();
229 } catch (final NotAvailableException ignore) {
230 // never happens
231 }
232
233 final var numberOfBins = internalEstimator.getNumberOfBins();
234
235 final var delta = (localMaxValue - localMinValue) / (numberOfBins - 1);
236
237 // pick two values around initial coarse solution
238 minEvalPoint = middleEvalPoint - delta;
239 maxEvalPoint = middleEvalPoint + delta;
240
241 if (!areMinMaxAvailable) {
242 this.minValue = localMinValue;
243 this.maxValue = localMaxValue;
244 areMinMaxAvailable = true;
245 }
246 } else {
247 if (!areMinMaxAvailable) {
248 computeMinMaxValues();
249 }
250
251 // use min/max values as a bracket to obtain optimal solution
252 minEvalPoint = minValue;
253 maxEvalPoint = maxValue;
254 middleEvalPoint = (minValue + maxValue) * 0.5;
255 }
256
257 if ((maxValue - minValue) < EPS) {
258 // min-max limits are almost equal, so we return it as the solution
259 locked = false;
260 return middleEvalPoint;
261 }
262
263 double solution;
264 try {
265 // Use an optimizer to find maximum value on histogram (PDF)
266 if (optimizer == null) {
267 optimizer = new BrentSingleOptimizer(new EvaluatorListener(),
268 BracketedSingleOptimizer.DEFAULT_MIN_EVAL_POINT,
269 BracketedSingleOptimizer.DEFAULT_MIDDLE_EVAL_POINT,
270 BracketedSingleOptimizer.DEFAULT_MAX_EVAL_POINT,
271 BrentSingleOptimizer.DEFAULT_TOLERANCE);
272 }
273
274 optimizer.setBracket(minEvalPoint, middleEvalPoint, maxEvalPoint);
275 optimizer.minimize();
276 solution = optimizer.getResult();
277 } catch (final Exception ignore) {
278 // if optimization fails, pick coarse solution if available
279 if (useHistogramInitialSolution) {
280 solution = middleEvalPoint;
281 } else {
282 // if coarse solution is not available, then compute it
283 internalEstimator.setInputData(inputData);
284 internalEstimator.setGaussianSigma(gaussianSigma);
285 internalEstimator.computeMinMaxValues();
286 solution = internalEstimator.estimate();
287 }
288 }
289
290 locked = false;
291 return solution;
292 }
293
294 /**
295 * Internal class used by the BrentSingleOptimizer in order to evaluate
296 * the aggregation of Gaussians for all the samples in input data array with
297 * a high degree of precision.
298 */
299 private class EvaluatorListener implements SingleDimensionFunctionEvaluatorListener {
300
301
302 /**
303 * Evaluates the aggregation of Gaussians for all the samples in input
304 * data array, by assuming that each sample has an associated small
305 * Gaussian centered at the sample value and with a small sigma value.
306 * The aggregation of Gaussians will generate the averaged PDF function
307 * of all the input values.
308 *
309 * @param point Point where the aggregation of Gaussians will be
310 * evaluated
311 * @return The value of the aggregation of samples at provided point.
312 */
313 @Override
314 public double evaluate(final double point) {
315 double out = 0.0;
316 double x;
317 for (final var data : inputData) {
318 x = point - data;
319 out += Math.exp(-x * x / (2.0 * gaussianSigma * gaussianSigma))
320 / (Math.sqrt(2.0 * Math.PI) * gaussianSigma);
321 }
322
323 // negate value because optimizers always attempt to
324 // minimize function
325 return -out;
326 }
327
328 }
329 }