A F Arteaga

K Badea

J Bekk

H Blümer

I M Bozdog

W D Apel

J C Arteaga

A F Badea

K Bekk

J Blümer

H Bozdog

I M Brancus

M Brüggemann

P Buchholz

F Cossavella

K Daumiller

V De Souza

P Doll

R Engel

J Engler

M Finger

H J Gils

R Glasstetter

C Grupen

A Haungs

D Heck

J R Hörandel

T Huege

P G Isar

K.-H Kampert

D Kang

D Kickelbick

H O Klages

Y Kolotaev

P Łuczak

H J Mathes

H J Mayer

J Milke

B Mitrica

S Nehls

J Oehlschläger

S Ostapchenko

S Over

M Petcu

T Pierog

H Rebel

M Roth

G Schatz

H Schieler

F Schröder

O Sima

M Stümpert

G Toma

H Ulrich

J Van Buren

W Walkowiak

A Weindl

J Wochele

M Wommer

J Zabierowski

D Fuhrmann

P Luczak

Energy Spectra of Elemental Groups of Cosmic Rays: Update on the KASCADE Unfolding Analysis

whether they are published or not. The documents may come

Introduction

Due to the rapidly falling intensity with increasing energy, cosmic rays of energies above 10 15 eV can be studied only indirectly by observations of extensive air showers (EAS) which are produced by the interactions of cosmic particles with nuclei of the Earth's atmosphere. The observation of a change of the power law slope [2] of the size spectrum of EAS and consequently of the all-particle energy spectrum at ∼ 3 • 10 15 eV 50 years ago has not yet been convincingly explained [3]. Several theories for the origin of the knee predict different knee positions for particles of different primary mass. Therefore the energy spectra of single elements or at least mass groups are of considerable interest. Indeed, recent analyses [4,1] find a steepening in the energy spectra of the light components in the knee region. Whereas the measurement method of detecting air showers alleviates statistical problems, one has to rely on the results of simulations and the description of hadronic interactions while reconstructing the properties of the primary particles. Since the required energy and important kinematic regions of these interactions are beyond the range of collider or fixed target experiments, the interaction models used are uncertain and differ in their predictions. On the other hand, a thorough analysis of EAS data offers the opportunity of testing [5] and improving the validity of these hadronic interaction models.

The KASCADE-Grande experiment [6], located on site of the Forschungszentrum Karlsruhe (Germany), is designed to measure EAS in the energy range between 0.5 PeV and 1 EeV. The installation consists of the original KAS-CADE [7] experiment and an extension by the Grande array, covering an effective area of 0.5 km2 .

The data analysis pursued by the KASCADE collaboration invokes an unfolding procedure of the two-dimensional shower size spectrum (total elec-tron number vs. muon number) into energy spectra of five individual mass groups [1]. Despite the success of this method for the reconstruction of the shape of the spectral forms, a strong dependence of the result for the elemental abundances on the interaction model underlying the analysis was found. Also, an insufficient description of the measured data by the employed simulations has been demonstrated. While in [1] the results for the analysis based on two different high energy interaction models were compared, we give in the present paper an update of the composition analysis concentrating on the influence of the low energy interaction model on the result. Furthermore, the analysis is repeated for measured data of different zenith angle intervals, thus testing the consistency of the analysis.

Composition analysis of the KASCADE shower size spectrum

Starting point of the analysis is the so-called two-dimensional shower size spectrum, i.e. the number of measured EAS depending on the electron number lg N e and the truncated muon number lg N tr µ (number of muons with shower core distances between 40 m and 200 m). In Fig. 1 this spectrum is given for showers inside the KASCADE array for three different ranges of inclination angle: EAS between 0 • and 18 • , between 18 • and 25.9 • , and between 25.9 • and 32.3 • , respectively.

The content N j of each histogram cell j can be written as

N j = C N A A=1 +∞ -∞ dJ A d lg E p A d lg E. ( 1
)
C is a normalizing constant (time, aperture), and the sum is carried out over all primary particle masses A. The functions p A = p A (lg N e,j , lg N tr µ,j | lg E) give the probability for an EAS with primary energy E and mass A to be measured and reconstructed with shower sizes N e,j and N tr µ,j . The probabilities p A include shower fluctuations, efficiencies, and reconstruction resolution. For reasons of clarity integration over solid angle and cell area is omitted in Eqn. 1. In case of KASCADE data p A is dominated by the shower fluctuations, whereas reconstruction systematics play an inferior role [1]. The data range was chosen in a way to minimize influences from inefficiencies.

Adopting this notation the two-dimensional shower size spectrum is regarded as a set of coupled integral equations. In the analysis the primary particles H, He, C, Si, and Fe are chosen as representatives for five mass groups. The corresponding probabilities p A are determined by Monte Carlo simulations using CORSIKA [8] and a detailed GEANT [START_REF]GEANT -Detector Description and Simulation Tool[END_REF] based simulation of the exper-ACCEPTED MANUSCRIPT iment. To solve the equation system for the mass group energy spectra, the unfolding algorithm proposed by Gold [START_REF] Gold | AEC Research and Development Report[END_REF] is applied. Details of the selection, reconstruction, and the analysis can be found in Ref. [1]. In particular, in that paper the determination and definition of the systematic uncertainties and inaccuracies of the analysis chain as well as the applicability of the chosen unfolding algorithm method is discussed in detail.

Influence of low energy hadronic interaction model

In the original analysis [1] the probabilities p A were determined using the high energy hadronic interaction models QGSJet [START_REF] Kalmykov | [END_REF] (2001 version) and SIBYLL [START_REF] Engel | Proc. 26 th Int. Cosmic Ray Conf[END_REF] (version 2.1) in the simulations. For both cases the GHEISHA [START_REF] Fesefeldt | [END_REF] code was used for interactions with laboratory energy of < 80 GeV. GHEISHA has been widely used as low energy hadronic model in air shower simulations in the last decade.

Though the high energy interaction models are believed to have a larger effect on the KASCADE observables, the influence of the low energy hadronic model ought to be investigated. In particular, the number of secondary muons at sea level could be affected because these are decay products of low energy charged mesons. However, only at large distances from the core (where KASCADE is not sensitive), the contribution from low energy pions attains importance [START_REF] Zabierowski | Proc. 29 th Int. Cosmic Ray Conf[END_REF].

ACCEPTED MANUSCRIPT

It is known [START_REF] Drescher | [END_REF]15] that the GHEISHA model suffers from deficiencies in handling the reaction kinematics properly leading to, e.g., a notably flatter lateral distribution than the FLUKA model. In addition, GHEISHA produces too many pions as compared to accelerator data, while the FLUKA (version 2002.4) [START_REF] Fassò | FLUKA: status and prospective for hadronic applications[END_REF] package is in line with the measured data.

Therefore, in the more recent analysis of the KASCADE data GHEISHA is replaced by the FLUKA model (preserving QGSJet version 2001 as high energy interaction model). Since FLUKA is known to describe the accelerator data to higher interaction energies, the change to the high energy interaction model is now at 200 GeV, i.e. FLUKA replaces in the intermediate energy range also the QGSJet model.

When performing the simulations including detector response and the stan-

ACCEPTED MANUSCRIPT

dard KASCADE reconstruction, the differences between these two models are found to be rather small, as can be seen in Fig. 2. As an example, the distributions of electron and truncated muon sizes for nearly vertical 1 PeV iron induced showers are displayed, one simulated with the combination QGSJet/GHEISHA, the other with QGSJet/FLUKA. In case of FLUKA simulations the shower size distribution is shifted by ∆ lg N e ≈ 0.01 (≈ 2% for both, proton and iron induced showers) towards larger electron numbers, and by ∆ lg N tr µ ≈ 0.02 towards (≈ 5% for iron induced EAS, but only 2% for proton induced EAS) smaller muon numbers, i.e. the differences are larger for the muon number than for the electron number. 4 For illustration in Fig. 3 the most probable values of the simulated electron-muon number correlation for the two models and for different primaries are displayed on top of the data distribution. The differences are much smaller than compared to those observed for different high-energy interaction models (see Fig. 22 in ref. [1]). This behavior is expected in view of the above mentioned differences between GHEISHA and FLUKA: GHEISHA produces a higher pion multiplicity but flatter muon lateral distributions, or in other words, FLUKA predicts fewer muons in the range of 40 m -200 m core distance.

After carrying out the complete unfolding analysis, the results for the FLUKA case differ only little from the GHEISHA case, as might have been expected from the small differences in particle numbers. In Fig. 4, left panel, the results for the energy spectra of H, He, and C using GHEISHA and FLUKA are compared with each other. It is a specific characteristic of the QGSJet model that the analysis of the KASCADE data always results in a Helium dominated composition. Differences between the two solution sets are small, especially when compared to the systematic uncertainties imposed by the unfolding procedure [1] (represented by the shaded bands in the figure, which are displayed for the QGSJet/GHEISHA solution only, but which are of the same order for all the discussed solutions). This also holds for the corresponding spectra of Si and Fe, which are compared to each other in the right panel of Fig. 4. Here, the influence of replacing the low energy interaction model on the result is larger than for the light elements. The relative abundance of the heavy elemental groups increases slightly when using FLUKA for the unfolding. This is understandable taking into account the shift of the simulated electron-muon numbers of primary iron towards the bulk of the data distribution (Fig. 3). Finally, the results for the all-particle energy spectrum (which is the sum of the individual mass group spectra) are displayed in Fig. 9 together with results discussed in the next section. No systematic difference in the all-particle spectrum could be observed by changing the low energy interaction model.

ACCEPTED MANUSCRIPT

