Cécile Durot

Hanna Jankowski

Behaviour of the monotone single index model under repeated measurements Fadoua Balabdaoui •

Keywords: Goodness-of-fit, isotonic regression, maximum likelihood, shapeconstrained estimation, single index model, set identifiability

published or not. The documents may come

Introduction

The class of generalized linear models was introduced as a unified, and hence more general, framework which encompasses many popular tools such as linear and logistic regression [START_REF] Nelder | Generalized linear models[END_REF]. Both these examples, which are widely used in practice, involve an increasing link function: the identity and logit function respectively. In fact, link functions in the usual generalized linear models are increasing. Thus, these generalized linear models can be themselves viewed as a sub-family of the even more general class of isotonic single index models. Furthermore, isotonic single index models offer flexibility and robustness against picking the wrong model. As rightly pointed out in Härdle et al. (1997, page 213) "the advantage of [the single index] approach is that an interpretable linear single index, a weighted sum of of the predictor variables, is still produced."

Estimation in the isotonic single index model has recently received increasing attention, see for example [START_REF] Kakade | Efficient learning of generalized linear and single index models with isotonic regression[END_REF]; [START_REF] Chen | Generalised additive and index models with shape constraints[END_REF], [START_REF] Foster | Variable selection in monotone single-index models via the adaptive lasso[END_REF] for the high-dimensional setting, and the recent work of [START_REF] Balabdaoui | Least squares estimation in the monotone single index model[END_REF]. The model is closely related to the semi-parametric single index model which was much popularized in econometrics [START_REF] Li | Nonparametric econometrics[END_REF]. Unlike the isotonic single index model, no shape assumptions on the link function are made. However, it is necessary to impose some smoothness conditions on the unknown link function, which is typically estimated via kernel-based nonparametric methods. The most popular estimation methods are due to [START_REF] Ichimura | Semiparametric least squares (SLS) and weighted SLS estimation of single-index models[END_REF] and [START_REF] Klein | An efficient semiparametric estimator for binary response models[END_REF], which are implemented in the R package np. Despite their popularity, one of the main downsides of such approaches is the fact that they require bandwidth selection, known to be a hard problem to solve. For a review of existing estimation algorithms we refer the reader to Li and Racine (2007, Chapter 8). Unlike kernel-based methods, isotonic estimation is fully automatic and no tuning parameter needs to be selected to achieve optimal performance. In addition, [START_REF] Foster | Variable selection in monotone single-index models via the adaptive lasso[END_REF] argue that, whenever appropriate, the assumption of monotonicity "noticeably improves performance". The advantages of the monotone single index model have also been investigated in [START_REF] Qin | Using a monotone singleindex model to stabilize the propensity score in missing data problems and causal inference[END_REF]; [START_REF] Wan | Monotonic single-index models to assess drug interactions[END_REF] where applications to, respectively, causal inference and drug interaction modelling were considered. [START_REF] Qin | Using a monotone singleindex model to stabilize the propensity score in missing data problems and causal inference[END_REF] considered a similar estimator as we do here assuming continuous design, while [START_REF] Wan | Monotonic single-index models to assess drug interactions[END_REF] studied a spline approach, which has the drawback of having to use cross-validation in order to select to optimal knot positions.

To the best of our knowledge, theoretical developments for the isotonic single index model are restricted to the setting of random and continuous design where it is assumed that the covariates are generated from a multivariate density. In the recent work of [START_REF] Balabdaoui | Least squares estimation in the monotone single index model[END_REF] it is shown that the rate of convergence of the least squares estimator (LSE) of the isotonic link function in this setting is n 1/3 under some standard regularity assumptions. The methods used there cannot be adapted to datasets where repeated measurements of predictor values are observed. Such discrete design settings occur for example in drug-interaction studies, cf. [START_REF] Carter | An asymptotic confidence region for the ed 100p from the logistic response surface for a combination of agents[END_REF]; [START_REF] Wan | Monotonic single-index models to assess drug interactions[END_REF]. We note that although [START_REF] Wan | Monotonic single-index models to assess drug interactions[END_REF] considered an estimation method which is not fully automatic like the one presented here, they stressed the benefit to be gained from fitting a monotone single index model. Hence, it is cru-cial to understand the behaviour of this estimator as well as any potential technical challenges associated with the setting with repeated predictors.

Our main contributions in this paper are to give a precise description of the asymptotic behavior of the maximum likelihood estimator (MLE) in the monotone single index model under discrete design and develop an asymptotic goodness-of-fit test which could be used before applying our approach. Similarly to [START_REF] Qin | Using a monotone singleindex model to stabilize the propensity score in missing data problems and causal inference[END_REF]; [START_REF] Chen | Generalised additive and index models with shape constraints[END_REF], we study the case where the sample size n is allowed to grow while the dimension d of the predictors is held fixed. One of our main findings is that, as opposed to the classical monotone single index models where the nonparametric LSE or MLE of the link function converges at the n 1/3 -rate, the asymptotics of the MLE under the discrete design exhibit the faster rate of n 1/2 . We also investigate the finite sample performance of the MLE using synthetic data, where the MLE is found to have competitive predictive performance when compared to the MLE in the true generalized linear model from which the data are drawn.

The manuscript is organized as follows. In Section 2 we give a description of the isotonic single index model under discrete design and prove some of the fundamental properties of the MLE such as existence and lack of uniqueness of the estimator. In Section 3 we prove consistency of the MLE of the link function and give sufficient conditions under which it is asymptotically normal. Furthermore, we propose two goodness-of-fit tests. In Section 4 we illustrate through simulations the performance of the proposed method. In Section 5, following [START_REF] Wan | Monotonic single-index models to assess drug interactions[END_REF], we give an example involving drug interactions in a Leukemia dataset, where we compare the monotone single model approach with kernel-based nonparametric single-index and fully parametric methods. To keep the reading more fluid, we defer all proofs and technical details to the supplementary material.

The model and maximum likelihood estimation

The model and corresponding log-likelihood

Let n be the total size of a sample of observed responses organized into m groups. Observations in the i-th group for i = 1, . . . , m are assumed to have expectation of the form

µ 0i = f 0 (α T 0 x i ) (1)
where x i ∈ IR d is a design point, α 0 ∈ IR d \{0} is an unknown regression parameter, also called the "index", and f 0 is an unknown increasing link function. We use the term increasing to mean non-decreasing. All values of x i , i = 1, . . . , m are assumed to be distinct. We assume that the density function of the observed response in the i-th group, with respect to a given dominating measure λ which is either Lebesgue measure or counting measure, is given by

c(y, φ) exp y (µ 0i ) -B • (µ 0i ) φ , ( 2
)
where c is a normalizing function, φ > 0 is the dispersion parameter, and is a given real valued function with first derivative > 0, and inverse -1 = B . Balabdaoui et al. (2019, Proposition 5.1) study identifiability of the model when the distribution of the predictors has a density with respect to Lebesgue measure. The situation considered here is very different in that we observe f 0 (α T 0 x) with noise only at the discrete design points x = x 1 , . . . , x m . Hence, the model is not necessarily identifiable in the classical sense. However, we can show that the index parameter satisfies a weaker form of identifiability. Let S d-1 denote the unit sphere in d-dimensional Euclidean space.

Proposition 1 Assume that there exist an increasing function f 0 and α 0 ∈ S d-1 such that observations in the i-th group have expectation µ 0i of the form (1). Define the set

S 0 = α ∈ S d-1 : for all i, j, α T x i ≤ α T x j implies µ 0i ≤ µ 0j .
(3)

Consider any increasing function f and index α ∈ S d-1 such that µ 0i = f (α T x i ), i = 1, . . . , m. Then, the index α is set identifiable in the sense that α ∈ S 0 .

Heuristically, the set S 0 identifies all α which, via the projected predictors α T x i , preserve the correct order in the components of µ 0 . This can be rephrased as "S 0 is the subset of S d-1 such that α T x i < α T x j whenever µ 0i < µ 0j for i, j = 1, . . . , m". Hence, if µ 0 is constant, then S 0 = S d-1 . At the other extreme, if µ 0 satisfies µ 01 < . . . < µ 0m , then S 0 contains all α ∈ S d-1 such that α T x 1 < . . . < α T x m .

Assume that we have n i independent observations Y ij (j = 1, . . . , n i ) from group i (i = 1, . . . , m) with density as in (2) with respect to λ, and where µ 0i = f 0 (α T 0 x i ) for some unknown increasing function f 0 and index α 0 ∈ S d-1 . The overall sample size is n = m i=1 n i , and the log-likelihood is

l n (f, α) = m i=1 n i Y i {f (α T x i )} -B • {f (α T x i )} . (4)
Here,

Y i = n -1 i ni j=1
Y ij is the average of the observations from group i.

