Keywords: Linear inverse problem, Force reconstruction, Bayesian inference, Generalized Gaussian priors, Gibbs sampler, Hamiltonian Monte Carlo

In a previous paper, the authors introduced a flexible methodology for reconstructing mechanical sources in the frequency domain from prior local information on both their nature and location over a linear and time invariant structure. The proposed approach was derived from Bayesian statistics, because of its ability in mathematically accounting for experimenter's prior knowledge. However, since only the Maximum a Posteriori estimate was computed, the posterior uncertainty about the regularized solution given the measured vibration field, the mechanical model and the regularization parameter was not assessed. To answer this legitimate question, this paper fully exploits the Bayesian framework to provide, from a Markov Chain Monte Carlo algorithm, credible intervals and other statistical measures (mean, median, mode) for all the parameters of the force reconstruction problem.

On a full Bayesian inference for force reconstruction problems

Introduction

In structural dynamics, most of the research papers are generally focused on the modeling of the dynamic response of structures subjected to perfectly determined excitation sources. However, the latter are only roughly or partially known. As a result, an error in the definition of the excitation vector is propagated to the dynamic response through the model and can have a significant impact on the subsequent mechanical analysis. A possible solution to this problem consists in directly measuring the excitation sources. The main drawback of such an intrusive approach is twofold. First, it is potentially time-consuming to implement and calibrate. Second, direct measurements can be impossible to perform on industrial structures. To circumvent these undesirable features, a possible alternative is to perform indirect measurements using accessible quantities such as displacement or acceleration fields and a model of the dynamic behavior of the studied structure. Unfortunately, the reconstruction of mechanical sources from vibration measurements is an ill-posed inverse problem. A classical approach to bypass this difficulty consists in constraining the space of solutions by using prior information on the noise and the sources to reconstruct. A convenient and efficient way to deal with such prior information is the Bayesian framework, because it allows combining both probabilistic and mechanical data. The most widespread and popular approaches deriving from Bayesian statistics are certainly additive regularizations, such as the Tikhonov regularization [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF] or the Lasso regularization [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. Formally, they correspond to the maximum a posteriori (MAP) solution when normal distributions [START_REF] Busby | Optimal regularization of an inverse dynamics problem[END_REF][START_REF] Thite | The quantification of structure-borne transmission paths by inverse methods. Part 2 : Use of regularization techniques[END_REF][START_REF] Liu | Reducing the impact of measurement errors when reconstructing dynamic forces[END_REF][START_REF] Renzi | Vibratory source identification by using the finite element model of a subdomain of a flexural beam[END_REF] or generalized Gaussian distributions [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF][START_REF] Aucejo | Bayesian source identification using local priors[END_REF] are used to mathematically express one's prior knowledge of the reconstruction problem. However, the major drawback of additive regularizations is that only a point estimate is generally computed, even if, in the case of the standard Tikhonov regularization (normal distributions), the covariance matrix associated to regularized solution can be easily obtained. Consequently, the experimenter has often no information about the uncertainty associated to the identified solution given the measured vibration field, the chosen mechanical model and the estimated regularization parameter. To fully exploit the Bayesian framework and obtain the credibility margins as well as other statistical measures on the identified solutions, the reconstruction problem is generally solved from Markov Chain Monte Carlo (MCMC) algorithms [START_REF] Kroese | Handbook of Monte Carlo Methods[END_REF][START_REF] Brooks | Handbook of Markov Chain Monte Carlo[END_REF][START_REF] Gelman | Bayesian Data Analysis[END_REF]. This has the great advantage of exploring the posterior probability distribution by drawing random samples and thus determine statistical indicators for all the parameters of the inverse problem. To the authors' knowledge, only a few studies have fully exploited the Bayesian framework in the field of force identification. One can nevertheless cite the work of Zhang et al. [START_REF] Zhang | Bayesian force reconstruction with an uncertain model[END_REF], in which they apply MCMC algorithm to solve the reconstruction problem considering an uncertain model as well as Gaussian priors for the noise and the sources and gamma distributions for the related hyperparameters. In the same vein, Faure et al. [START_REF] Faure | Vibration source characterization using force analysis technique and a Bayesian regularization[END_REF] present a full Bayesian approach of the inverse problem based on a model deriving from the equation of motion of the structure as in the Force Analysis Technique [START_REF] Pezerat | Two inverse methods for localization of external sources exciting a beam[END_REF]. In that contribution, two formulations are introduced. The first one is rather classical since Gaussian and gamma priors are employed, while the second one makes use of a Bernoulli-Gaussian prior in order to identify sparse source distribution.

The present paper presents a full Bayesian treatment applied to the formulation introduced in Ref. [START_REF] Aucejo | Bayesian source identification using local priors[END_REF], in which a MAP estimate has been obtained by considering generalized Gaussian distributions to exploit prior local information (nature and location) on the mechanical sources acting on a linear and time-invariant structure. More specifically, the proposed extension of our previous work aims at not only assessing the credibility of the identified solution given experimental data and a mechanical model, but also estimating the posterior uncertainties of all the a priori unknown parameters appearing in the formulation. To this end, the article is divided into four parts.

In section 2, the Bayesian formulation of the reconstruction problem is first introduced. This section allows highlighting the original features of the proposed formulation. Then, section 3 details the MCMC algorithm used to solve the problem. Practically, it consists of a Gibbs sampler [START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF][START_REF] Gelfand | Sampling-Based Approach to Calculating Marginal Densities[END_REF] including a Hamiltonian Monte Carlo [START_REF] Duane | Hybrid monte carlo[END_REF][START_REF] Neal | Hanbook of Markov Chain Monte Carlo, chapter 5 -MCMC Using Hamiltonian Dynamics[END_REF] step to draw samples for the parameters having a non-standard full conditional probability distribution.

Finally, the ability of the proposed approach in identifying all the parameters of the problem in the frequency domain is illustrated using synthetic and experimental data in sections 4 and 5, where partly sparse sources are reconstructed. In particular, the analysis of the credible intervals calculated for the identified excitation field reveals that the proposed formulation allows obtaining credible solutions.

Bayesian formulation of the reconstruction problem

This section aims at introducing the formulation of the Bayesian force reconstruction problem as well as the related parameters. To render this section more didactic, the formulation is introduced step by step in order to better highlight the main features of the proposed approach. It should be noticed that section 2.1 is partially extracted from Ref. [START_REF] Aucejo | Bayesian source identification using local priors[END_REF]. It is included here in order to make the article self-contained and to better explain the need for defining a more complete formulation.

Standard Bayesian formulation

Let us consider the practical situation where the vibration field X, measured over the surface of a structure, is caused by an unknown excitation field F. If the structure is linear and time invariant, its dynamic behavior is completely determined by the transfer functions matrix H, relating the vibration field X to the excitation field F, so that:

X = HF. (1)
Practically, this relation cannot be satisfied for two main reasons. First, the vibration field X is necessarily corrupted by measurement errors, because the transducers used to performed measurements are imperfect. Second, the transfer functions matrix H can be biased either by measurement errors if it is measured or by modeling errors if it is computed. Assuming that modeling errors are small enough, it is reasonable to model these two biases as an additive noise N. As a result, the reconstruction model is given by the following direct formulation:

X = HF + N. (2)
However, the force reconstruction problem consists in estimating the unknown excitation field F acting on a structure from the knowledge of the measured vibration field X on its surface and the transfer functions matrix H. In other words, the noise vector N is another unknown of the reconstruction problem. To solve this problem efficiently, the Bayesian framework is adopted. Formally, the Bayesian paradigm consists in considering all the parameters of the problem as random variables. Consequently, the uncertainty about each parameter is modeled by a probability distribution, describing the state of knowledge or the prior on this parameter. From a mathematical standpoint, the Bayesian reconstruction problem relies on the Bayes' rule:

p(F|X) ∝ p(X|F) p(F), (3) where:
• p(F|X) is the posterior probability distribution, representing the probability of observing F given a vibration field X. In other words, it defines what it is known about the excitation field F after making vibration measurements;

• p(X|F) is the likelihood function, representing the probability of measuring X given an excitation field F. It reflects the uncertainty related to the measurement of the vibration field X;

• p(F) is the prior probability distribution, representing our knowledge of the unknown excitation field F before measuring the vibration field X.

Generally, the quality of the force reconstruction strongly depends on the choice of the likelihood function and the prior probability distribution. That is why, the choices made in this paper and in Ref. [START_REF] Aucejo | Bayesian source identification using local priors[END_REF] have to be carefully explained.

Choice of the likelihood function

As mentioned previously, the likelihood function reflects the uncertainty related to the measurement of the vibrations field X. By definition, this uncertainty is mainly related to the noise vector N. Consequently, the likelihood function p(X|F) can be rewritten under the following form:

p(X|F) = p(X -HF|N), (4)
representing the probability of obtaining X-HF = 0 given a noise vector N.

In other words the likelihood function measures the fidelity of the dynamic model to the measured data. Accordingly, the likelihood function can be also written as

p(X|F) = p(N), (5)
where p(N) is the prior probability distribution of the noise vector N, i.e.

the a priori of the experimenter on the nature of the noise corrupting the data.

If the noise is supposed spatially white and due to multiple independent causes, then the prior probability distribution of the noise vector N can be represented by a complex multivariate normal distribution with zero mean and precision parameter τ n :

p(N|τ n ) = τ n π N exp -τ n N 2 2 , (6)
where • 2 is the 2 -norm and N is the number of measurement points.

From Eq. ( 2), it can be inferred that the likelihood function is written as:

p(X|F, τ n ) = τ n π N exp -τ n X -HF 2 2 . (7)
Finally, it is worth mentioning that in practical situations the noise may not be strictly white or even Gaussian. However, assuming a spatially Gaussian white noise allows limiting the number of hyperparameters of the Bayesian model and obtaining satisfying reconstructions as shown in sections 4 and 5.

Choice of the prior probability distribution

The prior probability distribution reflects the uncertainty related to the unknown excitation field F. Actually, it can be seen as a measure of the a priori knowledge of the experimenter on the sources to identify. If one supposes that the structure is excited in R different regions by uncorrelated excitations of various types (localized or distributed), then local excitation fields F r can be considered as independent random vectors. As a result, the prior probability distribution can be written as the product of local prior probability distributions p(F r ), that is:

p(F) = R r=1 p(F r ), (8)
where p(F r ) reflects the prior knowledge of the nature of the sources in the region r.

Furthermore, it is assumed, for practical reasons, that the local excitation vectors F r are real random vectors, whose components are supposed a priori to be independent and identically distributed random variables following a generalized Gaussian distribution with zero mean [START_REF] Nadarajah | A generalized normal distribution[END_REF]. As a consequence, each local excitation field follows a multivariate generalized Gaussian distribution with zero mean. From the mathematical standpoint, the local prior probability distributions are thus written:

p(F r |τ sr , q r ) = q r 2 Γ(1/q r ) Mr τ Mr qr sr exp -τ sr F r qr qr , (9)
where:

• q r is the shape parameter of the distribution in the region r. Its value is defined in the interval ]0, +∞[; • • qr is the q r -norm or quasi-norm, if q r ≥ 1 and q r < 1 respectively;

• τ sr is the scale parameter of the distribution, which can be viewed as a generalized measure of the precision of the distribution;

• M r is the number of reconstruction points in the region r;

• Γ(x) = +∞ 0 t x-1 e -t dt is the gamma function.

It should be noted that the choice of a multivariate generalized Gaussian distribution allows a great flexibility for describing prior knowledge of the sources to identify. However, the hypothesis of independence used to derive the prior probability distribution can be discussed, insofar as one can argue that assuming a priori a small spatial correlation can sometimes improved the results. Actually, this hypothesis has been introduced here because the spatial correlation between each component of the force vector is difficult to assess a priori. Nonetheless, as shown in Ref. [START_REF] Aucejo | Bayesian source identification using local priors[END_REF] and in the next of the paper, the a priori hypothesis of independence allows obtaining consistent reconstructions, especially when the sources are spatially independent or functionally dependent. At the present time, the reconstruction of spatially correlated sources has not been tested yet.

Summary and comments

From the explanations given above, the standard Bayesian formulation of the reconstruction problem finally writes:

p(F|X, τ n , τ sr , q r ) ∝ p(X|F, τ n ) R r=1 p(F r |τ sr , q r ). ( 10
)
This formulation is said standard, because it gives rise to additive regularizations (Tikhonov when q r = 2 and Lasso when q r = 1), which corresponds to the MAP estimate of Eq. [START_REF] Brooks | Handbook of Markov Chain Monte Carlo[END_REF]. Furthermore, samples can be drawn [see section 3.3] to explore the posterior probability distribution and determine its statistical properties such as the mean, the median, the mode and credible intervals. In doing so, one has information on the credibility of the identified excitation given the measured vibration field X, the precision parameters τ n and τ sr and the shape parameters q r . Actually, this is the main drawback of the standard formulation, because the quality of the inference is conditioned to the knowledge of the shape and precision parameters. If their values are poorly chosen, then the resulting inference won't be representative of the actual target distribution. As a consequence, it is compulsory to determine near-optimal values of q r , τ n and τ sr if one wants to perform a relevant statistical inference.

To properly choose reasonable values of the shape parameters q r in each region r, it can be noticed from Fig. 1b that the smaller q r is, the larger is the weight on small values of F r . On the contrary, for large values of F r , the smaller q r is, the smaller is the weight on these values. This difference in the weighting behavior for small and large values of F r can be directly observed in the solutions obtained, since for q r ≤ 1 the solution vector F r will tend to have only a few non-zero values, while for q r = 2 it will tend to have only a few very small values. From this analysis, it can be inferred that distributed sources are promoted for q r = 2 [START_REF] Boyd | Convex optimization[END_REF], while localized sources are favored for q r ≤ 1 [START_REF] Grasmair | Non-convex sparse regularization[END_REF]. Consequently, the possible values of q r are bounded and practically lies in the interval ]0, 2] for force reconstruction problems.

