L Raillon
email: loic.raillon1@insa-lyon.fr

C Ghiaus

An Efficient Bayesian Experimental Calibration of Dynamic Thermal Models

Keywords: Bayesian calibration, model selection and validation, dynamic thermal models, real house experiment, Metropolis-Hastings algorithm, robust gradient and Hessian computation, change of variables, prior distribution selection, identifiability

An efficient Bayesian experimental calibration of dynamic thermal models

L. Raillon, Christian

Introduction

The existing methods for characterizing building energy performance and energy saving provided by retrofitting are not relevant [START_REF] Turner | Energy Performance of LEED ® for New Construction Buildings[END_REF][START_REF] De Wilde | The gap between predicted and measured energy performance of buildings: A framework for investigation[END_REF]. The energy performance estimation of buildings and associated systems must be independent of weather conditions and user behavior. From this assessment, the Efficiency Valuation Organization has developed the International Performance Measurement and Verification Protocol (IPMVP) (EVO 2014). The idea is to construct a physical model which characterizes the building intrinsic thermal dynamic and relate inputs to outputs measured on-site. Hence, the gap between the energy use given by the pre-retrofit and post-retrofit models represents the energy gained by the refurbishment.

Minimizing heat losses from buildings is the most obvious solution to reduce the heating and cooling demand but the efficiency and sustainability of the energy chain, from the production to the HVAC systems, must also be improved. Nowadays, the dominant paradigm is that the energy sources need to respond to all requests at any moment. The complexity of this strategy will be augmented with the increasing of the share of renewable energy sources in the energy mix. Therefore, supply and demand must become more flexible by using demand response mechanisms and energy storage (European Commission 2016). In order to adapt the demand to the production, the energy demand must be known. Physical models characterizing the thermal dynamic of buildings associated with model predictive control can be used to forecast the energy demand while maintaining indoor comfort (Hazyuk et al. 2012a, Hazyuk et al. 2012b[START_REF] Ghiaus | Calculation of optimal thermal load of intermittently heated buildings[END_REF].

Two important societal needs are identified, the estimation of building energy demand and the estimation of energy savings brought by energy conservation measures. These two societal needs have the same scientific deadlock: the experimental estimation of the physical parameters of the dynamic thermal behavior of buildings. Such models can be obtained considering the energy balance between buildings and their surroundings and energy balance in buildings can be modelled by using thermal networks [START_REF] Naveros | Modelling Heat Transfer for Energy Efficiency Assessment of Buildings: Identification of Physical Parameters[END_REF][START_REF] Ghiaus | Causality issue in the heat balance method for calculating the design heating and cooling load[END_REF]. Stochastic state-space models are obtained by first transforming thermal networks in deterministic state-space, and then noise terms are added to represent the deviations between the differential algebraic equations and the true variations of the states. Stochastic state-space models relate inputs to outputs, where the dynamic of the states is given by the parameters; hence by knowing inputs and parameters, the output of the system can be simulated. However, the direct problem requires to know the parameters, so the inverse problem of parameter estimation must be solved first. The interest in parameter estimation for dynamic thermal models is not new [START_REF] Nielsen | A dynamic test method for the thermal performance of small houses[END_REF]) and experiments of various scales have been used to test the validity of different approaches (Bloem 1994, Baker & van [START_REF] Baker | PASLINK and dynamic outdoor testing of building components[END_REF][START_REF] Jiménez | Bayesian Calibration Of Residential Building Clusters Using A Single Geometric Building Representation Department of Engineering[END_REF].

It is essential when making prediction or decision from an identified model to assess the uncertainties in the parameter estimates; it should be done by taking all the information available. From this perspective, Bayesian estimation gets more and more consideration, for instance, in quantification of energy saving from retrofit [START_REF] Heo | Calibration of building energy models for retrofit analysis under uncertainty[END_REF][START_REF] Heo | Scalable methodology for large scale building energy improvement: Relevance of calibration in model-based retrofit analysis[END_REF][START_REF] Tian | Identifying informative energy data in Bayesian calibration of building energy models[END_REF][START_REF] Li | Assessment of linear emulators in lightweight Bayesian calibration of dynamic building energy models for parameter estimation and performance prediction[END_REF], in calibration of energy models [START_REF] Chong | Uncertainty analysis and parameter estimation of HVAC systems in building energy models[END_REF][START_REF] Chong | A Comparison of MCMC Algorithms for the Bayesian Calibration of Building Energy Models for Building Simulation[END_REF][START_REF] Zayane | Identification d 'un modèle de comportement thermique de bâtiment à partir de sa courbe de charge[END_REF], in estimation of thermal characteristic of a wall [START_REF] Berger | Bayesian inference for estimating thermal properties of a historic building wall[END_REF] and in estimation of heating consumption (Kristensen et al. 2017). Bayesian methods, such as the Metropolis-Hastings (MH), are usually employed for low-dimensional problems because by increasing the number of parameters it becomes more and more difficult to tune properly the algorithms.

The paper compare the three phases of an experimental model identification (selection, calibration and validation) from a Bayesian and frequentist point of view. More precisely, the paper treats the problem of parameter estimation of dynamic thermal models when the model structure is known. First, the choice of a model structure is discussed and then an implementation of the second-order Metropolis-Hastings for linear and Gaussian state-space models is proposed. Tools and guidelines for tuning and diagnosing the algorithm are also presented.

Next, different criteria are presented to assess the performance of models and to guide the selection of a model structure in agreement with data. The whole procedure is tested on a real test case where the differences between the Bayesian and the frequentist approach are illustrated.

Twin houses experiment

Twin houses are a real outdoor experiment conducted by the Fraunhofer Institute near Munich during April and May 2014. It is an unoccupied single family house (twin house O5) with 100 m 2 ground floor, a cellar and an attic space; a full description of the house and the experiment is given by [START_REF] Strachan | Empirical Whole Model Validation Modelling Specification Validation of Building Energy Simulation Tools, Team Stan Development[END_REF]. The experiment was designed such that the south zone (Figure 1) has only two boundary conditions: the external temperature and the adjacent spaces (cellar, attic, north zone). The adjacent spaces were held at 22 °C with blinds closed to reduce the chance of overheating and the doors separating the north and south zone were sealed off. The electric heaters on the south zone were synchronized on a Randomly Ordered Logarithmic Binary Sequence (ROLBS) to maintain a similar temperature in the different rooms (800 W in the living room, 500 W in the south bedroom and 500 W in the bathroom). The ROLBS signal was designed for three reasons: maximize the temperature difference with the boundary conditions, excite the range of time constant from 1 hour to 90 hours and decorrelate the heating signal with the solar radiation. The heaters are lightweight, with a time response estimated around 1 or 2 minutes by the Fraunhofer Institute [START_REF] Strachan | Empirical Whole Model Validation Modelling Specification Validation of Building Energy Simulation Tools, Team Stan Development[END_REF] and a split coefficient between convective and radiative heat gains of 70 30 % ⁄ . Mechanical ventilation was set to supply a volume flow rate of 60 m 3 /h into the living room and extract 30 m 3 /h in the bathroom and the south bedroom.

Figure 1: Layout of the twin house O5

The stratification of the air in each room of the south zone is measured with temperature sensors at 10 cm, 110 cm and 170 cm from the ground. This level of accuracy is not required to fit a simple dynamic thermal model; therefore the temperature of the rooms are chosen as the average of the three heights and the south zone temperature 𝑇 𝑠 (°C) is a weighted average of the spaces (living room, south bedroom, corridor and bathroom) by their respective surfaces. The boundary temperature 𝑇 𝑛 (°C) is also a weighted average of the temperatures in the different spaces (kitchen, lobby and north bedroom). A weather station near the house provides the outside air temperature 𝑇 𝑜 (°C) and the global solar irradiance measured on a horizontal surface 𝑄 ̇𝑔ℎ (W/m 2 ). The heating in the south zone is done by three electric heaters; therefore 𝑄 ̇ℎ (W) represents the total heat input injected.

Buildings are dynamic systems which may be modelled by using thermal networks from where state-space models can be deduced [START_REF] Ghiaus | Causality issue in the heat balance method for calculating the design heating and cooling load[END_REF]. This procedure is used to find a suitable model of the twin house where the south zone (in green on Figure 1) is considered as the main thermal zone and the north zone as an adjacent space.

Experimental calibration process

The experimental calibration process is decomposed in 3 phases: selection, calibration and validation [START_REF] Ljung | System identification: theory for the user (second edition)[END_REF]. First, a set of likely model structures characterized by unknown physical parameter vector 𝛉 is chosen based on a-priori knowledge (physics, experiment details, etc). Then, the calibration assesses how these model structures relate to observed data and to physical considerations. The calibration consists of finding a set of parameters which best represents the input-output relationship of a model through observed data, 𝐮 1:𝑘 = [𝐮 1 , 𝐮 2 , … , 𝐮 𝑘 ] and 𝐲 1:𝑘 =

[𝐲 1 , 𝐲 2 , … , 𝐲 𝑘 ]. Finally, the performances of the calibrated models are evaluated to select the model which is the best suited for its intended use.

The experimental calibration process is usually treated either from a frequentist or a Bayesian point of view. In Bayesian estimation, the unknown parameters are treated as random variables with a certain prior distribution 𝑝(𝛉), which represents the prior belief before looking at the data [START_REF] Dahlin | Accelerating Monte Carlo methods for Bayesian inference in dynamical models Accelerating Monte Carlo methods for Bayesian inference in dynamical models[END_REF]. Then, all the information available in the data is summarized in the likelihood function 𝑝(𝐲 1:𝑘 |𝛉). The prior belief and the data information are combined in the Bayes' theorem to compute the posterior distribution:

𝑝(𝛉|𝐲 1:𝑘 ) = 𝑝(𝐲 1:𝑘 |𝛉)𝑝(𝛉) 𝑝(𝐲 1:𝑘 ) ∝ 𝑝(𝐲 1:𝑘 |𝛉)𝑝(𝛉) (1)
where 𝑝(𝐲 1:𝑘 ) is a normalization constant independent of the parameters.

The posterior distribution 𝑝(𝛉|𝐲 1:𝑘 ) contains all the statistical information about 𝛉; the most probable value of the posterior distribution gives the maximum a posteriori (MAP) estimate:

𝛉 𝐌𝐀𝐏 = argmax 𝛉 (𝑝(𝐲 1:𝑘 |𝛉)𝑝(𝛉)) (2)
If only the information in the data is considered, maximizing the likelihood function 𝑝(𝐲 1:𝑘 |𝛉) gives the maximum likelihood (ML) estimate:

𝛉 𝐌𝐋 = argmax 𝛉 (𝑝(𝐲 1:𝑘 |𝛉)) (3)
The ML estimate can be seen as a MAP estimate with uniform prior distribution, 𝑝(𝛉) ∝ 1 [START_REF] Sarkka | Bayesian Filtering and Smoothing[END_REF].

Two philosophies exist for computing these estimates and their uncertainties. The frequentist approach relies on the fact that, as the number of observations increases, the influence of the prior distribution becomes negligible compared to the likelihood and the posterior distribution can be approximated by a Gaussian distribution (Gelman, Carlin, et al. 2014). ML estimation is popular because it requires only point estimates of the posterior modes and their corresponding uncertainties are determined by asymptotic properties. ML estimates are usually found by optimization routines as in the CTSM-R package (CTSM-R Development Team 2015). This strategy has been proven to be efficient at numerous cases [START_REF] Naveros | Setting up and validating a complex model for a simple homogeneous wall[END_REF][START_REF] Himpe | Characterisation of the thernial performance of a test house based on dynamic measurements[END_REF][START_REF] Nespoli | Grey-Box System Identification of Building Thermal Dynamics Using only Smart Meter and Air Temperature Data[END_REF][START_REF] Andersen | Characterization of heat dynamics of an arctic low-energy house with floor heating[END_REF][START_REF] Bacher | Identifying suitable models for the heat dynamics of buildings[END_REF][START_REF] Váňa | Building semi-physical modeling : On selection of the model complexity ˇ a c a n[END_REF]. However, the asymptotic theory does not hold for small number of observations, which is often the case in real experiment. From a Bayesian point of view, all the statistical information is summarized in the posterior distribution, thus no assumption is made.

Bayesian and frequentist methods are compared in the three phases of the calibration process in order to illustrate the differences and to take advantage of both methods.

Model selection

The choice of an appropriate model structure is the most crucial part in experimental calibration according to [START_REF] Ljung | System identification: theory for the user (second edition)[END_REF]. The model structure should be representative of the real physical system and also in agreement with the measured data. Based on a-priori knowledge of the experiment, a set of models of increasing complexity is defined which allows a forward selection between nested models, i.e. the simplest model can be recovered by putting extra parameters to zero [START_REF] Bacher | Identifying suitable models for the heat dynamics of buildings[END_REF]. In this paper, the distinction between model selection and model comparison is done because the model comparison requires to discuss about calibration first; therefore it is presented later. To illustrate the Bayesian experimental calibration process only the two most promising model candidates are presented. These two models, illustrated in Figure 2, have been selected from a larger set of models which is not presented here due to lack of space. The smallest thermal network in Figure 2, the 3 states model, ℳ 3 , is obtained by not taking into account the surrounded dotted part whereas the 4 states model, ℳ 4 , is obtained by adding this part to ℳ 3 , such that, ℳ 3 ⊂ ℳ 4 . In Error! Reference source not found. the nodes represent:

• 𝑥 𝑤 : the temperature of the building envelope (°C),

• 𝑥 𝑖 : the indoor air temperature (°C),

• 𝑥 𝑚 : the temperature of the medium (internal walls, furniture, etc) (°C),

• 𝑥 𝑠 : the temperature of the sensor (°C).

When the node 𝑥 𝑠 is used, it is selected as the model output and is equivalent to the south zone air temperature measured 𝑇 𝑠 (°C), otherwise 𝑥 𝑖 is the model output.

The different heat transfers are characterized by the thermal resistances:

• 𝑅 𝑜 : between the outside and the middle of the building envelope (K/W),

• 𝑅 𝑖 : between half the building envelope and the south zone (K/W),

• 𝑅 𝑚 : between the south zone and the medium (K/W),

• 𝑅 𝑧 : between the air in the north zone and the south zone (K/W),

• 𝑅 𝑠 : between the south zone and the sensor (K/W).

The accumulation of energy is modeled by the thermal capacities:

• 𝐶 𝑤 : for the building envelope (J/K),

• 𝐶 𝑖 : for the air in the south zone (J/K),

• 𝐶 𝑚 : for the medium (internal walls, furniture, etc) (J/K),

• 𝐶 𝑠 : for the sensor.

The parameters 𝑎 𝑊 and 𝑎 𝐼 are respectively the effective area through which the solar radiation enters the building envelope and the effective window area of the building.

State-space models can be easily deduced from thermal networks by considering heat balance in each temperature nodes [START_REF] Ghiaus | Causality issue in the heat balance method for calculating the design heating and cooling load[END_REF]:

𝐱̇= 𝐀(𝛉)𝐱 + 𝐁(𝛉)𝐮 𝐲 𝑘 = 𝐂𝐱 𝑘 ( 4
)
where 𝐀 is the state matrix, 𝐁 the input matrix, 𝐂 the output matrix 𝛉 the parameter vector, 𝐮 the input vector and 𝐲 𝑘 is an output vector which can be measured at discrete time instant 𝑡 𝑘 (e. g. 𝐱 𝑘 = 𝐱(𝑡 𝑘 )).

Both models from Figure 2 If the input is assumed constant in the time interval ∆ 𝑡 (zero order hold), then 𝛂 = 𝟎, whereas if the input is assumed to vary linearly (first order hold), then [START_REF] Kristensen | Continuous Time Stochastic Modelling[END_REF])

