1 // -*- C++ -*-
2 /**
3 * @brief The MixedHistogramUserDefined distribution
4 *
5 * Copyright 2005-2021 Airbus-EDF-IMACS-ONERA-Phimeca
6 *
7 * This library is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser 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 * This library 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 Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public License
18 * along with this library. If not, see <http://www.gnu.org/licenses/>.
19 *
20 */
21
22 #include <cmath>
23 #include "openturns/MixedHistogramUserDefined.hxx"
24 #include "openturns/ComposedDistribution.hxx"
25 #include "openturns/Dirac.hxx"
26 #include "openturns/Histogram.hxx"
27 #include "openturns/Uniform.hxx"
28 #include "openturns/DistFunc.hxx"
29 #include "openturns/OSS.hxx"
30 #include "openturns/Tuples.hxx"
31 #include "openturns/PersistentObjectFactory.hxx"
32 #include "openturns/ResourceMap.hxx"
33 #include "openturns/RandomGenerator.hxx"
34
35
36 BEGIN_NAMESPACE_OPENTURNS
37
38 CLASSNAMEINIT(MixedHistogramUserDefined);
39
40 static const Factory<MixedHistogramUserDefined> Factory_MixedHistogramUserDefined;
41
42 /* Default constructor */
MixedHistogramUserDefined()43 MixedHistogramUserDefined::MixedHistogramUserDefined()
44 : DistributionImplementation()
45 , ticksCollection_(1, Point(1))
46 , kind_(1, DISCRETE)
47 , probabilityTable_(1, 1.0)
48 , discreteIndices_(1, 0)
49 , allIndices_(1, 1, Indices(1, 0))
50 , normalizedProbabilityTable_(1, 1.0)
51 {
52 setName("MixedHistogramUserDefined");
53 setDimension(1);
54 computeRange();
55 }
56
57 /* Parameters constructor */
MixedHistogramUserDefined(const PointCollection & ticksCollection,const Indices & kind,const Point & probabilityTable)58 MixedHistogramUserDefined::MixedHistogramUserDefined(const PointCollection & ticksCollection,
59 const Indices & kind,
60 const Point & probabilityTable)
61 : DistributionImplementation()
62 , ticksCollection_(ticksCollection)
63 , kind_(kind)
64 , probabilityTable_(probabilityTable)
65 {
66 setName("MixedHistogramUserDefined");
67 const UnsignedInteger dimension = kind.getSize();
68 // Check the ticks
69 if (ticksCollection.getSize() != dimension) throw InvalidArgumentException(HERE) << "Error: expected a collection of ticks of size=" << dimension << ", got size=" << ticksCollection.getSize();
70 // Check the probability table
71 // kind[i] == 0 -> the ith marginal is discrete
72 // kind[i] == 1 -> the ith marginal is continuous
73 UnsignedInteger totalSize = 1;
74 for (UnsignedInteger i = 0; i < dimension; ++i)
75 {
76 if (kind[i] > 1)
77 throw InvalidArgumentException(HERE) << "Kind must be in [[0, 1]]";
78 if (!(ticksCollection[i].getSize() >= 1))
79 throw InvalidArgumentException(HERE) << "Empty ticks";
80 if (ticksCollection[i].getSize() == kind[i])
81 throw InvalidArgumentException(HERE) << "Need at least 2 ticks for continuous variable";
82 totalSize *= ticksCollection[i].getSize() - kind[i];
83 }
84 if (probabilityTable.getSize() != totalSize) throw InvalidArgumentException(HERE) << "Error: expected a probability table of size=" << totalSize << ", got size=" << probabilityTable.getSize();
85
86 // cache some data
87 Indices discretization(dimension);
88 for (UnsignedInteger i = 0; i < dimension; ++i)
89 // Here, kind[i] == 0 <-> kind[i] is false <-> i is discrete
90 discretization[i] = ticksCollection_[i].getSize() - kind_[i];
91 allIndices_ = Tuples(discretization).generate();
92 for (UnsignedInteger j = 0; j < dimension; ++ j)
93 if (kind_[j] == DISCRETE)
94 discreteIndices_.add(j);
95 else
96 continuousIndices_.add(j);
97 Scalar weightSum = 0.0;
98 for (UnsignedInteger i = 0; i < probabilityTable_.getSize(); ++ i)
99 weightSum += probabilityTable_[i];
100 normalizedProbabilityTable_ = probabilityTable_ / weightSum;
101
102 setDimension(dimension);
103 computeRange();
104 }
105
106 /* Comparison operator */
operator ==(const MixedHistogramUserDefined & other) const107 Bool MixedHistogramUserDefined::operator ==(const MixedHistogramUserDefined & other) const
108 {
109 if (this == &other) return true;
110 return (ticksCollection_ == other.ticksCollection_) && (kind_ == other.kind_) && (probabilityTable_ == other.probabilityTable_);
111 }
112
equals(const DistributionImplementation & other) const113 Bool MixedHistogramUserDefined::equals(const DistributionImplementation & other) const
114 {
115 const MixedHistogramUserDefined* p_other = dynamic_cast<const MixedHistogramUserDefined*>(&other);
116 return p_other && (*this == *p_other);
117 }
118
119 /* String converter */
__repr__() const120 String MixedHistogramUserDefined::__repr__() const
121 {
122 OSS oss(true);
123 oss << "class=" << MixedHistogramUserDefined::GetClassName()
124 << " name=" << getName()
125 << " dimension=" << getDimension()
126 << " ticksCollection=" << ticksCollection_
127 << " kind=" << kind_
128 << " probabilityTable=" << probabilityTable_;
129 return oss;
130 }
131
__str__(const String & offset) const132 String MixedHistogramUserDefined::__str__(const String & offset) const
133 {
134 OSS oss(false);
135 oss << offset << getClassName() << "(ticksCollection = " << ticksCollection_ << ", kind = " << kind_ << ", probabilityTable = " << probabilityTable_ << ")";
136 return oss;
137 }
138
139 /* Virtual constructor */
clone() const140 MixedHistogramUserDefined * MixedHistogramUserDefined::clone() const
141 {
142 return new MixedHistogramUserDefined(*this);
143 }
144
145 /* Compute the numerical range of the distribution given the parameters values */
computeRange()146 void MixedHistogramUserDefined::computeRange()
147 {
148 const UnsignedInteger dimension = getDimension();
149 Point lowerBound(dimension);
150 Point upperBound(dimension);
151 for (UnsignedInteger j = 0; j < dimension; ++j)
152 {
153 const Point ticks(ticksCollection_[j]);
154 lowerBound[j] = ticks[0];
155 upperBound[j] = ticks[0];
156 for (UnsignedInteger k = 1; k < ticks.getSize(); ++ k)
157 {
158 if (ticks[k] < lowerBound[j])
159 lowerBound[j] = ticks[k];
160 if (ticks[k] > upperBound[j])
161 upperBound[j] = ticks[k];
162 }
163 }
164 setRange(Interval(lowerBound, upperBound));
165 }
166
167
168 /* Get one realization of the distribution */
getRealization() const169 Point MixedHistogramUserDefined::getRealization() const
170 {
171 const UnsignedInteger dimension = getDimension();
172 if (!base_.getSize())
173 (void) DistFunc::rDiscrete(normalizedProbabilityTable_, base_, alias_);
174 const UnsignedInteger index = DistFunc::rDiscrete(base_, alias_);
175 Point realization(dimension);
176 for (UnsignedInteger j = 0; j < discreteIndices_.getSize(); ++j)
177 {
178 const UnsignedInteger jDiscrete = discreteIndices_[j];
179 const UnsignedInteger k = allIndices_(index, jDiscrete);
180 const Point ticks(ticksCollection_[jDiscrete]);
181 realization[jDiscrete] = ticks[k];
182 }
183 for (UnsignedInteger j = 0; j < continuousIndices_.getSize(); ++j)
184 {
185 const UnsignedInteger jContinuous = continuousIndices_[j];
186 const UnsignedInteger k = allIndices_(index, jContinuous);
187 const Point ticks(ticksCollection_[jContinuous]);
188 realization[jContinuous] = ticks[k] + (ticks[k + 1] - ticks[k]) * RandomGenerator::Generate();
189 }
190 return realization;
191 }
192
193 /* Get a sample of the distribution */
getSample(const UnsignedInteger size) const194 Sample MixedHistogramUserDefined::getSample(const UnsignedInteger size) const
195 {
196 return DistributionImplementation::getSample(size);
197 }
198
199 /* Get the PDF of the distribution */
computePDF(const Point & point) const200 Scalar MixedHistogramUserDefined::computePDF(const Point & point) const
201 {
202 const UnsignedInteger dimension = getDimension();
203 if (point.getDimension() != dimension)
204 throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=" << dimension << ", here dimension=" << point.getDimension();
205
206 // build the list of discrete ticks indices, with early exit if no tick matches
207 Indices discreteTicksIndices(discreteIndices_.getSize());
208 for (UnsignedInteger j = 0; j < discreteIndices_.getSize(); ++ j)
209 {
210 const Point ticks(ticksCollection_[discreteIndices_[j]]);
211 const UnsignedInteger index = ticks.find(point[discreteIndices_[j]]);
212 if (index >= ticks.getSize())
213 return 0.0;
214 discreteTicksIndices[j] = index;
215 }
216
217 // loop over the cpt
218 Scalar pdfValue = 0.0;
219 const UnsignedInteger totalSize = probabilityTable_.getSize();
220 for (UnsignedInteger i = 0; i < totalSize; ++i)
221 {
222 Bool skip = false;
223 // first, loop over the discrete components
224 for (UnsignedInteger j = 0; j < discreteIndices_.getSize(); ++j)
225 {
226 const UnsignedInteger k = allIndices_(i, discreteIndices_[j]);
227 if (discreteTicksIndices[j] != k)
228 {
229 skip = true;
230 break;
231 }
232 }
233
234 // exclude non-matching discrete terms
235 if (skip)
236 continue;
237
238 skip = false;
239
240 // now compute the pdf over continuous components
241 Scalar pdfI = 1.0;
242 for (UnsignedInteger j = 0; j < continuousIndices_.getSize(); ++j)
243 {
244 const UnsignedInteger k = allIndices_(i, continuousIndices_[j]);
245 const Point ticks(ticksCollection_[continuousIndices_[j]]);
246 const Scalar x = point[continuousIndices_[j]];
247 if ((x <= ticks[k]) || (x > ticks[k + 1]))
248 {
249 skip = true;
250 break;
251 }
252 pdfI *= 1.0 / (ticks[k + 1] - ticks[k]);
253 }
254
255 // exclude non-matching continuous terms
256 if (skip)
257 continue;
258
259 pdfValue += normalizedProbabilityTable_[i] * pdfI;
260 }
261 return pdfValue;
262 }
263
264
265 /* Get the CDF of the distribution */
computeCDF(const Point & point) const266 Scalar MixedHistogramUserDefined::computeCDF(const Point & point) const
267 {
268 const UnsignedInteger dimension = getDimension();
269 if (point.getDimension() != dimension)
270 throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=" << dimension << ", here dimension=" << point.getDimension();
271
272 // build the list of discrete ticks, with early exit if no tick matches
273 Indices discreteTicksIndices(discreteIndices_.getSize());
274 for (UnsignedInteger j = 0; j < discreteIndices_.getSize(); ++ j)
275 {
276 const Point ticks(ticksCollection_[discreteIndices_[j]]);
277 const Scalar x = point[discreteIndices_[j]];
278 UnsignedInteger index = ticks.getSize();
279 for (UnsignedInteger i = 0; i < ticks.getSize(); ++i)
280 if (ticks[i] <= x)
281 index = i;
282 if (index >= ticks.getSize())
283 return 0.0;
284 discreteTicksIndices[j] = index;
285 }
286
287 // loop over the cpt
288 Scalar cdfValue = 0.0;
289 const UnsignedInteger totalSize = probabilityTable_.getSize();
290 for (UnsignedInteger i = 0; i < totalSize; ++i)
291 {
292 Bool skip = false;
293 // first, loop over the discrete components
294 for (UnsignedInteger j = 0; j < discreteIndices_.getSize(); ++j)
295 {
296 const UnsignedInteger k = allIndices_(i, discreteIndices_[j]);
297 if (k > discreteTicksIndices[j])
298 {
299 skip = true;
300 break;
301 }
302 }
303
304 // exclude non-matching discrete terms
305 if (skip)
306 continue;
307
308 skip = false;
309
310 // now compute the cdf over continuous components
311 Scalar cdfI = 1.0;
312 for (UnsignedInteger j = 0; j < continuousIndices_.getSize(); ++j)
313 {
314 const UnsignedInteger k = allIndices_(i, continuousIndices_[j]);
315 const Point ticks(ticksCollection_[continuousIndices_[j]]);
316 const Scalar x = point[continuousIndices_[j]];
317 if (x <= ticks[k])
318 {
319 skip = true;
320 break;
321 }
322 else if (x < ticks[k + 1])
323 cdfI *= (x - ticks[k]) / (ticks[k + 1] - ticks[k]);
324 //else (x>ticks[k + 1]: nothing to do
325 }
326
327 // exclude non-matching continuous terms
328 if (skip)
329 continue;
330
331 cdfValue += normalizedProbabilityTable_[i] * cdfI;
332 }
333 return cdfValue;
334 }
335
computeComplementaryCDF(const Point & point) const336 Scalar MixedHistogramUserDefined::computeComplementaryCDF(const Point & point) const
337 {
338 return DistributionImplementation::computeComplementaryCDF(point);
339 }
340
341 /* Get the characteristic function of the distribution, i.e. phi(u) = E(exp(I*u*X)) */
computeCharacteristicFunction(const Scalar x) const342 Complex MixedHistogramUserDefined::computeCharacteristicFunction(const Scalar x) const
343 {
344 return DistributionImplementation::computeCharacteristicFunction(x);
345 }
346
347 /* Get the quantile of the distribution */
computeQuantile(const Scalar prob,const Bool tail) const348 Point MixedHistogramUserDefined::computeQuantile(const Scalar prob,
349 const Bool tail) const
350 {
351 return DistributionImplementation::computeQuantile(prob, tail);
352 }
353
354 /* Get the i-th marginal distribution */
getMarginal(const UnsignedInteger index) const355 Distribution MixedHistogramUserDefined::getMarginal(const UnsignedInteger index) const
356 {
357 const UnsignedInteger dimension = getDimension();
358 if (index >= dimension) throw InvalidArgumentException(HERE) << "The index of a marginal distribution must be in the range [0, dim-1]";
359 if (dimension == 1) return clone();
360
361 // contract probability table
362 const Point ticks(ticksCollection_[index]);
363 const UnsignedInteger size = ticks.getSize();
364 Point marginalProbabilityTable((kind_[index] == DISCRETE) ? size : size - 1);
365 const UnsignedInteger totalSize = probabilityTable_.getSize();
366 for (UnsignedInteger i = 0; i < totalSize; ++i)
367 {
368 const UnsignedInteger k = allIndices_(i, index);
369 marginalProbabilityTable[k] += probabilityTable_[i];
370 }
371
372 Distribution marginal;
373 if (kind_[index] == DISCRETE)
374 {
375 SampleImplementation support(size, 1);
376 support.setData(ticks);
377 marginal = UserDefined(support, marginalProbabilityTable);
378 }
379 else
380 {
381 marginal = Histogram(ticks, marginalProbabilityTable);
382 }
383 marginal.setDescription(Description(1, getDescription()[index]));
384 return marginal;
385 }
386
387 /* Get the distribution of the marginal distribution corresponding to indices dimensions */
getMarginal(const Indices & indices) const388 Distribution MixedHistogramUserDefined::getMarginal(const Indices & indices) const
389 {
390 const UnsignedInteger dimension = getDimension();
391 if (!indices.check(dimension)) throw InvalidArgumentException(HERE) << "The indices of a marginal distribution must be in the range [0, dim-1] and must be different";
392
393 Indices full(dimension);
394 full.fill();
395 if (indices == full) return clone();
396
397 // contract probability table
398 Indices marginalKind;
399 PointCollection marginalTicksCollection;
400 UnsignedInteger marginalTotalSize = 1;
401 Description description(getDescription());
402 Description marginalDescription;
403 Indices discretization;
404 for (UnsignedInteger j = 0; j < indices.getSize(); ++j)
405 {
406 const UnsignedInteger index = indices[j];
407 marginalKind.add(kind_[index]);
408 marginalTicksCollection.add(ticksCollection_[index]);
409 const UnsignedInteger size = ticksCollection_[index].getSize();
410 discretization.add((kind_[index] == DISCRETE) ? size : size - 1);
411 marginalTotalSize *= discretization[j];
412 marginalDescription.add(description[index]);
413 }
414 IndicesCollection marginalAllIndices(Tuples(discretization).generate());
415 Point marginalProbabilityTable(marginalTotalSize);
416 const UnsignedInteger totalSize = probabilityTable_.getSize();
417
418 // compute the base of discretization to quickly retrieve the global marginal index
419 Indices productDiscretization(indices.getSize(), 1);
420 for (UnsignedInteger j = 1; j < indices.getSize(); ++j)
421 {
422 productDiscretization[j] = productDiscretization[j - 1] * discretization[j - 1];
423 }
424
425 for (UnsignedInteger i = 0; i < totalSize; ++i)
426 {
427 for (UnsignedInteger j = 0; j < indices.getSize(); ++j)
428 {
429 // find the global marginal index
430 UnsignedInteger marginalProbabilityTableIndex = 0;
431 for (UnsignedInteger k = 0; k < indices.getSize(); ++k)
432 {
433 marginalProbabilityTableIndex += allIndices_(i, indices[k]) * productDiscretization[k];
434 }
435 marginalProbabilityTable[marginalProbabilityTableIndex] += probabilityTable_[i];
436 }
437 }
438
439 MixedHistogramUserDefined marginal(marginalTicksCollection, marginalKind, marginalProbabilityTable);
440 marginal.setDescription(marginalDescription);
441 return marginal;
442 } // getMarginal(Indices)
443
444 /* Check if the distribution is continuous */
isContinuous() const445 Bool MixedHistogramUserDefined::isContinuous() const
446 {
447 const UnsignedInteger size = kind_.getSize();
448 for (UnsignedInteger i = 0; i < size; ++i)
449 if (kind_[i] == DISCRETE) return false;
450 return true;
451 }
452
453 /* Check if the distribution is discrete */
isDiscrete() const454 Bool MixedHistogramUserDefined::isDiscrete() const
455 {
456 const UnsignedInteger size = kind_.getSize();
457 for (UnsignedInteger i = 0; i < size; ++i)
458 if (kind_[i] == CONTINUOUS) return false;
459 return true;
460 }
461
462 /* Check if the distribution is integral */
isIntegral() const463 Bool MixedHistogramUserDefined::isIntegral() const
464 {
465 const Scalar epsilon = ResourceMap::GetAsScalar("DiscreteDistribution-SupportEpsilon");
466 const UnsignedInteger size = kind_.getSize();
467 for (UnsignedInteger i = 0; i < size; ++i)
468 {
469 if (kind_[i] == CONTINUOUS) return false;
470 const UnsignedInteger supportSize = ticksCollection_[i].getSize();
471 for (UnsignedInteger j = 0; j < supportSize; ++j)
472 {
473 const Scalar x = ticksCollection_[i][j];
474 if (std::abs(x - floor(x + 0.5)) > epsilon) return false;
475 }
476 } // i
477 return true;
478 }
479
480
481 /* Compute the mean of the distribution */
computeMean() const482 void MixedHistogramUserDefined::computeMean() const
483 {
484 const UnsignedInteger dimension = getDimension();
485 mean_ = Point(dimension);
486 const UnsignedInteger totalSize = probabilityTable_.getSize();
487 for (UnsignedInteger i = 0; i < totalSize; ++i)
488 {
489 Point meani(dimension);
490 for (UnsignedInteger j = 0; j < dimension; ++j)
491 {
492 const UnsignedInteger k = allIndices_(i, j);
493 const Point ticks(ticksCollection_[j]);
494 if (kind_[j] == DISCRETE)
495 meani[j] = ticks[k];
496 else
497 meani[j] = 0.5 * (ticks[k] + ticks[k + 1]);
498 }
499 mean_ += meani * normalizedProbabilityTable_[i];
500 }
501 isAlreadyComputedMean_ = true;
502 }
503
504 /* Get the standard deviation of the distribution */
getStandardDeviation() const505 Point MixedHistogramUserDefined::getStandardDeviation() const
506 {
507 const UnsignedInteger dimension = getDimension();
508 Point standardDeviation(dimension);
509 for (UnsignedInteger i = 0; i < dimension; ++i)
510 standardDeviation[i] = getMarginal(i).getStandardDeviation()[0];
511 return standardDeviation;
512 }
513
514 /* Get the skewness of the distribution */
getSkewness() const515 Point MixedHistogramUserDefined::getSkewness() const
516 {
517 const UnsignedInteger dimension = getDimension();
518 Point skewness(dimension);
519 for (UnsignedInteger i = 0; i < dimension; ++i)
520 skewness[i] = getMarginal(i).getSkewness()[0];
521 return skewness;
522 }
523
524 /* Get the kurtosis of the distribution */
getKurtosis() const525 Point MixedHistogramUserDefined::getKurtosis() const
526 {
527 const UnsignedInteger dimension = getDimension();
528 Point kurtosis(dimension);
529 for (UnsignedInteger i = 0; i < dimension; ++i)
530 kurtosis[i] = getMarginal(i).getKurtosis()[0];
531 return kurtosis;
532 }
533
534 /* Compute the covariance of the distribution */
computeCovariance() const535 void MixedHistogramUserDefined::computeCovariance() const
536 {
537 const UnsignedInteger dimension = getDimension();
538 covariance_ = CovarianceMatrix(dimension);
539 for (UnsignedInteger j = 0; j < dimension; ++j)
540 covariance_(j, j) = 0.0;
541 // First, compute E(X.X^t)
542 const UnsignedInteger totalSize = probabilityTable_.getSize();
543 for (UnsignedInteger i = 0; i < totalSize; ++i)
544 {
545 Point meanI(dimension);
546 Point varianceI(dimension);
547 for (UnsignedInteger j = 0; j < dimension; ++j)
548 {
549 const UnsignedInteger k = allIndices_(i, j);
550 const Point ticks(ticksCollection_[j]);
551 if (kind_[j] == DISCRETE)
552 meanI[j] = ticks[k];
553 else
554 {
555 meanI[j] = 0.5 * (ticks[k] + ticks[k + 1]);
556 const Scalar eta = ticks[k + 1] - ticks[k];
557 varianceI[j] = eta * eta / 12.0;
558 }
559 }
560 for(UnsignedInteger row = 0; row < dimension; ++row)
561 for(UnsignedInteger column = 0; column <= row; ++column)
562 covariance_(row, column) += normalizedProbabilityTable_[i] * (((row == column) ? varianceI[row] : 0.0) + meanI[row] * meanI[column]);
563 }
564 // Then, subtract E(X).E(X)^t
565 const Point mean(getMean());
566 for(UnsignedInteger row = 0; row < dimension; ++row)
567 for(UnsignedInteger column = 0; column <= row; ++column)
568 covariance_(row, column) -= mean[row] * mean[column];
569 isAlreadyComputedCovariance_ = true;
570 }
571
572 /* Get the moments of the standardized distribution */
getStandardMoment(const UnsignedInteger n) const573 Point MixedHistogramUserDefined::getStandardMoment(const UnsignedInteger n) const
574 {
575 const UnsignedInteger dimension = getDimension();
576 Point standardMoment(dimension);
577 for (UnsignedInteger i = 0; i < dimension; ++i)
578 standardMoment[i] = getMarginal(i).getStandardMoment(n)[0];
579 return standardMoment;
580 }
581
582 /* Get the standard representative in the parametric family, associated with the standard moments */
getStandardRepresentative() const583 Distribution MixedHistogramUserDefined::getStandardRepresentative() const
584 {
585 return clone();
586 }
587
588 /* Ticks collection accessor */
setTicksCollection(const Collection<Point> & ticksCollection)589 void MixedHistogramUserDefined::setTicksCollection(const Collection<Point> & ticksCollection)
590 {
591 *this = MixedHistogramUserDefined(ticksCollection, kind_, probabilityTable_);
592 computeRange();
593 }
594
getTicksCollection() const595 Collection<Point> MixedHistogramUserDefined::getTicksCollection() const
596 {
597 return ticksCollection_;
598 }
599
600 /* Kind accessor */
setKind(const Indices & kind)601 void MixedHistogramUserDefined::setKind(const Indices & kind)
602 {
603 *this = MixedHistogramUserDefined(ticksCollection_, kind, probabilityTable_);
604 computeRange();
605 }
606
getKind() const607 Indices MixedHistogramUserDefined::getKind() const
608 {
609 return kind_;
610 }
611
612 /* Probability table accessor */
setProbabilityTable(const Point & probabilityTable)613 void MixedHistogramUserDefined::setProbabilityTable(const Point & probabilityTable)
614 {
615 *this = MixedHistogramUserDefined(ticksCollection_, kind_, probabilityTable);
616 computeRange();
617 }
618
getProbabilityTable() const619 Point MixedHistogramUserDefined::getProbabilityTable() const
620 {
621 return probabilityTable_;
622 }
623
624 /* Conversion as a Mixture */
asMixture() const625 Mixture MixedHistogramUserDefined::asMixture() const
626 {
627 const UnsignedInteger dimension = getDimension();
628 const UnsignedInteger totalSize = probabilityTable_.getSize();
629 Mixture mixture;
630 // Special case: dimension 1
631 if (dimension == 1)
632 {
633 const Point ticks(ticksCollection_[0]);
634 if (kind_[0] == DISCRETE)
635 {
636 const UnsignedInteger size = ticks.getSize();
637 SampleImplementation support(size, 1);
638 support.setData(ticks);
639 mixture = Mixture(Mixture::DistributionCollection(1, UserDefined(support, probabilityTable_)));
640 }
641 // Continuous
642 else
643 mixture = Mixture(Mixture::DistributionCollection(1, Histogram(ticks, probabilityTable_)));
644 } // dimension == 1
645 else
646 {
647 Bool allDiscrete = (discreteIndices_.getSize() == dimension);
648 // Multivariate discrete
649 if (allDiscrete)
650 {
651 Sample support(totalSize, dimension);
652 for (UnsignedInteger i = 0; i < totalSize; ++i)
653 {
654 for (UnsignedInteger j = 0; j < dimension; ++j)
655 support(i, j) = ticksCollection_[j][allIndices_(i, j)];
656 }
657 mixture = Mixture(Mixture::DistributionCollection(1, UserDefined(support, probabilityTable_)));
658 } // allDiscrete
659 else
660 {
661 Mixture::DistributionCollection atoms(totalSize);
662 for (UnsignedInteger i = 0; i < totalSize; ++i)
663 {
664 Mixture::DistributionCollection subAtoms(dimension);
665 for (UnsignedInteger j = 0; j < dimension; ++j)
666 {
667 const UnsignedInteger k = allIndices_(i, j);
668 const Point ticks = ticksCollection_[j];
669 if (kind_[j] == DISCRETE)
670 subAtoms[j] = Dirac(Point(1, ticks[k]));
671 // Continuous
672 else
673 subAtoms[j] = Uniform(ticks[k], ticks[k + 1]);
674 } // j
675 atoms[i] = ComposedDistribution(subAtoms);
676 } // i
677 mixture = Mixture(atoms, probabilityTable_);
678 } // At least one continuous
679 } // dimension > 1
680 mixture.setDescription(getDescription());
681 return mixture;
682 }
683
save(Advocate & adv) const684 void MixedHistogramUserDefined::save(Advocate & adv) const
685 {
686 DistributionImplementation::save(adv);
687 adv.saveAttribute( "ticksCollection_", ticksCollection_ );
688 adv.saveAttribute( "kind_", kind_ );
689 adv.saveAttribute( "probabilityTable_", probabilityTable_ );
690 }
691
692 /* Method load() reloads the object from the StorageManager */
load(Advocate & adv)693 void MixedHistogramUserDefined::load(Advocate & adv)
694 {
695 DistributionImplementation::load(adv);
696 adv.loadAttribute( "ticksCollection_", ticksCollection_ );
697 adv.loadAttribute( "kind_", kind_ );
698 adv.loadAttribute( "probabilityTable_", probabilityTable_ );
699 computeRange();
700 }
701
702 /* Description accessor */
setDescription(const Description & description)703 void MixedHistogramUserDefined::setDescription(const Description & description)
704 {
705 DistributionImplementation::setDescription(description);
706 }
707
708 END_NAMESPACE_OPENTURNS
709