To test the quality of the solution the obtained mass spectra are folded forward and compared to the measured data by a χ 2 -test [1]. Noticeably, the χ 2 dof parameter improves considerably from 2.38 to 1.34 by going from GHEISHA to FLUKA. In Fig. 5 the individual χ 2 j distributions for both cases are displayed, suggesting that FLUKA describes the correlation of muon to electron number in air showers better than GHEISHA.

To summarize, using the FLUKA model instead of the GHEISHA model seems to have no significant effect on the overall picture of the solution, but a slightly better description of the data is achieved. This improvement in describing the KASCADE data is not surprising considering the fact that FLUKA has been tuned to provide a good description of recent accelerator data.

ACCEPTED MANUSCRIPT 4 Investigation of different KASCADE data sets

In the analysis described so far, only EAS reaching the detector with zenith angles below 18 • were included in the composition studies. The analysis of more inclined shower data can serve as consistency check. Due to the limited reproduction of the measured observable correlations by the models, one cannot expect to obtain identical results for the energy spectra as compared to the vertical data set. Nevertheless, large differences between the solution sets for different zenith angle ranges could indicate a severe problem in either the simulation chain or the performed analysis technique. In particular, characteristics of the predictions of the attenuation of the electromagnetic and muonic shower components in the atmosphere can be tested. For this reason the analysis based on the QGSJet/FLUKA simulations were repeated for two additional data sets, the first one containing EAS with zenith angles ranging from 18 • to 25.9 • , the second one from 25.9 • to 32.3 • ; thus covering the same acceptance on the sky. Data selection is identical to the one applied to the vertical data set, i.e. the same runs in the same period and the same cuts (except for reconstructed zenith angle). The two-dimensional measured shower size spectra are included in Fig. 1.

Again, first the shower fluctuations are compared to get an impression on the differences in the data sets. In Fig. 6, as an example the distributions of electron and truncated muon sizes for 1 PeV proton induced showers are displayed for all three angular ranges. Going from nearly vertical to ≈ 30 • inclined showers, the shower size distribution is shifted by ∆ lg N e ≈ 0.5 towards smaller electron numbers (corresponding to an attenuation of the electron size relative to vertical EAS by ≈ 70%), and by ∆ lg N tr µ ≈ 0.1, (attenuation of muons relative to vertical EAS by ≈ 20%) towards smaller muon numbers. The widths of the distributions vary only slightly for the different angular ranges. Qualitatively, these differences are reflected in the two-dimensional shower size spectra shown in Fig. 1, considering the coverage of the measured distribution in particle numbers. 5By applying the full unfolding analysis for the resulting mass group spectra only small differences can be detected (Fig. 7). The display of the systematic error bands is omitted in this figure for a visibility. They are of the same order as those shown in Fig. 4 and therefore cover all three proton spectra, in particular at high energies. In the case of the proton spectra, systematic differences are present for energies above the proton knee. As can be seen, the knee in the proton spectrum gets less pronounced, i.e. the change of index decreases with increasing zenith angle. This may be related to the fact that shower fluctuations increase significantly for inclined showers and are larger also for lower primary masses.

For Helium, the most abundant group in all QGSJet based analyses [1], the spectra derived from the three data sets coincide within their statistical uncertainties. For the three heavier mass groups no significant difference within systematic and statistical accuracy could be found either.

The observed systematic differences between the solution sets are small and can be understood by the interplay of the dependence of shower sizes on energy and primary particle type, increasing shower fluctuations with increasing zenith angle, and shifted energy threshold (caused by the fixed data range in N e and N tr µ ) due to attenuation effects. The latter is also the reason for ACCEPTED MANUSCRIPT the increasing minimum primary energy of the results with increasing inclination of the incident showers. The corresponding results for the all-particle energy spectrum are depicted in Fig. 9. These coincide very well within their statistical uncertainties.

For completeness, in Fig. 8 the χ 2 j distributions for the two inclined shower data sets are shown. The χ 2 dof parameter, which is 1.34 for vertical showers improves further to 1.17 for the range 18 • -25.9 • and to 0.92 for the third range. This could be expected following the arguments above, as the QGSJet model was shown to have problems in describing the data at lower energies. When going to larger inclinations, keeping the cuts on measured particle numbers, results in a higher energy threshold. The analysis in terms of energy spectra of individual mass groups as described in [1] has been applied to KASCADE data based on a different low energy hadronic interaction model and to data sets in different intervals of the zenith angle.

Using the FLUKA model instead of the GHEISHA model has no significant effect on the overall picture of the solution, but a slightly better description of the data is achieved, i.e. i.e. the correlation of electron to muon number in the simulations with FLUKA is in better agreement with the data which results in a lower χ 2 dof distribution of the solution. The results of the unfolding analyses of the two-dimensional shower size spectrum for different zenith angular ranges show no strong or unexplainable systematic differences. Thus, the results give no hint to any severe problem in the simulation or the analysis, and reaffirm the conclusions [1] drawn from the analysis of the nearly vertical shower set: The knee is observed at an energy around ≈ 5 PeV with a change of the index ∆γ ≈ 0.4 . Considering the results of the mass group spectra, in all analyses an appearance of knee-like features in the spectra of the light elements is ascertained. In all solutions the positions of the knees in these spectra is shifted to higher energy with increasing element number.

By applying the analysis to different data sets and based on different interaction models, it has been demonstrated that unfolding methods are capable to reconstruct energy spectra of individual mass groups from air shower data, in addition to the all-particle spectrum. But still, the limiting factor of the analysis are the properties of the hadronic interaction models used and not the quality or the understanding of the KASCADE data. Furthermore, the procedure of the KASCADE data analysis, and in future also the analysis of KASCADE-Grande data measuring higher primary energies and muons at larger distances [START_REF] Zabierowski | Proc. 29 th Int. Cosmic Ray Conf[END_REF], gives valuable hints for the improvement of hadronic interaction models. The data can be confidently used when improved interaction models, based on more and extended accelerator experiments, become available.

Fig. 1 .

 1 Fig.1. Two-dimensional shower size spectra as measured by KASCADE for three different zenith angular ranges, where the first one shows similar data as used for the analysis described in[1]. The effective time of the measurements amounts to 840 days.

Fig. 2 .

 2 Fig. 2. Distribution of electron (left) and truncated muon (right) shower size for iron induced showers of 1 PeV using GHEISHA and FLUKA as low energy interaction models and in both cases QGSJet as high energy interaction model.

Fig. 3 .

 3 Fig. 3. Measured two-dimensional shower size spectrum (-18 • ) together with lines of the most probable values for proton and iron induced showers for both simulations.

Fig. 4 .

 4 Fig. 4. Comparison between QGSJet/FLUKA based results and QGSJet/GHEISHA based results for the energy spectra of H, He, and C (left) and of Si and Fe (right). Shaded bands correspond to estimates of the systematic uncertainties for the QGSJet/GHEISHA solutions.

Fig. 5 .

 5 Fig. 5. Distribution of individual χ 2 j in the data range for the QGSJet solution.The low energy interaction model used is GHEISHA in the left panel (picture taken from[1]) and FLUKA in the right panel.

Fig. 6 .

 6 Fig. 6. Distribution of electron (left) and truncated muon (right) shower size for proton induced showers of 1 PeV for different zenith angular ranges. The EAS are simulated with the QGSJet/FLUKA model combination. The lines show the parameterization of the distributions used in the analysis [1].

Fig. 7 .

 7 Fig. 7. Results for the energy spectra of primary H and He, based on the analysis of EAS data originating from different zenith angle ranges (QGSJet/FLUKA). The display of the systematic uncertainties is omitted for reasons of clarity.

Fig. 8 .

 8 Fig. 8. Distribution of individual χ 2 j for the data of two different zenith angular ranges for the QGSJet/FLUKA solution.

Fig. 9 .

 9 Fig.9. Comparison between results for the all-particle energy spectrum, using GHEISHA and FLUKA in the simulations. In the FLUKA case the data sets originate from different zenith angle intervals. Error bars show the statistical uncertainty; the shaded band gives an estimate of the systematic uncertainties for the QGSJet/GHEISHA solution, which are of similar order for the other solutions.

corresponding author, E-mail address: finger@ik.fzk.de

now at: Universidad Michoacana, Morelia, Mexico

now at: Radboud University Nijmegen, The Netherlands

Investigations of these kinds by using SIBYLL instead of QGSJet as high energy interaction models lead for 1 PeV to a similar general behavior (≈ 1% higher electron number in case of FLUKA and a smaller muon number), but interestingly for SIBYLL there is a ≈ 4% effect in muons for primary protons and only ≈ 2% for primary iron.

The variation of the particle numbers with inclination is in first order independent of hadronic interaction models as it is due to attenuation in the atmosphere, which is handled by the CORSIKA code itself rather than by the interaction models.

Acknowledgements

This work is supported by the BMBF of Germany, the Ministry of Science and Education of Poland, and the Romanian Ministry of Education and Research (grant CEEX 05-D11-79/2005).

ACCEPTED MANUSCRIPT


Anestis Antoniadis

Jéremie Bigot

Irène Gijbels

Penalized wavelet monotone regression

Keywords: Besov spaces, Constrained curve fitting, Monotonicity, Splines, Wavelets, Wavelet nonparametric regression, Wavelet thresholding. AMS 1991 subject classifications: Primary 62G07, secondary 65Dxx