The MLE in the monotone single index model

The MLE of (f 0 , α 0 ) maximizes l n (f, α) in (4) over all increasing functions f and index vectors α ∈ S d-1 . To maximize l n we use the following profile approach: we first hold α fixed and maximize in f ; then we maximize the resulting likelihood in α.

Consider first maximization in f for a fixed α. Let M denote the class of increasing functions f : IR → IR. For k ∈ IN\{0}, vectors w ∈ IR k + and W ∈ IR k , let iso(w, W ) denote the weighted isotonic regression of W, which is the unique solution to

iso(w, W ) = argmin β∈M k k i=1 w i (W i -β i ) 2 , ( 5
)
where M k denotes the class of vectors β ∈ IR k such that β 1 ≤ . . . ≤ β k . Some well-known properties of isotonic regression, such as uniqueness, are reviewed in the supplementary material, Section A.

Proposition 2 Fix α ∈ S d-1 . Let z 1 , . . . , z k denote the unique values of α T x 1 , . . . , α T x m . The maximizer of l n (f, α) over f ∈ M exists and is unique at the points z 1 , . . . , z k . Moreover, if π denotes the unique permutation on {1, . . . , k} such that z π(1) < . . . < z π(k) , then the values of the maximizer at the points z π(1) , . . . , z π(k) are the components of iso(w, W ), where

w j = i:α T xi=z π(j)
n i , and

W j = (w j ) -1 i:α T xi=z π(j) n i Y i , for j = 1 . . . , k.
Let f n | α denote the maximizer of l n (f, α) over f ∈ M for fixed α. We now maximize l n ( f n | α , α) in α to find the MLE. The MLE of (f 0 , α 0 ) is any pair ( f n , α n ) where α n ∈ { α n }, where the set { α n } is defined in the next proposition, and

f n = f n | αn , whereas the MLE µ n of µ 0 has components µ ni = f n ( α T n x i ), i = 1, .
. . , m. The MLE exists, but is not unique, as the following propositions shows.

Proposition 3 Let { α n } denote the set of all maximizers of l n ( f n | α , α) over α ∈ S d-1 . We can find α n ∈ { α n } such that α T n x i = α T n x j for all i = j. Furthermore, { α n } has positive Lebesgue measure in S d-1 .
It is also important to note the ways in which the maximum likelihood estimator is not unique. First, the maximum likelihood estimator f n | α of the link function for a fixed α is uniquely defined only on the set of projections α T x i (i = 1, . . . , m), and any monotone interpolation of the values at these points gives a maximum likelihood estimator. More importantly, with α n as above, it can be seen from the proof of Proposition 3 that any α ∈ S d-1 such that α T x i < α T x j for all i, j satisfying α T n x i < α T n x j , maximizes the likelihood.

3 Asymptotics of the maximum likelihood estimator

Consistency

Assume that the design points x 1 , . . . , x m are fixed. Throughout, we make one of the two following assumptions about the design.

Fixed design: We assume that there are n i (i = 1, . . . , m) observations collected from group i at each design point x i (i = 1, . . . , m) with a non-random n i . We consider what happens as n increases in such a way that λ n,i = n i /n converges to some λ i ∈ (0, 1) for all i = 1, . . . , m. We call this assumption (A1).

Random design: We assume that at each design point x i there is a random number n i of observations Y ij where (n 1 , . . . , n m ) follows a multinomial distribution with parameters n and λ i (i = 1, . . . , m). This setting occurs if one observes independent copies of (X, Y ) where X is a discrete variable such that pr(X = x i ) = λ i ∈ (0, 1) for i = 1, . . . , m and Y is from the group i if and only if X = x i . Again, we consider what happens as n increases, and call this assumption (A2).

Recall the definition of S 0 given above in (3).

Proposition 4 Under either assumption (A1) or (A2), almost surely there exists a positive integer

N 0 = N 0 (ω) such that for all n ≥ N 0 , { α n } ⊂ S 0 .
This consistency result shows that any monotone single index model MLE is consistent for the set S 0 . The model is set identifiable for α 0 up to the set S 0 , see Proposition 1, and hence a stronger statement should not be expected. However, simulations in Section 4.1 indicate that the method estimates α 0 reasonably well for finite sample sizes, with efficiency similar to or only slightly worse than a fully parametric model. One important consequence of the above proposition is the following.

Proposition 5 Let µ n denote the MLE of µ 0 , and (f 0 , α 0 ) ∈ M × S d-1 such that µ 0j = f 0 (α T x j ) for j = 1, . . . , m. Then, under either assumption (A1) or (A2), almost surely, there exists a positive integer N 0 such that for all n ≥ N 0

m j=1 n j Ψ ( µ 0j -µ 0j ) ≤ m j=1 n j Ψ (Y j -µ 0j ),
for any convex function Ψ. In particular, we have

m j=1 n j ( µ 0j -µ 0j ) 2 ≤ m j=1 n j (Y j -µ 0j ) 2 .
The preceding result shows the benefits of the proposed single index model over the naive model where nothing is assumed about the true means µ 0i in (2). In the naive model, the natural estimator for µ 0i is Y i , and thus the proposition shows that our estimator has smaller risk than this naive estimator. Similar results were obtained for the discrete Grenander estimator [START_REF] Jankowski | Estimation of a discrete monotone distribution[END_REF]. Some simulations studying Proposition 5 are given in the supplementary material, Section C.3.

Limiting distribution

For a vector µ ∈ IR m we write part(µ) to denote the partition of the set {1, . . . , m} generated by the flat regions of µ. For example, suppose that m = 5 and µ 1 < µ 2 = µ 4 < µ 3 = µ 5 . Then, in this example, part(µ) = {{1}, {2, 4}, {3, 5}}. To avoid ambiguity, the indices in each J ∈ part(µ 0 ) are arranged in increasing order. Let Z = (Z 1 , . . . , Z m ) where Z i are independent N (0, 1/λ i ) random variables. (6) Also, for any vector β ∈ IR m and J ∈ part(µ 0 ), let β J = {β i , i ∈ J} be the restriction of β to J. With some abuse of notation we write µ J 0 for the constant value taken by µ 0 on a subset J ∈ part(µ 0 ). Our main theorem follows.

Theorem 1 For all α ∈ S 0 , let g| α be such that for all J ∈ part(µ 0 ), ( g| α ) J is the unique minimizer of i∈J λ i (g i -Z i ) 2 over the set of vectors g of length |J| such that g i ≤ g j whenever α T x i ≤ α T x j , for all (i, j) ∈ J × J. Then, the function

α → m i=1 λ i (( g| α ) i -Z i ) 2 (7)
achieves its minimum over S 0 . Let α denote an arbitrary minimizer. Then, under either assumption (A1) or (A2),

√ n( µ n -µ 0 ) ⇒ h| α where ( h| α ) i = (φ/ (µ 0i )) 1/2 ( g| α ) i (i = 1, . . . , m).
Despite its complex appearance, it is simple to sample from the limiting distribution described in Theorem 1, assuming that µ 0 is known.

It follows from Theorem 1 that in the important setting where all components of µ 0 differ from each other, the limiting distribution of 6). This holds in particular if the true link function f 0 is strictly monotone. Practitioners comfortable with this assumption may therefore use the Gaussian limiting distribution to obtain confidence intervals for µ 0 , as the variance can be easily estimated as φ n /λ n,i ( µ n,i ) (estimation of dispersion is reviewed in Section E.2 of the supplementary material, while comes from the exponential family assumption). We also note that the limiting distribution depends on the exponential family only through the normalization (φ/ (µ 0i )) 1/2 .