However, it can be argued that in real-life applications the proposed guidelines seems impractical, because the force distribution is unknown or there is no idea about the force to recover. We are prone to think that in real-life applications rough information about the force to recover is available. This information is provided by a careful analysis of the mechanical system. As a result, we think that a rough idea of the location and the type of the sources acting on the structure can be obtained. This obviously supposed that the user has some knowledge in mechanics.

The choice of τ n and τ sr is generally less obvious. To this end, the more straightforward approach is first to determine the optimal values of the regularization parameters λ r = τ sr /τ n from automatic selection procedure based, for instance, on a generalization of the L-curve framework [START_REF] Belge | Efficient determination of multiple regularization parameter in a generalized L-curve framework[END_REF] or the marginalized MAP approach [START_REF] Pereira | Empirical bayesian regularization of the inverse acoustic problem[END_REF]. In a second step, a MAP estimate and (b) Corresponding penalty function |x| qr for (-) q r = 2, (--) q r = 1 and (-• -)

q r = 0.5
of the precision parameters τ sr can be obtained from the data only using the approach proposed in Ref. [START_REF] Pereira | Empirical bayesian regularization of the inverse acoustic problem[END_REF]. Finally, from the two first steps, the value of τ n can be obtained since τ n = τ sr /λ r .

When following the previous guidelines, we can reasonably expect that the standard Bayesian formulation allows performing an acceptable inference.