In this paper we focus on nonparametric estimation of a constrained regression function using penalized wavelet regression techniques. This results into a convex optimization problem under linear constraints. Necessary and sufficient conditions for existence of a unique solution are discussed. The estimator is easily obtained via the dual formulation of the optimization problem. In particular we investigate a penalized wavelet monotone regression estimator. We establish the rate of convergence of this estimator, and illustrate its finite sample performance via a simulation study. We also compare its performance with that of a recently proposed constrained estimator. An illustration to some real data is given.

Introduction

Researchers in the physical and medical sciences are often interested in investigating an assumed monotonic relationship between an independent variable X and a dependent variable Y. Typical examples include the analysis of dose-response curves in pharmakinetics, growth curves in biology and many specific practical problems discussed in the literature cited below. Linear regression is usually too restrictive in these situations. Incorporating monotonicity constraints into the estimation of regression functions is then natural and dates back to the literature on isotonic regression. The isotonic regression model assumes that the expected value of Y is constant within disjoint regions of X and moreover that the mean levels within regions are nondecreasing in increasing X. An early exposition of this literature appears in Barlow, Bartholomew,

A c c e p t e d m a n u s c r i p t

Bremner & Brunk [START_REF] Barlow | Statistical Inference under Order Restrictions[END_REF] and later in Robertson, Wright & Dykstra [START_REF] Robertson | Order Restricted Statistical Inference[END_REF]. Consistency of monotonic regression is proved in Hanson, Pledger & Wright [START_REF] Hanson | On consistency in monotonic regression[END_REF].

Smoother estimators of monotone regression can be found in Ramsay [START_REF] Ramsay | Monotone regression splines in action (with comments)[END_REF], Kelly and Rice [START_REF] Kelly | Monotone smoothing with application to dose response curves and the assessment of synergism[END_REF], Mammen [START_REF] Mammen | Estimating a smooth monotone regression function[END_REF], Mammen and Thomas-Agnan [START_REF] Mammen | Smoothing splines and shape restrictions[END_REF], Hall and Huang [START_REF] Hall | Nonparametric kernel regression subject to monotonicity constraints[END_REF] and Mammen, Marron, Turlach and Wand [START_REF] Mammen | A general projection framework for constrained smoothing[END_REF].

The previous literature on isotonic regression is based on determining the fitted values of the estimator on a finite set of points (usually the observed covariates) and uses a set of inequality constraints to impose restrictions on the value of the regression function at these points. The algorithms used to compute these estimators can be computationally intensive and involve a large set of inequality restrictions and require a special structure of the support. Moreover, as pointed out by Gijbels [START_REF] Gijbels | Monotone regression, to appear in Encyclopedia of Statistical Sciences[END_REF] many of the constrained estimates proposed in the literature reduce the smoothness of the estimator with which they started by using isotonic regression techniques or by projecting an unconstrained curve estimate onto a constrained subspace of regression functions. As a consequence the monotone estimates appear less smooth as the unconstrained estimates and often have jump discontinuities.

Series estimators provide a convenient alternative to the isotonic regression literature. Gallant [START_REF] Gallant | On the bias inßexible functional forms and an essential unbiased form : The fourier ßexible form[END_REF] proposes the Fourier Flexible Form (FFF) estimator which is based on the trigonometric functions base. He identifies the set of restrictions on the coefficients of the FFF expansion that are sufficient to impose convexity. Monotonicity, however, cannot be easily imposed on the estimator. Especially convenient series estimators are those that are based on wavelets. The estimator used in this paper is a least-squares series estimator of the regression function based on wavelet basis functions. Given a grid of points on the support of the covariate the proposed estimator imposes restrictions on the values of the estimator at the grid points and then uses interpolation to compute the predicted values in any desired point on the support. This yields a finite number of constraints which are translated to linear inequality restrictions on the values of the coefficients of the wavelet functions. The regression function f is assumed to belong to a class of functions, F , that satisfies certain regularity conditions.

As it is often the case in nonparametric regression when the regression function is assumed to belong to a large class of functions, the wavelet based least-squares estimator is inconsistent and a penalized or regularization procedure is a good remedy to provide a consistent estimator. Regularization is also a convenient framework for taking into account shape restrictions on the estimator and this is the approach chosen in this paper.

The remainder of the paper is structured as follows. Section 2 lays the foundations for wavelet based penalized estimation of a constrained regression function, in particular of a monotone regression function. Section 3 describes the optimization procedures for computing the estimator. Asymptotic properties of the estimator are established in Section 4. Section 5 presents a short Monte Carlo study on the efficiency and rate of convergence of the estimator and on a comparison with another estimator developed in the recent literature. Some concluding remarks are given in Section 6. Appendix A recalls some facts from convex optimization while Appendix B provides the proofs of the results.

A c c e p t e d m a n u s c r i p t 2 Description of the estimator

The objective of this section is to describe a regression estimator that takes into account monotonicity constraints on the shape of the regression function but does not use functional form assumptions. This section focuses on the technical description of the estimator. A discussion on the asymptotic properties is left for the next section. For ease of presentation we restrict the discussion of the estimator to the one dimensional case. All results can be extended to a higher dimensional case but not without a considerable technical effort. In shape-restricted estimation, the curse of dimensionality has an additional effect: the number of constraints, needed to assure that the estimator satisfies certain restrictions, increases with dimensionality. An increase in the number of constraints in addition to the usual curse of dimensionality can make the problem intractable.

The following notations are used throughout the paper. For any n-dimensional vector v = (v 1 , . . . , v n ) T with components v i from [0, 1] and for any real-valued function f defined on [0, 1] define f = f (v) = ( f (v 1 ), . . . , f (v n )) T . Furthermore, for a vector v, v ≤ 0 means coordinate wise. Finally, for x, y ∈ R n , x, y denotes the usual inner product in R n .

Assume that we have noisy observations of an unknown function f : [0, 1] → R at discrete time positions x i = i n , i = 1, . . . , n:

y i = f (x i ) + σǫ i , i = 1, . . . , n, (2.1)
where ǫ i are i.i.d. normal variables with zero mean and variance 1, and σ is an unknown noise level parameter. The regression function f (•) is assumed to belong to a class of functions, F , that satisfies certain regularity conditions. These conditions are further discussed later.

In an unconstrained setting, a least-squares estimator of the regression function in (2.1) is based on the empirical analogue of the expected value of a squared loss function and is the solution to the following problem

min f ∈F 1 n n ∑ i=1 (y i -f (x i )) 2 .
(2.2)

Assuming that n is a power of 2, let W n be the discrete orthogonal wavelet transform matrix corresponding to some suitable wavelet basis of L 2 ([0, 1]). Letting θ be the vector of discrete wavelet coefficients of f , i.e. θ = W n f and d = W n y be the corresponding vector of noisy wavelet coefficients, an equivalent formulation of problem (2.2) is given by

min θ∈W 1 2 n ∑ i=1 ( di -θ i ) 2 = min θ∈W 1 2 θ -d 2 2 , (2.3)
where W denotes the set of discrete wavelets coefficients of f when f ∈ F . Recall that, for regular enough wavelets, if β i denote the continuous wavelet coefficients of f , then the corresponding discrete wavelet coefficients θ i of f are of the order √ nβ i (see e.g. Antoniadis and Fan [START_REF] Antoniadis | Regularization of Wavelets Approximations (with discussion)[END_REF]).

Often, in a nonparametric setup, the class F (and hence W) is large and the above least-squares estimator is inconsistent. In order to obtain a consistent estimator one therefore must use some kind of regularization or thresholding. Antoniadis and Fan

A c c e p t e d m

a n u s c r i p t [START_REF] Antoniadis | Regularization of Wavelets Approximations (with discussion)[END_REF] have shown that wavelet thresholding amounts to solve a penalized least-squares problem of the following form:

min θ∈R n 1 2 θ -d 2 2 + λ n ∑ i=i 0 ψ(|θ i |) , (2.4)
where ψ : R → R + is an appropriate penalty function depending on the method that is being used, and i 0 is the coarse level of smoothing. To simplify the notations, we will take hereafter i 0 = 1 without any loss of generality. We will concentrate on the following choice ψ(|x|) = |x| which yields the so-called soft thresholding estimator:

θi = 0 if | di | ≤ λ, θi = di -λ sign( di ) if | di | > λ, for i = 1, . . . , n.
Antoniadis and Fan [START_REF] Antoniadis | Regularization of Wavelets Approximations (with discussion)[END_REF] show that, for a large variety of function classes F , controlling the rate at which the penalty λ/n converges to 0 as n tends to infinity, leads to estimators that achieve an optimal rate of convergence. The regularization procedure above is also a convenient framework for taking into account monotone restrictions on the estimator. Our task in what follows is to build such a penalized wavelet estimator of the regression function that satisfies the desired monotonicity constraint.

It should be mentioned that the form of the penalty function, i.e. ψ(|x|) = |x| is not without importance. Indeed, this form will allow us to rely upon a simple dual formulation of the optimization problem. The equivalence between the original and dual optimization problem is carried by the convexity of the penalty function and by the fact that the absolute value of a wavelet coefficient is decomposed into the positive part and the minus-negative part. Extending the approach to a larger class of penalties is a challenge for future research.

Let f = W T n θ be the estimation of f at the design points. In a discrete setting, imposing shape constraints consists in determining a constrained estimate fc ∈ R n such that fc , l j ≤ 0 for j = 1, . . . , m, where the l j 's are appropriate vectors depending on the shape constraints. For instance with m = n -1, the choice l j i = 0 if i < j or i > j + 1, l j j = 1 and l j j+1 = -1 would yield an increasing estimate in the discrete sense. These restrictions, resulting in m inequality constraints (with m = (n -1)) to be imposed on the estimator can be written in matrix notation as D m fc (γ) ≤ 0, where D m denotes the m × n differentiation matrix given by