√ n( µ n -µ 0 ) is Gaussian since ( g| α ) i = Z i (i = 1, . . . , m), with Z i as in (

Goodness-of-fit

Define the null hypothesis H 0 : " µ 0 (x i ) = f 0 (α T 0 x i ), i = 1, . . . , m for some (f 0 , α 0 ) ∈ M × S d-1 ". The alternative hypothesis H 1 is that H 0 does not hold. Define the test statistic

D 2 n = n m i=1 λ n,i (Y i ) (Y i -µ ni ) 2 . ( 8
)
We consider a weighted sum of squares since the variables Y i have unequal variances.

Under assumption (A1), the variance of

Y i is φ[n i (µ 0i )] -1
, see for instance (Silvapulle and Sen, 2005, page 50). Under assumption (A2) the previous expression gives the conditional variance of Y i , (conditioning on n i , i = 1, . . . , m).

Proposition 6 Suppose that assumption (A1) or (A2) holds. Then, under H 0 the statistic D 2 n converges in distribution to

D 2 = φ m i=1 λ i (( g| α ) i -Z i ) 2 (9)
where g| α is taken from Theorem 1.

The distribution of the random variable D 2 , depends on several unknowns which require estimation. First, it involves the dispersion φ, though this can be estimated using standard techniques, see Section E.2 in the supplementary material. Second, it depends on the unknown true mean µ 0 in a nontrivial manner via S 0 and part(µ 0 ). This dependence makes it difficult to estimate in practice. To circumvent this issue, we calibrate the test using quantiles of a simpler random variable D 2 n which we conjecture to be asymptotically stochastically larger than D 2 .

To define D 2 n , consider independent mean zero Gaussian random variables Z n,1 , . . . , Z n,m with variances 1/λ n,i , i = 1, . . . , m, respectively. For all α ∈ S d-1 , let G α be the set of all vectors g ∈ IR m such that g i ≤ g j whenever α T x i ≤ α T x j , for all (i, j) ∈ {1, . . . , m} 2 , and ǧ| α be the unique minimizer of

n i=1 λ n,i (g i -Z n,i ) 2 over G α . Define D 2 n = φ n inf α∈S d-1 m i=1 λ n,i ((ǧ| α ) i -Z n,i ) 2 , ( 10
)
where φ n is any consistent estimator of the true dispersion parameter φ; see Section E.2 of the supplementary material. Clearly, as n → ∞, D 2 n converges in distribution to

D 2 = φ inf α∈S d-1 inf g∈ Gα m i=1 λ i (g i -Z i ) 2 . ( 11
)
If µ 0 is constant then part(µ 0 ) contains the single element {1, . . . , m} and S 0 = S d-1 . In this case, we have the equality D 2 = D 2 . Therefore, the asymptotic distribution of D 2 n is the same as D 2 when µ 0 is constant. In tests for monotonicity in the context of shape constrained estimation, the constant setting has previously been shown to be least favorable, see for example Ghosal et al. (2000, Inequality (2.7)) or Durot (2003, Theorem 1). We conjecture that this holds here as well; i.e. we conjecture that D 2 ≤ st D 2 .The conjecture clearly holds if µ 0 is constant. It also holds if µ 0 is such that µ 0i = µ 0j for all 1 ≤ i = j ≤ m (which is true if the link function is strictly increasing) since in this case D 2 = 0. Therefore, our conjecture holds in two key scenarios. To prove that the constant case is indeed least favorable is out of the scope of the paper.

Let q n,1-A denote the (1 -A)-quantile of D 2 n , for some fixed probability A ∈ (0, 1) The goodness-of-fit test rejects H 0 if D 2 n > q n,1-A . This guarantees a hypothesis test with asymptotic level no larger than A when µ 0 is constant and when µ 0i = µ 0j for i = j (i, j = 1, . . . , m). Moreover, the asymptotic level is exact when µ 0 is constant.

A second goodness-of-fit test is based on the following stochastic inequality.

1 2 3 4 5 6 7 8 -2 -1 0 1 2 (a) correct GLM 1 2 3 4 5 6 7 8 -2 -1 0 1 2 (b) monSIM
Fig. 1 The plot compares the estimates of α 0 for Example 1 for the correct GLM (a) and well-specified monotone single index method (b) over B = 200 simulations. Boxplots of α n,j -α 0j , j = 1, . . . , d = 2 are shown, with j appearing on the x-axis. Here, σ 2 0 = 1 and n i = 3.

Proposition 7 Let Z = (Z 1 , . . . , Z m ) T such that Z 1 , . . . , Z m are independent and Z i ∼ N (0, λ -1 i ). Also, let K 2 = Z T QZ denote a generalized chi-squared random variable, where Q = A m ΛA m , with Λ denoting the matrix with diagonal elements equal to λ i , i = 1, . . . , m and {A m } i,j = δ i,j -1/m. Then,

D 2 ≤ st φK 2 ,
where D 2 is the random variable defined in (9). When λ i = 1/m, i = 1, . . . , m, then K 2 is a chi-squared random variable with m -1 degrees of freedom.

The proof is straightforward and is provided in the supplementary material. Note that the same approach can also be used to show that D 2 ≤ st φK 2 , which is consistent with the conjecture that D 2 ≤ st D 2 . A histogram showing the extent of the difference between the distributions of D 2 , D 2 n , and φK 2 for a particular choice of design is given in Figure 2 in the the supplementary material. Section 4.2 studies the finite sample performance of the two proposed tests.

Finite sample performance

We next investigate the finite sample performance of the MLE in different scenarios. To compute the MLE, we use a modified stochastic search algorithm over α, and find f n | α (α T x j ), j = 1, . . . , m using PAVA. We use linear interpolation to obtain

f n | α (α T x).

Estimation and prediction

Example 1

This example is taken from [START_REF] Wan | Monotonic single-index models to assess drug interactions[END_REF], where it forms the basis of their simulation study. The true model is Y ij = µ 0 (x i ) + ε ij where ε ij are i.i.d. ∼ N (0, σ 2 0 ). In this example the different combinations of the levels of two drugs, d 1 and d 2 say, are used for prediction. The interaction of these drugs is modeled using the Emax model [START_REF] Wagner | Kinetics of pharmacologic response i. proposed relationships between response and drug concentration in the intact animal and man[END_REF][START_REF] Kong | A generalized response surface model with varying relative potency for assessing drug interaction[END_REF]. The dose levels for both drugs were set to (0, 0.1, 0.5, 1, 2, 4) and each dose level for drug 1 was combined with each dose level for drug 2, for a total number of interactions equal to m = 36. The true mean response curve is given by

µ 0 (x) = 100β 0 x 1 + β T 0 x
, where

x T = d 1 , d 2 , d 1 d 2 , d 1 d 2 , d 2 d 2 , d 3 1 d 2 , d 1 d 3 2 , ( 12
)
where again d i is the level of drug i, i = 1, 2. In this model, we can re-write µ 0 (x) = f 0 (α T 0 x) with α 0 = β 0 / β 0 and f 0 (t) = 100 β 0 t/(1 + β 0 t), which is a monotone increasing function. As in [START_REF] Wan | Monotonic single-index models to assess drug interactions[END_REF], we take β T 0 = (1, 1, 0.9, 0, 0, -0.3, -0.3, 0).

In our simulation study, we considered several estimation methods, as follows.

(a). Generalized linear model with inverse link function as in the true model. We call this the "correct GLM". (b). Monotone single index model with µ 0 (x) = f 0 (α T 0 x) where x is as given in ( 12). (c). Empirical approach where µ 0 (x i ) is estimated by Y i . This yields no estimates of α 0 and cannot be used in predicting the response for new predictor values. (d). A misspecified monotone single index model with µ 0 (x) = f 0 (α T 0 x) where (e). A misspecified single index model with µ 0 (x) = f 0 (α T 0 x) where x T = (d 1 , d 2 ) and no assumptions are made on f 0 . We use a kernel-based estimator calculated with the R package np [START_REF] Hayfield | Nonparametric econometrics: The np package[END_REF]. (f). A misspecified generalized linear model with inverse link function f 0 as in the true model but with µ 0 (x) = f 0 (α T 0 x) where x T = (d 1 , d 2 ). (g). A misspecified generalized linear model with link function equal to the identity and the predictor x as given in ( 12). (h). A well-specified single index model with µ 0 (x) = f 0 (α T 0 x) where the predictor x is as given in ( 12) but no assumptions are made on f 0 . We use a kernel-based estimator calculated with the R package np [START_REF] Hayfield | Nonparametric econometrics: The np package[END_REF]. (i). A misspecified generalized linear model with identity link function and

x T = (d 1 , d 2 ).
x T = (d 1 , d 2 ).
Table 1 gives the mean of the estimated α 0 values for approaches (a), (b), (f), and (h) -as these are the only approaches which enable its estimation. Although nothing beats the correct GLM approach, the monotone single index model (b) performs well. Figure 1 shows boxplots of the estimates for models (a) and (b), from which one can note the striking fact that the MLE in the monotone single index model does not exhibit a big loss efficiency when compared with the MLE in the true GLM model. Note that in Table 1 and Figure 1 we only report the case where σ 2 0 = 1 and n i = 3. The five other settings performed similarly and hence their results are not reported here. In Tables 2 and3, we consider estimation and prediction of the response surface: To study the estimation performance, we report in Table 2 the average sum of squared deviations of the different estimators from the true mean response surface at the observed predictors. The prediction performances are compared in Table 3 where we report the average sum of squared errors evaluated at 100 new and randomly selected predictor values. As expected, the correct GLM approach is best in all instances. In terms of prediction and estimation, the second best model is (b), the well-specified montone single index model -with the exception of model (c) which sometimes outperforms (b) in Table 2. Theoretically, we know from Proposition 5 that (b) will outperform (c) asymptotically. Recall also that (c) does not produce estimates of α 0 and cannot be used for prediction. We note that the kernel-based single-index approach (e) and (h) does not perform well throughout. The comparison here may not necessarily be fair, as one would not expect a nonparametric method perform well with small values of m. We include it here to highlight the importance of the monotonicity assumption in method (b). As mentioned in the introduction, Foster et al. ( 2013) also find that the assumption of monotonicity "noticeably improves performance". In fact, even the misspecified monotone model (d) outperforms methods (e)-(i) in Tables 2 and3.

Example 2

In our second example, we consider the setting of a logistic regression. We take Y ij to be independent Bernoulli variables, with probability of success given by µ 0 (x i ) = logit -1 h 0 (α T 0 x i ) and consider three different cases of h 0 . We will assume throughout that x = (x 1 , x 2 ) T , and take α T 0 = (1, 1)/ √ 2. For the first variable x 1 , the design levels are set to (0, 0.01, 0.1, 0.25) and (0, 0.01, 0.1, 1.0) for the second variable x 2 . Each level for x 1 was combined with each level for x 2 , for a total of m = 16. This is the same design as that in case study in Section 5, except that the levels have been rescaled. We take values n i = 3, 10, 100 as in Example 1 and consider the following four estimation methods. Note that the samples n i in the data set of Section 5 are close 100.

(A). Generalized linear model with µ 0 (x) = β 0 +β 1 x 1 +β 2 x 2 . This yields the estimate

µ n (x) = β n,0 + β n,1 x 1 + β n,2 x 2 = β n,0 + β T
n,1,2 x q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

Fig. 2 Plots of µ 0 (x i ), i = 1, . . . , m against α T 0 x i for cases (i)-(iii), left to right, in Example 4. Since the plot shows the response level mean, the function in case (i) is not linear after the logit transformation which we rewrite as

q 1 2 1 2 1 2 -1.0 -0.5 0.0 0.5 1.0 (A) GLM (B) monSIM (D) kerSIM q q q q q q q q 1 2 1 2 1 2 -1.0 -0.5 0.0 0.5 1.0 (A) GLM (B) monSIM (D) kerSIM q q q q q q q 1 2 1 2 1 2 -2 -1 0 1 2 (A) GLM (B) monSIM (D) kerSIM
µ n (x) = β n,0 + β n,1,2 α T n x to obtain the GLM estimate of α 0 , α n = ( β n,1 , β n,2 ) T / β n,1,2 , with β n,1,2 2 = β 2 n,1 + β 2 n,2 . (B). Monotone single index model with µ 0 (x) = f 0 (α T 0 x) . (C).
Empirical approach where µ 0 (x i ) is estimated by Y i . Again, this yields no estimates of α 0 and cannot be used for making prediction. surface at the observed predictors, while Table 6 reports the same, except evaluated at 100 randomly selected predictor values, chosen to be within the range of the fixed design. In case (i), as expected, GLM performs best in all metrics, followed by the monotone SIM. For a semi-parametric approach, the monotone SIM method performs remarkably well when compared with the correct parametric method. In case (ii), where approach (A) is misspecified, the monotone SIM performs better than GLM.

Both the GLM and monotone SIM methods also do not perform well for case (iii), as expected. As in Example 1, the kernel-based estimator of the single index model performs poorly for prediction of the response surface, despite being well-specified in all three cases.

Goodness-of-fit test

We perform a simulation study of the proposed goodness-of-fit test for all examples of the previous section, with results given in Table 7. For Example 1, Example 2 q q q q q q q q q q q q q q 0 10 25 0 10 100 0.25 0.5 0.75 q q q q q q q q q q q q q q q q 0 10 25 0 10 100 0.25 0.5 0.75 q q q q q q q q q q q q q q q q Fig. 4 From left to right: grey-scale heatmap of the response surface for logistic regression, isotonic single index model with linear contours, and isotonic single index model with elliptical contours. case (i) and (ii), the conservative nature of the test is clearly visible. Furthermore, in Example 1, Example 2 case (i), the true mean vector is strictly increasing, in which case asymptotically D 2 = 0, contributing to the effect. In Example 2 case (ii), the true mean has a region where it is flat. However, for the particular choice of covariates, this region is quite small. For this reason, we added an artificial Example 2 case (iv) where µ 0i = E[Y i |x i ] = 0.5, i = 1, . . . , m. In this example, the mean function is flat in which case we know that the asymptotic distributions D 2 and D 2 are the same. This fact that can be seen from the observed significance levels. Lastly, in Example 2 case (iii) the model is misspecified. Here, the power performs reasonably well, noting that when n i = 3, the overall sample size is n = 108. Additional simulations exploring the goodness-of-fit test for a variety of ridge functions in more depth are given in the Appendix.

When the exponential family is binomial, it is possible for small sample sizes to observe Y i equal to 0 or 1, which makes computation of D 2 n in (8) difficult, as (Y i ) = 1/0 in this case. Therefore, for the binomial family we replaced the term (Y i ) with ( µ ni ), as µ ni is smoother and hence has fewer such instances. This alteration does not change the asymptotic behaviour of the test statistic. For cases where ( µ ni ) = 1/0, we simply set the term artificially to zero (i.e. ( µ ni ) = 0). The effect of this only augments the conservative nature of the tests.

A drug interaction case study

As in [START_REF] Wan | Monotonic single-index models to assess drug interactions[END_REF], we consider a drug interaction data set. [START_REF] Carter | An asymptotic confidence region for the ed 100p from the logistic response surface for a combination of agents[END_REF] present an analysis of cytotoxicity data, where the effect of treatment combinations 0 10 25 0 10 100 0.25 0.5 0.75 q q q q q q q q q q q q q q q q 0 10 25 0 10 100 0.25 0.5 0.5 0.75 q q q q q q q q q q q q q q q q 0 10 25 0 10 100 0.25 0.5 0.75 q q q q q q q q q q q q q q q q Fig. 5 From left to right: grey-scale heatmap of the response surface µ 3 , µ 4 , µ 5 . 0 10 25 0 10 0.25 0. 5 0 .7 5 0 .9 5 q q q q q q q q q q q q 0 10 25 0 10 0.25 0.5 q q q q q q q q q q q q 0 10 25 0 10 0.25 0.5 0.75 q q q q q q q q q q q q 0 10 25 0 10 0.25 0.5 q q q q q q q q q q q q Fig. 6 Close up of the grey-scale heatmap of the response surface µ 1 (top left) µ 3 (top right), µ 4 (bottom left), and µ 5 (bottom right). q q q q q q q q q q q q q q q q 0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 1.0 death rate q q q q q q q q q q q q q q q q 0 5 10 method of [START_REF] Klein | An efficient semiparametric estimator for binary response models[END_REF] and let µ 2 (x) = f n,2 ( α T n,2 x) denote the nonparametric single index estimate. When f 0 is assumed to be increasing, we use the maximum likelihood estimator proposed in this work, and let µ 3 (x) = f n,3 ( α T n,3 x) denote the corresponding estimate. The nonparametric estimate was computed using the R package np [START_REF] Hayfield | Nonparametric econometrics: The np package[END_REF]. We also consider the monotone single index model with elliptical contours, that is, µ 4

(x) = f n,4 ( α T n,4 y) where y = (x 1 , x 2 , x 2 1 , x 1 x 2 , x 2
2 ) T , and the model µ 5 (x) = f n ( α T n,5 y), with y = (x 1 , x 2 , x 2 1 , x 2 2 ) T , where the ellipses have been constrained to have axes parallel to the coordinate axes. We do not include the intercept term as it is not identifiable. We choose the convention that f n,3 , f n,4 and f n,5 are piecewise linear.

Figure 4, from left to right, shows grey-scale heatmap plots of µ 1 (x), µ 3 (x), and µ 5 (x). The phorbol 12-myristrate, 13-acetate variable appears on the vertical axis and the transformed methylmethanesulfonate variable appears on the horizontal axis. Black dots indicate the locations of the design points, and the size of the dot is proportional to n i . While all three approaches appear to identify the phorbol 12-myristrate, 13-acetate axis as the more important direction, the heatmaps appear quite different. We also implemented the proposed goodness of fit test for models µ 3 , µ 4 , µ 5 . In all three cases, the test fails to reject at the 5% level the hypothesis that the monotone single index model is valid for the datatset. The residual sum of squares is 264. 04, 264.39, 261.82, and 261.82 for µ 1 , µ 3 , µ 4 , µ 5 respectively. The residual sum of squares is very similar for µ 1 and µ 3 , however, the contour plots away from the design points are more reasonable for µ 3 . The residual sum of squares for µ 4 and µ 5 are the same, as are the fitted values for both models. Additional plots of the models are given in Figures 5 and6.

An unusual characteristic of the parametric logistic regression is that it selects a polynomial model where a decrease in the response surface is seen for large values of phorbol 12-myristrate, 13-acetate. [START_REF] Carter | An asymptotic confidence region for the ed 100p from the logistic response surface for a combination of agents[END_REF] discuss this and suggest that this is due to lysis; High concentrations of the drug render the cells uncountable. Comparing µ 1 with µ 3 explains that this is an artefact of the parametric model, and not a true phenomenon. Figure 7, left, plots the monotone estimate f n,3 (t), dotted line. For comparison, the solid line denoting µ 1 (x) projected, with smoothing, onto the line t = α T n,3 x has been added. The observed values Y i are shown as grey dots.

We immediately see that the parabolic effect of the parametric regression is artificial and not due to lysis. Figure 7, right, similarly compares the parametric regression with the nonparametric single index regression. Here, we see f n,2 (t), dashed line, plotted with µ 1 (x) in solid, projected, with smoothing, onto the line t = α T n,2 x. The observed values Y i are again shown as grey dots. The nonparametric method identifies the methylmethanesulfonate axis as the important direction, which contradicts both the parametric and isotonic results, as well as a three dimensional visual inspection of the data set. We believe that this anomalous result is due to the small value of m. In fact, the comparison of the kernel-based is not necessarily a fair one, as this approach would not be expected to perform well with such a small design. We include it to, yet again, emphasize the benefits of the montonicity assumption in nonparametric single index models. R scripts with the analysis and additional details can be found in the supplementary material.

Conclusions

From the results obtained in the simulation study, one can observe that the main advantage of the approach proposed in this work is robustness and protection against model specification. The assumption of monotonicity of the link function is a natural one, and does not require the correct link function to be known which is a clear limitation with the parametric GLM. The only prior knowledge required from the user is the exponential family (e.g. Gaussian, binomial, Poisson) to which the distribution of the response belongs. Our simulation results show also competitiveness of the MLE. Indeed, of all estimation procedures considered, it is rivaled only by the correct parametric GLM approach, which in practice may not be known. Although the index vector in the monotone single index model with discrete design is only setidentifiable, our estimator shows remarkably small loss of efficiency against the MLE in the correct GLM. This is a very important feature which gives much incentive to apply our method in case the user is not ready to impose some parametric form on the link function.

Fig. 3

 3 Fig. 3 Estimates of α 0 for Example 2 for methods (A)-(C) over B = 200 simulations. Boxplots of α n,i -α 0i are shown, with i appearing on the x-axis. From top to bottom, the plots show cases (i)-(iii) for n i = 100.

Fig. 7

 7 Fig.7From left to right: monotone and nonparametric estimates of the ridge function compared with a projected and smoothed parametric response surface onto the line given by the estimate of the index. In both cases the single index model has linear contours.

Table 1

 1 Estimation of α 0 in Example 1 with σ 0 = 1 and n i = 3: mean values of estimated α across B = 200 simulations.

		α 1	α 2	α 3	α 4	α 5	α 6	α 7	α 8
	true α 0	0.659 0.659 0.593	0	0	-0.198 -0.918	0
	(a)	0.575 0.573 0.514 -0.010	0.016	-0.168 -0.177 -0.001
	(b)	0.568 0.567 0.305	0.105	0.130	-0.176 -0.183 -0.049
	(f)	0.132 0.132 0.775 -0.407 -0.404 -0.007 -0.008	0.184
	(h)	0.015 0.068 0.288 -0.192 -0.228	0.226	0.178	-0.216

Table 2

 2 Estimation of the response surface in Example 1: mean values of sum of squared errors at observed predictors across B = 200 simulations.

	σ 2 0	n i	(a)	(b)	(c)	(d)	(e)	(f)	(g)	(h)	(i)
		3	0.026 0.290 0.123 174.97 624.34 862.97 1826.27 3119.24 7860.14
	0.1	10	0.008 0.207 0.036 174.92 620.34 862.97 1826.25 3304.51 7860.13
		100 0.001 0.173 0.004 174.91 623.10 862.96 1826.24 3190.56 7860.13
		3	2.585 10.70 11.94 181.35 641.30 863.97 1829.15 3549.25 7861.09
	1	10	0.765	3.34	3.62	176.82 619.86 863.25 1827.11 3663.78 7860.41
		100 0.081	0.40	0.35	175.09 621.65 863.00 1826.33 3252.70 7860.16

Table 3

 3 Prediction of the response surface in Example 1: mean values of sum of squared errors at 100 randomly sampled new predictor values across B = 200 simulations.

	σ 2 0	n i	(a)	(b)	(d)	(g)	(f)	(i)	(h)	(e)
		3	0.066 20.62 268.30 498.65 3226.68 11998.56 40884.34 269698.1
	0.1	10	0.020 21.94 270.54 502.19 3214.53 12003.50 39089.48 273460.0
		100 0.002 18.33 267.26 499.00 3243.34 12062.79	36247.1	276533.0
		3	6.40	56.60 274.03 517.20 3164.58 11995.58 36222.50 255228.5
	1	10	1.94	39.46 274.13 512.07 3248.82 12124.38 39956.67 253959.1
		100	0.21	20.18 260.82 505.66 3177.71 11905.38 39552.99 270996.2

Table 5

 5 Estimation of the response surface in Example 2: mean values of sum of squared errors at observed predictors across B = 200 simulations.

	case	n i	(A)	(B)	(C)	(D)
		3	0.161 0.350 0.935 0.230
	(i)	10	0.055 0.120 0.305 0.091
		100	0.005 0.015 0.027 0.014
		3	0.496 0.397 1.026 0.540
	(ii)	10	0.364 0.140 0.303 0.202
		100	0.320 0.017 0.029 0.018
		3	0.194 0.382 0.926 0.263
	(iii)	10	0.091 0.141 0.308 0.111
		100	0.040 0.024 0.028 0.013

Table 6

 6 Prediction of the response surface in Example 2: mean values of sum of squared errors at 100 randomly sampled new predictor values across B = 200 simulations.

	case	n i	(A)	(B)	(D)
		3	0.655	1.508	4.046
	(i)	10	0.284	0.532	4.102
		100	0.022	0.221	5.493
		3	23.015 14.698	47.919
	(ii)	10	16.502	5.172	59.482
		100	14.032	1.570	119.579
		3	1.845	2.553	4.899
	(iii)	10	1.279	1.477	4.012
		100	0.938	1.167	4.065

Acknowledgements

We thank Georges Monettes and Konstatninos Ntentes from York University for helpful discussions. The second author is supported in part by MME-DII (ANR11-LBX-0023-01) and by the FP2M federation (CNRS FR 2036).

 D). A well-specified single index model with µ 0 (x) = f 0 (α T 0 x) but with no assumptions on f 0 . We use a kernel-based estimator calculated with the R package np [START_REF] Hayfield | Nonparametric econometrics: The np package[END_REF].