However, the standard formulation is limited in some respects, essentially because the parameters q r , τ n and τ sr are fixed. As a result, a poor choice of these parameters leads to a poor reconstruction. Nevertheless, this formulation will be used in section 3 to initialize the chains of the proposed MCMC algorithm in a reasonably high propability region of the target posterior distribution thanks to the prior estimation of q r , τ n and τ sr .

Extended Bayesian formulation

In the previous section, it has been shown that the inference ability of the standard formulation is clearly limited, because the parameters q r , τ n and τ sr need to be determined beforehand. To alleviate this limitation, the approach generally adopted in the literature consists in considering the precision parameters τ n and τ sr as random variables, while letting fixed the shape parameters q r . If we further consider the precision parameters as independent variables, the following extended formulation is obtained:

p(F, τ n , τ sr |X, q r ) ∝ p(X|F, τ n ) p(τ n ) R r=1 p(F r |τ sr , q r ) p(τ sr ), (11)
where p(τ n ) and p(τ sr ) are the prior probability distribution of the precision parameters τ n and τ sr respectively.

Choice of the prior probability distribution of the precision parameters

τ n and τ sr

The choice of the priori probability distributions p(τ n ) and p(τ sr ) is first limited to distribution having a strictly positive support, because the precision parameters τ n and τ sr are real positive numbers. The common choice, made in the literature, is the Gamma distribution. The reason for this is rather clear, since the conjugate prior for the precision of a generalized Gaussian distribution is a Gamma distribution [START_REF] Bishop | Pattern recognition and machine learning[END_REF]. The main interest of conjugate priors is to simplify the inference because it corresponds in general to classical probability distribution for which random number generators have already been implemented. This makes efficient the implementation of MCMC algorithms like the Gibbs sampler.

Practically, the Gamma distibution is defined by:

G(x|α, β) = β α Γ(α) x α-1 exp(-β x) with α > 0, β > 0, (12)
where α and β are respectively the scale parameter and the rate parameter of the distribution. The knowledge of these hyperparameters allows characterizing the probability distribution, since:

E[x] = α β and Var[x] = α β 2 . ( 13
)

Summary and comments

Considering the previous choice, the extended Bayesian formulation of the reconstruction problem is given by:

p(F, τ n , τ sr |X, q r ) ∝ p(X|F, τ n ) p(τ n |α n , β n ) R r=1 p(F r |τ sr , q r ) p(τ sr |α sr , β sr ), (14)
where (α n , β n ) are the hyperparameters related to the precision τ n , while (α sr , β sr ) are the hyperparameters related to the precision parameters τ sr .

To the authors' knowledge, this extended formulation is generally restricted to only one region and to Gaussian priors (i.e. q r = 2) [START_REF] Bardsley | Mcmc-based image reconstruction with uncertainty quantification[END_REF] and has given rise to the augmented Tikhonov regularization [START_REF] Bin | Augmented Tikhonov regularization[END_REF]. This method provides a point estimate corresponding to a critical point of the opposite of the logarithm of the posterior probability distribution. The main advantage of this approach is to determine the regularized solution and the precision simultaneously using an iterative process.

On the other hand, the use of a Gamma distribution is questionable, since it has been chosen for mathematical convenience and does not reflect any real prior information on the precision parameters, except their positiveness.

That is why, the prior distribution on τ n and τ sr should be as minimally informative as possible [START_REF] Hoff | A First Course in Bayesian Statistical Methods[END_REF]. To this end, one generally sets α n = α sr = 1

and

β n = β sr → 0.
Finally, the quality of the inference is strongly dependent on a proper choice of the shape parameters q r . Consequently, when setting q r = 2, as classically done in the literature to perform the inference using a Gibbs sampler [START_REF] Zhang | Bayesian force reconstruction with an uncertain model[END_REF][START_REF] Faure | Vibration source characterization using force analysis technique and a Bayesian regularization[END_REF],

one takes the risk to draw erroneous conclusions if the structure is actually excited by localized sources.

Complete Bayesian formulation

The derivation of the standard and extended formulations points out the need for properly defining the value of the shape parameters q r . However, choosing a priori relevant value is far from an easy task for non-experienced user. For this particular reason, it is interesting to infer the shape parameters from a Bayesian analysis. Practically, this is done by considering these parameters as independent and identically distributed variables random variables. In doing so, one obtains the complete Bayesian formulation:

p(F, τ n , τ sr , q r |X) ∝ p(X|F, τ n ) p(τ n |α n , β n ) R r=1 p(F r |τ sr , q r ) p(τ sr |α sr , β sr ) p(q r ), (15)
where p(q r ) is the prior probability distribution of the shape parameters q r .

Choice of the prior probability distribution of the shape parameters q r

The only available information is that the value of the shape parameters is bounded and positive. In absence of more precise knowledge on the parameters, the probability distribution is not only chosen to reflect the available information but also for its mathematical tractability. A probability distribution that meets these requirements is the truncated Gamma distribution defined by:

G T (x|α, β, l b , u b ) = Γ(α) γ(α, βu b ) -γ(α, βl b ) G(x|α, β) I [l b ,u b ] (x), (16)
where:

• G(x|α, β) is the Gamma distribution defined in Eq. ( 12);

• I [l b ,u b ] (x)
is the truncation function defined between the lower bound l b and the upper bound u b . More precisely, this function simply writes:

I [l b ,u b ] (x) =      1 if x ∈ [l b , u b ] 0 otherwise ; (17)
• γ(s, x) = x 0 t s-1 exp(-t) dt is the lower incomplete Gamma function.

Summary and comments

From the above considerations, one obtains the proposed complete Bayesian formulation, namely:

p(F, τ n , τ sr , q r |X) ∝ p(X|F, τ n ) p(τ n |α n , β n ) × R r=1 p(F r |τ sr , q r ) p(τ sr |α sr , β sr ) p(q r |α r , β r , l b , u b ), (18)
where (α r , β r ) are the hyperparameters related to the shape parameters q r .

As evoked previously, the choice of the truncated Gamma distribution has been made for mathematical convenience, because other continuous truncated distribution could have theoretically been used. To avoid biasing the inference, the shape of the prior distribution needs to be weakly informative. Here, this means that the hyperparameters should be defined such that α r = 1 and β r → 0.

At last, even if the lower and upper bounds l b and u b can theoretically take any positive value, one knows from section 2.1.3 that the value of q r practically lies in the interval ]0, 2]. For this particular reason, we set l b = 0.05 and u b = 2.05 for all the validations presented in the next of the paper.

Resolution of the reconstruction problem -MCMC algorithm

Classically, the Bayesian framework is adopted because it allows quantifying the uncertainties about the unknown parameters. To some extent, this problem consists in determining the shape of the posterior probability distribution p(F, τ n , τ sr , q r |X) in high probability regions, since it integrates all the available information about the parameters to identify. In the present situation, this task is analytically intractable, because of the dimension of the parameter space. An alternative solution is to discretize the parameter space to explore the posterior probability distribution. An efficient way to perform this exploration consists in sampling from the posterior distribution using MCMC algorithms.

In the present paper, the prior distributions have been chosen with conjugacy relations in mind. Here, this implies that most of the full conditional probability distributions have the same form as the corresponding prior distribution and can then be easily sampled using standard statistical packages.

For the complete Bayesian formulation, the full conditional probabilities are for:

• the shape parameters q r : p(q r |X, F, τ n , τ sr ) ∝ p(F r |τ sr , q r ) p(q r |α r , β r , l b , u b )

∝ τ Mr qr sr Γ(1/q r ) Mr q αr+Mr-1 r exp -β r q r -τ sr F r qr qr I [l b ,u b ] (q); (19)
• the precision parameters τ sr :

p(τ sr |X, F, τ n , q r ) ∝ p(F r |τ sr , q r ) p(τ sr |α sr , β sr ) ∝ τ αsr+ Mr qr -1 sr exp -τ sr β sr + F r qr qr ∝ G τ sr α sr + M r q r , β sr + F r qr qr ; (20)
• the precision parameters τ n :

p(τ n |X, F, τ sr , q r ) ∝ p(X|F, τ n ) p(τ n |α n , β n ) ∝ τ αn+N -1 n exp -τ n β n + X -HF 2 2 ∝ G τ n α n + N, β n + X -HF 2 2 ; (21)
• the force vector F:

p(F|X, τ n , τ sr , q r ) ∝ p(X|F, τ n ) R r=1 p(F r |τ sr , q r ) ∝ exp -τ n X -HF 2 2 - R r=1 τ sr F r qr qr , (22)
which corresponds to the posterior distribution of the standard Bayesian formulation introduced in section 2.1.

Regarding MCMC methods, the previous full conditional probability distributions can be efficiently used to sample the target posterior probability distributions thanks to the Gibbs sampler, which is a special case of the Metropolis-Hastings algorithm [START_REF] Metropolis | Equation of State Calculations by Fast Computing Machines[END_REF][START_REF] Hastings | Monte Carlo sampling methods using Markov chain and their applications[END_REF]. The idea of the Gibbs sampler is to draw samples form the posterior distribution by sampling each variable from its full conditional distribution given the current values of the remaining variables. Technically, the Gibbs sampler is divided into three main steps:

1. Set k = 0 and initialize q

(0) r , τ (0)
sr , τ

n and F (0) ; 2. Draw N s samples from full conditional distributions

for k = 1 : N s a. for each region r, draw q (k) r ∼ p q r |X, F (k-1) , τ (k-1) n , τ (k-1) sr ;
b. for each region r, draw τ

(k) sr ∼ p τ sr |X, F (k-1) , τ (k-1) n , q (k) r ; c. draw τ (k) n ∼ p τ n |X, F (k-1) , τ (k)
sr , q

(k) r ; d. draw F (k) ∼ p F|X, τ (k) n , τ (k)
sr , q (k) r

; end for

Monitor the convergence of the Markov chains

The implementation of the previous Gibbs sampler requires some comments. First of all, the parameters can be theoretically randomly initialized.

However, this potentially leads to a long burn-in period because the chains are slowly mixing, meaning that the chains slowly explore the support of the posterior distribution. To limit the burn-in period, a possible solution is to start at a point known to have a reasonably high probability, such as the mode [START_REF] Geyer | Handbook of Markov Chain Monte Carlo, chapter 1 -Introduction to Markov Chain Monte Carlo[END_REF]. This approach is adopted here and the calculation of the starting point is detailed in section 3.1.

Then, a particular attention has to be paid when drawing samples from full conditional distributions. While the samples for the precision parameters

τ (k) n and τ (k)
sr are easily obtained, drawing samples for the shape parameters q (k) k and the force vector F (k) requires special treatments, that are thoroughly explained in sections 3.2 and 3.3.

Finally, inferences performed from MCMC samples are based on the assumption that the chains have properly mixed, i.e. that the chains have converged to the stationary distribution. To this end, several methods have been proposed in the literature such as the Gelman-Rubin statistic based on multiple parallel chains [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF]. Here, all the inferences are performed from one single chain for each parameter. Hence, adapted procedures have to be used to assess the convergence of the Markov chains. The procedure used in the proposed validations to check the convergence of the Gibbs sampler is explained in section 3.4.

Initialization of the sampler

As explained previously, it has been chosen to initialize the sampler from a starting point having a reasonably high probability. In the present case, the initial force vector F (0) is obtained from the MAP estimate of the standard Bayesian formulation, which is defined by [see Ref. [START_REF] Aucejo | Bayesian source identification using local priors[END_REF] for more details]:

F (0) = argmax F p F|X, τ (0) n , τ (0) sr , q (0) r = argmin F X -HF 2 2 + R r=1 λ (0) r F r q (0) r q (0) r , (23)
where λ

(0) r = τ (0) sr /τ (0)
n is the regularization parameter in the zone r.

To simplify somewhat the calculation, it is supposed that the precision parameters τ sr are equal to a unique constant value τ

s . Accordingly, Eq. ( 23) becomes:

F (0) = argmin F X -HF 2 2 + λ (0) R r=1 F r q (0) r q (0) r , (24)
where

λ (0) = τ (0) s /τ (0) n .
To obtain a relevant initial force vector, it is necessary to determine reasonable values of the shape parameters q (0) r and the precision parameters τ (0) s and τ

