{'fraction_non_alphanumeric': 0.04277974087161366, 'fraction_numerical': 0.010388692579505301, 'mean_word_length': 4.149320718098011, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Robust Bayesian inference in complex models with possibility theory

We propose a general solution to the problem of robust Bayesian inference in complex settings where outliers may be present. In practice, the automation of robust Bayesian analyses is important in the many applications involving large and complex datasets. The proposed solution relies on a reformulation of Bayesian inference based on possibility theory, and leverages the observation that, in this context, the marginal likelihood of the data assesses the consistency between prior and likelihood rather than model fitness. Our approach does not require additional parameters in its simplest form and has a limited impact on the computational complexity when compared to non-robust solutions. The generality of our solution is demonstrated via applications on simulated and real data including matrix estimation and change-point detection.


Introduction
Robustness is a crucial element for bridging the gap between simulation and real data: i) as summarised by George Box's famous aphorism "all models are wrong, but some are useful" [5], there will always be deviations between the model and the true data-generating mechanisms, and ii) most real data sets are contaminated by outliers which become increasingly difficult to remove manually as the dimension of the data or the size of the data set becomes larger. The need for robustness has long been acknowledged in the literature, starting with M-estimator [18], and continues to be the central motivation of a large body of work, e.g. [25,20,24,12,7,4,21]. Particularly relevant are the methods of [20] and [24], which rely on different forms of likelihood discounting, either via the use of a suitable loss function [20] or as a consequence of "coarsening" the conditioning in Bayes' theorem [24]. Related to robust inference is the problem of identifying inconsistencies between prior and data [2,27,26] where it is generally the prior that is deemed unsuitable in the presence of a conflict. The objective of this work is to study the opportunities to address robust inference which are arising when viewing Bayesian inference through the lens of possibility theory [10].
Bayesian inference is well known for the performance it naturally yields in terms of model selection, via the notion of evidence or marginal likelihood. The evidence assesses the tightness of a model, which combines a measure of the consistency between the model and the data with a measure of the simplicity of the model. When following the Bayesian approach in the context of possibility theory, the analogue of the notion of evidence only assesses the consistency between the model and the data; although this aspect implies that model selection cannot be carried out in the same way as in the standard Bayesian framework, it creates an opportunity to discount the data based on this notion of consistency and devise robust inference algorithms. The main advantages of this approach are that it is free from tuning parameter and that it can be applied in closed-form under conjugacy.
The relevant background in possibility theory will be covered in Section 2 before introducing the proposed approach to robust inference in Section 3. Limitations of our approach are then discussed in Section 4 and simulations on simulated and real data are presented in Section 5.
2 Possibility theory 2.1 General setting Possibility theory [10] can be interpreted as modelling deterministic uncertainty [35] as opposed to uncertainty stemming from random phenomena. As in probability theory, we consider a sample space Ω which contains all the possible states of nature and define a (deterministic) uncertain variable as a mapping θ from Ω to a parameter set Θ; yet, instead of endowing Ω with a probabilistic structure, we simply define ω * ∈ Ω as the true state of nature from which it follows that θ(ω * ) is the true value of the parameter of interest. The available information about θ can be described by a possibility function f θ on Θ, i.e. f θ is non-negative and verifies sup θ∈Θ f θ (θ) = 1. The law of large numbers and central limit theorem for uncertain variables [16] motivate the following notions of expected value and variance and with the variance being defined when E * (θ) is a singleton and when f θ is twice differentiable at E * (θ). The expected value verifies that E * (T (θ)) = T (E * (θ)) for any mapping T . These notions are local, which corresponds to the fact that we are only interested in a single point θ(ω * ) rather than in a whole distribution as in probability theory. In particular, E * (θ) and V * (θ) match with the mean and variance in the Laplace approximation. One important difference between possibility functions and probability density functions (p.d.f.s) is that the former are not densities; it follows that the change of variable formula for uncertain variables does not contain a Jacobian term: if θ is an uncertain variable in Θ described by the possibility function f θ and if ψ = T (θ) for some mapping T on Θ then with the assumption that sup ∅ = 0. The credibility of the event θ ∈ A is defined via the set functionP θ (A) = sup θ∈A f θ (θ) for any A ⊆ Θ; it follows thatP θ is an outer measure additionally verifyingP θ (Θ) = 1, so that we refer to it as an outer probability measure (o.p.m.). The scalarP θ (A) ∈ [0, 1] can be interpreted as the maximum subjective probability that we are ready to assign to the event θ ∈ A. If, for a fixed θ ∈ Θ, Y is a random variable in Y distributed according to a p.d.f. p Y (· | θ), then, omitting measure-theoretic details, a more general o.p.m. [14] can be defined as where 1 B is the indicator function of the set B. Defining conditioning withP θ,Y in the same way as in probability theory yields the posterior possibility function [6,36] f θ|Y (θ | y) where y is a given observation such that the evidence e(y) . = sup θ∈Θ p Y (y | θ)f θ (θ) is positive, and where " . =" emphasises that the l.h.s. is a notation for the r.h.s. The main differences between (1) and the standard Bayes theorem are that the prior is a possibility function and the denominator is based on maximisation rather than integration. These differences are small enough to allow for most of the intuition about Bayesian inference to translate to this approach and substantial enough to create opportunities at the methodological and practical levels. In particular, it is always possible to model the total absence of information a priori by considering f θ = 1, with 1 the function equal to 1 everywhere; as opposed to the probabilistic case, this prior is proper as a possibility function and therefore avoids the pitfalls of improper prior p.d.f.s [8]. The posterior expected value defined as is the maximum a posterior (MAP), which transforms coherently under re-parametrisation: if ψ = T (θ) is a new parametrisation then E * (ψ | Y = y) = T (E * (θ | Y = y)); this does not hold in general in a probabilistic setting [9], which makes the probabilistic MAP more subjective.
The considered framework allows for likelihoods defined as possibility functions, e.g. as the exponential of negative loss functions [3] or negative energy functions [30,29], which often model deterministic uncertainty. Moreover, since these functions naturally have a minimum equal to zero, the exponential of their negative counterparts have maximum one. We focus however on the case of a probabilistic likelihood in order to ease the comparison with standard Bayesian inference.