The three different cases for h 0 are chosen to consider a variety of scenarios.

(i). h 0 (x) = 5x. With this choice of h 0 , models (A)-(C) are well-specified. (ii). h 0 (x) = 7.5 -max{(5x -3) 2 , 6.25}. This function has both increasing and flat regions. With this choice of h 0 , models (B)-(C) are well-specified, but (A) is not.

Plots of µ 0 (x i ) under the three scenarios are given in Figure 2.

The results of our simulations follow. Table 4 gives the mean of the estimated α 0 values for approaches (A), (B), (D), and Figure 3 shows boxplots of the estimates for n i = 100 (this is the setting most similar to the real data set). Table 5 reports the mean sum of squared errors for the different estimators from the true mean response from multiple drug agents on the binary response is studied. The response counts the number of surviving human promyelocytic Leukemia cells after being exposed to a combination of doses of methylmethanesulfonate and phorbol 12-myristrate, 13-acetate. In the analysis, the dose-response surface is estimated using multivariate parametric logistic regression. The data consists of a total of 16 design points, four per agent, with anywhere from 83 -98 observations per point and an overall sample size of n = 1436 (Carter et al., 1986, Table 1, page 125).

The goal of the data analysis in [START_REF] Carter | An asymptotic confidence region for the ed 100p from the logistic response surface for a combination of agents[END_REF] was to evaluate the doseresponse surface for the combination of two drug agents. Following [START_REF] Carter | An asymptotic confidence region for the ed 100p from the logistic response surface for a combination of agents[END_REF], we similarly transform the methylmethanesulfonate variable. In their paper, the authors perform a second order polynomial logistic regression, resulting in the fitted parametric mean function µ 1 (x) = -1.33 -0.84x 1 + 0.159x 2 + 0.003883x 2 1 -0.001308x 2 2 , after backward stepwise model selection. Here, x = (x 1 , x 2 ), where x 1 is the transformed methylmethanesulfonate variable, and x 2 is the phorbol 12myristrate, 13-acetate variable.