𝛂 = 𝒖 𝑘+1 -𝒖 𝑘 ∆ 𝑡 (7)
In order to make the notation lighter, the parameter dependence in 𝐀 𝐝 , 𝐁 𝐝𝟎 , 𝐁 𝐝𝟏 , 𝚺 𝐰 and 𝚺 𝐯 is omitted such that 𝐀 𝐝 = 𝐀 𝐝 (𝛉).

Some parameters of the state-space model (6) may not be known and, then, they have to be estimated based on measured data.

Bayesian calibration with Markov Chain Monte Carlo

Markov chain Monte Carlo (MCMC) is a general method for constructing posterior distributions. The main idea is to simulate a Markov chain which has been constructed such that it has the posterior distribution as its stationary distribution [START_REF] Sarkka | Bayesian Filtering and Smoothing[END_REF]. The Metropolis-Hastings (MH) algorithm is the most common type of MCMC method due to its simplicity. MH is an iterative scheme, where a new candidate 𝛉 * is suggested from a proposed distribution 𝑞(𝛉 * |𝛉 𝑖-1 ) given the previous one 𝛉 𝑖-1 . The candidate is then accepted or rejected according to some acceptance probability

𝛼 𝑖 = min {1, 𝑝(𝐲 1:𝑘 |𝛉 * )𝑝(𝛉 * ) 𝑝(𝐲 1:𝑘 |𝛉 𝑖-1 )𝑝(𝛉 𝑖-1 ) 𝑞(𝛉 𝑖-1 |𝛉 * ) 𝑞(𝛉 * |𝛉 𝑖-1 ) } (8)
where the ratio 𝑞(𝛉 𝑖-1 |𝛉 * ) 𝑞(𝛉 * |𝛉 𝑖-1 ) ⁄ corrects the asymmetry in the proposed distribution.

If the candidate 𝛉 * increases significantly the posterior probability, the candidate is always accepted. However, a candidate which decrease the posterior probability can still be accepted as opposed to optimization algorithms; it allows the MH algorithm to explore regions of high posterior probability. Hence, by its stochastic nature, the MH algorithm may escape from local extrema which is a problem for many optimization algorithms used for ML estimation [START_REF] Dahlin | Accelerating Monte Carlo methods for Bayesian inference in dynamical models Accelerating Monte Carlo methods for Bayesian inference in dynamical models[END_REF].

The performance of the MH algorithm is highly dependent on the choice of the proposed distribution. A commonly used choice is the Gaussian random walk,

𝑞(𝛉 * |𝛉 𝑖-1 ) = 𝒩(𝛉 * |𝛉 𝑖-1 , 𝚺 𝛉 𝑖-1 ) (9)
where 𝒩(𝛉 * |𝛉 𝑖-1 , 𝚺 𝛉 𝑖-1 ) is a Gaussian probability density function of a random variable 𝛉 * ∈ ℝ 𝑁 𝑝 with mean