D m =        1 -1 0 0 • • • 0 0 1 -1 0 • • • 0 . . . . . . . . . . . . 0 • • • 0 1 -1 0 0 • • • 0 0 1 -1       
and where γ is the equidistant grid vector ( 1 n , 2 n , . . . , 1) T . For θ ∈ R n , let g j (θ) j=1,...,n-1 be the components of the vector defined by D m W T n θ. With the choice ψ(|x|) = |x|, wavelet regression under shape constraints can then be formulated as a convex optimization problem under linear constraints: Proof: see Appendix B.

min θ∈R n h(θ) = min θ∈R n 1 2 θ -d 2 2 + λ n ∑ i=1 |θ i | (2.5) s.t g i (θ) ≤ 0, i = 1, . . . , n -1
The next section explains how to construct the estimator.

Practical Implementation of the estimator

Note that the convex optimization problem that we consider in this paper is of the following form:

min θ∈R n 1 2 θ -d 2 2 + λ n ∑ i=1 |θ i | , (3.1) s.t Φθ ≤ 0,
for some d ∈ R n , and where Φ is an m × n real matrix. To retrieve our previous setting, take d = d and Φ = D m W T n . Since the objective function for problem (3.1) is not differentiable, the optimization results recalled in Appendix A, while useful, cannot be used directly. In this section, we show that the non-differentiable problem (3.1) falls into the category of constrained non-smooth optimization problems (CNSO) as described in Section 14.6 of Fletcher [START_REF] Fletcher | Practical methods of optimization[END_REF], and that necessary and sufficient conditions similar to the KKT equations of Appendix A can be derived.

Necessary and sufficient conditions for optimality

Following the notations in Fletcher [START_REF] Fletcher | Practical methods of optimization[END_REF], problem (3.1) can be written as a penalized constrained convex optimization problem:

min θ∈R n {h(θ) + k(c(θ))} , (3.2) s.t t(r(θ)) ≤ 0, with r(θ) = (g 1 (θ), . . . , g m (θ)) = Φθ ∈ R m , t(r) = max i=1,...,m g i (θ), h(θ) = 1 2 θ -d 2 2 , c(θ) = λθ, k(c(θ)) = c(θ) 1 .
Now, define the Lagrangian function as:

L(θ, µ, u, π) = h(θ) + µ T c(θ) + πu T r(θ) .
For non-smooth convex functions, the notion of differentiability can be extended via the introduction of the subdifferential:

∂h(θ) = {g : h(θ + δ) ≥ h(θ) + δ T g; ∀θ + δ ∈ K} .
From Fletcher [START_REF] Fletcher | Practical methods of optimization[END_REF] p. 363, we have that for θ ∈ R n , ∂ max i=1,...,n

θ i = {µ ∈ R n ; n ∑ i=1 µ i = 1, µ i ≥ 0 and θ i < max i=1,...,n θ i ⇒ µ i = 0} ∂ θ 1 = {µ ∈ R n ; |µ i | ≤ 1, and θ i = 0 ⇒ µ i = sign(θ i )} .
From Fletcher [START_REF] Fletcher | Practical methods of optimization[END_REF], Theorem 14.6.1, we have that if θ * is a local minimizer of (3.1), then there exist µ ∈ ∂k(θ * ), u ∈ ∂t(θ * ) and a real π ≥ 0 such that:

t(r(θ * )) ≤ 0, πt(r(θ)) = 0, ∇ θ L(θ * , µ, u, π) = ∇h(θ * ) + µ T ∇c(θ * ) + πu T ∇r(θ * ).
In our setting, the above necessary conditions can be written as: there exists µ ∈ R n and ν ∈ R m such that

θ * -d + λµ + Φ T ν = 0, g j (θ * ) ≤ 0, j = 1, . . . , m ν ≥ 0, g j (θ * ) < 0 ⇒ ν j = 0, j = 1, . . . , m |µ i | ≤ 1 θ * i = 0 ⇒ µ i = sign(θ * i ), i = 1, . . . , n .
Since ∇ 2 θ L(θ, µ, u, π) = I is positive definite for all θ ∈ R n , we have that the above conditions are also sufficient (see Theorem 14.6.3 in Fletcher [START_REF] Fletcher | Practical methods of optimization[END_REF]). Note that if we impose no shape constraints (i.e. ν = 0), the soft thresholding estimator satisfies the above necessary and sufficient conditions of optimality.

Dual formulation

For θ ∈ R n , let θ + i = max(θ i , 0) and θ - i = max(-θ i , 0). Note that, θ + i -θ i ≥ 0, θ - i + θ i ≥ 0 and θ + i + θ - i = |θ i |.
Let e ∈ R n be the unit vector with all entries equal to one, and consider the following constrained convex optimization problem: Note that the condition ∇ θ L(θ, µ, ν) = 0 implies that θ = dλµ -Φ T ν. By eliminating θ, the dual problem (3.4) can then be written as:

min θ,θ + ,θ -∈R n 1 2 θ -d 2 2 + λe T (θ + + θ -) , (3.3)
s.t Φθ ≤ 0, θ + -θ ≥ 0, θ + ≥ 0, θ -+ θ ≥ 0, θ -≥ 0 Proposition
max µ,ν∈R n L(µ, ν) = max µ,ν∈R n - 1 2 λ 2 µ T µ - 1 2 ν T ΦΦ T ν + λµ T d -λµ T Φ T ν + ν T Φd , s.t -e ≤ µ ≤ e, ν ≥ 0 Let B = λI Φ T (note that B is an n × (n + m)
real matrix which may be singular). Hence, solving the dual problem (3.4) is equivalent to minimizing the following leastsquares problem with box constraints:

min x∈R n+m 1 2 Bx -d 2 2 , (3.5)
s.t -e ≤ µ ≤ e, ν ≥ 0 , for x T = [µ T ν T ].
Hence, in our discrete setting, wavelet regression under shape constraints amounts to solving a least-squares minimization problem with bound constraints for which many algorithms are available.

Asymptotic properties

A large body of literature is devoted to establishing the consistency of wavelet thresholding estimators and to the rate at which they converge to the true parameter. Most of this literature treats the case where no shape restrictions are imposed on the estimator. Therefore, before using any of the results from this literature, one has to make sure that they remain intact under shape restrictions. The results presented in this section show that the monotonicity restricted estimator described in the previous section is consistent and that it achieves the same optimal rate of convergence as in the unrestricted estimation problem.

Let

R( fλ , f ) = n -1 n ∑ i=1 { fλ (x i ) -f (x i )} 2 ,
be the risk function of the unrestricted regularized wavelet estimator fλ . By translating the problem from the function space into the wavelet domain, Antoniadis and Fan [START_REF] Antoniadis | Regularization of Wavelets Approximations (with discussion)[END_REF] A c c e p t e d m a n u s c r i p t

have shown that the above estimator, with λ = σ 2 log n nearly achieves the optimal rate of convergence for a large variety of function classes F , a typical example of which are the Besov spaces. Assume that the unknown signal f is in a Besov ball. Because of simple characterization of this space via the wavelet coefficients of its members, for C > 0, the Besov space ball B r p,q (C) (r > 0, 1 ≤ p, q ≤ ∞) can be defined as

B r p,q = f ∈ L p : ∑ j 2 j(r+1/2-1/p) β j• p q < C , (4.1)
where β j• is the vector of wavelet coefficients at the resolution level j. Here, r indicates the degree of smoothness of the underlying signal f . Note again that the wavelet coefficients β in the definition of the Besov space are continuous wavelet coefficients.

They are approximately a factor of n 1/2 larger than the discrete wavelet coefficients W n f. This is equivalent to assuming that the noise level is of order 1/n.

The following theorem evaluates the rate of convergence of the maximum risk of our estimator. The proof of this theorem is provided in Appendix B. The proof consists in computing the difference between the wavelet coefficients of the constrained and the unconstrained estimate and showing that it shrinks to zero at a faster rate than the optimal convergence rate of the unconstrained estimator.

Remark.

For ease of presentation we stated the theorem for functions belonging to the Besov ball B r 1,1 (C). However, using various embedding theorems for Besov spaces (see for example Corollary 9.2, page 124 of Härdle et al. [START_REF] Härdle | Wavelets, Approximation, and Statistical Applications[END_REF]), under appropriate conditions on s, p and q the result is also true for more general Besov spaces B s p,q . We know, for example, that for any ǫ < 1/2 and q ≥ 1 we have B r-ǫ 1/(1-ǫ),q ⊆ B r 1,1 and the inclusion is continuous. From the proof of Theorem 4.1, it follows that the ∑ i |β i | is still uniformly bounded (the embedding is continuous) and the result follows again by applying results on the behavior of R( f , f ) given by Antoniadis and Fan [START_REF] Antoniadis | Regularization of Wavelets Approximations (with discussion)[END_REF] but this time for the Besov ball B r-ǫ 1/(1-ǫ),q (C).

Monte-Carlo Simulations and a Real Example

We have designed a small scale simulation study for illustrating the performance of our penalized wavelet monotone smoother and for comparing it with another simple smoothing splines based monotone smoother developed recently by Zhang [START_REF] Zhang | A simple and efficient monotone smoother using smoothing splines[END_REF]. This estimator is also monotone in a discrete sense. The main idea of the smoothing splines smoother is to shift the monotonicity constraint on the underlying regression function to the positiveness or negativeness constraint on the associated derivative curve. The interested reader is referred to the paper by Zhang [START_REF] Zhang | A simple and efficient monotone smoother using smoothing splines[END_REF], where the smoothing parameter is chosen by generalized cross validation (GCV) and closed form formulas for the estimation are derived.