Combining information
One operation that is defined in general for possibility functions is the combination of information: if θ and ψ are uncertain variables on Θ jointly described by the possibility function f θ,ψ on Θ × Θ and if we are given the information that θ and ψ represent in fact the same unknown quantity then we can compute the conditional possibility function describing θ (equiv. ψ) given that θ = ψ as The normalising constant c = sup θ∈Θ f θ,ψ (θ, θ), which is assumed to be positive, corresponds to the credibility of the event θ = ψ, i.e. to the credibility that θ and ψ do represent the same unknown value; we will refer to c as the consistency. The conditioning used in (2) will be crucial for the proposed approach to robust inference but would lead to paradoxes [19,Chapter 15.7] when applied in the context of probability theory. We will be particularly interested in the case where θ and ψ are independently described, that is f θ,ψ (θ, ψ) = f θ (θ)f ψ (ψ) for any θ, ψ ∈ Θ. This notion of independence corresponds to the case where the available information about θ is unrelated to the one about ψ.

Conjugate priors
The concept of conjugate prior, which is central in Bayesian statistics, naturally extends to possibility functions [16] via (1): if f θ and f θ|Y (· | y) take the same parametric form then f θ is said to be a conjugate prior (possibility function) for the likelihood p Y (y | ·). In fact, each conjugate prior p.d.f. has an analogue as a possibility function up to shifts in the parameter set. For instance, for some µ ∈ R d and some d × d positive semidefinite matrix P , is the Gaussian/normal possibility function parametrised by its expected value µ and its precision matrix P and is a conjugate prior for the multivariate normal likelihood N(y; θ, Σ) for any positive definite matrix Σ. The absence of normalising constant implies that P does not have to be positive definite and, in fact, can be set to the d × d zero matrix 0 d,d in which case N(θ; µ, 0 d,d ) = 1. Another example is the beta possibility function defined for some parameters α ≥ 0 and β ≥ 0 as with expected value α/(α + β) and variance αβ/(α + β) 3 , which is a conjugate prior for Bernoulli and binomial likelihood like its probabilistic counterpart. Although the expected value matches between the probabilistic and possibilistic versions of the beta form, there is a shift by one unit in both parameters, as becomes apparent when noticing that B(0, 0) = B(1, 1) = 1, with B(α, β) the beta p.d.f. with parameters α and β. The inverse-Wishart possibility function will also be introduced in the simulations. As opposed to p.d.f.s, possibility functions are closed under power, i.e. f γ θ remains a possibility function for any γ ≥ 0. Also, if θ is described by f θ and if ψ is described by f ψ = f γ θ then E * (θ) = E * (ψ) and V * (θ) = V * (ψ)/γ. Remarkably, conjugate prior families are also closed under power, with for instance N(µ, P ) γ = N(µ, γP ) and B(θ; α, β) γ = B(θ; γα, γβ). They also all contain the uninformative prior 1 which corresponds to the case where γ = 0.
Although convenient, conjugate priors are not always applicable; in this case, synergies between probability and possibility theories can be leveraged to solve possibilistic inference problems using Markov chain Monte Carlo as in [17] or sequential Monte Carlo as in [31,32]. This is similar to the use of random exploration in optimisation problems, but with the definition of the underlying p.d.f.s being directly motivated, e.g. by seeing o.p.m.s as upper bounds for p.d.f.s.