𝛉 𝑖-1 ∈ ℝ 𝑁 𝑝 and covariance 𝚺 𝛉 𝑖-1 ∈ ℝ 𝑁 𝑝 x 𝑁 𝑝 , 𝒩(𝛉 * |𝛉 𝑖-1 , 𝚺 𝛉 𝑖-1 ) = 1 (2𝜋) 𝑛 2 ⁄ |𝚺 𝛉 𝑖-1 | 1 2 ⁄ exp (- 1 2 (𝛉 * -𝛉 𝑖-1 ) T (𝚺 𝛉 𝑖-1 ) -1 (𝛉 * -𝛉 𝑖-1 )) (10)
Finding a suitable covariance matrix 𝚺 is a hard task which involves many trials and becomes unrealistic for highdimensional problems [START_REF] Sarkka | Bayesian Filtering and Smoothing[END_REF]. The Markov chain should converge to the stationary distribution in a reasonable time (burn-in phase) and the Markov chain should not be highly autocorrelated, such that the number of iteration for exploring the stationary distribution is minimized.

The performance and the tuning of the MH algorithm can be respectively improved and simplified by using the gradient and Hessian of the posterior distribution in order to construct a better proposed distribution [START_REF] Dahlin | Accelerating Monte Carlo methods for Bayesian inference in dynamical models Accelerating Monte Carlo methods for Bayesian inference in dynamical models[END_REF].

The next section shows a robust and accurate method for computing the likelihood, gradient and Hessian for linear Gaussian state-space models (6).

Construction of the posterior and proposal distribution

The construction of the posterior distribution (1) requires the evaluation of the likelihood 𝑝(𝐲 1:𝑘 |𝛉) and the prior distribution 𝑝(𝛉). The challenging part is the evaluation of the likelihood because the prior distribution is usually chosen such that it is easy to evaluate [START_REF] Sarkka | Bayesian Filtering and Smoothing[END_REF]. For a state-space model, the likelihood can be computed by using the prediction error decomposition

𝑝(𝐲 1:𝑘 |𝛉) = 𝑝(𝐲 1 |𝛉) ∏ 𝑝(𝐲 𝑘 |𝐲 1:𝑘-1 , 𝛉) 𝑁 𝑘=2 (11)
where the predictive likelihood can be computed recursively by

𝑝(𝐲 𝑘 |𝐲 1:𝑘-1 , 𝛉) = ∫ 𝑝(𝐲 𝑘 |𝐱 𝑘 , 𝛉) 𝑝(𝐱 𝑘 |𝐲 1:𝑘-1 , 𝛉) d𝐱 𝑘 (11.a)
with 𝑝(𝐲 𝑘 |𝐱 𝑘 , 𝛉) and 𝑝(𝐱 𝑘 |𝐲 1:𝑘-1 , 𝛉) representing respectively the measurement model and the predictive distribution of the state.

To avoid computational inaccuracy and instability, the logarithm of the unnormalized posterior distribution (righthand side of ( 1)), named log-posterior, is computed instead.

Since not all the states are not observed, the parameter estimation problem also requires to solve the state estimation problem [START_REF] Dahlin | Accelerating Monte Carlo methods for Bayesian inference in dynamical models Accelerating Monte Carlo methods for Bayesian inference in dynamical models[END_REF]. For linear Gaussian state-space model, the integral in (11.a) can be computed in closed form by the Kalman filter. The gradient and an approximation of the Hessian are obtained by differentiation of the Kalman filter equations with respect to the unknown parameters, referred to as the sensitivity equations. Therefore, computation of 𝑁 𝑝 sensitivity equations is required in parallel to the Kalman filter, where 𝑁 𝑝 is the number of unknown parameters. However, this strategy is numerically unstable due to rounding errors, i.e., the state covariance matrix 𝐏 may cease to be symmetric and positive definite, which leads to the failure of the computational process. This problem can be solved by using a robust square root implementation [START_REF] Kulikova | A unified square-root approach for the score and Fisher information matrix computation in linear dynamic systems[END_REF][START_REF] Tsyganova | On efficient parametric identification methods for linear discrete stochastic systems[END_REF], where only the square root factor 𝐏 1/2 is propagated instead of the full state covariance matrix. The upper triangular matrix 𝐏 1/2 is obtained by Cholesky decomposition, such that 𝐏 = 𝐏 T/2 𝐏 1/2 .

The recursion starts by updating the mean value of the prior state 𝐱 𝑘|𝑘-1 and prior square root factor 𝐏 𝑘|𝑘-1 1 2 ⁄ with the current measurements

𝐱 𝑘|𝑘 = 𝐱 𝑘|𝑘-1 + 𝐊 ̅ 𝑘 𝐞 ̅ 𝑘 (12)
with the standardized residuals

𝐞 ̅ 𝑘 = 𝐒 𝑘 -1 2 ⁄ (𝐲 𝑘 -𝐂 𝐝 𝐱 𝑘|𝑘-1 ) (13)
The normalized Kalman gain 𝐊 ̅ 𝑘 and the square root factor of the residual covariance, 𝐒 𝑘 1 2 ⁄ , are directly read from the post-array, in the right-hand side of equation ( 14)

𝐐 𝑘 MU [ 𝚺 𝐯 1 2 ⁄ 0 𝐏 𝑘|𝑘-1 1 2 ⁄ 𝐂 𝐝 T 𝐏 𝑘|𝑘-1 1 2 ⁄ ] ⏟ Pre-array = [ 𝐒 𝑘 1 2 ⁄ 𝐊 ̅ 𝑘 T 0 𝐏 𝑘|𝑘 1 2 ⁄ ] ⏟ Post-array (14)
where 𝐐 𝑘 MU is an orthogonal rotation matrix obtained by QR decomposition of the pre-array, such that the postarray is upper triangular; this notation is used throughout the paper.

The partial derivatives of the quantities in the post-array ( 14) are computed with

𝜕𝐒 𝑘 1 2 ⁄ 𝜕𝜃 𝑖 = ((𝓛 𝑖 ‡ ) T + 𝓓 𝑖 ‡ + 𝓤 𝑖 ‡ ) 𝐒 𝑘 1 2 ⁄ (15) 𝜕𝐊 ̅ 𝑘 𝜕𝜃 𝑖 = 𝐘 𝑖 T + ((𝓛 𝑖 ‡ ) T -𝓛 𝑖 ‡ ) 𝐊 ̅ 𝑘 T + (𝐕 𝑖 𝐒 𝑘 -1 2 ⁄ ) T 𝐏 𝑘|𝑘 1 2 ⁄ (16) 𝜕𝐏 𝑘|𝑘 1 2 ⁄ 𝜕𝜃 𝑖 = ((𝓛 𝑖 † ) T + 𝓓 𝑖 † + 𝓤 𝑖 † ) 𝐏 𝑘|𝑘 1 2 ⁄ ( 17
)
Where 𝓛 † , 𝓓 † , 𝓤 † , 𝓛 ‡ , 𝓓 ‡ and 𝓤 ‡ are obtained by first multiplying the orthogonal rotation matrix 𝐐 𝑘 MU to the partial derivative of the pre-array in ( 14)

𝐐 𝑘 MU [ 𝜕𝚺 𝐯 1 2 ⁄ 𝜕𝜃 𝑖 0 𝜕(𝐏 𝑘|𝑘-1 1 2 ⁄ 𝐂 𝐝 T ) 𝜕𝜃 𝑖 𝜕𝐏 𝑘|𝑘-1 1 2 ⁄ 𝜕𝜃 𝑖 ] = [ 𝐗 𝑖 𝐘 𝑖 𝐕 𝑖 𝐖 𝑖 ] (18)
and then the post-array ( 18) is multiplied by the inverse of the pre-array ( 14)

𝐆 = [ 𝐗 𝑖 𝐘 𝑖 𝐕 𝑖 𝐖 𝑖 ] [ 𝐒 𝑘 1 2 ⁄ 𝐊 ̅ 𝑘 T 0 𝐏 𝑘|𝑘 1 2 ⁄ ] -1 (19)
The matrices 𝓛 † , 𝓓 † and 𝓤 † are respectively the lower triangular, diagonal, and upper triangular parts of the submatrix 𝐆, from row and column 𝑁 𝑦 + 1 to row and column 𝑁 𝑦 + 𝑁 𝑥 . Matrices 𝓛 ‡ , 𝓓 ‡ and 𝓤 ‡ are respectively the lower triangular, diagonal, and upper triangular parts of the submatrix 𝐆, from the first row and column to row and column 𝑁 𝑦 [START_REF] Tsyganova | On efficient parametric identification methods for linear discrete stochastic systems[END_REF].

The partial derivative of the prior state mean update is computed by

𝜕𝐱 𝑘|𝑘 𝜕𝜃 𝑖 = 𝜕𝐱 𝑘|𝑘-1 𝜕𝜃 𝑖 + 𝜕𝐊 ̅ 𝑘 𝜕𝜃 𝑖 𝐞 ̅ 𝑘 + 𝐊 ̅ 𝑘 𝜕𝐞 ̅ 𝑘 𝜕𝜃 𝑖 (20)
with

𝜕𝐞 ̅ 𝑘 𝜕𝜃 𝑖 = -𝐒 𝑘 -1 2 ⁄ ( 𝜕𝐒 𝑘 1 2 ⁄ 𝜕𝜃 𝑖 𝐞 ̅ 𝑘 + 𝜕𝐂 𝐝 𝜕𝜃 𝑖 𝐱 𝑘|𝑘-1 + 𝐂 𝐝 𝜕𝐱 𝑘|𝑘-1 𝜕𝜃 𝑖 ) (21)
The posterior state mean and square root factor are propagated forward in time by the state equation similar to ( 6)

𝐱 𝑘+1|𝑘 = 𝐀 𝐝 𝐱 𝑘|𝑘 + 𝐁 𝐝𝟎 (𝐮 𝒌 + 𝛂∆ 𝑡 ) -𝐁 𝐝𝟏 𝛂 (22)
and by