A c c e p t e d m a n u s c r i p t

We investigate the regression model with an equidistant design on [0, 1], normally distributed errors, sample sizes n = 128 and n = 256 and signal to noise ratios (SNR) of 3 and 5. Signal-to-noise ratios are measured as sd(m(x))/σ, where sd(m(x)) is the estimated standard deviation of the regression function, m(x i ) over the sample i = 1, . . . , n and σ is the true standard deviation of the noise in the data. The monotone regression functions that we consider are

m 1 (x) = exp(20(x -1/2)) 1 + exp(20(x -1/2)) m 2 (x) = x + 1 6π sin(6πx), m 3 (x) = u(x) + 0.1I {x>0.5} (x),
where

u(x) = x 0 1 2 exp(- (t -0.2) 2 2(0.05 2 ) ) + 1.5 exp(- (t -0.6) 2 2 * (0.1 2 ) ) dt and m 4 (x) = xI {x<0.3} (x) + 1.5xI {x>0.6} (x) + 0.3I {x>0.3} (x) + 2x * I {x>0.8} (x) .
These functions correspond to, respectively, a function which changes several times from a strongly increasing part to a flat part, a function with a continuous "jump", a strictly increasing curve with a discontinuous jump; and an increasing function with linear parts and some flat parts. The different functions are displayed in Figure 1. We have used 100 simulations runs, for each regression function, each sample size and each signal to noise ratio. For our penalized wavelet smoother estimates we used symlets of order 8, i 0 = 3 and the universal threshold λ = σ 2 log(n) with a ro- bust estimate of σ based on the median absolute deviation of the wavelet transform

A c c e p t e d m a n u s c r i p t

at the finest scale. Note that imposing monotonicity constraints tends to averaging the oscillating parts (pseudo-Gibbs phenomema) of the unconstrained estimate. For the monotone spline smoother we have used the matlab procedure implementing the smoother, available at http://www.stat.nus.edu.sg/˜zhangjt. We have compared the constrained estimates for both methods (smoothing splines and penalized wavelet smoother). For the 100 simulations and each setting of the simulation design we calculate the pointwise mean squared error (MSE) for the two estimates mw (wavelet estimator) and ms (splines smoother). In the following, we present curves for the estimates and their MSE, where only results for the sample size n = 128, SNR=5 are displayed (see Figure 2 and Figure 3). The results corresponding to the cases n = 256 and SNR=3 are quite similar and available from the authors (see also Table 1 for results of the integrated MSE for all cases of the simulation design). The dashed curves in the figures correspond to the monotone splines smoother, whereas our wavelet estimator is plotted as a solid curve. As one can see, the monotone spline smoother has a tendency to oversmooth, leading to a large squared bias. The MSE comparison shows a substantial difference between the two estimates, mainly caused by the amount of oversmoothing of the spline smoother. The superiority of the wavelet smoother is also confirmed by the results displayed in Table 1.

As one can see from Table 1, the order of the integrated mean squared error for each estimator depends on the regression model under consideration, but whatever simulation setup is used the wavelet estimator outperforms the spline GCV-based monotone smoother. Table 1: Simulated integrated mean squared error of the estimators mw and ms for the four regression functions, based on samples of size n = 128.

To end this section, we apply our methodology to a real example concerning monotone regression, namely the fuel consumption data that relate fuel efficiency (in miles per gallon) to engine output. These data have also been used by Mammen et al. [START_REF] Mammen | A general projection framework for constrained smoothing[END_REF] to illustrate a projection-based constrained smoothing procedure. A scatter plot of the data is displayed in Figure 4. The data are available at the Statlib internet repository at Displayed in Figure 4 are also a spline smooth monotone fit (dashed curve) and a wavelet monotone fit (solid curve). The number of data points is 392 but with only 93 distincts x-coordinates. In order to apply our wavelet monotone estimate, we have used, preliminary to the estimation, the interpolation method to a fine regular grid of 128 points by Kovac and Silverman [START_REF] Kovac | Extending the scope of wavelet regression methods by coefficient-dependent thresholding[END_REF] (their makegrid procedure). The wavelet esti- mator was then computed with symlets of order 8, i 0 = 3 and the universal threshold λ = σ 2 log(n) with a robust estimate of σ based on the median absolute deviation of the wavelet transform at the finest scale. The spline monotone fit was produced by Zhang's monotone smoother using a smoothing parameter 4.06E -011 chosen by GCV. It can be seen that, overall, the difference between the two fits is small, except at the left end of the support, where many observations are available.

Concluding Remarks

A nonparametric monotone restricted estimator based on a penalized wavelet leastsquares method is described in this paper. The benefits of the penalized wavelet framework are as follows. First, assumptions like monotonicity are easily incorporated into the estimator. Second, the estimator is computed using convex programming with linear inequalities constraints. Therefore, this estimator is easy to implement. Finally, the estimator behaves well even for functions that may present jumps or discontinuities.

A drawback of our estimation procedure is that it requires that observations are sampled on an equidistant grid and that the sample size is a power of 2. However, for a deterministic design, one may use several procedures that have been developed to relax these requirements without affecting the results obtained in this paper; such as, for example, the interpolation method of Hall and Turlach [START_REF] Hall | Interpolation methods for nonlinear wavelet regression with irregularly spaced design[END_REF], the binning method of Antoniadis et al. [START_REF] Antoniadis | Random design wavelet curve smoothing[END_REF], the transformation method of Cai and Brown [START_REF] Cai | Wavelet Shrinkage for nonequispaced samples[END_REF], the isometric method of Sardy et al. [START_REF] Sardy | Wavelet shrinkage for unequally spaced data[END_REF], the interpolation method to a fine regular grid of Kovac and Silverman [START_REF] Kovac | Extending the scope of wavelet regression methods by coefficient-dependent thresholding[END_REF] and the penalized wavelet method of Antoniadis and Fan [START_REF] Antoniadis | Regularization of Wavelets Approximations (with discussion)[END_REF]. In principle, any of these methods can fit into the framework of this paper, with each interpolation method inducing a different choice of a function base for the series expansion.

For a random design one could use warped wavelets recently investigated by Kerky-

A c c e p t e d m a n u s c r i p t

acharian and Picard [START_REF] Kerkyacharian | Regression in random design and warped wavelets[END_REF]. It is shown there that for designs having a property of Muckenhoupt type, these new bases have a behavior quite similar to a regular wavelet basis, leading to estimation algorithms that mimick exactly the equi-spaced case. This would in principle enable us to prove that the associated constrained penalized procedure achieves rates of convergence which have been proved to be optimal in the uniform design case. However, it poses some challenges to extend this method to our setting and we intend to pursue this in a future work.

It should be mentioned that the general setup of Section 2 allows to treat more general constraint problems than the monotonicity constraint problem. Indeed, any restriction on the regression function that involves signing of any discrete functional of the wavelet coefficients of any order can be treated by our method. It is also possible to allow a number of such restrictions to apply simultaneously. For example, a regression function can be constrained to be both monotone and concave using the same technique.

The developed procedure produces an estimate that is monotone on the grid of design points. When the design is dense, i.e the grid can be as fine as one wishes (asymptotically), the approach is still feasible, but the algorithm can be computationally intensive since it involves then a very large set of inequality restrictions. But with a limited number of constraints, we still get a reasonable estimator. In order to get an estimated curve that satisfies the monotonicity constraints on the whole interval one would need to use a monotonicity preserving wavelet basis so that the estimator will satisfy the required restriction globally on the interval. An appropriate wavelet basis to insure this is a wavelet B-spline basis of fixed order, but such a basis is not anymore orthogonal and most of the derivations in the paper should be changed accordingly. While we believe that such a bi-orthogonal wavelet basis can fit into the framework of the paper and will lead to the same asymptotic rates, it involves hard theoretical work (insuring monotonicity, decomposing on this basis but reconstructing with the dual basis, extending the asymptotic results of Antoniadis and Fan [START_REF] Antoniadis | Regularization of Wavelets Approximations (with discussion)[END_REF] to this case). This is beyond the scope of the paper and subject of future research work. Hence, if f is in the Besov ball B r 1,1 (C) then ∑ n i=1 |β i | is uniformly bounded and we obtain that for the classical universal threshold λ = σ 2 log(n):

R( f c , f ) = 1 n θc -θ 2 2 ≤ R( f , f ) + 0( log(n) √ n )
The result now follows using Theorem 4 of Antoniadis and Fan [START_REF] Antoniadis | Regularization of Wavelets Approximations (with discussion)[END_REF].

A c c e p t e d m a n u s c r i p t

A c c e p t e d m a n u s c r i p t Proposition 2 . 1

 21 The optimization problem (2.5) has a unique solution.

A c c e p t e d m a n u s c r i p t Definition 3 . 1

 31 If f is a convex function defined on a convex set K, then the subdifferential of f at θ is defined by:

3 . 1

 31 The optimization problems(3.1) and(3.3) are equivalent in the sense that they yield the same unique solution θc ∈ R n .Proof: see Appendix B.

A c c e p t e d m a n u s c r i p t Proposition 3 . 2 2 θ -d 2 2 +

 3222 Problem (3.3) is a smooth convex optimization problem whose dual can be written as max θ,µ,ν∈R n L(θ, µ, ν) = max θ,µ,ν∈R n 1 λµ T θ + ν T Φθ , (3.4) s.t ∇ θ L(θ, µ, ν) = 0, -e ≤ µ ≤ e, ν ≥ 0 Proof: see Appendix B.