From evidence to consistency
The evidence e(y) behaves differently from the evidence e (y) = p Y (y | θ)π(θ)dθ that would be obtained in the standard approach with a prior probability distribution π on Θ. Considering for instance the case where Θ = R, it holds that one of the prior possibility functions f θ maximising e(y) is the uninformative prior 1; this is in stark contrast with e (y) which tends to 0 if π(θ) ∝ 1 B (θ) with B = [−s, s] and s → ∞. Improper prior distributions cannot be used to evaluate the evidence as e (y) would essentially be arbitrary in this case. While the behaviour of e (y) is ideal for model selection where the fitness of the prior is key, the evidence e(y) appears to be suitable for evaluating the consistency between the observation and the prior. However, the value of e(y) depends on the choice of the reference measure when defining the likelihood p Y (y | θ), which introduces some arbitrariness. To circumvent this dependence, we first characterise the information in the likelihood via the posterior possibility function where we have assumed that the likelihood p Y (y | ·) is bounded. Indeed, the possibility function f ψ|Y (· | y), which is related to the frequentist notion of likelihood ratio test, only contains information stemming from the likelihood and from the received observation y. On the other hand, we have some prior information described by f θ , and we assume as is usual that the information encoded in f θ is unrelated to either the likelihood or the observation. It follows that θ and ψ|Y are independently described and the underlying information can be combined as Although the resulting possibility function f θ (· | θ = ψ, Y = y) is equal to the posterior possibility function f θ|Y (· | y), the corresponding consistency c(y) = sup θ∈Θ f ψ|Y (θ | y)f θ (θ) differs from e(y).
The advantages of using c(y) rather than e(y) are as follows: the scalar c(y) does not depend on the choice of reference measure when defining the likelihood and it is calibrated, i.e. c(y) ∈ [0, 1], with c(y) = 1 being achieved when f θ = 1. Henceforth, for the sake of simplicity, we will write L(θ | y) instead of f ψ|Y (θ | y) for the information in the likelihood and we will suppress the explicit conditioning on the observations in the indices, e.g. we will write f θ (· | y) instead of f θ|Y (· | y) for the posterior possibility function given Y = y.

Discounting the information in the likelihood
Discounting or tempering some or all of the components in statistical inference is a standard procedure, from the classical composite likelihood [34] to the more recent variational tempering [23], and it is particularly meaningful for possibility functions which are closed under powers. The objective in this section is to introduce a suitable discounting mechanism for the information in the likelihood in order to gain robustness properties. We first consider an example in order to better understand how the consistency behaves, especially w.r.t. the dimension of the parameter set.
The expression of the consistency c(y) as a product of the component-wise consistencies c (y i ) ∈ [0, 1], i ∈ {1, . . . , d} shows that c(y) will tend to decrease as the dimension d increases. We could however consider the geometric average d c(y) to compensate for this dependence on the dimension.
A candidate for a robust posterior possibility function follows from the observation made in Example 1 as where the discount γ is equal to d c(y), with the consistency c(y) defined as If the prior and observation are fully consistent, i.e. c(y) = 1, then there is no discount and the robust posterior f r θ (· | y) is equal to the standard posterior f θ (· | y). On the other hand, it is no longer necessary to assume that c(y) > 0 since c(y) = 0 leads to a well-defined robust posterior which is in fact equal to the prior f θ . The idea of discounting the likelihood to obtain robustness has been explored in the standard Bayesian context in [24] and in a generalised Bayesian context in [20,4] although, in both cases, the value of the discount has to be separately optimised in general. In the context of prior-data checks [2,27,26], one would generally apply the discount γ to the prior f θ , which would be equally well-defined in our formulation.
We now consider the case where several observations y 1 , . . . , y n are available as realisations of i.i.d. copies of Y . One important property of posterior probability distributions that extends to posterior possibility functions is that the way in which the observations are taken into account does not affect the overall posterior: updating with (y 1 , y 2 ) is the same as updating with y 1 and y 2 separately. This property no longer holds with the robust posterior since the consistency of y 1 and y 2 will be jointly assessed when updating with (y 1 , y 2 ); for instance, if y 1 is an outlier then both y 1 and y 2 will be mostly ignored when jointly updating with (y 1 , y 2 ). It follows that independent observations should be taken into account one by one in order to retain as much information as possible. However, this statement only applies to situations where each observation is sufficiently informative to yield a meaningful assessment of the consistency. This is particularly crucial for parameters that are not observed directly, unlike most location parameters for instance. This aspect is illustrated in the following example.
Example 2. Consider the observations y 1 , . . . , y n ∈ {0, 1} corresponding to Bernoulli experiments with unknown probability of success θ ∈ Θ = [0, 1] together with a prior beta possibility function. In this case, the information in a single observation y i , i ∈ {1, . . . , n}, is so limited that the ensuing assessment of the consistency bears little meaning. For instance, if we have two differently biased coins, then one cannot determine which coin is being flipped by looking at the data from a single flip. If n is sufficiently large, one can however consider the binomial likelihood corresponding to the observation y = y 1 + · · · + y n . Taking n large is not an issue as long as the observations y 1 , . . . , y n are either all inliers or all outliers.