𝐐 𝑘 TU [ 𝐏 𝑘|𝑘 1 2 ⁄ 𝐀 𝐝 T 𝚺 𝐰 1 2 ⁄ ] = [ 𝐏 𝑘+1|𝑘 1 2 ⁄ 0 ] (23)
The partial derivative of the state equation ( 22) is

𝜕𝐱 𝑘+1|𝑘 𝜕𝜃 𝑖 = 𝜕𝐀 𝐝 𝜕𝜃 𝑖 𝐱 𝑘|𝑘 + 𝐀 𝐝 𝜕𝐱 𝑘|𝑘 𝜕𝜃 𝑖 + 𝜕𝐁 𝐝𝟎 𝜕𝜃 𝑖 (𝐮 𝑘-1 + 𝛂∆ 𝑡 ) - 𝜕𝐁 𝐝𝟏 𝜕𝜃 𝑖 𝛂 (24)
where the partial derivatives of the state and input discrete matrices with respect to the continuous parameters are computed by [START_REF] Mbalawata | Parameter estimation in stochastic differential equations with Markov chain Monte Carlo and non-linear Kalman filtering[END_REF])

[ 𝐀 𝐝 0 𝜕𝐀 𝐝 𝜕𝜃 𝑖 𝐀 𝐝 ] ⏟ 𝐀 𝐈𝐌 = exp ( [ 𝐀 0 𝜕𝐀 𝜕𝜃 𝑖 𝐀 ] ⏟ 𝐀 𝐌 ∆𝑡 ) (24.a) [ 𝐁 𝐝𝟎 𝜕𝐁 𝐝𝟎 𝜕𝜃 𝑖 ] = 𝐀 𝐌 -1 (𝐀 𝐈𝐌 -𝐈) [ 𝐁 𝜕𝐁 𝜕𝜃 𝑖 ] (24.b) [ 𝐁 𝐝𝟏 𝜕𝐁 𝐝𝟏 𝜕𝜃 𝑖 ] = 𝐀 𝐌 -1 (-𝐀 𝐌 -1 (𝐀 𝐈𝐌 -𝐈) + 𝐀 𝐈𝐌 ∆ 𝑡 ) [ 𝐁 𝜕𝐁 𝜕𝜃 𝑖 ] (24.c)
The partial derivative of the square root factor in the post-array ( 23)

𝜕𝐏 𝑘+1|𝑘 1 2 ⁄ 𝜕𝜃 𝑖 = (𝓛 𝑖 T +𝓓 𝑖 + 𝓤 𝑖 )𝐏 𝑘+1|𝑘 1 2 ⁄ (25)
requires to first multiply the orthogonal rotation matrix 𝐐 𝑘 TU by the partial derivative of the pre-array ( 23)

𝐐 𝑘 TU [ 𝜕(𝐏 𝑘|𝑘 1 2 ⁄ 𝐀 𝐝 T ) 𝜕𝜃 𝑖 𝜕𝚺 𝐰 1 2 ⁄ 𝜕𝜃 𝑖 ] = [ 𝐀 𝑖 𝐂 𝑖 ] (26)
where the matrices 𝓛, 𝓓 and 𝓤 are respectively the lower triangular, diagonal, and upper triangular parts of the

matrix product 𝐀 𝑖 𝐏 𝑘|𝑘 -T 2 ⁄ .
The log-likelihood is recursively computed with

ln 𝑝(𝐲 1:𝑘 |𝛉) = 1 2 ∑ - 𝑁 𝑦 2 𝑁 𝑘=1 ln(2𝜋) -ln (det(𝐒 𝑘 T 2 ⁄ 𝐒 𝑘 1 2 ⁄ )) -𝐞 𝑘 T 𝐞 ̅ 𝑘 (27)
where the standardized innovations 𝐞 ̅ 𝑘 and the square-root of the innovation covariance matrix 𝐒 𝑘 1 2

⁄ are computed respectively by ( 13) and ( 14).

The gradient and the Hessian approximation of ( 27 (31) and are used in the MH algorithm to construct an efficient proposal distribution [START_REF] Dahlin | Accelerating Monte Carlo methods for Bayesian inference in dynamical models Accelerating Monte Carlo methods for Bayesian inference in dynamical models[END_REF])

𝑞(𝛉 * |𝛉 𝑖-1 ) = 𝒩 (𝛉 * |𝛉 𝑖-1 + 1 2 𝚺𝐇 ̂-1 (𝛉 𝑖-1 )𝐠(𝛉 𝑖-1 ), 𝚺𝐇 ̂-1 (𝛉 𝑖-1 )) (32)
where 𝚺 is a diagonal matrix which control the step length of the proposal distribution.

Using the geometric information of the posterior distribution has the advantage of steering the Markov chain towards areas of high posterior probability [START_REF] Nemeth | Parameter Estimation for State Space Models using Sequential Monte Carlo Algorithms[END_REF], which reduces the burn-in phase since the Markov chain takes larger steps when the Markov chain is far from the posterior mode and smaller steps as it gets closer.

This approach allows to save user time and also computational time. Firstly, because the covariance matrix is given by the inverse of the Hessian approximation, thus only the step length has to be tuned. Secondly, because it increases the mixing of the Markov chain, so the MH algorithm needs less iterations [START_REF] Dahlin | Accelerating Monte Carlo methods for Bayesian inference in dynamical models Accelerating Monte Carlo methods for Bayesian inference in dynamical models[END_REF][START_REF] Nemeth | Parameter Estimation for State Space Models using Sequential Monte Carlo Algorithms[END_REF].

The proposal distribution ( 32) is based on Newton-type optimization; consequently, the suggested candidates 𝛉 * are unconstrained and can violate the physical meaning of the system. A simple solution would be to reject candidates outside specified bounds but it could increase the autocorrelation of the Markov chain if too many candidates are rejected. A better solution is to reparametrize the model [START_REF] Dahlin | Accelerating Monte Carlo methods for Bayesian inference in dynamical models Accelerating Monte Carlo methods for Bayesian inference in dynamical models[END_REF].

Reparametrization of the model

The idea is to use non-linear functions to transform a constrained problem into an unconstrained one. In this way, The Jacobian matrix is triangular if each transformed parameter only depends on a single untransformed parameter, which simplifies the determinant computation to the product of the diagonal elements.

The gradient (30) and the Hessian (31) are with respect to 𝛉, so they need to be multiplied by the Jacobian 𝐉 in order to be with respect to 𝛈 (chain rule). The partial derivative of the new term in (35) needs also to be taken into account. The gradient and the Hessian in the unconstrained space are obtained by

𝐠(𝛈) = 𝐉 T 𝐠(𝛉) + 𝜕 ln|det(𝐉)| 𝜕𝛈 (36) 𝐇 ̂(𝛈) = 𝐉 T 𝐇 ̂(𝛉)𝐉 + 𝜕 2 ln|det(𝐉)| 𝜕 2 𝛈 (37)
A problem arises with this type of reparametrization when 𝜃 𝑖 gets close to a bound: its corresponding Jacobian term goes towards zero, so 𝐠(𝛈) and 𝐇 ̂(𝛈) are unreliable. Moreover, the Hessian estimate 𝐇 ̂(𝛈) could become illconditioned which prevents from its inversion. This issue is solved in optimization by adding a penalty function [START_REF] Kristensen | Continuous Time Stochastic Modelling[END_REF] which increases the gradient near the bounds, but this strategy is not suitable here because the proposal distribution (32) has a stochastic part, which means that candidates can still be projected towards the bounds. Instead, a prior distribution 𝑝(𝛉) is used to assign a low probability near the bounds. The details of the MH algorithm with gradient and Hessian information for constrained parameters are given in Algorithm 1.

Algorithm 1: Second order Metropolis-Hastings [START_REF] Dahlin | Accelerating Monte Carlo methods for Bayesian inference in dynamical models Accelerating Monte Carlo methods for Bayesian inference in dynamical models[END_REF] Inputs: 𝑁 (number of iterations), 𝛉 0 (initial parameters), 𝚺 (step length)

Output: 𝛈 1:𝑁 (samples from the posterior distribution)

1. Transformation to the unconstrained parameter space 𝛈 0 = 𝐟(𝛉 0 ) with ( 33) and ( 34)

2. Compute: ln 𝑝(𝛈 0 |𝐲 1:𝑘 ), 𝐠(𝛈 0 ) and 𝐇 ̂(𝛈 0 ) with ( 35), ( 36) and ( 37)

3. 𝑓𝑜𝑟 𝑖 = 1 𝑡𝑜 𝑁 4. Suggest a new candidate 𝛈 * ~ 𝒩 (𝛈 𝑖-1 + 1 2 𝚺𝐇 ̂-1 (𝛈 𝑖-1
)𝐠(𝛈 𝑖-1 ), 𝚺𝐇 ̂-1 (𝛈 𝑖-1 ))

5.

Compute: ln 𝑝(𝛈 * |𝐲 1:𝑘 ), 𝐠(𝛈 * ) and 𝐇 ̂(𝛈 * ) with ( 35), ( 36) and (37)

6.

Compute the acceptance probability 𝛼 𝑖 with (8)

7.

Generate a uniform random variable 𝑢 ~ 𝒰(0,1) and set

Choice of prior distribution

The knowledge of possible parameter values before anything has been observed is represented probabilistically by the prior distribution 𝑝(𝛉). A prior distribution, which is relevant with the experiment, the physical nature of the problem, or for other reasons, has to be specified by the user. Three categories of prior information are considered:

non-informative, weakly informative and informative (see [START_REF] Gelman | Understanding predictive information criteria for Bayesian models[END_REF], for a complete discussion on the subject).