Theorem 4 . 1

 41 For the penalty function ψ(x) = |x| and for r > 1/2, the maximum risk of the restricted monotone penalized least-squares estimator f c over the Besov ball B r 1,1 (C) is of rate O(n -2r/(2r+1) log n) when the universal thresholding σ 2 log n is used.

Figure 1 :

 1 Figure 1: Signals (a): m 1 , (b): m 2 , (c): m 3 and (d): m 4 used in the simulations.

Figure 2 :

 2 Figure 2: Average estimates over 100 simulations for the 4 signals for a sample size n = 128 and a signal to noise ratio of 5. The true regresssion function is represented as a dotted curve, the estimate mw as a solid curve and the estimate ms as a dashed curve.

Figure 3 :

 3 Figure 3: Simulated mean squared error of the estimators mw (solid curves) and ms (dashed curves) for the four regression functions. The sample size is n = 128 and the signal to noise ratio is 5.

Figure 4 :

 4 Figure 4: Fuel efficiency versus Engine output.

θ 2 2 ≤ θ -θ 2 2 +

 22 following inequality holds: θcθ2 2 ≤ θθ 2 2 + 2 θ, θθc If we assume that the function f satisfies the constraints of problem (3.1) namely Φθ ≤ 0, then θ, Φ T ν = (Φθ) T ν ≤ 0 since ν ≥ 0. Hence, if f verifies the constraints that we want to impose, the following inequality holds θc -2λ θ, μcμ Given that -e ≤ μ ≤ e and -e ≤ μc ≤ e, we have that| θ, μcμ | ≤ 2 ∑ n i=1 |θ i |.Recall that the empirical wavelet coefficients are such that θ i ≈ √ nβ i where β i is the corresponding continuous wavelet coefficient.

Acknowledgement. The authors wish to thank Dr Berwin Turlach for fruitful and helpful discussions on the topic. The authors are grateful to the reviewers for their remarks which led to an improvement of the paper.

Appendix A

We recall here some classical results for convex optimization problems that we have used in Section 3. Further details can be found in the books by Rockafeller [START_REF] Rockafellar | Convex Analysis[END_REF], Fletcher [START_REF] Fletcher | Practical methods of optimization[END_REF], and Boyd and Vandenberghe [START_REF] Boyd | Convex optimization[END_REF].

Consider the following optimization problem:

For µ ∈ R m , define the Lagrangian function as:

A c c e p t e d m a n u s c r i p t

• Karush-Kuhn-Tucker conditions

If there exists a strictly feasible point x ∈ C (Slater's condition), then there exist Lagrange multipliers µ * ∈ R m , such that x * , µ * satisfy the following system of equations, described as the Karush-Kuhn-Tucker (KKT) conditions:

Note that the KKT conditions also hold if the functions c i (x) ≤ 0, i = 1, . . . , m are all linear constraints (see Theorem 9.1.1 in Fletcher [START_REF] Fletcher | Practical methods of optimization[END_REF]) without assuming Slater's condition. These conditions are also sufficient if Slater's condition holds.

• Duality

Define the Lagrange dual function as:

Note that g is a convex function. Then, for any µ ≥ 0 which is dual feasible i.e. Note that if x * is a solution of the primal problem (A.1) and if x * , µ * satisfy the above KKT conditions, then under appropriate assumptions (e.g. if the functions c i (x) ≤ 0, i = 1, . . . , m are all linear constraints, see Theorem 9.5.1 in Fletcher [START_REF] Fletcher | Practical methods of optimization[END_REF]) x * , µ * solves the dual problem: max

and p * = L(x * , µ * ).

Appendix B. Proofs

Proof of Proposition 2.1:

be the unique minimizers of problems (3.1) and (3.3) respectively. Let θ

which finally implies that θ 1 = θ 2 and completes the proof.

Proof of Proposition 3.2: by definition, the dual of problem (3.3) is:


Marie-Francoise Barme-Delcroix

Ursula Gather

Limit laws for multidimensional extremes

Keywords: isobars, order statistics, extreme values, stability, limit laws, concomitant classification : 62 G 30, 60 D 05

A definition of the ordering of a multivariate sample -based on the isobar surfaces -is used in order to obtain limit laws for the extreme values of a multidimensional sample.

A c c e p t e d m a n u s c r i p t 1 Introduction

Today many papers on extreme values deal with multivariate extreme values, see [Nadarajah,Anderson,Tawn 98], [Schlatler, Tawn 03], [Fougères 04], [Heffernan,Tawn 04], [Falk,Reiss 05]; see also for recent developments about max-stable processes and multivariate extremal index, [Smith, Weissman 96],

[Zhang, Smith 04]. Unfortunately, ordering multivariate data can be done in various ways and many definitions have been proposed e.g. [Barnett 76], [Maller 90]. Some papers of [Abdous,Theodorescu 92], [Einmahl,Mason 92], [Massé,Theodorescu 94]

Preliminaries

We consider random variables defined on a probability space (Ω, A, P) and with values in the Euclidean space R k .

For every

x in R k \ {0} we define a pair ( x , x x ) = (r, θ) in R + * × S k-1
, where • is the Euclidean norm and R + * is the set of the positive real numbers. The unit sphere S k-1 in R k is endowed with the induced topology of R k .

For each random variable X = (R, Θ) = ( X , X X ) ∈ R k , we assume that the distribution of Θ, and for all θ, the distribution of R given Θ = θ, have a continuous density. F θ denotes the continuous conditional distribution function of R, given Θ = θ and F -1 θ the generalized inverse of F θ .

For each 0 < u < 1, we call the mapping θ F -1 θ (u) a u-level isobar from the distribution of R given Θ = θ . We suppose that this mapping is continuous and strictly positive; the surface given by ρ = F -1 θ (u) is also called an isobar, see Figure 1.

Let E n = (X 1 , . . . , X n ) be a sample of independent random variables in R k each with the same distribution as X. Clearly for each 1 ≤ k ≤ n there is almost surely a unique u k -level isobar from the distribution of R given

Θ = θ which contains (R k , Θ k ) = ( X k , X k X k ).
X n,1 , • • • , X n,i , • • • , X n,n-1 , X n,n = X * n ... ... X 0 X * n = X n,n X n,i X n,n-1 x x x
P (R * n ≤ r/Θ * n = θ) = (F θ (r)) n . Proof : see [Delcroix, Jacob 91].
The second assertion of the previous theorem shows that this multidimensional problem can be reduced to a unidimensional one . That will be done in the next section to find limit laws for the pair (R * n , Θ * n ). In a third section, we will focus on limit laws for the variables R * n . From the previous theorem, we know that the distributions of (R * n , Θ * n ) and (R, Θ) possess the same set of isobars. So it has been possible in [Delcroix,Jacob 91] to give the following definition.

Definition 1 (X * n ) n is called stable in probability (resp. a.s.) if and only if there is a sequence (g n ) n of isobars satisfying

(1) R * n -g n (Θ * n ) P ---→ (a.s.) 0 .
We are well aware that the above definition of multivariate maxima depends on the underlying distribution. This may be crucial in practical situations when given data are to be ordered 1 . To find out more about the limit distribution of this maximum however it is not a severe additional obstacle. Note, that even to check the type of extreme value distribution in one dimension affords at least the knowledge of the tail of the underlying distribution.

1 Another work is in progress on the estimation of isobars.

A c c e p t e d m a n u s c r i p t

3 Limit laws for some functionals of (R * n , Θ * n )

We know the three classical types of limit distributions for unidimensional extreme values and their domains of attraction, [Gnedenko 43], [Leadbetter 83],

[Galambos 87]. Recall the three corresponding distribution functions with the Galambos notations :

i)G 3,0 (x) = exp(-e -x ) for all x. ii)G 1,α (x) = exp(-x -α ) for x > 0, with α > 0. iii)G 2,α (x) = exp(-(-x) α )1 (x≤0) + 1 (x>0) , with α > 0.
Similarly, in this section, the goal is to find a nondegenerate distribution function H and sequences a n (.) > 0 and b n (.) such that for all x

P {a n (Θ * n )(R * n -b n (Θ * n )) ≤ x} w → H(x)
where w → denotes convergence at continuity points of H.

Let g n (θ) denote the (1

-1 n )-level isobar; that is g n (θ) = F -1 θ (1 -1 n ).
Under given conditions, [Delcroix, Jacob 91], we have proved that these isobars stabilize the sequence (R * n , Θ * n ), according to the definition given in the first section.

Theorem 2 Suppose that for all θ, sup

x {F θ (x) < 1} = +∞. If for all θ the distribution function F θ satisfies lim t→∞ 1 -F θ (tx) 1 -F θ (t) = x -α(θ) f or all x > 0,
where α is a strictly positive function, then for all x > 0

P {R * n ≤ x g n (Θ * n )} → S k-1 exp(-x -α(θ) )P Θ (dθ).

A c c e p t e d m a n u s c r i p t

Proof : From Theorem 1.6.2 of [Leadbetter 83], for all θ there exist sequences a n and b n depending on θ such that

F θ n x a n (θ) + b n (θ) → exp(-x -α(θ) ).
Then from the Dominated Convergence Theorem we have

S k-1 F θ n x a n (θ) + b n (θ) P Θ (dθ) → S k-1 exp(-x -α(θ) )P Θ (dθ).
But, we know that

F n θ is the conditional distribution function of R * n given Θ * n = θ, so S k-1 P (R * n ≤ x a n (θ) + b n (θ)/Θ * n = θ)P Θ (dθ) → S k-1 exp(-x -α(θ) )P Θ (dθ).
And since Θ * n and Θ have the same distribution,

P (R * n ≤ x a n (Θ * n ) + b n (Θ * n )) → S k-1 exp(-x -α(θ) )P Θ (dθ).
Hence