Limitations and extensions 4.1 Ordering
A property of standard Bayesian inference that does not extend to the proposed approach to robust inference is the independence w.r.t. the order in which the observations are taken into account. This is not an issue in models where there is a natural order between observations such as time series, however, it is an aspect to be considered in the absence of such an order. We will show in the simulations that the ordering has little impact on the results obtained with the proposed approach when assuming that the first observation is an inlier. Regardless of the existence of a natural order between observations, a generalisation of the proposed approach would be needed in order to allow for the first observation to be an outlier. Indeed, in this case, the information in both the likelihood and the prior should be allowed to be discounted in the presence of evidence against one or the other. In particular, prior and likelihood should be similarly discounted if the information in the prior is of the same order as the one in the likelihood; this will be the topic of future work.

Information loss
The two main advantages of the approach based on (4) are the absence of parameters in the way outliers are dealt with as well as the fact that conjugate prior families can still be used. Yet, the fact that the likelihood for most observations will be discounted to some extent means that there is an information loss between the robust posterior defined in (4) and the standard posterior based on inliers only. To reduce this information loss one can introduce a threshold τ ∈ (0, 1] such that no discount is applied when γ > τ and the usual discount of gamma is applied otherwise. The corresponding posterior possibility function is where ρ(τ ) equals γ if γ ≤ τ and 1 otherwise. The value of the parameter τ does not depend strongly on the considered statistical model, and should be chosen so as to separate the values that d c(y) takes for inliers from the ones it takes for outliers. It would be possible to consider ρ(τ ) = 0 when γ ≤ τ , however, this would prevent the algorithm from recovering from early outliers that might have been assimilated due to the lack of information.

Unbounded likelihood
The information in the likelihood L(ψ | y) = f ψ (ψ | y) ∝ p Y (y | ψ), ψ ∈ Θ, might be undefined if the likelihood is unbounded. In this situation, one can make the prior possibility function f ψ weakly informative in order to ensure that ψ → p Y (y | ψ)f ψ (ψ) is bounded. The prior f ψ should however be independent from f θ so as to avoid any redundancy when combining the latter with L(· | y). Yet, when multiple observations are available, the same prior f ψ should not be used for each observation as this will artificially inflate the available information. At iteration t ∈ {1, . . . , T }, one can split the information in the prior f θ (· | y 1 , . . . , y t−1 ) into two independent bits of information as follows f r θ (θ | y 1 , . . . , y t ) From the viewpoint of possibility theory, f ω θ and f 1−ω θ are independent for any possibility function f θ and any ω ∈ [0, 1], so that there is no unwanted reuse of information in the posterior f r θ (· | y 1 , . . . , y t ). The interpretation of this recursion is as follows: 100(1 − ω)% of the posterior information given by y 1 , . . . , y t−1 is used to ensure that L t (θ | y t ) is well defined and the rest is used as a prior. The parameter ω should be kept close to 1 since the amount of information invested in L t (θ | y t ) might be lost if the observation y t is an outlier. A numerical example including an unbounded likelihood can be found in Appendix A.

