in lib/maths/time_series/CTimeSeriesModel.cc [1080:1218]
bool CUnivariateTimeSeriesModel::correlatedProbability(
const common::CModelProbabilityParams& params,
const TTime2Vec1Vec& time,
const TDouble2Vec1Vec& value,
common::SModelProbabilityResult& result) const {
TSize1Vec correlated;
TSize2Vec1Vec variables;
TMultivariatePriorCPtrSizePr1Vec correlationModels;
TModelCPtr1Vec correlatedTimeSeriesModels;
if (this->correlationModels(correlated, variables, correlationModels,
correlatedTimeSeriesModels) == false) {
return false;
}
double neff{effectiveCount(variables.size())};
common::CProbabilityOfExtremeSample aggregator;
common::CBasicStatistics::COrderStatisticsStack<double, 1> minProbability;
// Declared outside the loop to minimize the number of times they are created.
maths_t::EProbabilityCalculation calculation{params.calculation(0)};
TDouble10Vec1Vec sample{TDouble10Vec(2)};
maths_t::TDouble10VecWeightsAry1Vec weights{
maths_t::CUnitWeights::unit<TDouble10Vec>(2)};
TDouble10Vec2Vec pli;
TDouble10Vec2Vec pui;
TTail10Vec ti;
core_t::TTime mostAnomalousTime{0};
double mostAnomalousSample{0.0};
TPriorPtr mostAnomalousCorrelationModel;
TTail2Vec tail;
bool conditional{false};
TSize1Vec mostAnomalousCorrelate;
TSize1Vec correlateIndices;
if (params.mostAnomalousCorrelate() != std::nullopt) {
if (*params.mostAnomalousCorrelate() >= variables.size()) {
LOG_ERROR(<< "Unexpected correlate " << *params.mostAnomalousCorrelate());
return false;
}
correlateIndices.push_back(*params.mostAnomalousCorrelate());
} else {
correlateIndices.resize(variables.size());
std::iota(correlateIndices.begin(), correlateIndices.end(), 0);
}
for (std::size_t i = 0; i < correlateIndices.size(); ++i) {
if (value[i].empty()) {
aggregator.add(1.0, neff);
} else {
std::size_t correlateIndex{correlateIndices[i]};
std::size_t v0{variables[correlateIndex][0]};
std::size_t v1{variables[correlateIndex][1]};
TDecompositionPtr trendModels[2];
trendModels[v0] = m_TrendModel;
trendModels[v1] = correlatedTimeSeriesModels[correlateIndex]->m_TrendModel;
const auto& correlationModel = correlationModels[correlateIndex].first;
bool isNonNegative[2];
isNonNegative[v0] = m_IsNonNegative;
isNonNegative[v1] = correlatedTimeSeriesModels[correlateIndex]->m_IsNonNegative;
sample[0][0] = trendModels[0]->detrend(time[i][0], value[i][0],
params.seasonalConfidenceInterval(),
isNonNegative[0]);
sample[0][1] = trendModels[1]->detrend(time[i][1], value[i][1],
params.seasonalConfidenceInterval(),
isNonNegative[1]);
weights[0] = CMultivariateTimeSeriesModel::unpack(params.weights()[i]);
if (correlationModel->probabilityOfLessLikelySamples(
calculation, sample, weights, {v0}, pli, pui, ti)) {
LOG_TRACE(<< "Marginal P(" << sample
<< " | weight = " << weights << ", coordinate = " << v0
<< ") = " << (pli[0][0] + pui[0][0]) / 2.0);
LOG_TRACE(<< "Conditional P(" << sample
<< " | weight = " << weights << ", coordinate = " << v0
<< ") = " << (pli[1][0] + pui[1][0]) / 2.0);
} else {
LOG_ERROR(<< "Failed to compute P(" << sample << " | weight = " << weights
<< ", coordinate = " << v0 << ")");
continue;
}
double pl{std::sqrt(pli[0][0] * pli[1][0])};
double pu{std::sqrt(pui[0][0] * pui[1][0])};
double pi{(pl + pu) / 2.0};
aggregator.add(pi, neff);
if (minProbability.add(pi)) {
tail.assign(1, ti[0]);
conditional = ((pli[1][0] + pui[1][0]) < (pli[0][0] + pui[0][0]));
mostAnomalousCorrelate.assign(1, correlateIndex);
mostAnomalousTime = time[i][v0];
mostAnomalousSample = sample[0][v0];
mostAnomalousCorrelationModel =
conditional ? correlationModel
->univariate(NOTHING_TO_MARGINALIZE,
{{v1, sample[0][v1]}})
.first
: correlationModel
->univariate({v1}, NOTHING_TO_CONDITION)
.first;
}
}
}
double pSingleBucket;
aggregator.calculate(pSingleBucket);
TDouble4Vec probabilities{pSingleBucket};
common::SModelProbabilityResult::TFeatureProbability4Vec featureProbabilities;
featureProbabilities.emplace_back(
common::SModelProbabilityResult::E_SingleBucketProbability, pSingleBucket);
double pOverall{pSingleBucket};
if (m_AnomalyModel != nullptr && params.useAnomalyModel()) {
TDouble2Vec seasonalWeight;
this->seasonalWeight(0.0, mostAnomalousTime, seasonalWeight);
double residual{(mostAnomalousSample - mostAnomalousCorrelationModel->nearestMarginalLikelihoodMean(
mostAnomalousSample)) /
std::max(std::sqrt(seasonalWeight[0]), 1.0)};
m_AnomalyModel->sample(params, mostAnomalousTime, residual, pSingleBucket, pOverall);
double pAnomaly;
std::tie(pOverall, pAnomaly) = m_AnomalyModel->probability(pSingleBucket, pOverall);
probabilities.push_back(pAnomaly);
featureProbabilities.emplace_back(
common::SModelProbabilityResult::E_AnomalyModelProbability, pAnomaly);
}
result.s_Probability = pOverall;
result.s_Conditional = conditional;
result.s_FeatureProbabilities = std::move(featureProbabilities);
result.s_Tail = std::move(tail);
result.s_MostAnomalousCorrelate = std::move(mostAnomalousCorrelate);
return true;
}