Non-informative prior distributions attempt to not affect the posterior distribution, such that only the information in the data are relevant; this is the idea behind the ML estimation. But, these flat or almost flat prior distributions put more probability mass outside the expected range of values than inside, which can have unforeseen effect on the posterior distribution [START_REF] Dahlin | Accelerating Monte Carlo methods for Bayesian inference in dynamical models Accelerating Monte Carlo methods for Bayesian inference in dynamical models[END_REF], especially for small data set. Moreover, non-informative prior distributions, such as 𝒰(-∞, +∞) may be improper (they do not integrate to one), thus they cannot be expressed as a probability density function. In some cases, proper posterior distribution can be obtained with an improper prior distribution, but the result must be interpreted with care (Gelman, Carlin, et al. 2014).

Weakly informative prior distributions provide sufficient information to keep the parameters in a reasonable range and unlike informative prior distributions, they are not likely to outweigh the likelihood. For cases where the data set is too short or not enough informative, weakly informative prior distribution contains enough information to regularize the posterior distribution and prevent from identifiability issues; the curvature around the expected solution is increased (Team Stan Development 2015).

For the parameters transformed with the logarithm (33), the prior distributions are 𝜃 𝑖 ~ 𝒢(𝑎, 𝑏), where 𝒢(𝑎, 𝑏) denotes a Gamma distribution with shape 𝑎 and expected value 𝑏. The hyper-parameters 𝑎 and 𝑏 are chosen such that the probability near zero is low and that the distribution covers the expected range of values (Figure 3). For the transformed parameter with lower and upper bounds, the prior distributions are The choice of the prior distribution is an important decision which can strongly influence the posterior distribution [START_REF] Dahlin | Accelerating Monte Carlo methods for Bayesian inference in dynamical models Accelerating Monte Carlo methods for Bayesian inference in dynamical models[END_REF]) but the exploration of the posterior distribution depends on the tuning of the algorithm. The use of the Hessian approximation reduces the tuning of the proposal distribution to the choice of the step length matrix 𝚺 (32). Separate step lengths for each parameter can be used, but to simplify the tuning, when a single step length is used, such that 𝚺 = ε𝐈 𝑁𝑝 . The step length affects directly the acceptance rate of the MH algorithm (the percentage of accepted candidates at stationarity). A too large 𝚺 produces broad jumps which are more likely to be rejected, which increases the autocorrelation of the Markov chains and give a low acceptance rate. On the contrary, if the step length is too small, short jumps are likely to be accepted, which gives a higher acceptance rate. However, it limits the exploration of the posterior distribution to a small neighborhood, which also increases the autocorrelation. Consequently, the acceptance rate alone is not a correct indicator of the algorithm performance.

A better solution is to look at the mixing of the Markov chains at stationarity, which can be quantified by the integrated autocorrelation time (IACT) [START_REF] Dahlin | Accelerating Monte Carlo methods for Bayesian inference in dynamical models Accelerating Monte Carlo methods for Bayesian inference in dynamical models[END_REF])

IACT ̂(𝜃 𝑗 𝑁 𝑏 :𝑁 ) = 1 + 2 ∑ 𝜌 ̂𝑙 (𝜃 𝑗 𝑁 𝑏 :𝑁 ) 𝐿 𝑙=1 (38)
where 𝜌 𝑙 denotes the autocorrelation coefficient at lag 𝑙 of 𝜃 𝑗 𝑁 𝑏 :𝑁 , and 𝜃 𝑗 𝑁 𝑏 :𝑁 is the Markov chain of 𝜃 𝑗 from the burn-in time 𝑁 𝑏 to the last iteration 𝑁. The number of lags 𝐿 is determined as the first index for which |𝜌 ̂𝑙(𝜃 𝑗 𝑁 𝑏 :𝑁 )| < 2 √𝑁 -𝑁 𝑏 ⁄ , when the autocorrelation coefficient becomes statistically insignificant.

The IACT represents the number of iterations between two uncorrelated samples; therefore, the step length should be chosen such that it minimizes the IACT. The number of iteration 𝑁 should be chosen sufficiently large, such that, once the burn-in phase has been removed, the number of samples left are sufficient to represent the posterior distribution. Tools for diagnosing convergence are discussed in the next section.

Convergence diagnosis

The procedure of [START_REF] Gelman | Understanding predictive information criteria for Bayesian models[END_REF] is used here and presented briefly. The procedure consists of simulating 𝑀 Markov chains of 𝑁 samples, where the starting points of the 𝑀 Markov chains are randomly sampled from the prior distributions. The first step consists of inspecting visually the trace plots of the different Markov chains to determine the burn-in time 𝑁 𝑏 and to check if they converge to the same posterior distributions.

The 𝑀 Markov chains, with the burn-in phase removed, are split in two to give 2𝑀 chains of length (𝑁 -𝑁 𝑏 ) 2 ⁄ ;

then the variations between and within the 2𝑀 chains are compared (Gelman, Carlin, et al. 2014). The stationarity implies that the first and the second half of each sequence come from the same distribution. A good mixing requires that the variance inside chains should be closed to the variance between chains; this is quantified by the potential scale reduction 𝑅 ̂ (see [START_REF] Gelman | Understanding predictive information criteria for Bayesian models[END_REF] for computational details). The number of iterations 𝑁 should be increased until 𝑅 ̂ is near one or at least 𝑅 ̂< 1.1. The mixing of the Markov chains can also be quantified by the effective sample size (ESS) which approximates the number of independent samples in the 2𝑀 sequences. [START_REF] Gelman | Understanding predictive information criteria for Bayesian models[END_REF] suggest that the ESS should be at least superior to 5 × 2𝑀. If the aforementioned criteria are satisfied, the samples from the 2𝑀 sequences can be used to estimate the posterior distribution.

Model validation

After having calibrated different models, how to assess their reliability? The purpose of a model is to reproduce an input-output relationship; an intuitive way to start is by looking at what the model is not able to reproduce, the residuals 𝐞 ̅ (13) [START_REF] Ljung | System identification: theory for the user (second edition)[END_REF]. A plot of the residuals and the data allows to understand which features are not properly described by the model and who might be responsible for. The noise terms in model ( 6) are assumed to be white noise sequences which implies that it should also be the case for the residuals. The white noise sequence is uncorrelated, normally distributed with a zero mean and is uniformly distributed on all frequencies [START_REF] Madsen | Introduction to General and Generalized Linear Models[END_REF]. These properties are assessed by plotting the autocorrelation function (ACF) and the cumulated periodogram (CP) with their respective 95% confidence intervals. Furthermore, the residuals should be independent of the past inputs, which is tested by plotting the cross-correlation function (CCF).

The reliability of the model is also tested on a data set which has not been used for the calibration (validation data set). New values for the inputs are introduced in the model and the simulated output is compared to the measured one. If the identification data set is informative enough, i.e. the different dynamics of the system are observable in the data, the model should be representative of the system and therefore the simulation should be close to the measurement.

The model validation assesses the reliability of a model and gives insight on the model order selection and directions for improvement; but how to select the best model?

Model comparison

In section 3.1, the selection of a model structure based on insights of the experiment has been presented. This section discusses the agreement of calibrated models with measured data and how to compare different models in order to select the most appropriate one. The model fit to the data is summarized by the log-likelihood ln 𝑝(𝐲 1:𝑘 |𝛉)

(27) and the prior distribution is not relevant for assessing the accuracy of a model. The best model is not necessarily characterized by the highest log-likelihood value because as the complexity of a model increases, the number of degrees of freedom increases, the parameters adjust themselves to fit a particular realization of the noise (overfitting) [START_REF] Ljung | System identification: theory for the user (second edition)[END_REF]. In order to adjust for overfitting, the Akaike's Information Criterion (AIC) and

Bayesian Information Criterion (BIC) penalize the log-likelihood in function of the complexity of the model:

AIC = -
where ln 𝑝(𝐲 𝑘 |𝛉) ∈ ℝ (𝑁-𝑁 𝑏 ) x 𝑁 𝑠 is the log-likelihood at time instant 𝑘 and (𝑁 -𝑁 𝑏 ) is the number of sample used to approximate the posterior distribution ln 𝑝(𝛉|𝐲 1:𝑘 ). Differently of the AIC and BIC, the WAIC is penalized by the dispersion of the log-likelihood.

Evaluating these criteria on same data set used for the calibration introduced a bias in the model selection process and therefore it is advised to evaluate them with the validation data set instead; this point is illustrated in the next section for a real test case.

Application to the Twin houses experiment

Model comparison

The capabilities of the second order MH (Algorithm 1) are now tested on the twin houses experiment presented in section 2. The purpose is to calibrate the models ℳ 3 and ℳ 4 (Figure 2) where the south zone temperature (green zone Figure 1) is the output and two boundary conditions are considered, the outdoor temperature 𝑇 𝑜 (°C) and the north zone temperature 𝑇 𝑛 (°C). The vector of unknown parameters 𝛉 for ℳ 3 and ℳ 4 are respectively given in 2015) used a model with four states which was correlated with the HVAC system and the solar radiations. They improved their model by scaling the system noise with respect to the heater and solar radiation signal. A validation data set has not been used to demonstrate the effectiveness of the model; only the improvement of the residuals is shown. With a zero order-hold assumption, the model is not able to understand the fast dynamic of the heating signal. Indeed, the time response of the electric heaters is estimated between 1 and 2 minutes, which is faster than the sampling time of the data (10 minutes). This issue is solved by considering that the inputs vary linearly between two samples (first order-hold).