Likelihood with observation-dependent support
The reliance of the proposed approach on discounting means that the likelihood cannot be tempered if it is equal to 0 on parts of the parameter space. This might be problematic as soon as the likelihood has an observation-dependent support, such as with the uniform likelihood p Y (y | θ) = θ −1 1 [y,∞) (θ). Indeed, in this case, the likelihood will overwrite the prior regardless of the consistency. One solution is to change the likelihood to make sure that p Y (y | ·) has the same support for all y, e.g. by setting exponential decays where the likelihood was originally equal to 0. An example of this is given in Appendix B. The main shortcoming of such an approach is that additional parameters have to be introduced in order to define how fast the likelihood will tend to 0. The same limitation would apply to standard probabilistic methods; however, it is possible to work in the space of probability distributions instead [25,7], in which case a likelihood of 0 can be overcome by defining a suitable notion of distance, e.g. the Wasserstein metric.

Simulations
We present three challenging applications for our method, the first and second on simulated data and the third one on real data. These applications show the generality of our approach as well as the more general strengths of possibilistic inference.

Extended feature estimation
We consider the estimation of the location and extent of a Gaussian-shaped feature in dimension d. The data is simulated so that the performance of the approach can be better assessed. The extent of the feature is modelled by a positive definite matrix Σ, sampled at random from the inverse-Wishart distribution IW(I d , d √ d), with I d the identity matrix of dimension d. There are T ∈ N iterations and, at each iteration t ∈ {1, . . . , T }, we receive n observations y t,i , i ∈ {1, . . . , n}, which are sampled independently from a normal distribution with mean zero and covariance Σ. We define Y t as the set {y t,1 , . . . , y t,n }. From an inference viewpoint, the possibility function NIW(µ, Σ; µ 0 , λ, Ψ, ν) = N(µ; µ 0 , λΣ −1 )IW(Σ; Ψ, ν) is a conjugate prior for the likelihood n i=1 N(y t,i ; µ, Σ), where IW(Ψ, ν) is the inverse Wishart possibility function with Ψ ∈ R d×d positive semi-definite and ν ≥ 0, which is defined as with | · | and tr(·) respectively denoting the determinant and the trace. The matrix Ψ is allowed not to be positive definite when ν = 0. Outliers are generated with probability and are sampled from the same data generating process as the inliers but with an inflated covariance of αΣ with α > 1.
We model the unknown mean µ and covariance Σ by the respective uncertain variables µ and Σ and define θ = (µ, Σ). We define the prior possibility function f θ as uninformative, i.e. f θ = 1, which corresponds to the NIW(µ 0 , λ 0 , Ψ 0 , ν 0 ) possibility function when λ 0 = ν 0 = 0 and Ψ 0 = 0 d,d . The updated parameters in the discounted normal inverse-Wishart model are where γ t is the discount, where S t = n i=1 (y i −ȳ t ) (y i −ȳ t ), and whereȳ t is the mean of the observations y t,1 , . . . , y t,n . The non-discounted parameters used in the expression of c t in Section 5.1 can be easily recovered by setting γ t = 1.
The consistency c t . = c t (Y t ) verifies , withȳ t the mean of the observations y t,1 , . . . , y t,n . We consider in particular the case where the dimension is d = 10 and where there are n = 25 observations at each iteration t ∈ {1, . . . , T } with T = 500. The probability of outliers is = 0.02 and we consider an inflation of the variance by a coefficient α = 3. The discount at iteration t is defined as γ t = d √ c t with d = d + d(d + 1)/2 the effective dimension of θ which follows from the fact that Σ is positive definite. We compare the original method of Section 3.2 with the alternative approach described in Section 4.2, the latter being implemented with a threshold τ = 0.25. The estimation error is defined as the Frobenius norm Σ − E * (Σ|Y 1:T ) F , with E * (Σ|Y 1:T ) = Ψ T /ν T the posterior expected value of Σ.
We first show that the results depend weakly on the arbitrary order in which the observations are considered. For this purpose, we fix a realisation of Σ and of the observations Y 1:T and consider 1000 random permutations of the indices {1, . . . , T } with the constraint that the observation at the first iteration after permutation must not be an outlier. The resulting standard deviation of the error is less than 2 × 10 −4 for both methods, which confirms the weak dependence of the error on the ordering of the observations. Next, we compute the evolution of the error as a function of the iteration t ∈ {1, . . . , T } in 4 different cases: i) standard Bayesian inference without outliers, ii) standard Bayesian inference with all observations, iii) the original method, and iv) the alternative method with a threshold τ = 0.25. The results displayed in Figure 1a show that the use of the consistency  as a discount in the likelihood, as advocated in Section 3.2, allows for largely compensating for the presence of outliers when compared to the standard Bayesian recursion. The use of a threshold allows for gaining further accuracy at the cost of adding a tuning parameter. Figure 1b illustrates the effect of the dimension on the consistency, with a sample size of n = 5d/2; it shows that, with such a linear dependence between dimension and sample size, the respective consistencies of inliers and outliers tend to stabilise when the dimension increases. The characterisation of this tendency is the topic of future work.