P (a n (Θ * n )(R * n -b n (Θ * n )) ≤ x) → S k-1
exp(-x -α(θ) )P Θ (dθ).

Moreover from Corollary 1.6.3. of [Leadbetter 83] we know that for all θ b n = 0 and (a n (θ)

) -1 = F -1 θ (1 -1 n ) = g n (θ). Example 1 If F θ (x) = (1 -K(θ)x -α(θ) )1 " x>(K(θ)) 1 α(θ)
« where α and K are strictly and continuous functions, then for all x > 0

P {R * n ≤ x(nK(Θ n,n )) 1/ α(Θ * n ) } → S k-1 exp(-x -α(θ) )P Θ (dθ).
For the two others domains of attraction it is possible to give similar limit laws according to the results of [Leadbetter 83]. For all θ we use the notation

x F θ = sup(x : F θ (x) < 1).

A c c e p t e d m a n u s c r i p t

Theorem 3 If for all θ x F θ < +∞ and if for all x > 0 and for all θ there exists a strictly positive function α such that

lim h→0 1 -F θ (x F θ -xh) 1 -F θ (x F θ -h) = x α(θ) then P { x F Θ * n -g n (Θ * n ) -1 (R * n -x F Θ * n ) ≤ x} w → S k-1 exp(-(-x) α(θ) )1 (x≤0) + 1 (x>0) P Θ (dθ).
Theorem 4 If for all θ there exists a strictly positive function l θ such that for all θ and for all x

lim t→x F θ 1 -F θ (t + xl θ (t)) 1 -F θ (t) = e -x
then for all x

P { l Θ * n (g n (Θ * n )) -1 (R * n -g n (Θ * n )) ≤ x} → S k-1 exp(-e -x )P Θ (dθ) = exp(-e -x ). Example 2 If F θ (x) = (1 -e -α(θ)x )1 (x>0)
where α is a strictly positive continuous function, then for all

x P α(Θ * n ) R * n - Log(n) α(Θ * n ) ≤ x → exp(-e -x ).
Example 3 For a bivariate Gaussian sample with covariance matrix

⎛ ⎝ σ 2 0 0 τ 2 ⎞ ⎠ , the density over R 2 is f (x, y) = 1 2πστ e -x 2 2σ 2 -y 2 2τ 2 ; the conditional distribu- tion of R given Θ = θ is then given by F θ (x) = (1 -e -x 2 Φ(θ) )1 (x>0)
, where

Φ(θ) = cos 2 (θ) 2σ 2 + sin 2 (θ) 2τ 2 .
Note that in this case, the isobars are also the level curves of the bivariate density. The condition of the previous theorem is satisfied and for all x, By definition,

P 2 Φ(Θ * n )Log(n) R * n - Log(n) Φ(Θ * n ) ≤ x → exp(-e -x ).
R * n is the concomitant of U n,n = Max(U 1 , . . . , U n ) where for all 1 ≤ k ≤ n, U k = F Θ k (R k ); see e.
Theorem 5 a) (U 1 = F Θ 1 (R 1 ), . . . , U n = F Θn (R n )) is a sample from the uni- form distribution on [0, 1]. b) U k is independant of Θ k for all 1 ≤ k ≤ n. Proof : see [Delcroix, Jacob 91].
Note that the uniform distribution belongs to the type III -according to the notation of [Leadbetter 83]-of extreme values with α = 1. It will be useful for the following. From Galambos results on concomitants it is possible to state this limit law for R * n .

Theorem 6 If there exist A n , B n > 0 such that lim n→∞

P {R ≤ A n + B n u/U = 1 + z/n} = T (u, z) (for -n ≤ z ≤ 0 and u > 0) is a non degenerate distribution function, then lim n→∞ P {R * n < A n + B n u} = T (u)
where

T (u) = 0 -∞ T (u, z)e z dz.
Proof : We have

P (R * n ≤ r) = [0,1] P (R * n ≤ r/U n,n = u)P Un,n (du). But (R * n , U n,n ) is a point of the sample ((R 1 , U 1 ), • • • , (R n , U n )) and the (R k , U k ) are independent, so P (R * n ≤ r) = [0,1] P (R ≤ r/U = u)nF 1 n-1 (u)f 1 (u)du

A c c e p t e d m a n u s c r i p t

where F 1 and f 1 denote the uniform distribution function over [0, 1] and the uniform density function respectively. Now, put u = 1 + z n . We obtain

P (R * n ≤ r) = [-n,0] P R ≤ r/U = 1 + z n 1 + z n n-1 dz.
And finally put r = A n + B n u, hence

P (R * n ≤ r) = [-n,0] P R ≤ A n + B n u/U = 1 + z n 1 + z n n-1
dz.

We conclude with the Dominated Convergence Theorem.

Remark 2

The condition of the previous theorem may be written more easily with the conditional law of R given U = u if we suppose F θ one-to-one . Since U = F Θ (R), we have

P (R ≤ r/U = u) = P (U ≤ F Θ (r)/U = u).
And since the variables U and Θ are independent, we obtain

P (R ≤ r/U = u) = P (u ≤ F Θ (r)).
So the condition of the theorem is

P 1 + z n ≤ F Θ (A n + B n u) → T (u, z).
Of course this condition depends heavily on the form of the conditional dis- Math. 44,

tribution of R given Θ = θ. Example 4 If F θ (x) = (1 -e α(θ)x )1 (

  , [ De Haan, Huang 95], and more recently, of [Berlinet, Gannoun, Matzner-Lober 01], [Serfling 02] develop the notion of multivariate quantiles. In the classical scheme (cartesian coordinates) the multivariate variables are ordered coordinate by coordinate -see for example [Galambos 87] and the references therein. And in this way the maximum value thus obtained is not a sample point.In this paper we will use the definition of the maximum value of a multidimensional sample, given in[Delcroix, Jacob 91]. This definition is based on surfaces we have called isobars. Our approach is more geometric and as recalled just below, we use the level-surfaces of the conditional distribution function of the radius R given the angle Θ. So the maximum value is a sample point. The first motivation, [Barme 93], was to find the asymptotic location of a multidimensional sample without using the convex hull of the sample as it is done classically in [Geffroy 61], [Geffroy 58]. By a unidimensional approach, stability, strong behaviour and outlier-proneness properties have been explored for this new notion of extreme value in[Barme, Brito 01] and[Barme, Gather 02]. In this paper, the goal is to find limit laws for such mul-Through the use of the isobar surfaces we will see that the classical results -[Gnedenko 43] -about limit laws for unidimensional extremes are still efficient in this multidimensional context.

  Figure 1: u-level isobar

Figure

  Figure 2: isobar-maximum

  It has been shown in [Delcroix, Jacob 91] that the conditional distribution of R * n given Θ * n is F n θ . More exactly we have the following theorem. Theorem 1 a)The variables Θ * n and Θ have the same distribution. b)The conditional distribution function of R * n given Θ * n = θ is given by

  consider the pairs (R k , U k ). In this part we use the following remark.

  g. [Galambos 87], [David 81]. The two properties of the next theorem will be used in this section.

  x>0) and α(Θ) is uniformly distributed over [0, 1] then the previous theorem is true with A n = Log n, B n = Log n, the limit law is Pareto.

  It would be possible to consider the pairs (U k , Θ k ) as done in the last section for the variables U and R, but we obtain only the asymptotic independence between Θ and U and the fact that asymptotically the variables Θ * n and Θ have the same distribution, which is weaker than the results stated in Theorems 1 and 5.

  Gannoun, Matzner-Lober 01] BERLINET A., GANNOUN A., MATZNER-LOBER E. (2001) -Asymptotic normality of convergent estimates of conditional quantiles. Statistics 35

  d'appui d'un échantillon laplacien à k dimension, Publ. Inst. Stat. Univ. Paris, 1961, no59. [Gnedenko 43] GNEDENKO B.V. (1943) -Sur la distribution limite du terme maximum d'une série aléatoire. Ann.

[Zhang, Smith 04] ZHANG Z. and SMITH R.L. ( 2004) -The behavior of multivariate maxima of moving maxima processes. J. Appl. Probab. 41, 1113-1123.


Johannes Forkman
email: johannes.forkman@bt.slu.se

Steve Verrill

The distribution of McKay's approximation for the coeffcient of variation

Keywords: Coefficient of variation, McKay's approximation, Noncentral beta distribution

McKay's chi-square approximation for the coefficient of variation is type II noncentral beta distributed and asymptotically normal with mean n -1 and variance smaller than 2(n -1).

A c c e p t e d m a n u s c r i p t 1 Introduction

The coefficient of variation is defined as the standard deviation divided by the mean.

This measure, which is commonly expressed as a percentage, is widely used since it is often necessary to relate the size of the variation to the level of the observations. McKay (1932) introduced a χ 2 approximation for the coefficient of variation calculated on normally distributed observations. It can be defined in the following way.

Definition 1. Let y j , j = 1, ..., n, be n independent observations from a normal distribution with expected value µ and variance σ 2 . Let γ denote the population coefficient of variation, i.e. γ = σ/µ, and let c denote the sample coefficient of variation, i.e.

c = 1 m 1 n -1 n j=1 (y j -m) 2 , m = 1 n n j=1 y j .
McKay's approximation K n is defined as

K n = 1 + 1 γ 2 (n -1) c 2 1 + (n-1)c 2 n . ( 1
)
As pointed out by [START_REF] Umphrey | A Comment on McKay's Approximation for the Coefficient of Variation[END_REF]

The distribution of McKay's approximation