n . Practically, the values of the shape parameters q (0) r can be chosen from subjective considerations using the guidelines given in section 2.1.3. On the contrary, near-optimal values of the precision parameters are difficult to assess a priori, i.e. without any calculation, because they are strongly related to the optimization problem given in Eq. ( 24) through the regularization parameter λ (0) , whose optimal value is partly conditioned to q (0) r . In addition, since the shape parameter q (0) r can take any value in the range ]0,2], the solution of the optimization problem has generally no closed-form expression.

Incidentally, the optimal value of the regularization parameter λ (0) can not be directly estimated from the marginalized MAP [START_REF] Pereira | Empirical bayesian regularization of the inverse acoustic problem[END_REF], the L-curve principle [START_REF] Hansen | Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion[END_REF] or the Generalized Cross Validation [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF] before solving the minimization problem [see Ref. [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF] for details]. To this end, the optimization problem is solved iteratively using an Iteratively Reweighted Least Squares (IRLS) algorithm [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF][START_REF] Rodriguez | An iteratively weighted norm algorithm for total variation regularization[END_REF], which allows determining in an iterative manner the regularized force vector as well as the optimal regularization parameter associated to the minimization problem.

The core idea of the IRLS algorithm consists in replacing the direct resolution of the minimization problem [see Eq. ( 24)] by an equivalent iterative process having an explicit solution at each iteration. Technically, the algorithm is based on the fact that:

∀ q (0) r , R r=1 F r q (0) r q (0) r = R r=1 W 1/2 r F r 2 2 = W 1/2 F 2 2 , ( 25
)
where W is a diagonal global weighting matrix depending explicitly on F r and q (0) r and defined from local weighting matrices W r such that [see Refs. [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF][START_REF] Aucejo | Bayesian source identification using local priors[END_REF][START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF][START_REF] Rodriguez | An iteratively weighted norm algorithm for total variation regularization[END_REF] for details]:

W = diag[W 1 , . . . , W r , . . . , W R ]. (26)
In doing so, one obtains, after convergence of the iterative process, the optimal force vector F (0) , the global weighting matrix W, as well as the optimal value of the regularization parameter λ (0) , that has been updated at each iteration using an automatic selection procedure such as the L-curve principle [START_REF] Hansen | Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion[END_REF] or the Bayesian estimator [START_REF] Antoni | A Bayesian approach to sound source reconstruction: Optimal basis, regularization, and focusing[END_REF].

From a Bayesian standpoint, the IRLS procedure can be viewed as a way of calculating iteratively the MAP estimate by approximating Generalized Gaussian priors by a series of multivariate Gaussian-like priors. As a result, at the final iteration of the IRLS algorithm, Eq. ( 25) should hold, meaning that:

p(F|τ (0) s , q (0) r ) = R r=1 p(F r |τ (0) s , q (0) r ) ∝ p(F|W, τ (0) s ) ∝ exp -τ (0) s W 1/2 F 2 2 . ( 27
)
At this stage of the initialization procedure, it remains to determine the values of τ (0) n and τ

s . For this purpose, we follow the approach proposed by Pereira et al. [START_REF] Pereira | Empirical bayesian regularization of the inverse acoustic problem[END_REF] consisting in finding the optimal values of τ . Formally, this condition writes:

τ s (0) , τ n (0) = argmax τ (0) s ,τ (0) n p(τ (0) n , τ (0) s |X). (28)
It follows that the optimal values of τ

n and τ (0) s are obtained after some simple calculations detailed in Appendix A and are expressed as:

τ (0) s = N X H λ (0) I + HW -1 H H -1 X and τ (0) n = τ (0) s λ (0) . (29)
3.2. Drawing samples from p q r |X,

F (k-1) r , τ (k-1) n τ (k-1) sr
The Gibbs sampler is generally the first choice for conditionally conjugate models, where samples can be drawn directly from each conditional probability distribution. Unfortunately, this is not the case for the shape parameters q r , for which the conditional distribution is not standard. In this situation, a usual practice is to sample this variable using a Metropolis-Hastings update. Such a combination of Gibbs sampler and Metropolis-Hastings update is generally called Metropolis-within-Gibbs algorithm [START_REF] Tierney | Markov Chains for Exploring Posterior Distributions[END_REF].

The main drawback of the Metropolis-Hastings algorithm is its random walk behavior when using a normal proposal distribution (or any random walk proposal). Although its popularity, such proposal distributions can lead to a somewhat inefficient exploration of the probability space implying slow mixing of the chains. To avoid such a behavior, a Hamiltonian Monte Carlo update is used instead.

Hamiltonian Monte Carlo (HMC), a.k.a. Hybrid Monte Carlo, is a MCMC algorithm based on Hamiltonian dynamics [START_REF] Duane | Hybrid monte carlo[END_REF]. In HMC, the Hamiltonian function can be written as follows:

H(q r , s) = U (q r ) + K(s), (30)
where U (q r ) is called the potential energy and is defined as:

U (q r ) = -log p q r |X, F (k-1) r , τ (k-1) n , τ (k-1) sr , (31)
while K(s) is called the kinetic energy and is chosen for mathematical convenience. Usually, it is defined such that:

K(s) = -log [N (s|0, 1)] = 1 2 s T s + const, (32)
where N (s|0, 1) is the standard univariate normal distribution.

In this framework, q r and s define the state space of the dynamical systems and can be viewed as the position and the momentum respectively. The evolution of the state (q r , s) over the time t is governed by the following system of equations:

     dq r dt = ∂H(q r , s) ∂s ds dt = - ∂H(q r , s) ∂q r . (33)
In practice, Eq. ( 33) is difficult to solve analytically. That is why, Hamilton's equations are approximated by discretizing time and solved using the leapfrog method, that requires the computation of the gradient of the potential energy U (q r ) [see Appendix B]. The latter aspect is, actually what partly avoids the random-walk behavior of the sampler, because the gradient information provides a direction of exploration. Practically, a constrained version of the HMC sampler has been implemented to take the bounded character of the shape parameters q r into account [see Ref. [START_REF] Neal | Hanbook of Markov Chain Monte Carlo, chapter 5 -MCMC Using Hamiltonian Dynamics[END_REF] for details].

Drawing samples from p F|X, τ

sr , q

(k) r
The full conditional distribution p F|X, τ

(k) s , τ (k)
sr , q

(k) r is not standard.

Consequently, a natural approach for drawing samples from this probability distribution consists in performing a Metropolis-Hastings or a HMC update. This updating stage can be implemented component-wise or block-wise.

However, this can lead to convergence and computational efficiency issues, because of a low acceptance rate, which is likely to happen when the number of identification points is large. Another possibility, adopted in this work, is to approximate the full conditional probability distribution by noting that, for an appropriate choice or estimation of the local weighting matrices W r , one can write:

p F|X, τ (k) n , τ (k) sr , q (k) r ∝ exp -τ (k) n X -HF 2 2 - R r=1 τ (k) sr F r q (k) r q (k) r ∝ exp -τ (k) n X -HF 2 2 - R r=1 τ (k) sr W 1/2 r F r 2 2 ∝ N c (F|µ F , Σ F ), (34)
where N c (F|µ F , Σ F ) is the (circurlarly-symmetric) complex multivariate Gaussian distribution with complex mean vector µ F and complex covariance ma-trix Σ F . Formally, the multivariate Gaussian distribution is expressed as:

N c (F|µ F , Σ F ) = 1 |πΣ F | exp -(F -µ F ) H Σ -1 F (F -µ F ) , (35)
where

|Σ F | is the determinant of Σ F .
Here, the complex mean vector µ F and the complex covariance matrix Σ F are expressed as:

µ F = τ (k) n Σ F H H X and Σ F = τ (k) n H H H + T (k) s W -1 . (36)
In the previous equation, W is the global weighting matrix defined form the local weighting matrices W r after Eq. ( 26), while T (k) s is the global precision matrix related to the local precision parameters τ (k) sr . By construction, T (k) s is defined such that:

T (k) s = diag τ (k) s1 I M 1 , . . . , τ (k) sr I Mr , . . . , τ (k) sR I M R , (37)
where I Mr is the identity matrix of dimension M r .

Practically, µ F is computed from the IRLS algorithm described in Refs. [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF][START_REF] Aucejo | Bayesian source identification using local priors[END_REF], since the IRLS solution converges to the conditional mean (i.e. µ F ). As explained earlier, at the last iteration of the IRLS algorithm, one obtains the complex mean vector µ F , as well as the related local weighting matrices W r enabling the calculation of the complex covariance matrix Σ F .

Once the complex mean vector µ F and covariance matrix Σ F of the complex multivariate Gaussian distribution are obtained, samples of the force vector F (k) can be easily drawn [START_REF] Gentle | Computational Statistics[END_REF]. At this stage, it must be mentioned that, from a technical standpoint, the proposed sampling procedure of the force vector can be seen as a normal approximation around the mode of the full conditional probability distribution. Regarding the overall sampling procedure, this results in an approximate Bayesian inference when the shape parameters are different from 2. However, in the light of the results presented in the next of the paper on plate-like structures and in Ref. [START_REF] Aucejo | Bayesian formulations for force reconstruction problems[END_REF] on a beam, the proposed approximation seems to lead to relevant Bayesian inferences.

Convergence diagnostics

From Markov chain theory, one can expect the Markov chains to converge to their equilibrium distribution, corresponding to the target posterior probability distribution. However, there is no guarantee that the convergence can be achieved. Worse, one can never be sure that the chains have properly mixed. Nevertheless, several tests can be performed to determine whether the chains appear to converge.

Convergence diagnostics remains an open topic, but there is essentially two paradigms for establishing convergence tests. The first one consists in considering one long chain, while the second one is based on the exploitation of multiple chains running in parallel. Each standpoint has its pros and cons. In the present paper, all the inferences have been performed from one single chain for each parameter. More precisely, we have monitored the convergence of the precision parameters (τ n , τ sr ) and the shape parameters q r .

The procedure used for convergence monitoring is based on the combination of two complementary diagnostics. The proposed procedure is conse-quently divided into three steps:

1. Apply the Raftery-Lewis diagnostic [START_REF] Raftery | One long run with diagnostics: Implementation strategies for Markov chain Monte Carlo[END_REF] to each individual chain to monitor for estimating the total run length and the burn-in period.

Apply the Geweke diagnostic [41] to the chains resulting from step 1

to test the convergence of the Markov chains.

3. If the convergence is not achieved after step 2, repeat steps 1 and 2 on the remaining parts of the chains.

General comment

The proposed MCMC algorithm is relatively general and can be used to derive the MCMC algorithms corresponding to the standard and the extended Bayesian formulations. The MCMC algorithm related to the extended Bayesian formulation presented in section 2.2 is obtained by not considering the step 2a of the Gibbs sampler (i.e. by fixing the values of q r ). On the other hand, the MCMC algorithm associated to the standard Bayesian formulation is derived from the general Gibbs sampler by computing steps 1 and 2c only (i.e. by considering fixed values for τ n , τ sr and q r ).

Numerical validation

The present numerical study intends to demonstrate how the proposed full Bayesian approach operates to obtain consistent reconstructions as well as relevant and useful information on the posterior uncertainties on the model parameters. This last aspect is certainly the most important from an industrial point of view, because having credible intervals allows assessing not only the quality of identified solutions but also the pertinence of the proposed Bayesian model. To properly simulate experimental measurements, the exact vibration displacement field X exact is first computed from a FE mesh of the plate made up with 187 shell elements, assuming that only bending motions are measurable. Then, the exact displacement field is corrupted by an additive Gaussian white noise with a signal-to-noise ratio equal to 34 dB. It should be added that a structural damping has been introduced in the calculation to avoid infinite displacement amplitudes at resonance frequencies. Here, the structural damping ratio is equal to 0.01.

Regarding, finally, the transfer functions matrix H, a FE model of the plate with free boundary conditions is used, assuming that only bending motions are measured. In other words, the computed transfer functions matrix H is dynamically condensed over the measurable degrees of freedom only [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF][START_REF] Ibrahim | Structural force identification at unmeasured locations[END_REF]. The main interest in using free boundary conditions to model the dynamic behavior of the plate is to allow the identification of external excitations acting on the structure as well as reaction forces at boundaries [START_REF] Renzi | Vibratory source identification by using the finite element model of a subdomain of a flexural beam[END_REF].

Application

To numerically validate any force reconstruction strategy, it is first necessary to define the reference force vector F ref that could serve as a proper benchmark. This reference force vector is computed from the transfer functions matrix H and the exact displacement field X exact thanks to the following relation:

F ref = H -1 X exact . (38)
This numerical study is focused on the identification of the excitation field at 350 Hz, i.e. outside the resonance frequencies of the plate. As shown in Fig. 2, the reference force vector corresponds to the description of the test case given in the previous section, since it exhibits smooth reaction forces at boundaries of the plate as well as a unit point force F 0 at (x 0 , y 0 ) = (0.42 m, 0.25 m).

Before applying the complete Bayesian formulation to the proposed numerical test case, it is important to stress the interest of defining several identification region and inferring the shape parameters through a Bayesian analysis. To this end, the extended Bayesian formulation presented, in section 2.2, is first applied by considering only one identification region. In this situation, a single shape parameter q needs to be defined. In this situation, when setting q to typical values, namely q = 2 and q = 1, one clearly observes from the analysis of Fig. 3 and Tables 1 and2, obtained after initially drawing 7000 samples, that defining only one identification region does not On the contrary, when defining several identification region, the quality of the inference can be significantly improved. For the present test case, two reconstruction regions plotted in Fig. 4 are defined: (i) a central region associated to the shape parameter q 1 and containing the point force only and