Kalman filtering
We now consider a truly dynamical problem where preserving the closed-form recursion of existing solutions is particularly crucial. One of the most important tools for dynamical systems is the Kalman filter which applies in the linear-Gaussian case. The Kalman filter can be shown to hold with an identical recursion in the context of possibility theory [15] and the objective of this section is to illustrate how it can easily be made more robust with the proposed approach. There has been a lot of research on robust versions of the Kalman filter, see e.g. [37,33,13] as well as more generally for online inference with hidden Markov models [22,4]. For instance, the method introduced in [13] requires the use of a pre-whitening method based on the sample median and MAD, followed by an application of the iteratively-reweighted least squares algorithm; in comparison the approach we describe is a simple modification of the Kalman filter with negligible effects on the computational time.
Specifically, we consider a scenario where the state follows a nearly-constant velocity model over T time steps of duration ∆, i.e. θ t is described by N(F θ t−1 , Q) given that θ t−1 = θ t−1 , with F the transition matrix and Q the covariance matrix of the dynamical noise. We assume that the observation y t is the realisation of a normal random variable with mean Hθ t , with H the observation matrix, and with covariance matrix R. In this context, we have L(θ t | y t ) = N(θ t ; y t , R) and the consistency is c t (y t ) = N(y t ; Hθ t , HΣ t H + R), with µ t the predicted expected value at time t and Σ t the corresponding covariance matrix. This result is the direct analogue of the marginal likelihood in a probabilistic context. We consider in particular a nearly-constant velocity model over T = 250 time steps of duration ∆ = 1, which is defined via where σ a = 0.05 is the standard deviation of the acceleration noise. The observation matrix is H = 1 0 and we consider R = 1. The true state is initialised at θ * 1 = 0 0.1 and is propagated according to θ * t = F θ * t−1 + u t with u t ∼ N(0, Q) independently from all other variables. Outliers are generated with probability = 0.1 and are sampled, at time t, from a normal distribution with   mean Hθ * t and standard deviation 10. A realisation of the state and observation processes for the considered scenario is displayed in Figure 2a as an illustration. The performance of the proposed approach is shown in Figure 2b which includes box plots of the time-averaged absolute error for 500 repeats for each of the considered methods: "std inliers" for the standard Bayesian approach with inliers only, "std all" for the standard Bayesian approach with all observations and "discount" for the proposed approach.