We compare their parametric model with a single index model where µ(x) = f 0 (α T 0 x) with x = (x 1 , x 2 ) T as above, under two assumptions: f 0 is either entirely unknown or increasing. When no assumptions are made on f 0 , we use the fitting

Supplementary material

Supplementary material contains proofs and additional technical details. A collection of R scripts for the Leukemia data example of Section 5 is also available there.


Bertrand Mercier

+ ሺ1 -ݕ

ሻ ݏ

A study of the Riemann problem for diphasic water.

A study of the Riemann problem for diphasic water.

Bertrand Mercier March 2, 2021 Introduction

To study the Riemann problem we need an equation of state for water in the pure liquid case, in the diphasic case and in the pure steam case. To get simple computations, many people are considering two different stiffened gas equations of state (EOS) : one for the liquid phase and one for the steam phase, and use thermodynamic rules to combine them in the diphasic domain (see e.g. [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF]). In the present paper, we rather use a table, for the diphasic domain, giving ܶ, ሺܶሻ, ߬ ሺܶሻ, ߬ ௩ ሺܶሻ, ߝ ሺܶሻ, ߝ ௩ ሺܶሻ, ݏ ሺܶሻ, ݏ ௩ ሺܶሻ

where ሺܶሻ is the saturation pressure, ߬ ሺܶሻ, .ݏ݁ݎሺ ߬ ௩ ሺܶሻሻ is the specific volume at saturation in the liquid (resp. steam) phase, ߝ ሺܶሻ, ߝ ௩ ሺܶሻ (resp. ݏ ሺܶሻ, ݏ ௩ ሺܶሻ)are similarly the specific energy (resp. entropy) at saturation. We show how to use our table to derive an EOS in the diphasic domain. ( §1).