If U and V are independent central χ 2 distributed random variables with u and v degrees of freedom respectively, the ratio R = V/ (U + V ) is beta distributed with v/2 and u/2 degrees of freedom respectively. If V is instead a noncentral χ 2 distributed random variable the ratio R is noncentral beta distributed [START_REF] Johnson | Distributions in Statistics: Continuous Univariate Distributions-2[END_REF]. In this case [START_REF] Chattamvelli | A Note on the Noncentral Beta Distribution Function[END_REF] calls the distribution of R the type I noncentral beta distribution and the distribution of 1 -R the type II noncentral beta distribution. We shall in agreement with [START_REF] Chattamvelli | A Note on the Noncentral Beta Distribution Function[END_REF]

U/ (U + V ).
The following theorem states that the random variable K n , claimed by [START_REF] Mckay | Distribution of the Coefficient of Variation and the Extended "t" Distribution[END_REF] to be approximately χ 2 distributed, is type II noncentral beta distributed.

Theorem 3. The distribution of McKay's approximation K n , as defined in Definition 1,

is n 1 + 1 γ 2 Beta II n -1 2 , 1 2 , n γ 2 . ( 2
)
Proof. Let s denote the standard deviation, i.e. s = cm. Then the second factor in ( 1)

can be written (n -1) c 2 1 + (n-1)c 2 n = n j=1 (y j -m) 2 m 2 + 1 n n j=1 (y j -m) 2 = n n j=1 (y j -m) 2 n j=1 (y j -m) 2 + n j=1 m 2 = nU U + V , where U = n j=1 (y j -m) 2 /σ 2 and V = n j=1 m 2 /σ 2 .
Here U is central χ 2 distributed with n-1 degrees of freedom. The average m is normally distributed with expected value µ and variance σ 2 /n. Consequently nm 2 /σ 2 , i.e. V , is χ 2 distributed with 1 degree of freedom and noncentrality parameter nµ 2 /σ 2 = n/γ 2 . Since the sums of squares n j=1 (y j -m) 2 and n j=1 m 2 are independent the theorem follows.

It is well known that

√ n/c is noncentral t distributed with n -1 degrees of freedom and noncentrality parameter √ n/γ (e.g. [START_REF] Owen | A Survey of Properties and Applications of the Noncentral t-Distribution[END_REF]). Theorem 3 is easily proven from this starting point as well.

We also note that the factor n(1 + 1/γ 2 ) in ( 2) is the expected value of U + V as defined in the proof of Theorem 3. This observation suggests application of the law of large numbers when investigating the convergence of McKay's approximation.

A c c e p t e d m a n u s c r i p t

Theorem 4. The distribution of McKay's approximation K n as defined in Definition 1, equals the distribution of U n W n , where U n is a central χ 2 distributed random variable with n -1 degrees of freedom and W n is a random variable that converges in probability to 1.

Proof. Let Z k , k = 1, 2, . . . , n -1, be independent standardized normal random variables.

Then

U n n -1 d = 1 n -1 n-1 k=1 Z 2 k ,
which converges almost surely to 1. Let also Z denote a standardized normal random variable, and let

V n n = 1 n Z + √ n γ 2 = Z 2 n + 2Z γ √ n + 1 γ 2 .
V n /n converges in probability to 1/γ 2 . Thus

U n + V n n p → 1 + 1 γ 2 . ( 3
)
By Theorem 3

K n d = n 1 + 1 γ 2 U n U n + V n = W n U n ,
where 3), converges in probability to 1.

W n = n(1 + 1/γ 2 )/ (U n + V n ), by (
Given Theorem 4 one might assume that McKay's approximation K n is asymptotically normal with mean n -1 and variance 2(n -1). Instead the following result holds.

Theorem 5. Let K n be McKay's approximation and γ the coefficient of variation as defined in Definition 1. Then

K n -(n -1) 2(n -1) d → N 0, 1 + 2γ 2 1 + 2γ 2 + γ 4 .
(4)

Proof. Let Z denote a standardized normal random variable, and let

V n = (Z + √ n/γ) 2 .
Let U n be a central χ 2 distributed random variable with n -1 degrees of freedom, independent of V n . Then, by Theorem 3,

K n -(n -1) 2(n -1) d = 1 2(n -1) n(1 + 1/γ 2 )U n U n + V n -(n -1) = A n B n , (5)

A c c e p t e d m a n u s c r i p t

where, by (3),

A n = n U n + V n p → γ 2 1 + γ 2 (6)
and

B n = 1 2(n -1) U n (γ 2 + 1) γ 2 - (n -1)(U n + V n ) n .
We obtain

B n = C n + D n + E n + F n (7)
where

C n = U n -(n -1) γ 2 2(n -1) d → N 0, 1 γ 4 , ( 8
)
D n = - √ 2(n -1)Z γ n(n -1) d → N 0, 2 γ 2 , ( 9
)
E n = U n n 2(n -1) p → 0, (10)
F n = -(n -1)Z 2 n 2(n -1) p → 0. ( 11
)
Since C n is independent of D n , results ( 7) -( 11) imply that

B n d → N 0, 1 γ 4 + 2 γ 2 . ( 12
)
Results ( 5), (6) (12) yield the theorem.

Discussion

We have seen that McKay's χ 2 approximation for the coefficient of variation is exactly type II noncentral beta distributed. This observation provides insight into the approximation, originally derived by complex analysis. We showed that McKay's χ 2 approximation in distribution equals the product of a χ 2 distributed random variable and a variable that converges in probability to 1. Nevertheless, according to Theorem 5, McKay's χ 2 approximation is asymptotically normal with mean n -1 and variance 2(n -1)(1 + 2γ Though the inverse of the coefficient of variation is noncentral t distributed and algorithms for calculating the cumulative density function of this distribution nowadays exist [START_REF] Lenth | Algorithm AS 243: Cumulative Distribution Function of the noncentral t Distribution[END_REF]), McKay's approximation is still adequate and may be useful for various composite inferential problems on the coefficient of variation in normally distributed data.

Algorithms for computing the cumulative distribution function of the noncentral beta distribution were reviewed by [START_REF] Chattamvelli | A Note on the Noncentral Beta Distribution Function[END_REF]. The open source software R makes use of algorithms given by [START_REF] Lenth | Algorithm AS 226: Computing Noncentral Beta Probabilities[END_REF] and [START_REF] Frick | Algorithm AS R84: A Remark on Algorithm AS 226: Computing Non-central Beta Probabilities[END_REF].

  Vangel, M.G.(1996), Confidence Intervals for a Normal Coefficient of Variation. Amer.Statist.15, 21-26. Warren, W.G. (1982), On the Adequacy of the Chi-Squared Approximation for the Coefficient of[659][660][661][662][663][664][665][666]

c c e p t e d m a n u s c r i p t

	it is "quite adequate for any practical purpose." Also Pearson (1932) examined the new approximation is accurate. Umphrey (1983) corrected a similar study made by Warren (1982) and concluded that McKay's approximation is adequate. Vangel (1996) analytically compared the cumulative density function of McKay's approximation with the cumulative density function of the "naïve" χ 2 approximation N n = (n -1) c 2 γ 2 A compared results according to the approximation with exact results and found that the and showed that McKay's approximation is substantially more accurate. Vangel also pro-
	approximation and found it "very satisfactory." Later Iglewicz & Myers (1970) studied
	the usefulness of McKay's approximation for calculating quantiles of the distribution of
	the sample coefficient of variation c when the underlying distribution is normal. They

formula (1) appears slightly different in the original paper by

[START_REF] Mckay | Distribution of the Coefficient of Variation and the Extended "t" Distribution[END_REF]

since McKay used the maximum likelihood estimator of σ 2 , with denominator n, instead of the unbiased estimator with denominator n -1.

[START_REF] Mckay | Distribution of the Coefficient of Variation and the Extended "t" Distribution[END_REF]

claimed that K n is approximately central χ 2 distributed with n -1 degrees of freedom provided that γ is small (γ < 1/3). This result was established by expressing the probability density function of c as a contour integral and making an approximation in the complex plane. McKay did not theoretically express the size of the error of the approximation. For this reason

[START_REF] Fieller | A Numerical Test of the Adequacy of A. T. McKay's Approximation[END_REF]

, in immediate connection to McKay's paper, investigated McKay's approximation K n numerically and concluded that posed a small modification of McKay's approximation useful for calculating approximate confidence intervals for the coefficient of variation. Forkman (2006) suggested McKay's approximation for testing the hypothesis that two coefficients of variation are equal. Another test for the hypothesis of equal coefficients of variation, also based on McKay's approximation, was proposed by Bennett (1976). It is thus well documented that McKay's approximation is approximately central χ 2 distributed with n -1 degrees of freedom, and useful applications have been suggested. In this paper it is shown that McKay's approximation is type II noncentral beta distributed, and its asymptotic behavior is investigated.

  use the following definition.

A c c e p t e d m a n u s c r i p t

Definition 2. Let U be a central χ 2 distributed random variable with u degrees of freedom, and let V be a noncentral χ 2 distributed random variable, independent of U , with v degrees of freedom and noncentrality parameter λ. The type II noncentral beta distribution with parameters u/2, v/2 and λ, denoted Beta II (u/2, v/2, λ) is defined as the distribution of

Acknowledgements

We thank Prof. Dietrich von Rosen for ideas and discussions. The Centre of Biostochastics, Swedish University of Agricultural Sciences and Pharmacia Diagnostics AB financed the research.

Bennett, B.M. (1976), On an Approximate Test for Homogeneity of Coefficients of Variation, in: W.J. Ziegler, eds. Contributions to Applied Statistics dedicated to A.