Change-point detection
We formulate Bayesian change-point detection [1,11,20] in the context of possibility theory as follows: a sequence of observations y 1 , . . . , y T is divided into intervals of unknown run lengths, which induces a partition Π of {1, . . . , T }. The objective is to detect the interfaces between these intervals which are referred to as change points. For the sake of simplicity, we assume that the observations within one subset S of the partition are conditionally i.i.d. according to p Y (· | θ S ), θ S ∈ Θ, with θ S = θ S for any S, S ∈ Π such that S = S . We introduce an integer-valued uncertain variable r t modelling the unknown run-length of the current interval as well as a Θ-valued uncertain variable θ modelling the unknown parameter for the current interval. The run-length r t is described by the possibility function f rt|rt−1 defined as where h(r) ∈ [0, 1] models the credibility that the total run-length of the current interval is equal to r.
Since possibility functions can be seen as upper bounds for p.d.f.s (via o.p.m.s), h(r) can be seen as the maximum subjective probability for a change point to occur after r iterations. Although it would be possible to set h(r) to 1 for all r, this would model the fact that there is no prior information on durations between change points, which would make the detection of the latter impossible; indeed, in this situation, the most likely outcome would be that a change point occurs at every iteration. Yet, one can model the limited information about interval lengths by setting h(r) to a sufficient large value. The predicted possibility function as iteration t can now be expressed as where y t :t stands for the sequence (y t , . . . , y t ) when t ≥ t and for the empty sequence when t > t and wheref rt−1,θ (· | y 1:t−1 ) is the posterior at iteration t − 1. The predicted possibility function describing θ given r t satisfies f θ|rt (θ | r, y 1:t−1 ) = f θ (θ | y t−r:t−1 ), that is, θ only depends on the last r t observations and does not otherwise depend on r t ; in particular, f θ|rt (θ | 0, y 1:t ) is equal to the prior possibility function f θ describing θ. The posterior possibility function at iteration t is then characterised bŷ As in the probabilistic algorithmic solutions to this problem, the tail of the posterior f rt (· | y 1:t ) has to be truncated in practice in order to control the computational cost; however, this truncation does not affect the possibility of non-truncated values of r t since no normalisation is needed. Change-point detection in the presence of outliers is particularly challenging in a filtering setting since only future observations can help distinguish between the two. In order to devise a robust solution to the change-point detection problem, we introduce the consistency c t (y t ) at iteration t as Although the likelihood only depends on θ, the consistency c t (y t ) also take into account the run-length r t via the conditional f θ|rt (· | r, y 1:t−1 ); e.g., if a change most likely just happened (E * (r t ) = 0) and if f θ is uninformative then c t (y t ) = 1. The robust posterior can then be defined recursively aŝ where L(θ | y t ) ∝ p Y (y t | θ) and where γ t = c t (y t ) since this model is two-dimensional.
To assess the performance of the proposed solution, we consider the classical well-log data set [28] which is a univariate time series including both change points and outliers. Within a given subset S of Π, the observations are assumed to be independently sampled from p Y (· | θ S ) = N(θ S , σ 2 ) with σ = 2500. The uncertain variable θ on Θ = R is therefore the mean of the observation process and we consider the uninformative prior f θ = N(0, 0) = 1, the normal possibility function N(µ t , λ t ) being conjugate for the normal likelihood. In order to take into account the deviations from the constant mean model considered for the behaviour between change points, a discount of 0.9 is implemented for the precision λ t , i.e. we consider that the precision evolves as λ t = 0.9λ t−1 + σ −2 ; this discounting can be defined formally for possibility functions via a power. The possibility of change point h(r) is assumed constant and equal to 2.5 × 10 −3 , which models that there is no information about specific run-lengths and that, overall, the maximum probability of change point is smaller than 1 in 400 iterations. Our method considers that four of the rapid changes in the time series are change points, as opposed to the method of [20] which considers these events as outliers. Details of two of these events in Figure 3b show that both interpretations can be deemed correct. If the size of the jump in these events indicates that the corresponding observations are indeed outliers, then a possibility function describing the likely magnitude of the jumps can be added in the model. In terms of computational time, our method takes less than 12 milliseconds per observation using a 2.3 GHz Intel Core i5, whereas the standard Bayesian procedure takes about 7 milliseconds per observation; the computational complexity is however of the same order.