Around 24 days of data are used, where the first 14 days are the identification data set (ROLBS sequence, Figure 4) and the last 10 days are the validation data set (Figure 5); the detail of the inputs and outputs is given in section 2. such that it minimizes the IACT. This tuning gives an acceptance rate of approximately 30% for ℳ 3 and 25% for ℳ 4 . The diagnosis procedure presented in section 3.2.5 has been applied with 𝑀 = 6 Markov chains with initial parameter values randomly sampled from their respective prior distributions (Table 1 andTable 2); this is illustrated by the trace plot of the first thousand iterations in Figure 6. The first 500 samples of 𝑁 = 5500 are discarded as burn-in for the model ℳ 3 whereas the first 1500 samples of 𝑁 = 6500 are discarded for ℳ 4 .

The chains at stationarity for both models are split in two to give 12 chains of 2500 samples which are used to quantify the mixing of the MH algorithm. The results are summarized in Table 1 and Table 1, with the worst values highlighted in red and the best in green. The worst potential scale reductions 𝑅 ̂ are below the threshold of 1.1 and the worst effective samples sizes (ESS) are easily above 5 × 12. Consequently the 12 × 2500 samples can be used to approximate the posterior distributions of the parameters. The posterior distributions of the 6 simulated Markov chains are represented in Figure 7 for ℳ 3 and in Figure 8 for ℳ 4 with different colors; the black line is the global approximation of the posterior distributions using all samples. The reliability of the model is also tested by comparing the measured south zone temperature of the validation data set with the simulated output. A clear advantage of Bayesian estimation is that it is possible to simulate directly from the posterior distribution, which gives a simulated output with all the uncertainties. This is very useful for model predictive control; the weather forecast is introduced in the estimated model to predict the indoor temperature. Afterwards, the trade-off between comfort and energy saving is chosen by taking either the lowest temperature prediction, such that the HVAC system is sure to maintain the comfort or by taking the highest predicted temperature such that the HVAC system uses the minimal amount of energy. The simulation from the posterior distribution is plotted in Figure 11. The measured south zone temperature is always included in the simulated output distribution for both models, but the dispersion for the model ℳ 4 is more important.

In order to select the most appropriate model, the performances of both models (section 3.4) are summarized in Table 3. For the identification data set, the model ℳ 4 should be accepted against ℳ 3 ; the AIC, BIC and WAIC are smaller for ℳ 4 than for ℳ 3 , and the 𝑝 𝑉𝑎𝑙𝑢𝑒 of the LRT confirms this choice. However, the criteria evaluated on the validation data set indicate that the model ℳ 3 should be preferred; the log-likelihood of the model ℳ 3 is higher and less dispersed than the log-likelihood of the model ℳ 4 , as shown in Figure 11. The performance gap between both models is significant for the identification data set in comparison to the validation data set. In the identification data set, the ROLBS introduces an unconventional dynamic which is not representative of the intended use of the building. A more complex model is required to fit the fast variations of the south zone air temperature, but these fast variations are not present in a conventional use and therefore the extra complexity of the model ℳ 4 is not required. As mentioned in section 3.4, using the same information from the data for calibration and for selection may be misleading (Gelman, Hwang, et al. 2014). In this case, it is also illustrated by the whiteness improvement of the residuals for the model ℳ 4 . It can be concluded that the model ℳ 3 is more representative of the south zone and consequently, only the model ℳ 3 is considered in the following of the paper.

Physical interpretation of the results

The posterior distributions are compared against the building characteristics which are available in [START_REF] Strachan | Empirical Whole Model Validation Modelling Specification Validation of Building Energy Simulation Tools, Team Stan Development[END_REF]. The envelope thermal resistance is given by

𝑅 ̅ 𝑤 = (∑ 𝑈 𝑗 𝑆 𝑗 3 𝑗=1 + 𝑈 𝑤𝑖 𝑆 𝑤𝑖 ) -1 = 4.47 • 10 -2 K/W (43)
where 𝑈 𝑗 and 𝑆 𝑗 are the U-values and surfaces of the different walls (south, east and west) and 𝑈 𝑤𝑖 and 𝑆 𝑤𝑖 are the U-value and surfaces of the windows.

The envelope thermal resistance is estimated with the posterior distributions of 𝑅 𝑜 and 𝑅 𝑖 such that

𝑅 𝑤 = 𝑅 𝑜 + 𝑅 𝑖 ∈ [ 4.11 • 10 -2 9.38 • 10 -2 ] (43.a)
The resistance 𝑅 ̅ 𝑤 belongs to the estimated interval (43.a), but this interval is large which means that the uncertainties are important.

The thermal resistance between the south and the north zone is given by

𝑅 ̅ 𝑧 = (𝑈 𝑧1 𝑆 𝑧1 + 𝑈 𝑧2 𝑆 𝑧2 + 3𝑈 𝑑𝑜𝑜𝑟 𝑆 𝑑𝑜𝑜𝑟 ) -1 = 1.53 • 10 -2 K/W (44)
where 𝑈 𝑧1 and 𝑆 𝑧1 are the U-value and the surface of the north wall of the living room, 𝑈 𝑧2 and 𝑆 𝑧2 are the Uvalue and the surface of the north wall of the bathroom and the corridor, and 𝑈 𝑑𝑜𝑜𝑟 and 𝑆 𝑑𝑜𝑜𝑟 are the U-value and the surface of the doors.

The estimated range from the posterior distribution of 𝑅 𝑧 is [ 4.59 • 10 -3 5.52 • 10 -3 ] which is far smaller than the value computed in (44). This gap could mean that the infiltration between the two zones are significant.

The envelope thermal capacity is defined by

𝐶 ̅ 𝑤 = ∑ 𝐶 𝑗 𝑆 𝑗 3 𝑗=1 = 1.30 • 10 7 J/K (45)
with 𝐶 𝑗 the heat capacity. The thermal capacity of the windows is negligible as compared to 𝐶 𝑤 .

The estimated posterior distribution of 𝐶 𝑤 covers the following range of values [ 2.21 • 10 6 4.47 • 10 6 ], which is more than half less than the expected value (45).

The thermal capacity of the medium 𝐶 𝑚 consists of the inner walls of the south zone but also of parts of the ceiling and ground floor, such as

𝐶 ̅ 𝑚 = 𝐶 𝑖𝑤 1 + 𝐶 𝑖𝑤 2 + 𝐶 𝑔𝑟𝑜𝑢𝑛𝑑 + 𝐶 𝑐𝑒𝑖𝑙𝑖𝑛𝑔 = 8.00 • 10 7 J/K (46)
where the subscripts 𝑖𝑤 1 and 𝑖𝑤 2 denote respectively the east wall of the living room and the other light walls.

The estimated range for 𝐶 𝑚 is [ 1.22 • 10 7 1.54 • 10 7 ] which is the expected order of magnitude since ( 46) is overestimated by taking into account all the volume of the ground floor and the ceiling.

The thermal capacity of the indoor air is simply given by 𝐶 ̅ 𝑖 = 𝜌 𝑎 𝑐 𝑎 𝑉 = 1.79 • 10 5 J/K (47)

with 𝑐 𝑎 and 𝜌 𝑎 the specific heat and the density of the air, and 𝑉 the volume of the south zone.

In this case as well, the estimated posterior distribution of 𝐶 𝑖 is consistent with the order of magnitude of the expected value, where 𝐶 𝑖 ∈ [ 6.60 • 10 5 6.90 • 10 5 ].

Determining prior knowledge on the convective resistance 𝑅 𝑚 between the indoor air and the medium is not an easy task and it is not of main interest in this study. The parameters 𝑎 𝐼 and 𝑎 𝑊 are interpreted as effective areas because 𝑄 ̇𝑔ℎ is measured on a horizontal surface. Nevertheless, it is interesting to see that 𝑎 𝐼 is superior to 𝑎 𝑊 which shows the importance of direct solar radiations into the south zone.

The time constants 𝛕 of the continuous system (4) are computed by

𝛕 = - 1 𝛌 ( 48
)
where 𝛌 are the eigenvalues of the state matrix and 𝛕, 𝛌 ∈ ℝ ,.

The time constants of the model ℳ 3 , given in Table 4, are consistent in the range of the fast dynamics of the air and the slow accumulation of energy in the medium.

The performances of the second-order MH are compared to a ML estimation in the next section. Furthermore, the regularization effect of the prior distribution is illustrated by identifiability analysis. The performance of the second-order MH algorithm is now tested against a ML estimation with a quasi-Newton optimization. The log-likelihood and its gradient are supplied to the unconstrained MATLAB's function fminunc 1 (optimality tolerance and step tolerance fixed to 10 -8 ) which use a BFGS approximation of the Hessian. The same parameter transformations are used except for the standard deviations; they are bounded between 10 -8 and 5 because it has been observed that with the logarithm transformation (33), the standard deviations can be too close to zero which cause numerical instabilities. The penalty function given by [START_REF] Kristensen | Continuous Time Stochastic Modelling[END_REF] is used to repulse the parameters near the bounds. The ML estimation is repeated 6 times with the same initial parameters as for MH algorithm; the ML parameter estimates are given in Table 5. The results are highly dependent on the initial conditions for most of the parameters; some of them are different of several orders of magnitude. It can be concluded in this case that the second-order MH algorithm has a better global convergence than the ML optimization.