In §2, we show how to combine our diphasic EOS with a stiffened gas EOS in the pure liquid domain and a perfect gas EOS in the steam domain.

In §3 we address the solution of the Riemann problem with our combined equations of state. We show that the Hugoniot curves we obtain are convex, and provide a general method for solving the Riemann problem. Finally we give some specific examples in connection with depressurization.

Method A : to compute , ܶ and ,ݏ when ߬ and ߝ are given :

Let ݕ ఛ ሺܶሻ = ൫߬ -߬ ௩ ሺܶሻ൯ ቀ߬ ሺܶሻ -߬ ௩ ሺܶሻቁ ൗ ݕ ఌ ሺܶሻ = ൫ߝ -ߝ ௩ ሺܶሻ൯ ቀߝ ሺܶሻ -ߝ ௩ ሺܶሻቁ ൗ
To compute ܶ we just have to solve the equation ݕ ఛ ሺܶሻ = ݕ ఌ ሺܶሻ. This is a non linear equation with one unknown ܶ which can be easily solved by

• finding ݅ such that ݕ ఛ ሺܶ ሻ > ݕ ఌ ሺܶ ሻ and ݕ ఛ ሺܶ ାଵ ሻ < ݕ ఌ ሺܶ ାଵ ሻ ܶ ≤ ߠ ≤ ܶ ାଵ
• solving a second degree equation to find ߠ such that ܶ = ሺ1 -ߠሻܶ ାଵ + ߠܶ and then .

(Indeed, we have to solve

൫߬ -߬ ௩ ሺߠሻ൯ ቀߝ ሺߠሻ -ߝ ௩ ሺߠሻቁ = ൫ߝ -ߝ ௩ ሺߠሻ൯ ቀ߬ ሺߠሻ -߬ ௩ ሺߠሻቁ
where the functions ߬ ሺߠሻ, ߬ ௩ ሺߠሻ, ߝ ሺߠሻ and ߝ ௩ ሺߠሻ are all linear in ߠ ) • Let ݕ * denote the common value of ݕ ఛ ሺܶሻ and ݕ ఌ ሺܶሻ we let

ݏ = ݕ * ݏ ሺܶ ሻ + ሺ1 -ݕ * ሻ ݏ ௩ ሺ ܶሻ∎
Method B : to compute , ܶ and ߝ, when ߬ and ݏ are given :

In the same way, we solve

ݕ ఛ ሺܶሻ = ݕ ௦ ሺܶሻ = ݕ * where ݕ ఛ ሺܶሻ = ൫߬ -߬ ௩ ሺܶሻ൯ ቀ߬ ሺܶሻ -߬ ௩ ሺܶሻቁ ൗ ݕ ௦ ሺܶሻ = ൫ݏ -ݏ ௩ ሺܶሻ൯ ቀݏ ሺܶሻ -ݏ ௩ ሺܶሻቁ ൗ
The details are left to the reader. ∎ Method C : to compute ߝ, when ߬ and  are given : This is still easier :  being given, first we evaluate ܶ and then compute

ݕ * = ൫߬ -߬ ௩ ሺܶሻ൯ ቀ߬ ሺܶሻ -߬ ௩ ሺܶሻቁ ൗ So that ߝ = ݕ * ߝ ሺܶሻ + ሺ1 -ݕ * ሻ ߝ ௩ ሺ ܶሻ. ∎

Test of our equation of state.