Conclusion
We have introduced a general method for robust Bayesian estimation in the context of possibility theory. Our method leverages the fundamental difference in the meaning of the evidence between  possibility theory and probability theory. The obtained robustness complements the strengths of possibilistic inference in practical settings such as the availability of proper uninformative priors and the ability to take into account limited knowledge about some aspects of the model without introducing additional levels of hierarchy. This last point was illustrated in the change-point detection problem but applies more generally to complex models. Future work include the theoretical analysis of our method in terms of robustness and asymptotic properties.
Based on the change of variable formula for possibility functions, a prior on the variance s = λ −1 could be equivalently defined via the inverse-gamma possibility function characterised by IG(s; α, β) = G(s −1 ; α, β) with the added advantage that the expected value transforms coherently between the two, i.e. E * (λ) = α/β = E * (s) −1 .
The prior being defined, we can describe the general recursion at iteration t ∈ {1, . . . , T }. Instead of considering a split of the overall posterior information about θ into two parts, as suggested for the general case in Section 4.3, we consider a more careful approach and only split the information about the precision by defining for some ω ∈ (0, 1), where µ t−1 , k t−1 , α t−1 and β t−1 are the updated parameters at iteration t − 1. We then proceed with the recursion as where the dependence of L t (· | y t ) on y 1:t−1 is omitted for the sake of simplicity. We emphasise once more that, although both the prior and the likelihood depend on the previous observations, there is no double counting of information thanks to the two-way splitting of the posterior at the previous iteration. The considered model leads to the consistency Since there are two parameters in the model, we define the discount as γ t = c(y t ). In particular, one can check that if ω = k t−1 = 0 then γ t = 1; indeed, these prior parameters imply that the prior f ω θ (· | y 1:t−1 ) is uninformative. In practice, ω should not be set to 0 since this would yield the standard (non-robust) posterior. The update for the parameters of the model is as follows: The derivation of these formulas follows from the usual calculations and from G(α, β) ω = G(ωα, ωβ).
To illustrate the behaviour of the proposed approach on this model, we consider the case where µ * = 2, λ * = 1 /4, β 0 = 1 and ω = 0.99. The probability of outlier is denoted and each outlier is sampled from a standard Cauchy distribution. The alternative approach of Section 4.2 is considered with a threshold τ equal to 0.1. We focus on this approach since discounting each observation would lead to a positive bias in the precision due to inliers in the tail of the distribution being more strongly discounted than other inliers. This effect arises because of the limited information about the precision in a single observation, similarly to the case of Example 2. Figure 4 shows the evolution of the root mean square error (RMSE) in the estimation of the two parameters of interest as a function of . In addition to the standard results obtained by applying Bayesian inference to all the observations or to the inliers only, we display one of the standard ways of obtaining robust estimate of location and scale, i.e., the median and the median absolute deviation (MAD) respectively. In spite of being averaged over 1000 repeats, the error in location for Bayesian inference based all the observations displays a very erratic behaviour due to the Cauchy distribution being heavy-tailed; the corresponding graph is partially cropped in the left panel of Figure 4 in order to highlight the difference between the other methods. The performance of the proposed approach is particularly good for the precision parameter.
In order to compare with methods based on maximum mean discrepancy (MMD) [7] 1 , we also consider the simplified model where the precision of the observations is known. Under the same conditions as above, we compute the RMSE for all methods and display it in Figure 5 as a function of the probability of outliers. While the "discount" approach can be applied in this case, its performance is lesser than the other robust methods. The alternative based on a threshold τ = 0.1 allows for improving the performance. While, in this particular scenario, the MMD-based approach displays a similar performance to the median, it can be applied in more general setting including uniform likelihoods as in the following section. The average run-times per repeat on an 1.8 GHz Intel Core i7 were of the order of 1 milliseconds for the median and standard Bayesian approach, of the order of 10 milliseconds for the proposed approach and of the order of 10 seconds for the MMD-based approach.

B Likelihood with observation-dependent support
As discussed in Section 4.4, the proposed approach cannot take into account outliers that yield a credibility of 0 in the region of the true parameter, since a positive discount will have no effect on such a likelihood. In order to illustrate this aspect, we consider the uniform likelihood p Y (· | θ) = θ −1 1 [·,∞) (θ) and multiple i.i.d. observations y 1 , . . . , y T . We denote the unknown parameter by θ and assume it to The expression of c(y t ) when y t−1 ≥ s t−1 is the most useful for our purpose since it evaluates how likely is the jump from s t−1 to y t−1 given the number of observations already received. We therefore assume that the observations are ordered, i.e. y 1 ≤ · · · ≤ y T . However, raising L(θ | y t ) = Pa(θ; 1, y t ) to the power c(y t ) still leads to a posterior that is supported by [y t , ∞), which is not desirable if y t is an outlier; this is because the indicator function in L(θ | y t ) cannot be discounted.
To allow for the proposed approach to be applied, we introduce a "soft" version of L(θ | y t ) as L(θ | y t ) = y t θ ∨ y t exp c(θ − θ ∨ y t ) with c the coefficient of the exponential decay and ∨ the maximum as a binary operator (assumed to have lower precedence than multiplication but higher precedence than addition). The possibility functionL(· | y t ) is plotted in Figure 6a for different values of c.
The considered prior/posterior possibility functions cannot be considered as conjugate since their complexity increases with the number of observations. Yet, they can be computed exactly with a complexity of the order of t, which is easily achievable even for large values of t. The main limitation of this model is that it introduces an additional tuning parameter c, which must be suitably chosen.
Finding automatic ways to set c depending on prior information on the problem at hand will be the topic of future work. The proposed solution is illustrated in a scenario where the true parameter is θ * = 100 and where there are T = 100 observations. The probability of outlier is = 0.25 and outliers are sampled from a normal distribution with mean 200 and standard deviation 50. The results for one specific realisation of the observations are displayed in Figure 6b, where both the robust posterior and the standard posterior based on inliers only are included. This is a challenging case where the outliers are close to the inliers. The coefficient c of the exponential decay is set to 10 in this scenario; however this parameter depends on both θ * and T in ways that remain to be explored. The difficulty of the scenario means that, in some cases, the estimate E * (θ | y 1:T ) will be either much smaller or much larger than the true parameter θ * . This can be seen by comparing the RMSE and the median absolute error over 1000 repeats, which were found to be equal to 17.7 and 1.27 respectively.