1 https://fr.mathworks.com/help/optim/ug/fminunc.html A closer look is given to the parameter 𝜎 𝑤 22 which seems to be unidentifiable compared to the parameter 𝐶 𝑖 . The profiles of the log-likelihood and log-posterior are plotted in Figure 12. These profiles are obtained by maximizing the log-likelihood and the log-posterior with respect to all parameters except 𝜎 𝑤 22 and 𝐶 𝑖 , which are fixed for each optimization to different values (x-axis, Figure 12). Confidence intervals for the profile likelihood can be computed by [START_REF] Madsen | Introduction to General and Generalized Linear Models[END_REF] ln 𝑝(𝐲 1:𝑘 |𝛉 ~𝑖) -ln 𝑝(𝐲 1:𝑘 |𝛉) > -1 2 𝜒 1-𝛼 2 (𝑝)

where ln 𝑝(𝐲 1:𝑘 |𝛉 ~𝑖) is the maximum log-likelihood with respect to all parameters except 𝜃 𝑖 , ln 𝑝(𝐲 1:𝑘 |𝛉) is the maximum log-likelihood and 𝜒 1-𝛼 2 (𝑝) is a chi squared distribution with 𝑝 degrees of freedom and confidence level 𝛼.

The 95% confidence intervals given by -1 2 𝜒 0.95 2 (1) = -1.92 are represented in Figure 12 by the dotted black lines. The profile of the log-posterior is similarly computed by ln 𝑝(𝛉 ~𝑖|𝐲 1:𝑘 ) -ln 𝑝(𝛉|𝐲 1:𝑘 ). These profiles are computed around the posterior modes given in Table 1; the idea is to visualize the quantity of information given by the data and the regularization brought by the prior distribution. The profile of the log-likelihood of 𝜎 𝑤 22 is asymmetric and almost flat around its maximum, which means that any value on the flat region has negligible effects on the log-likelihood; it explains the dispersion of the ML estimation in Table 5. The profile of the logposterior shows how the prior distribution (𝒢(2, 0.03), Figure 3) increases the curvature, especially towards zero, and regularizes the identifiability problem. For the parameter 𝐶 𝑖 , the information from the data overweight the prior distribution which means that the prior distribution has only a small effect on the posterior distribution. Two models of respectively 15 and 18 unknown parameters have been easily calibrated with the improved Metropolis-Hastings algorithms where a unique step-length has been used, which illustrates the gain of this method over a classical Metropolis-Hastings with random walk. Furthermore, it is shown that the second-order Metropolis-Hastings algorithm has a better robustness against the initial conditions than a maximum likelihood estimation with a quasi-Newton algorithm, and it is illustrated through an identifiability analysis, that the prior distributions act as regularization when the data are not informative enough.

It is highlighted that model selection criteria should be computed on a different data set than the one used for the calibration to avoid a biased selection. Indeed, in this experiment the unconventional excitation generated by the heaters implies that a more complex model should be selected, but this extra complexity is not required for a more conventional use of the HVAC system.