We let ߝ = ݂ሺ߬, ݏሻ : a well-known result in thermodynamics (see e.g. [START_REF] Faccanoni | Étude d'un modèle fin de changement de phase liquide-vapeur[END_REF]) is that we should have To check that this is the case, we have selected ݏ = 4.4957 and 5.81494≤ ߬ ≤ 10.46689. We compute the derivative of ߝ w.r.t ߬ both by forward and backward difference. The results given in Fig. 1 show a rather good agreement that make us confident with the validity of our equation of state.

Sound speed

We notice on Fig. 1 that, if we write  = ݂ሺ߬, ݏሻ, for fixed ,ݏ ݂ሺ߬, ݏሻ is a decreasing function of ߬, so that the sound speed exists :

We have

(2) ܿ = ߬ ඥ-߲ ߲߬ ⁄
provided we use international units for each variable.

When we select ߬, ߝ as the primitive thermodynamic variables, we use that

(3) ܿ = ߬ ඥ. ߲ ߲ߝ ⁄ -߲ ߲߬ ⁄
also in international units.

In the following test we replace partial derivatives by finite differences, and we get the results given in We notice that the sound speed in a diphasic mixture is much lower than in the liquid phase, where it is of the order of 800 to 1200 m/s. This result is well known. §2 Equation of state for the liquid phase For the pure liquid phase, we shall use a stiffened gas EOS. [START_REF] Faccanoni | Étude d'un modèle fin de changement de phase liquide-vapeur[END_REF] (4)

 = ߛ- ஶ + ሺߛ -1ሻሺߝ -ݍሻ/߬ From(3), we have ܿ = ߬ ඥ. ߲ ߲ߝ ⁄ -߲ ߲߬ ⁄
From (4) we get

߲ ߲ߝ ⁄ = ሺߛ -1ሻ/߬ ߲ ߲߬ ⁄ = -ሺߛ -1ሻሺߝ -ݍሻ/߬ ଶ ܿ ଶ = ߬ ଶ ሺ. ߲ ߲ߝ ⁄ -߲ ߲߬ ⁄ ሻ = ߬ ଶ ሺ-ߛሺߛ -1ሻ  ஶ ߬ ⁄ + ߛሺߛ -1ሻሺߝ -ݍሻ/߬ ଶ ሻ = ߛ߬ሺ-ߛ ஶ +  ஶ + ሺߛ -1ሻሺߝ -ݍሻ/߬ሻ = ߛሺ +  ஶ ሻ߬ So that we get ܿ = ඥߛሺ+ ஶ ሻ߬

Remark :

In [START_REF] Faccanoni | Étude d'un modèle fin de changement de phase liquide-vapeur[END_REF] Faccanoni shows that, when ߬ and ݏ are selected as the primitive variables, the SG equation of state can be written

 = - ஶ + ሺߛ -1ሻ ߬ ିఊ exp ሺሺݏ -݉ሻ ܥ ௩ ሻ ⁄
Along the isentrope, we have

+ ஶ = ሺ +  ஶ ሻ ቀ ఛ బ ఛ ቁ ఊ so that డ డఛ  = డ డఛ ሺ +  ஶ ሻ ቀ ఛ బ ఛ ቁ ఊ = ሺ-ߛሻ߬ ିଵ ሺ +  ஶ ሻ ሺ߬ ሻ ఊ ߬ ିఊ = ሺ-ߛሻ߬ ିଵ +ሺ ஶ ሻ
Another way to compute the sound speed is to use (2) which gives

ܿ = ߬ ඥ-߲ ߲߬ ⁄ = ඥߛሺ+ ஶ ሻ߬. ∎
In what follows, we shall select  = 5.664 ܽܲܯ and ߬ = 1.3083 ݃݇/ܮ which correspond to saturated liquid water at ܶ = 545 ,ܭ and ݏ = 2.9935 kJ/kg/K . To define our equation of state we just have to select  ஶ and ߛ.

We have selected  ஶ = 186 ܽܲܯ and ߛ = 2.79, but other choices are possible.

Our sub-saturated fluid will be initially at specific volume ߬ ଵ = 1.30098 .݃݇/ܮ We complement our isentrope in the two-phase mixture domain by using the second method described in §1. We get the result shown on Fig. 3. Obviously the isentrope is continuous but there is a strong slope discontinuity between both parts. This corresponds to a strong discontinuity of the sound speed ܿ (see Fig 4). We then have ݑ ோ = ݑ = 0 and we shall assume that  ோ >  . We anticipate that we shall have a 1-shock (propagating to the left) and a 3-rarefaction wave propagating to the right. For t > 0 we shall have an intermediate constant state ݑ * ,  * , itself subdivided in 2 parts separated by a contact discontinuity. On the left (resp. on the right) of the contact discontinuity, we shall have ߬ = ߬ ଵ (resp. ߬ = ߬ ଶ ).

We have 4 unknowns ݑ * ,  * , ߬ ଶ , ߬ ଵ , and we need 4 scalar equations. First we shall use the fact that the following Riemann invariant is constant along a 3-rarefaction wave. We remind the reader that in Eulerian coordinates

߲ ݐ߲ ቆ ߬ ݑ ߝ ቇ + ൭ ݑ -߬ 0 ߬  ఛ ݑ ߬  ఌ 0 ߬ ݑ ൱ ߲ ݔ߲ ቆ ߬ ݑ ߝ ቇ = 0
Let us call λ ଵ , λ ଶ and λ ଷ the 3 eigenvalues of the matrix of this hyperbolic sytem they satisfy

ሺݑ -λሻ ଷ + ߬ ଶ .  ఛ . ሺݑ -λሻ -߬ ଶ . .  ఌ . ሺݑ -λሻ = 0 So that ሺݑ -λሻሾሺݑ -λሻ ଶ -߬ ଶ ሺ.  ఌ - ఛ ሻሿ = ሺݑ -λሻሾሺݑ -λሻ ଶ -ܿ ଶ ሿ
And we get the well-known result that λ ଵ = ݑ -ܿ, λ ଶ = ,ݑ λ ଷ = ݑ + ܿ.

Riemann invariants :

uL, pL, τL uR, pR, τR

m/s L/kg
We check that

ݎ ଷ = ൭ -߬ ܿ  ߬ ൱ is the eigenvector associated to λ ଷ indeed ൭ -ܿ -߬ 0 ߬  ఛ -ܿ ߬  ఌ 0 ߬ -ܿ ൱ ൭ -߬ ܿ  ߬ ൱ = ൭ 0 -߬ ଶ  ఛ -ܿ ଶ +   ఌ ߬ ଶ 0 ൱ = ൭ 0 0 0 ൱ A function ܴ = ܴሺ߬, ,ݑ ߝሻ is a 3-Riemann invariant iff ∇ܴ. ݎ ଷ = 0 i-e -ܴ߬ ఛ + ܿ ܴ ௨ + ߬ ܴ ఌ = 0 Then ܴ = ݑ -݃ሺ߬ሻ is a 3-Riemann invariant iff ܿ = ߬ ݃ ᇱ ሺ߬ሻ or ݃ ᇱ ሺ߬ሻ = ܿ ߬ ⁄
As a second Riemann invariant we can choose the entropy ݏ which is constant in a rarefaction wave.

Let ݏ ோ denote the entropy of the right state, we let

ܿ ோ ሺ߬ሻ = ܿሺ߬, ݏ ோ ሻ We can choose ݃ሺ߬ሻ =  ܿ ோ ሺߪሻ ߪ ⁄ ݀ߪ ఛ ఛ బ
We now get our first two equations :

(5)

ݑ * -݃ሺ߬ ଶ ሻ -൫ݑ ோ -݃ሺ߬ ோ ሻ൯ = 0 (6)  * -݂ሺ߬ ଶ ሻ = 0
Remark : 1°/ By using Method B introduced in §1, we can tabulate the isentrope associated to ݏ ோ . More precisely, we compute a 5-column table such that we find ߬, , ,ݕ ܿ and ݃ in the 5 columns. So that we have tabulated values for ݃ሺ߬ሻ, but also for  = ݂ሺ߬ሻ By assuming linear interpolation, we can also evaluate ݃ ᇱ ሺ߬ሻ and ݂ ᇱ ሺ߬ሻ. 2°/ Knowing the velocity ݑ ோ ,for all values of ߬, we can compute ݑ = ݑሺ߬ሻ by using ݑሺ߬ሻ -݃ሺ߬ሻ -൫ݑ ோ -݃ሺ߬ ோ ሻ൯ = 0.

We still call "isentrope" the so obtained curve ߬ → ሼݑሺ߬ሻ, ሺ߬ሻሽ. ∎

Hugoniot curves. Now what happens along the 1-shock ?

We have the Rankine-Hugoniot relations. Let ߪ denote the speed of the shock, we should have ( 7)

ߩ ଵ ݑ ଵ -ߩ ݑ = ߪሺߩ ଵ -ߩ ሻ (8) ሺߩ ଵ ݑ ଵ ଶ +  * ሻ -ሺߩ ݑ ଶ +  ሻ = ߪሺߩ ଵ ݑ ଵ -ߩ ݑ ሻ (9) ሺߩ ଵ ܧ ଵ +  * ሻݑ * -ሺߩ ܧ +  ሻݑ = ߪሺߩ ଵ ܧ ଵ -ߩ ܧ ሻ
where (noting

ݑ ଵ = ݑ * ሻ ܧ ଵ = ߝ ଵ + ଵ ଶ ݑ ଵ ଶ ܧ = ߝ + ଵ ଶ ݑ ଶ
Proceeding as DESPRÉS B. Numerical Methods for Eulerian and Lagrangian Conservation Laws Springer International Publishing, 2017, p.155, we introduce :

݆ = ߩ ሺߪ -ݑ ሻ = ߩ ଵ ሺߪ -ݑ ଵ ሻ So that ݑ = ߪ -݆߬ and ݑ ଵ = ߪ -݆߬ ଵ ݑ ଵ -ݑ = ݆ሺ߬ -߬ ଵ ሻ or ݆ = -ሾݑሿ/ሾ߬ሿ
Then from (8) we get

ሺߩ ଵ ݑ ଵ ଶ +  ଵ ሻ -ሺߩ ݑ ଶ +  ሻ = ߩ ଵ ߪ ݑ ଵ -ߩ ݑߪ = ሺߩ ଵ ݑ ଵ + ݆ሻ ݑ ଵ -ሺߩ ݑ + ݆ሻ ݑ = ߩ ଵ ݑ ଵ ଶ + ݑ݆ ଵ -ߩ ݑ ଶ -jݑ
So that we get :

 ଵ - = ݆ሺݑ ଵ -ݑ ሻ or ݆ = ሾሿ/ሾݑሿ
We also have

(10) ሺ ଵ - ሻ ௨ భ ା௨ ಽ ଶ = ݆ ௨ భ మ ି௨ ಽ మ ଶ
Finally from (9) we get

ሺߩ ଵ ܧ ଵ +  ଵ ሻݑ ଵ -ሺߩ ܧ +  ሻݑ = ߪሺߩ ଵ ܧ ଵ -ߩ ܧ ሻ = ሺߩ ଵ ݑ ଵ + ݆ሻܧ ଵ -ሺߩ ݑ + ݆ሻܧ  ଵ ݑ ଵ - ݑ = ݆ሺܧ ଵ -ܧ ሻ Finally ܧ ଵ -ܧ = ߝ ଵ -ߝ + ଵ ଶ ሺݑ ଵ ଶ -ݑ ଶ ሻ So that  ଵ ݑ ଵ - ݑ = ݆ሺܧ ଵ -ܧ ሻ = ݆ሺߝ ଵ -ߝ ሻ + ݆ ௨ భ మ ି௨ ಽ మ ଶ = ݆ሺߝ ଵ -ߝ ሻ + ሺ ଵ - ሻ ௨ భ ା௨ ಽ ଶ ݆ሺߝ ଵ -ߝ ሻ =  ଵ ݑ ଵ - ݑ -ሺ ଵ - ሻ ௨ భ ା௨ ಽ ଶ = ଵ ଶ ሾ ଵ ݑ ଵ - ݑ - ଵ ݑ +  ݑ ଵ ሿ = ଵ ଶ ሺݑ ଵ -ݑ ሻሺ ଵ +  ሻ = ଶ ሺ߬ -߬ ଵ ሻሺ ଵ +  ሻ Hence (11) ሺߝ ଵ -ߝ ሻ + ଵ ଶ ሺ ଵ +  ሻሺ߬ ଵ -߬ ሻ = 0
Since ߝ ଵ = ݂ሺ߬ ଵ ,  ଵ ሻ equation (11) defines a (so called Hugoniot) curve in the plane (߬, ሻ.

We denote by

 ଵ =  ுை ሺ߬ ଵ ሻ
the relation so obtained between  ଵ et ߬ ଵ . If ݆ < 0 we have a 1-shock, whereas with ݆ > 0 a 3-shock (and for ݆ = 0 we have a contact discontinuity).

We have ݆ ଶ = -ሾሿ/ሾ߬ሿ which proves that the Hugoniot curve is decreasing.∎ Remark : Lax's entropy condition for a 1-shock requires that ߪ < ݑ -ܿ ܽ݊݀ ݑ ଵ -ܿ ଵ < ߪ < ݑ ଵ For a 3-shock, we have

ݑ ଵ + ܿ ଵ < ߪ ܽ݊݀ ݑ < ߪ < ݑ + ܿ . ∎

Convexity of the Hugoniot curves

An example is shown in [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF] where the Hugoniot curves both in the ሼ߬, ሽ plane and the ሼ,ݑ ሽ plane are not convex. This is not what we find here, but we are not in the same situation, since mass exchange between phases are possible. We shall consider 2 examples :

-The first one is a case where the Hugoniot curve is crossing the saturation curve on the steam side -The second one is a case where the Hugoniot curve is crossing the saturation curve on the liquid side ݈݁݉ܽݔܧ 1 ∶ On Fig. 6, we represent the Hugoniot curve crossing the saturation line in ሼ34.53 ,݃݇/ܮ 5.664 ܽܲܯሽ.

L/kg

MPa

The steam saturation line is in grey. The Hugoniot curve is made of a blue part (in the diphasic domain) and a red part (in the pure steam domain. Both seem to be tangent. In any case the curve is convex. We have assumed that steam is a perfect gas with ߛ = 1.21 but it does not change significantly if we take ߛ = 1.4. For the pure liquid phase we have selected ߛ = 2 and  ஶ = 215,863 .ܽܲܯ The Hugoniot curve in the diphasic (resp. liquid) domain is in blue (resp. grey). The saturation curve is in red.

We see that the Hugoniot curve is continuous, but it has a slope discontinuity when it crosses the saturation curve. However it is convex. We proceed as follows :

1. We build the isentrope starting from ߬ ோ ; ܲ ோ ; ݑ ோ 2. We build the Hugoniot curve starting from ߬ ; ܲ ; ݑ

Fig. 1

 1 Fig.1 pressure p vs -߲ߝ ߲߬ ⁄

Fig 2

 2 Fig 2 Sound speed evaluated either with (2) or (3) as a function of the steam mass fraction x (there are 2 superposed curves)

Fig 3

 3 Fig 3 Isentrope in a ߬,  diagram. The SG part is shown in blue. The diphasic part in red.

Fig 4

 4 Fig 4 Sound speed ܿ along the isentrope ݏ = 2.9935 kJ/kg/K. §3 : Solution of the Riemann problem with real equation of state We shall first consider the case where we have the same diphasic fluid with two different states separated by a diaphragm which is to be removed at time t=0.

Fig. 5

 5 Fig.5 Hugoniot curve and Isentrope starting from the same point ሼ߬ ,  ሽ. On Fig. 5 we compare the isentrope passing at ߬ = 313.7083 ݃݇/ܮ ; ܲ = 0.15 ܽܲܯ ; ݑ = 0 and the Hugoniot curve starting at the same point. (ܲ ଵ in MPa and ߬ ଵ in L/kg). We notice that both curves are very close to each other around the point ሼ߬ , ܲ ሽ, but this is a well known result. Note that on Fig 5 the isentrope starts from a saturated liquid state ሼ1.3083 ,݃݇/ܮ 5.664 ܽܲܯሽ Remark : From relations ݆ = -ሾݑሿ/ሾ߬ሿ and ݆ = ሾሿ/ሾݑሿ we get that ሾሿሾ߬ሿ + ሾݑሿ ଶ = 0. The parameter ݆ is called the Lagrangian velocity of the shock.

Fig 6

 6 Fig 6 Hugoniot curve crossing the steam saturation line in ሼ 34.53 ,݃݇/ܮ 5.664 ܽܲܯሽ in ሼ߬, ሽ axes.

Fig 7 2 ∶

 72 Fig 7 Hugoniot curve in ሼ,ݑ ሽ axes for example 1. ݈݁݉ܽݔܧ 2 ∶ On Fig. 8, we represent the Hugoniot curve crossing the saturation line in ሼ1,3083 ,݃݇/ܮ 5.664 ܽܲܯሽ.

Fig 9 :

 9 Fig 9 : Graphical solution to the Riemann Problem in a diagram ሼ,ݑ ሽ Computer solution of the Riemann Problem.We proceed as follows :1. We build the isentrope starting from ߬ ோ ; ܲ ோ ; ݑ ோ 2. We build the Hugoniot curve starting from ߬ ; ܲ ; ݑ 3. We define a function ߬ ଶ → ݂ሺ߬ ଶ ሻ such that ݂ሺ߬ ଶ ሻ = ݑ ீ -ݑ where a. ሼ߬ ଶ , ܲ * ሽ is on the same isentrope as ሼ߬ ோ , ܲ ோ ሽ b. ݑ = ݃ሺ߬ ଶ ሻ + ൫ݑ ோ -݃ሺ߬ ோ ሻ൯ = 0 c. ሼ߬ ଵ , ܲ * ሽ is on the same Hugoniot curve as ሼ߬ , ܲ ሽ

3 .

 3 We define a function ߬ ଶ → ݂ሺ߬ ଶ ሻ such that ݂ሺ߬ ଶ ሻ = ݑ ீ -ݑ where a. ሼ߬ ଶ , ܲ * ሽ is on the same isentrope as ሼ߬ ோ , ܲ ோ ሽ b. ݑ = ݃ሺ߬ ଶ ሻ + ൫ݑ ோ -݃ሺ߬ ோ ሻ൯ = 0 c. ሼ߬ ଵ , ܲ * ሽ is on the same Hugoniot curve as ሼ߬ , ܲ ሽ d. ݑ ீ = ݑ -ඥሺ߬ -߬ ଵ ሻሺܲ * -ܲ ሻ 4. We use the dichotomy method to solve ݂ሺ߬ ଶ ሻ = 0 Example 3: With ߬ ோ = 1.3083 ݃݇/ܮ ; ܲ ோ = 5.664 ܽܲܯ ; ݑ ோ = 0 ߬ = 313.7083 ݃݇/ܮ ; ܲ = 0.15 ܽܲܯ ; ݑ = 0 we get |݂ሺ߬ ଶ ሻ| < 10 ିଵ in 40 steps. We get : ߬ ଶ = 49.5734 L/kg ߬ ଵ = 54.9847 L/kg ܲ * = 0.80977 MPa ݑ * = ݑ ீ = ݑ = -413.137m/sWe give below a plot of the solution of this Riemann problem at t = 2.5ms.

(

  Fig 10 Solution to the Riemann Problem at t=2.5 ms. Specific volume L/kg wrt x (m).

Fig 11

 11 Fig 11 Solution to the Riemann Problem at t=2.5 ms. Pressure in MPa wrt x (m)..

To graphically solve the Riemann problem, we just have to find the intersection in the plane ሼ,ݑ ሽ of the "isentrope" starting from the state ሼ߬ ோ ,  ோ , ݑ ோ ሽ and the "Hugoniot" starting from the state ሼ߬ ,  , ݑ ሽ.

Here is an example : We start from ሼ1.4746 ܮ ݇݃ ⁄ , ,ܽܲܯ51 0 ݉/ݏሽ on the right and ሼ600. We note that the rarefaction wave propagates relatively slowly (~50 m/s) to the right.

We also note that, on this specific case, ߬ has a weak jump at the contact discontinuity. This is due to the fact that we have selected the right state and the left state with the same entropy. We can see that the rarefaction wave is made of 2 parts :

-a fast wave propagating at about 800 m/s in the pure liquid phase which reduces its pressure from 15 MPa to the saturation pressure (11.4 MPa) on the same isentrope. -a slow wave propagating at about 50 m/s in the high pressure domain. This phenomenon explains some facts about the Chernobyl accident (see [START_REF] Mercier | A simplified analysis of the Chernobyl accident[END_REF] ).

Conclusion :

We have shown that even though we need a table with 8 columns and 99 lines (only), we obtain a sufficiently accurate EOS to reproduce the depressurization process.