(ii) a region associated to the shape parameter q 2 and corresponding to the boundaries of the plate. In this situation, when applying the extended Bayesian formulation for q 1 = 0.5 and q 2 = 2 [see guidelines proposed in section 2.1.3], one obtains, after drawing initially 7000 samples, the results presented in Figure 5 and Table 3. Here, it is clear that the definition of different identification regions associated to a relevant choice of the corresponding shape parameters allows greatly improving the quality of the resulting inference. This is actually the major shortcoming of the extended Bayesian formulation, since only partial information on the excitation field are practically available. More specifically, this means that the value of the precise value of the shape parameters are not known precisely, even if the guidelines given previously allow determining reasonable values. Actually, this is this conclusion that motivated the development of the complete Bayesian formulation. To perform the inference from the complete Bayesian formulation, the initial values of the shape parameters have been set to q (0) 1 = 0.5 and q (0) 2 = 2 in order to start the sampler in a reasonably high posterior probability region.

After initially drawing 7000 samples for each chain, the results of the Bayesian inference are summarized in Table 4 and Figs. 6 and7.

Obtained results clearly show that the proposed inference provides consistent parameters estimates as well as a quantification of the posterior un- certainty about those parameters after the numerical experiment was carried out. In particular, it can be stressed that the median of the reconstructed force vector agrees well with the reference one. On the other hand, the analysis of the inference results indicates that the width of the related 95% credible intervals on the estimated excitation field are larger at boundaries than in the central region, indicating that our knowledge of the reaction forces is more uncertain than that of the point force. This observation is confirmed by the comparison of the precision τ s1 and τ s2 , insofar as the median and the mode of τ s1 are greater than those of τ s2 . Furthermore, as a verification of the proposed statistical model, it is noticed that the true noise precision τ n = 1.433 × 10 16 lies in the 95% credible interval [see Table 4]. Finally, it is also worth noting that the estimations of the shape parameters q 1 and q 2 are in a very good agreement with the general guidelines proposed in section 2.1.3. Consequently, the proposed Bayesian model seems suitable for force reconstruction problems.

The previous application of the complete Bayesian formulation has been carried out from a particular set of initial shape parameters (q

(0) 1 , q (0)
2 ) = (0.5, 2) because it allows starting the sampler at a point having a reasonably high probability, i.e. at a point close to the mode of the distribution.

However, as explained earlier, it could be difficult to assess in a real-life application the precise values of these parameters. In practical situation and after a thorough analysis of the mechanical problem, the guidelines provided previously allow having only a rough idea of the values of the shape parameters. As a consequence, if the proposed sampler is sensitive to a correct initialization of the shape parameters, the proposed approach would be seriously invalidated. To have better insights into the sampler behavior, the trace plots of the Markov chains of each parameters of the complete Bayesian formulation are presented in Fig. 8 for different initial values of the shape parameters q (0) 1 and q (0) 2 . The analysis of the trace plots shows that, regardless of the initial values of the shape parameters, the Markov chains converge to the same stationary state. Consequently, a poor prior estimation of the shape parameters seems not detrimental to the sampler convergence, but it generally prolongs the burn-in period. That is why, it is beneficial to start in a high probability region of the posterior distribution.

Finally, the analysis presented above has been performed outside the resonance frequencies of the plate. Indeed, force reconstruction at natural frequencies of a lightly damped structure is a far more challenging task (especially at low frequencies), because, H is close to be singular and the vibration response is mainly driven by one structural mode (reaction forces) and not , (e) q 1 and (f) q 2 for different initial values of shape parameters (q

(0) 1 , q (0)
2 ) -(-) (q

(0) 1 , q (0)
2 ) = (0.5, 2) , (--) (q

(0) 1 , q (0)
2 ) = (2, 0.5) and (-• -) (q (0) 1 , q (0)

2 ) = (1, 1). The burn-in period is included in the trace plots 39 by the point force itself. At these particular frequencies, a special treatment or formulation should be implemented, but, to the author's knowledge, this problem still remains an open question.

Experimental validation

This section aims at confirming the main conclusions drawn in the previous section by extending the analysis to a real-world application. Since the proposed method is a natural extension of the work presented in Ref. [START_REF] Aucejo | Bayesian source identification using local priors[END_REF], we have decided to perform this validation using the same experimental set-up.

Description of the experimental set-up

The structure under test is a steel parallelepiped box, excited on one of its faces by a shaker fed by a white noise signal and equipped with a force sensor [see Fig. 9]. The parameters this experimental validation are given in Table 5.

Application

A careful analysis of the experimental set-up suggests the definition of two identification regions. The first region associated to the shape parameter q 1 contains the point force only, while the second region, associated to the shape parameter q 2 , corresponds to the boundaries of the plate [see Fig. 11].

According to section 2.1.3, the Gibbs sampler is thus initialized using q (0) 1 = 0.2 and q (0) 2 = 2. As in the numerical validation, the Bayesian inference is first applied at 525 Hz, which lies outside the resonance frequencies of the plate. After initially drawing 7000 samples for each chains, obtained results are presented in Figs. 12, 13 and 14 and are summarized in Table 6. As observed in the numerical validation, the proposed Bayesian formula- terval indicates that given the data and the model, there is 95% chance (or probability) that the value of interest lie in that interval. The width of the CI is a way of measuring the posterior uncertainty about the parameters estimates, or, equivalently, our current knowledge of the parameters after the experience was performed. Accordingly, for the identified point force, it does not mean that the reconstructed value is necessarily close to 1, but, that it is credible to think that it lies in the calculated interval given all the information available a priori (dynamic model + specific contextual information from prior distribution) and/or measured (data). Fortunately, the estimate obtained in Fig. 15a shows that the proposed Bayesian model is relevant for coping with force reconstruction problems outside resonance frequencies.

Around resonance frequencies, the inference leads to the solution given the larger posterior uncertainty, as suggested by the median value of the shape parameter q 1 which tends to 2. Even if this result is not the best one, it is, after all, reassuring, because at resonance frequencies the response of the structure is governed by the vibration modes, implying that the identified excitation field is necessarily not unique.

Conclusion

In structural dynamics, regularization approaches, such as additive regularizations, can be used to deal with source identification problems. Although widely used and deeply studied, only point estimates are generally computed, even if the variance information can be easily obtained in some particular cases. In other words, there is often no quantification of the uncertainty about the regularized solution given the measured data and the mechanical model. However, such information is of primary interest for industrial applications, in which it is essential to guarantee the quality of obtained results. In the present paper, a full Bayesian formulation has been proposed to tackle this issue and an original hybrid Gibbs sampler has been implemented

to perform the inference. One of its merits is to provide estimates of all the parameters of the model, including the shape parameters of the generalized Gaussian used to express our prior knowledge on the sources to identify. The proposed numerical and experimental validations clearly highlight the practical interest of the proposed approach in terms of parameters estimations and posterior uncertainty quantification.

Since the potential energy is defined up to a constant, one obtains:

U (q r ) = M r log Γ(1/q r ) -M r q r log τ sr -[α r + M r -1] log q r + β r q r + τ sr F r qr qr .

(B.3)

From the previous relation, the gradient is easily obtained:

∂U (q r ) ∂q r = β r +τ sr Mr i=1 (|F ri | qr log |F ri |)- α r + M r -1 q r + 1 q 2 r
[M r log τ sr -M r ψ(1/q r )] ,

(B.4)
where F ri is the i th component of the force vector F r and ψ(x) is the digamma function.

Figure 1 :

 1 Figure 1: Selection of the shape parameter q r -(a) Probability distribution exp [-|x| qr ]

  vibration field X by calculating the MAP estimate of

4. 1 .

 1 Description of the test case In the present numerical validation, one seeks to identify a point force of unit amplitude (i.e. F 0 = 1 N) acting on a thin simply supported steel plate with dimensions 0.6 m × 0.4 m × 0.005 m. The coordinates of the point force, measured from the lower left corner of the plate, are (x 0 , y 0 ) = (0.42 m, 0.25 m). Practically, this configuration allows studying the influence of the definition of local regularization terms, since the present excitation field exhibits two types of spatial distribution over the structure, namely a smooth distribution of the reaction forces at boundaries and a singular distribution around the location of the point force.

Figure 2 :

 2 Figure 2: Reference force vector F ref at 350 Hz

Figure 3 :

 3 Figure 3: Surface plot of the real part of the reconstructed force vector (2.5 th , 50 th and 97.5 th percentiles) at 350 Hz obtained after applying the extended Bayesian formulation for (a) q = 2 and (b) q = 1

Figure 4 :

 4 Figure 4: Definition of the reconstruction regions -(•) region 1 (Point force), (×) region 2 (Reaction forces) and (♦) location of the point force

Figure 5 :

 5 Figure5: Surface plot of the real part of the reconstructed force vector (2.5 th , 50 th and 97.5 th percentiles) at 350 Hz obtained after applying the extended Bayesian formulation for (q 1 , q 2 ) = (0.5, 2)

Figure 6 :Figure 7 :

 67 Figure 6: Real part of the reconstructed force vector at 350 Hz. (a) Surface plot -2.5 th , 50 th (median) and 97.5 th percentiles and (b) Section view at y 0 = 0.25 m -(-) Reference, (--) Reconstruction and ( a ) 95% credible interval

Figure 8 :

 8 Figure 8: Trace plots of (a) F 0 , (b) τ n , (c) τ s1 , (d) τ s2 , (e) q 1 and (f) q 2 for different initial

Figure 9 : 1 ]Figure 10 :

 9110 Figure 9: Experimental set-up

Figure 11 :Figure 12 :Figure 13 :Figure 14 :

 11121314 Figure 11: Definition of the selected zones -(•) zone 1 (sparse), (×) zone 2 (smooth) and (♦) location of the point force

Fig. 15 .Figure 15 :

 1515 Fig. 15. The results presented in this figure have been computed between 300 Hz and 1200 Hz (frequency resolution: 5 Hz) after initially drawing 7000 samples for each parameters at each frequency.

Table 1 :

 1 Summary of inference results at 350 Hz obtained after applying the extended Bayesian formulation for q = 2

	Parameter	Median	Mode	95% CI
	F 0 [N]	0.349	0.349	[0.136, 0.561]
	τ n	1.62×10 16 1.61×10 16 [1.34, 1.96]×10 16
	τ s	27.47	26.94	[21.84, 33.54]

Table 2 :

 2 Summary of inference results obtained at 350 Hz after applying the extended

	Bayesian formulation for q = 1
	Parameter	Median	Mode	95% CI
	F 0 [N]	0.799	0.799	[0.711, 0.886]
	τ n	1.37×10 16 1.35×10 16 [1.17, 1.59]×10 16
	τ s	15.74	15.53	[13.66, 18.03]

Table 3 :

 3 Summary of inference results obtained at 350 Hz after applying the extended Bayesian formulation for (q 1 , q 2 ) = (0.5, 2)

	Parameter	Median	Mode	95% CI
	F 0 [N]	0.976	0.976	[0.957, 0.994]
	τ n	1.35×10 16 1.34×10 16 [1.15, 1.57]×10 16
	τ s1	38.27	38.05	[33.45, 43.17]
	τ s2	15.57	14.99	[11.51, 20.81]

Table 4 :

 4 Summary of inference results at 350 Hz

	Parameter	Median	Mode	95% CI
	F 0 [N]	0.998	0.998	[0.98, 1.016]
	τ n	1.34×10 16 1.33×10 16 [1.15, 1.57]×10 16
	τ s1	35.03	34.61	[27.97, 44.87]
	τ s2	6.28	6.18	[4.10, 9.24]
	q 1	0.460	0.459	[0.405, 0.538]
	q 2	1.919	2.012	[1.564, 2.044]

Table 6 :

 6 Summary of inference results at 525 Hz

	Parameter Median	Mode	95% CI
	F 0 [N]	0.999	0.999	[0.981, 1.016]
	τ n	9.04×10 8 9.01×10 8 [7.97, 10.19]×10 8
	τ s1	56.94	55.78	[49.14, 66.96]
	τ s2	10.61	6.11	[3.22, 42.72]
	q 1	0.403	0.401	[0.381, 0.429]
	q 2	1.935	2.007	[1.578, 2.046]

Appendix A. Calculation of the optimal value of τ (0) s and τ (0) n The optimal values of τ

s |X). Formally, this condition writes:

By using the Bayes' rule, the conditional probability distribution p(τ

is written as:

Assuming that p(τ

s ) ∝ 1, one has:

After some calculations, one establishes that [see Ref. [START_REF] Bishop | Pattern recognition and machine learning[END_REF] for details]:

Practically, it is generally easier to find the solution of the dual minimization problem:

The previous minimization problem is equivalent to the minimization of the functional:

where

This functional can be simplified by introducing in the previous relation the regularization parameter λ (0) defined such that:

In doing so, the functional to minimize becomes:

To find the optimal value of τ (0) s , it remains to find a stationary point of the functional J(τ (0) s , λ (0) ) with respect to τ (0) s . Mathematically, this is performed by applying the first-order optimality condition, namely:

After rearranging the previous equation, one obtains:

Finally, the optimal value of τ

n is directly obtained from Eq. (A.7), since:

λ (0) . (A.11)

Appendix B. Calculation of the potential energy and its gradient

The implementation of the Hamiltonian Monte Carlo requires the calculation of the potential energy U (q r ) and its gradient ∂U (qr) ∂qr . By definition, the potential energy is given by:

Practically, the conditional probability distribution p (q r |X, F r , τ n , τ sr ) is written (the superscripts (k -1) have been removed for the sake of clarity):

p (q r |X, F r , τ n , τ sr ) = Cτ

Mr qr s Γ(1/q r ) Mr q αr+Mr-1 r exp -β r q r -τ sr F r qr qr , (B.2

where C is a constant.