Figure 2 :

 2 Figure 2: 3 states model, ℳ 3 and 4 states model, ℳ 4

  the proposal distribution cannot suggest candidates which are outside the bounds. The constrained parameters 𝛉 are transformed to unconstrained parameters 𝛈 by a one-to-one invertible functions 𝛈 = 𝐟(𝛉). Two parametrizations are used: 1) the log transform𝜂 𝑖 = ln(𝜃 𝑖 )(33)constraints 𝜃 𝑖 between the open interval ]0, +∞[ and 2) the following transformation (probability (8), the unnormalized posterior distribution is computed in the constrained space whereas the proposal distribution is evaluated in the unconstrained one. To homogenize the acceptance probability, the log-posterior distribution is transformed in the unconstrained space by using the Jacobian adjustment(Gelman, Carlin, et al. 2014)ln 𝑝(𝛈|𝐲 1:𝑘 ) = ln 𝑝(𝛉|𝐲 1:𝑘 ) + ln|det(𝐉)|(35)where |det(𝐉)| is the absolute value of the determinant of the Jacobian matrix 𝐉; |det(𝐉)| adjusts for the distortion caused by the non-linear transformation; 𝐉 is the Jacobian matrix of the inverse transform 𝛉 = 𝐟 -𝟏 (𝛈), such that

  candidate{𝛈 𝑖 , ln 𝑝(𝛈 𝒊 |𝐲 1:𝑘 ) , 𝐠(𝛈 𝒊 ), 𝐇 ̂(𝛈 𝒊 ) } ← {𝛈 * , ln 𝑝(𝛈 * |𝐲 1:𝑘 ) , 𝐠(𝛈 * ), 𝐇 ̂(𝛈 * candidate {𝛈 𝑖 , ln 𝑝(𝛈 𝒊 |𝐲 1:𝑘 ) , 𝐠(𝛈 𝒊 ), 𝐇 ̂(𝛈 𝒊 ) } ← {𝛈 𝑖-1 , ln 𝑝(𝛈 𝑖-1 |𝐲 1:𝑘 ) , 𝐠(𝛈 𝑖-1 ), 𝐇 ̂(𝛈𝑖-1

  𝜃 𝑗 ~ 𝛽(2, 2, 𝜃 𝑗 min , 𝜃 𝑗 max ), where 𝛽(𝑎, 𝑏, 𝜃 𝑗 min , 𝜃 𝑗 max ) is a Beta distribution with shape hyper-parameters 𝑎 and 𝑏, lower and upper bounds 𝜃 𝑗 min and 𝜃 𝑗 max . The Beta distribution with 𝑎 = 𝑏 = 2 is symmetric and assigns lowprobabilities for values near the bounds (Figure3).

Figure 3 :

 3 Figure 3: Prior pdf, 𝒢(2, 0.03) in blue and 𝛽(2 ,2, 10 -4 , 2 • 10 -1 ) in red

Figure 4 :

 4 Figure 4: Identification data set, output, standardized residuals of ℳ 3 (blue) and ℳ 4 (orange), and inputs

Figure 6 :Figure 7 :

 67 Figure 6: Trace plot of the 6 Markov chains in the first thousand iterations (ℳ 3 )

Figure 9 :

 9 Figure 9: ACF, CP and CCF of the standardized residuals, ℳ 3

Figure 11 :

 11 Figure 11: Left: measured (black) and simulated south zone temperature with the validation data set; right: loglikelihoods for the validation data set, (ℳ 3 : blue, ℳ 4 : orange)

Figure 12 :

 12 Figure 12: Profile log-likelihood (black) and profile log-posterior (blue). The dashed lines represent the 95% confidence intervals of the log-likelihood and the red dots the respective maximums of the curves

  It has been shown in this section that the parameter estimation problem is also a state estimation problem. For a linear Gaussian state-space model, the log-likelihood with its gradient and Hessian approximation can be computed by a square root version of the Kalman filter. This numerically stable strategy requires only to run the square root Kalman filer and 𝑁 𝑝 sensitivity equations forward in time.

		𝜕 ln 𝑝(𝐲 1:𝑘 |𝛉) 𝜕𝜃 𝑖	= -∑ tr (𝐒 𝑘 -1 2 ⁄ 𝜕𝐒 𝑘 1 2 ⁄ 𝜕𝜃 𝑖 𝑁 𝑘=1	) +	𝑇 𝜕𝐞 ̅ 𝑘 𝜕𝜃 𝑖	𝐞 ̅ 𝑘	(28)
	-	𝜕 2 ln 𝑝(𝐲 1:𝑘 |𝛉) 𝜕𝜃 𝑖 𝜕𝜃 𝑗	≈ ∑ tr ( 𝑁 𝑘=1	𝜕𝐞 ̅ 𝑘 𝜕𝜃 𝑖	𝜕𝐞 ̅ 𝑘 𝑇 𝜕𝜃 𝑗	) + tr (𝐒 𝑘 -1 2 ⁄ 𝜕𝐒 𝑘 1 2 ⁄ 𝜕𝜃 𝑖	𝐒 𝑘 -1 2 ⁄ 𝜕𝐒 𝑘 1 2 ⁄ 𝜕𝜃 𝑗	)	(29)
	The gradient and the Hessian of the log-posterior distribution are easily computed with
		𝐠(𝛉) = 𝐇 ̂(𝛉) = -	𝜕 ln 𝑝(𝛉|𝐲 1:𝑘 ) 𝜕𝛉 𝜕 2 ln 𝑝(𝛉|𝐲 1:𝑘 ) + 𝜕 2 𝛉 -𝜕 ln 𝑝(𝛉) 𝜕𝛉 𝜕 2 ln 𝑝(𝛉) 𝜕 2 𝛉	(30)
							) are respectively obtained with

  𝛉 𝐌𝐋 is the ML estimate, 𝑁 𝑝 the number of parameters and 𝑁 𝑠 the sample size.The smallest AIC or BIC between different models indicates the most appropriate model.For nested models, like ℳ 3 and ℳ 4 , the likelihood ratio test (LRT) can be used[START_REF] Bacher | Identifying suitable models for the heat dynamics of buildings[END_REF] the ML estimate of model ℳ 3 and ℳ 4 .As the number of samples 𝑁 𝑠 goes to infinity, the LRT converges to 𝜒 2 distributed variable with (𝑁 𝑝 ℳ 4 -𝑁 𝑝 ℳ 3 ) degrees of freedom. Usually, a 𝑝 𝑉𝑎𝑙𝑢𝑒 of the LRT below 0.05, indicates that the improvement of the larger model ℳ 4 over ℳ 3 is significant and consequently the model ℳ 4 should be preferred.

	2 ln 𝑝(𝐲 1:𝑘 |𝛉 𝐌𝐋 ) + 2𝑁 𝑝	(39)
	BIC = -2 ln 𝑝(𝐲 1:𝑘 |𝛉 𝐌𝐋 ) + 𝑁 𝑝 ln 𝑁 𝑠	(40)
	where LRT = -2(ln 𝑝(𝐲 1:𝑘 |𝛉 𝐌𝐋 ℳ 3 ) -ln 𝑝(𝐲 1:𝑘 |𝛉 𝐌𝐋 ℳ 4 ))	(41)
	with 𝛉 𝐌𝐋 ℳ 3 and 𝛉 𝐌𝐋 ℳ 4 𝑁 𝑠
	𝑘=1

These criteria are based on point estimate 𝛉 𝐌𝐋 and not on the posterior distribution 𝑝(𝛉|𝐲 1:𝑘 ); a more Bayesian criterion is given by the Watanabe-Akaike Information Criterio (WAIC)

(Gelman, Hwang, et al. 2014

) WAIC = -2 ∑ mean(ln 𝑝(𝐲 𝑘 |𝛉)) -var(ln 𝑝(𝐲 𝑘 |𝛉))

Table 1 and

 1 Table 2. The process noise covariance is defined by 𝚺 𝐰 (𝛉) = diag(𝜎 𝑤 11 The standard deviation 𝜎 𝑤 44 of the state 𝑥 𝑠 in ℳ 4 has been fixed to 10 -6 instead of putting an informative prior distribution. This problem has been investigated by different authors who used relatively similar optimization strategies but different model structures. De Coninck et al. (2015) and Rehab & André (2015) used second order thermal models but with different structure and inputs. Himpe & Janssens (

	2	𝜎 𝑤 22 2	𝜎 𝑤 33 2 ) and the

Table 1 :

 1 Prior distributions, posterior modes and diagnosis tests (worst: red and best: green) of ℳ 3

		Prior distributions	Posterior modes	𝑅 ̂	ESS	Min/Max IACT ̂
	𝑅 𝑜 𝑅 𝑖 𝑅 𝑚 𝑅 𝑧 𝐶 𝑤 10 8 ⁄ ⁄ 𝐶 𝑖 10 8 ⁄ 𝐶 𝑚 10 8 𝑎 𝑊 𝑎 𝐼 𝜎 𝑤 11 𝜎 𝑤 22 𝜎 𝑤 33 𝜎 𝑣 𝑥 𝑤 0 𝑥 𝑚 0	𝛽(2, 2, 10 -4 , 2 • 10 -1 ) 𝛽(2, 2, 10 -4 , 2 • 10 -1 ) 𝛽(2, 2, 10 -4 , 10 -1 ) 𝛽(2, 2, 10 -4 , 10 -1 ) 𝛽(2, 2, 10 -3 , 5 • 10 -1 ) 𝛽(2, 2, 10 -4 , 10 -1 ) 𝛽(2, 2, 10 -2 , 5) 𝛽(2, 2, 10 -1 , 5) 𝛽(2, 2, 10 -1 , 5) 𝒢(2, 0.03) 𝒢(2, 0.03) 𝒢(2, 0.03) 𝒢(2, 0.03) 𝛽(2, 2, 15, 45) 𝛽(2, 2, 15, 45)	5.53 • 10 -2 2.09 • 10 -3 2.31 • 10 -3 4.98 • 10 -3 3.12 • 10 -2 6.73 • 10 -3 1.38 • 10 -1 1.06 1.24 1.02 • 10 -1 1.38 • 10 -2 1.82 • 10 -2 1.69 • 10 -2 29.20 29.45	1.0100 𝟏. 𝟎𝟎𝟑𝟗 1.29 • 10 3 1.22 • 10 3 1.0045 1.25 • 10 3 1.0086 1.16 • 10 3 1.0049 1.31 • 10 3 1.0080 1.28 • 10 3 1.0061 𝟏. 𝟒𝟗 • 𝟏𝟎 𝟑 1.0132 1.19 • 10 3 1.0143 1.39 • 10 3 1.0065 1.26 • 10 3 𝟏. 𝟎𝟐𝟐𝟗 𝟓. 𝟐𝟖 • 𝟏𝟎 𝟐 𝟑𝟔. 𝟐𝟗 𝟓𝟒. 𝟗𝟎 18.86 29.47 18.25 29.50 17.81 30.81 18.79 28.16 19.46 27.42 17.49 32.01 17.18 24.58 18.90 24.12 18.40 𝟐𝟐. 𝟗𝟗 19.06 30.88 1.0084 1.15 • 10 3 22.12 33.57 1.0179 6.26 • 10 2 29.39 50.03 1.0146 1.15 • 10 3 𝟏𝟔. 𝟒𝟑 24.49 1.0099 1.33 • 10 3 18.09 27.53

Table 2 :

 2 Prior distributions, posterior modes and diagnosis tests (worst: red and best: green) of ℳ 4

		Prior distributions	Posterior modes	𝑅 ̂	ESS	Min/Max IACT ̂
	𝑅 𝑜 𝑅 𝑖 𝑅 𝑚 𝑅 𝑠 𝑅 𝑧	𝛽(2, 2, 10 -4 , 10 -1 ) 𝛽(2, 2, 10 -4 , 10 -1 ) 𝛽(2, 2, 10 -4 , 10 -1 ) 𝛽(2, 2, 10 -3 , 5 • 10 -2 ) 𝛽(2, 2, 10 -4 , 10 -1 )	4.93 • 10 -2 1.83 • 10 -3 2.13 • 10 -3 5.35 • 10 -3 5.27 • 10 -3	1.0073 1.0127 1.0136 1.0099 1.0094	1.06 • 10 3 7.26 • 10 2 6.04 • 10 2 7.41 • 10 2 9.89 • 10 2	[24.17 41.44] [30.34 58.48] [34.27 66.75] [31.15 53.36] [24.93 35.84]

Table 4 :

 4 Time constants of model ℳ 3

			Time constant [hours]	ℳ 3
			𝜏 1	[ 1.57 • 10 -1 1.66 • 10 -1 ]
			𝜏 2		[2.19 3.75]
			𝜏 3		[25.84 34.34]
	4.3. Performance comparison with maximum likelihood estimation
	Table 5: ML estimation with the same random initial parameters as the MH algorithm, the values in bold
			represent the parameters closed to their boundaries
		MLE 1	MLE 2	MLE 3	MLE 4	MLE 5	MLE 6
	𝑅 𝑜 𝑅 𝑖 𝑅 𝑚 𝑅 𝑧 ⁄ 𝐶 𝑤 10 8 ⁄ 𝐶 𝑖 10 8 𝐶 𝑚 10 8 ⁄	𝟏. 𝟗𝟖 • 𝟏𝟎 -𝟏 𝟏. 𝟎𝟐 • 𝟏𝟎 -𝟒 𝟏. 𝟗𝟖 • 𝟏𝟎 -𝟏 𝟏. 𝟗𝟖 • 𝟏𝟎 -𝟏 𝟏. 𝟔𝟐 • 𝟏𝟎 -𝟏 𝟏. 𝟗𝟖 • 𝟏𝟎 -𝟏 1.40 • 10 -3 4.78 • 10 -2 1.40 • 10 -3 1.18 • 10 -3 1.35 • 10 -3 1.42 • 10 -3 2.24 • 10 -2 1.12 • 10 -3 2.19 • 10 -2 3.32 • 10 -3 2.91 • 10 -3 2.08 • 10 -2 2.94 • 10 -3 5.76 • 10 -3 2.92 • 10 -3 4.00 • 10 -3 𝟖. 𝟔𝟏 • 𝟏𝟎 -𝟐 2.87 • 10 -3 6.50 • 10 -2 𝟏. 𝟎𝟏 • 𝟏𝟎 -𝟑 6.46 • 10 -2 1.19 • 10 -1 6.94 • 10 -2 6.39 • 10 -2 6.75 • 10 -3 6.86 • 10 -3 6.76 • 10 -3 6.72 • 10 -3 6.74 • 10 -3 6.81 • 10 -3 3.03 • 10 -1 1.24 • 10 -1 3.12 • 10 -1 𝟏. 𝟑𝟑 • 𝟏𝟎 -𝟐 2.3 3.29 • 10 -1
	𝑎 𝑊	1.01	3.15	1.02	1.74	1.08	1.03
	𝑎 𝐼 𝜎 𝑤 11 𝜎 𝑤 22 𝜎 𝑤 33 𝜎 𝑣	1.23 5.83 • 10 -2 3.28 • 10 -5 𝟏. 𝟏𝟓 • 𝟏𝟎 -𝟖 1.30 • 10 -2 1.27 1.24 1.14 • 10 -2 5.50 • 10 -2 6.45 • 10 -1 6.58 • 10 -2 6.41 • 10 -1 1.90 • 10 -2 1.70 • 10 -2 1.63 • 10 -2	1.27 8.71 • 10 -2 1.32 • 10 -2 1.57 • 10 -5 1.63 • 10 -2	1.23 5.09 • 10 -2 2.90 • 10 -5 9.79 • 10 -2 1.90 • 10 -2 𝟏. 𝟑𝟖 • 𝟏𝟎 -𝟖 1.25 4.67 • 10 -2 2.51 • 10 -2 6.38 • 10 -1
	𝑥 𝑤 0	29.60	36.68	29.59	29.54	29.70	29.60
	𝑥 𝑚 0	𝟒𝟒. 𝟒𝟕	29.61	𝟒𝟒. 𝟓𝟎	29.95	25.00	𝟒𝟒. 𝟓𝟕

1. Acknowledgements This work was financially supported by BPI France in the FUI Project COMETE. Thanks you to Dr. Paul Strachan for making the Twin houses data available online; to the Dynastee network, particularly to Henrik Madsen, Peder Bacher and Rune Juhl for their statistical guidelines; as well as Dr. Johan Dahlin and Dr. Maria Kulikova for their time and expertise.