Léo Bacon
email: bacon.leo@gmail.com

Alexandre Robert

Yves Hingrat

Long lasting breeding performance differences between wild-born and released females in a reinforced North African Houbara bustard (Chlamydotis undulata undulata) population: a matter of release strategy

Keywords: Captive-breeding, post-release effect, reinforcement, reproduction, translocation

The success of translocation programmes is reflected by the ability of translocated individuals to survive and reproduce in their new environment. However, it has previously been reported that translocated individuals have lower demographic performance than their wild-born conspecifics, due to management and individual factors (such as release conditions or age).

Here, we study six breeding parameters in free-ranging females of the North African Houbara bustard (Chlamydotis undulata undulata) and compare these parameters between captive-bred released (n=204) and wild-born (n=101) birds, considering the age of individuals and the period of release (autumn versus spring). Our results indicate that (1) captive-bred released females successfully breed in the wild; (2) for three out of the six breeding parameters studied, released females show lower performances than wild-born females; but, (3) Although we observed consistently reduced breeding performances in one year old females relative to older females, we did not uncover any interaction between age and the origin of females, suggesting that the impairment of breeding parameters in released females is long lasting; and, (4) interestingly, this impairment of breeding parameters depends on the period of release, with lower breeding performances for spring releases compared to autumn releases. Overall, our study highlights the capacity of captive-bred females to reproduce in the wild, contributing to the dynamics of the population beyond their individual history. Our results also uncover complex variations of breeding parameters in translocated birds, but suggest that these differences can be minimized through an appropriate translocation strategy.

Introduction

The success of conservation translocation programmes is commonly defined as the capacity of translocated populations to persist without further intervention [START_REF] Ewen | Reintroduction Biology: Integrating Science and Management[END_REF], IUCN 2013).

Assessment of success thus requires an understanding of the long-term dynamics and viability of translocated populations (Robert et al. 2015a), which first and foremost requires a long-term demographic assessment [START_REF] Sutherland | Standards for documenting and monitoring bird reintroduction projects[END_REF]. However, such assessments are rarely achieved, mainly due to a lack of resources to support appropriate monitoring and difficulties associated with data analysis [START_REF] Fischer | An assessment of the published results of animal relocations[END_REF][START_REF] Sutherland | Standards for documenting and monitoring bird reintroduction projects[END_REF].

Translocated populations are complex, heterogeneous and unbalanced systems [START_REF] Robert | Integration of demography and genetics in population restorations[END_REF], and the assessment of demographic performances is required to encompass this complexity. In particular, it is necessary to assess potential demographic variations among translocated and wild-born organisms.

Previous knowledge indicates that the demography of translocated populations is likely to be affected by a variety of genetic and non-genetic issues. Genetic issues include outbreeding depression (Huff et al. 2011) and ill adaptation [START_REF] Montalvo | Nonlocal transplantation and outbreeding depression in the subshrub Lotus scoparius (Fabaceae)[END_REF] as well as inbreeding and drift loads [START_REF] Robert | Captive breeding genetics and reintroduction success[END_REF]. In cases of translocations based on captive breeding, adaptation to captivity [START_REF] Frankham | Genetic adaptation to captivity in species conservation programs[END_REF]) might be an important additional issue. In animals, non-genetic issues are associated with captive-rearing and release strategies, leading to various factors potentially affecting demographic rates in the wild. These include individual condition, such as health status, physiology and behaviour [START_REF] Champagnon | Low survival after release into the wild: assessing "the burden of captivity" on Mallard physiology and behaviour[END_REF][START_REF] Dickens | Stress: An inevitable component of animal translocation[END_REF], Hardouin et al. 2014[START_REF] Tavecchia | Maximizing re-introduction success: Assessing the immediate cost of release in a threatened waterfowl[END_REF], proximity of the release site with other populations [START_REF] Mihoub | Post-Release Dispersal in Animal Translocations: Social Attraction and the "Vacuum Effect[END_REF], period of release (Hardouin et al. 2014), age at release [START_REF] Sarrazin | Breeding biology during establishment of a reintroduced Griffon Vulture Gyps fulvus population[END_REF] or the size and composition of the animal group (Hardouin et al. 2014(Hardouin et al. , 2015a)).

Additionally, the potential impairment of demographic performances of translocated individuals might be explained by the interaction of their phenotype with their new environment and release conditions.

Post-release effects have been defined as impairments of vital rates caused by translocation conditions [START_REF] Ewen | Reintroduction Biology: Integrating Science and Management[END_REF]. Such effects can be estimated either by assessing changes in vital rates over time in a release cohort (e.g., [START_REF] Armstrong | Assessing the value of follow-up translocations: a case study using New Zealand robins[END_REF][START_REF] Tavecchia | Maximizing re-introduction success: Assessing the immediate cost of release in a threatened waterfowl[END_REF] or by comparing vital rates of translocated and resident animals over the same time period (e.g., [START_REF] Brown | Wild-reared aplomado falcons survive and recruit at higher rates than hacked falcons in a common environment[END_REF]. A number of studies in translocated vertebrates have documented these post-release effects, mostly in terms of survival probabilities [START_REF] Armstrong | Using Bayesian mark-recapture modelling to quantify the strength and duration of post-release effects in reintroduced populations[END_REF][START_REF] Bertolero | Conservation diagnosis of reintroducing Mediterranean pond turtles: what is wrong?[END_REF], Hardouin et al. 2014[START_REF] Sarrazin | High Survival Estimates of Griffon Vultures (Gyps fulvus fulvus) in a Reintroduced Population[END_REF] but also in reproduction rates [START_REF] Bertolero | Conservation diagnosis of reintroducing Mediterranean pond turtles: what is wrong?[END_REF][START_REF] Converse | Demographics of reintroduced populations: Estimation, modelling, and decision analysis[END_REF][START_REF] Sarrazin | Breeding biology during establishment of a reintroduced Griffon Vulture Gyps fulvus population[END_REF][START_REF] Tavecchia | Maximizing re-introduction success: Assessing the immediate cost of release in a threatened waterfowl[END_REF]) and even in dispersal behaviour [START_REF] Le Gouar | Roles of Survival and Dispersal in Reintroduction Success of Griffon Vulture (Gyps fulvus)[END_REF][START_REF] Mihoub | Post-Release Dispersal in Animal Translocations: Social Attraction and the "Vacuum Effect[END_REF]). In the particular cases of translocated individuals originating from captive-breeding, any impairments of vital rates in translocated individuals can result either from translocation conditions (i.e., post-release effects per se), captivity conditions (including genetic and non-genetic issues), or the interaction of translocation and captivity conditions. However, disentangling the various sources of demographic variation in translocated populations (and in particular those related to the translocation protocol from other sources of variation) is challenging because [START_REF] Maho | Rovers minimize human disturbance in research on wild animals[END_REF] the comparison between translocated and wild-born animals requires the monitoring of wild-born individuals; a difficult task often hindered by several methodological and ethical issues, and ( 2) changes of vital rates over time in translocated individuals can be confounded by other longitudinal sources of demographic heterogeneity, such as age effects (Bacon et al. 2017a), or by inter-individual heterogeneity (e.g., apparent improvement of survival due to the selection of the best survivors).

The North African Houbara bustard (Chlamydotis undulata undulata, hereafter Houbara -Fig. 1) is a medium-sized bird historically distributed from Northern Mauritania to Egypt. It is threatened by over-hunting, poaching and habitat degradation, which have led to the decline of the species during the second half of the twentieth century [START_REF] Goriup | The world status of the houbara bustard Chlamydotis undulata[END_REF]. This decline led to the establishment of a reinforcement programme (the Emirates Center for Wildlife Propagation, ECWP, see Material and Methods) which employs captive breeding and regular releases of captive-bred individuals, combined with ecological research and hunting management, in order to increase the population size of the threatened Houbara throughout its range and to maintain viable populations in adequacy with traditional Arab falconry [START_REF] Lacroix | The Emirates Center for Wildlife Propagation: comprehensive strategy to secure self-sustaining wild populations of houbara bustard (Chlamydotis undulata undulata) in Eastern Morocco[END_REF].

Importantly, captive-bred Houbara are released as juveniles (within their first year of life) either before or during the breeding season (in autumn and in spring, respectively). Previous researches on translocated Houbara indicate that (1) captive-bred juveniles exhibit reduction in their short-term survival (the magnitude of which varies depending on the release season and meteorological conditions), compared to long-term survival (Hardouin et al. 2014), and (2) oneyear-old released females have reduced breeding performance relative to older released females (Bacon et al. 2017a). These results suggest that translocated individuals may experience changes of vital rates over time. However, in these previous assessments of demographic rates, potential changes caused by translocation conditions, captivity conditions and their interactions may have been confounded with age effects. Thus, a rigorous assessment of these potential captivity and translocation effects requires formal comparison between captive-bred released and wild-born individuals of different ages. Hereafter, demographic effects occurring only in the first year following release (i.e., in one-year-old individuals) are referred to as short-term differences and effects occurring independently of the time since release (i.e., in individuals of any age) are referred to as long lasting differences.

Here, we investigate potential differences in breeding performances between wild-born and released captive-bred females, while controlling for age, period of release, and temporal and individual effects. The Houbara is a ground-nesting, gyneparental incubating species. In Eastern Morocco, females breed from mid-February to mid-June. They lay generally two to three eggs, at intervals of 2.5 days, which they incubate for an average of 23 days [START_REF] Gaucher | Breeding biology of the houbara bustard Chlamydotis undulata in Algeria[END_REF]. In case of nest failure, females generally initiate a replacement clutch. Houbara chicks are nidifugous but still rely on their mother for food during the first 10 days after hatching (Saint Jalme & van Heezik 1996) and fledge at a mean age of 60 days (Hardouin et al. 2012). Maximum longevity observed in the wild was at least 13 years old for a wild-born female and 10 years old for captive-bred released females (Bacon et al. 2017a). The diet of the Houbara bustard is generalist and opportunistic and present strong seasonal variation in animal (e.g. Coleoptera, Hymenoptera) and vegetal (e.g. Asteraceae, Brassicaeae, Chenopodiaceae -Bourass et al.

2012) proportion.

We used data collected from a longitudinal nest survey conducted between 2002 and 2016 in the reinforced Houbara population in Eastern Morocco.

We focused on six breeding parameters: nesting effort (number of nesting attempts per breeding season per female), nest initiation date, clutch size, egg volume, daily nest survival and brood survival. We hypothesized that potential short-term reduction of performance in released Houbara are related to their potentially low condition and limited experience with their new environment [START_REF] Bertolero | Assessing the efficacy of reintroduction programmes by modelling adult survival: the example of Hermann's tortoise[END_REF][START_REF] Bertolero | Conservation diagnosis of reintroducing Mediterranean pond turtles: what is wrong?[END_REF]). Thus, potential differences between released and wild-born individuals should decrease with increasing age, reflecting both individual (i.e., gain in experience within an individual, [START_REF] Mauck | Experience versus effort: what explains dynamic heterogeneity with respect to age?[END_REF]) and populationscale processes (i.e., death of lower-quality individuals, assuming that survival and breeding performance may positively covary, Robert et al. 2015b). We thus hypothesize that the difference in breeding performances of wild-born and captive-bred released females will decrease with their age (which is confounded with the time since release for released females).

In other words, we expect to find an interaction between the origin of females (released vs. wild-born) and their age.

Material and methods

Study area

The study was conducted in Eastern Morocco, in the ECWP intervention area (Fig. 2). The study area encompasses approximately 50 000 km² and is characterized by an arid climate marked by irregular rainfalls (less than 200 mm/year). Mean temperatures vary from 6.80°C in winter (December to February) to 26.64°C in summer (June to August). The habitat in the study area is characterized by sparse, shrubby vegetation. Vegetation cover in the plains mainly comprises Chenopodiaceae such as Salsola spp., Hamada spp. and Compositae such as Artemisia herba-alba. In clayand silt-rich areas created by seasonal runoff waters, Salsola sp. is often associated with Atriplex spp. Drainage courses and wadis are characterized by Zizyphus lotus (Ramnaceae) and Retama sp. (Leguminosae). On the high plateaux and slopes, from the extreme east to the Oran region of Algeria, the vegetation is dominated by Stipa tenacissima (Graminea). The core of the study area is used for nomadic pastoral activities, herding and grazing. Mixed herds are generally composed of 100-300 heads of sheep and goats. Nomads settle camps within the study area, whereas some shepherds come from surrounding villages (Le Cuziat et al. 2005).

Emirates Center for Wildlife Propagation and reinforcement programme

The EWCP was established in 1995 to mitigate the severe decline in Houbara populations occurring in the second half of the 20 th century. Two ECWP breeding stations were built for captive breeding (in Missour in 1995 and in Enjil in 2005, see online appendix A). Detailed descriptions of the captive breeding programme can be found in [START_REF] Lesobre | Conservation genetics of Houbara Bustard (Chlamydotis undulata undulata): population structure and its implications for the reinforcement of wild populations[END_REF] and [START_REF] Chargé | Male health status, signalled by courtship display, reveals ejaculate quality and hatching success in a lekking species[END_REF][START_REF] Chargé | Does recognized genetic management in supportive breeding prevent genetic changes in life-history traits?[END_REF]. From 1996 to 2016, 108 486 birds were released in North Africa (94 374 in the study area). Locally, these intensive releases may induce density dependent process leading to deleterious biotic interactions (intraspecific competition, disturbance, predation) and negatively impact life history parameters [START_REF] Azar | Survival of reintroduced Asian houbara in United Arab Emirates' reserves[END_REF], Bacon et al. 2017b).

Within the study area, the overall number of birds released progressively increased over the years (from 28 in 1996 to 9 084 in 2016, see online appendix B), with an increase in the number of release sites and a decrease in the group size per release (Hardouin et al. 2014). The release group size varied between sites (from four to 498 individuals, median = 18 individuals).

Houbara were released at an average of 6 ± 3 (standard deviation) months of age in autumn (August to December), and an average of 9 ± 1 months of age in spring (February to May).

Autumn releases were managed to avoid direct mortality and disturbance due to hunting. In addition, since 2003, some birds were released in summer (June) at an average of 3 ± 0.6 months of age, but they represent only 8% of total releases in the study area. As such, in this study, we only focused on captive-bred females released in spring and autumn.

Hunting was banned between 2000 and spring 2005, and subsequently restricted to the autumn/winter period (October-January) confined within an area comprising 60% of the intervention area. Since 2014, it is estimated that on average 2000 birds are harvested every year, with captive-bred individuals representing 85% of the hunting bags (ECWP unpublished data). The satellite and GSM transmitters used in this study were backpack solar-powered transmitters: PTT-100 30 or 45 g Solar GPS PTT and 30 g solar GSM (Microwave Telemetry, Inc., Columbia, Maryland, USA). To identify breeding events of individuals equipped with satellite transmitters, individuals' movements were remotely monitored. When stationary locations were recorded, the last location was checked in the field to confirm whether the individual was dead or nesting.

Monitoring of captive-bred released females

Trapping and identification of wild-born nesting females

Trapping of wild-born females was performed during egg-collection campaigns in the breeding season of the periods 1997-1998, 2001-2009 and 2015-2016 (see details in online appendix C).

Egg collection campaigns were aimed at building the founder population of the ECWP captive-breeding program [START_REF] Lesobre | Conservation genetics of Houbara Bustard (Chlamydotis undulata undulata): population structure and its implications for the reinforcement of wild populations[END_REF][START_REF] Chargé | Does recognized genetic management in supportive breeding prevent genetic changes in life-history traits?[END_REF]. Females were caught using nylon snares around their nest or around their recently hatched chicks. In total, 143 wild-born females were caught while nesting and were equipped with a transmitter or transmitter replacement. In addition, 10 wild-born females captured and equipped with VHF/PTT transmitters as juveniles were found nesting in later years (see Hardouin et al. 2012 for more details of the capture and monitoring procedures).

Data collection

Breeding parameters

Six breeding parameters were assessed from breeding survey data collected from 2002-2016: the nesting effort (number of nesting attempts per breeding season per female), the clutch size, the egg volume, the nest initiation date, the daily nest survival rate and the daily brood survival rate. The number of years for which survey data was available varied slightly between breeding parameters, this is detailed in online appendix D.

-Nest survey Nests were located from a collaborative survey with shepherds and from individually tracked females. Because of the bias inherent to the trapping and egg-collection campaigns (all wildborn non-equipped nesting females had their eggs collected if not already hatched and were caught) and to avoid confounding factors between transmitter presence and female origin (wildborn and captive-bred), we only selected those nests for which the nesting female was equipped with a transmitter at nest initiation. In total, 703 nests from 305 females were selected for this study. For wild-born females, 300 nests were from 91 wild-born females first captured on their nest, and 21 were from 10 females captured and marked at the juvenile stage. For captive-bred females (n nests = 382), we selected nests from cohorts of females that were released in their first year of life from August to November (in autumn, before the breeding season, n nests = 122, n females = 74) and from February to May (in spring, during the breeding season, n nests = 260, n females = 130). Nests were visited on average every 6 ± 3.60 days. In addition, since 2013, camera traps were deployed at some of the monitored nests (n=45), providing a continuous survey effort and more accurate assessments of nest fate. Nests equipped with camera traps were visited on average every 8 ± 5.07 days. Preliminary unpublished results did not highlight any deleterious effects of the camera trap monitoring on nest survival in our study area (Bacon 2017a). The use of camera traps gave great details on the guild of species predating on Houbara nests. It is composed of the red fox (Vulpes vulpes), the golden jackal (Canis aureus), stray dogs (Canis lupus familiaris), desert hedgehog (Paraechiinus aetiopicus), common raven (Corvus corax) and brown-necked raven (Corvus ruficollis). Nest failure was recorded for a visit interval if the clutch disappeared without signs of hatching [START_REF] Mabee | Using Eggshell Evidence to Determine Nest Fate of Shorebirds[END_REF] or was abandoned. We considered a nest to have survived if at least one egg hatched [START_REF] Mayfield | Suggestions for Calculating Nest Success[END_REF]).

The clutch size was defined as the maximum number of eggs observed in the nest. Since 2003, the length, width (to the nearest hundredth of a centimetre) and mass (to the nearest tenth of a gram) of each egg were measured. The width and length of eggs were then used to calculate the egg volume (Vol) in cm 3 (see details in Bacon et al. 2017a). The nest initiation date was estimated from the calculation of the eggs' incubation stage (in days) at the time of measurement (based on the egg weight loss equation from Hoyt 1979, see details in Bacon et al. 2017a).

The nesting effort was calculated as the number of nests initiated per breeding season per female, including females initiating no nest. We selected females that were continuously monitored during the breeding season. Even if all females were fitted with monitoring devices, technical limitations could cause gaps in their monitoring, which could cause the non-detection of nesting attempts. To increase the detectability of nesting attempts, we only selected females that were fitted with VHF and monitored from February to the end of June, with at least one individual sighted every two weeks (resulting in 173 monitoring histories from 141 females from 2002 to 2014).

Trapping of nesting females causes systematic nest abandonment; thus, female trapping and egg collection have a positive effect on overall nesting effort, as they promote replacement clutches. However, replacement clutches are initiated later in the breeding season, affecting the distribution of nest initiation dates. Therefore, for those two parameters, we excluded females from the years where their eggs were collected (see online appendix D).

-Brood survey Assessing breeding success after hatching and before fledging in species with nidifugous chicks is difficult. It requires the marking and close monitoring of all individuals; a delicate operation that may affect chick survival at this very sensitive stage of life [START_REF] Barron | Meta-analysis of transmitter effects on avian behaviour and ecology[END_REF]). As such, we chose to not mark chicks and instead assessed the fate of broods through the monitoring of brooding females. Therefore, between 2002 and 2014, nesting females fitted with VHF transmitters that successfully hatched at least one egg were then searched for at least once per week. Brood failure was considered for a female at a given monitoring interval if the entire brood disappeared before 60 days of age (Hardouin et al. 2012). In total, we monitored 201 broods from 133 females. A few females (n=12) were captured while rearing chicks, at the nest or after leaving the nest, either to take measurements or to replace the transmitter.

Statistical analyses

Nest and brood survival analyses

We used the logistic exposure approach of [START_REF] Shaffer | A unified approach to analyzing nest success[END_REF] to estimate the daily nest survival rate (defined as the probability of a living nest to not fail in one day). Egg collections and female trapping cause the end of the nest history and produce an unknown "natural" final fate.

Theoretically, in occasional cases where they represent a subset of the sample of nests monitored, nests with an unknown final fate are right truncated at the last date they are known to be active for daily nest survival analysis. In our study, the number of clutches collected was relatively important and biased towards wild-born females. Therefore, nests with an unknown final fate were not representative of the sample of the nests monitored and could lead to an overestimation of the daily nest survival of wild-born females relatively to released females. Thus, we excluded all nests with an unknown final fate from the analysis (following Manolis

Clutch size, egg volume, nest initiation date and nesting effort

The nesting effort was analysed using a generalized linear mixed effects model with a Poisson distribution and a log link function. The clutch size was analysed based on ordinal regression using cumulative linked mixed-effects models (CLMMs). Egg volumes and nest initiation dates were analysed using linear mixed-effects models (LMMs) fitted by maximum likelihood with a Gaussian distribution and an identity link function.

Origin and age factor

We grouped nests into categories based on female origin (wild-born, captive-bred released in autumn and captive-bred released in spring). In addition, we assigned to each nest the age of the nesting female: one-year-old or at least two years old (see online appendix E). For wildborn females, only those captured as juveniles could be aged precisely. A previous study focussing only on released females indicated that one-year-old females had lower breeding performance than older individuals (Bacon et al. 2017a). In addition, the percentage of oneyear-old released females initiating at least one nest is low (8%, ECWP unpublished data).

Based on this information, we assumed that all wild-born females with unknown age were at least two years old at their first identification on a nest.

Adjustment variables

In monitored animals, events such as trapping and marking are stressful and may have some effect on subsequent behavioural and demographic patterns (e.g. unusual movements, increased mortality, cessation of breeding) once released [START_REF] Casas | Assessing the short-term effects of capture, handling and tagging of sandgrouse[END_REF][START_REF] Ponjoan | Adverse Effects of Capture and Handling Little Bustard[END_REF]. For nesting female Houbara, trapping causes nest abandonment (direct effect) but may also have a carry-over effect on subsequent breeding performance. For females trapped on live chicks, trapping may have a direct effect on brood survival by affecting either the female or the chicks. Therefore, we considered an intra-seasonal carry-over effect of trapping for the following parameters: clutch size, egg volume, daily nest survival and daily brood survival rate, when measured on the replacement clutch consecutive to the trapping of the female (carry-over trapping). In addition, specifically for the daily brood survival rate, we considered a direct effect of trapping (direct trapping) on the brood when trapping of the females occurred while rearing chicks. Carry-over trapping and direct trapping effects were implemented as binary adjustment variables. To account for the temporal variation of breeding performance values within the breeding season, we included a date covariate (in Julian days, 1 st of January = 1, implemented as the nest initiation date for egg volume and clutch size, as the median date between two nest visits for the daily nest survival rate, and as the hatching date of the nest for the daily brood survival rate, [START_REF] Grant | Time-specific patterns of nest survival for ducks and passerines breeding in North Dakota[END_REF]. Finally, we considered temporal variation of the daily survival rate of the nest and brood by including, respectively, a binary factor describing whether the nest was in laying or incubation stage and a continuous covariate indicating the age of the brood since hatching [START_REF] Grant | Time-specific patterns of nest survival for ducks and passerines breeding in North Dakota[END_REF].

Model selection

For each breeding parameter analysis, covariates were rescaled, centralized and standardized by two times the standard deviation [START_REF] Gelman | Scaling regression inputs by dividing by two standard deviations[END_REF], and factors with two levels were transformed into binary dummy variables (0-1).

We examined the relationship between the breeding parameters and their corresponding variables by developing sets of candidate models that included additive terms of explanatory variables, quadratic terms of continuous temporal covariates (date² and brood age²) and the interaction between the age of the nesting female and its origin. To adjust for potential effects of trapping, we kept the trapping carry-over effect and the direct trapping effect in all candidate models of the concerned breeding parameters. Additionally, current year, female identity and nest identity nested under female identity (only for egg volume analyses) were incorporated as random intercept effects in all models. Random intercept effects were included to avoid pseudoreplication and to account for potential heterogeneity among years, females and nests [START_REF] Zuur | Mixed effects models and extensions in ecology with R, Statistics for Biology and Health[END_REF]).

Thus, the starting (most complex) statistical models were of the following form:

𝑌 = 𝐹𝑒𝑚𝑎𝑙𝑒 𝐴𝑔𝑒 + 𝐹𝑒𝑚𝑎𝑙𝑒 𝑂𝑟𝑖𝑔𝑖𝑛 + 𝐹𝑒𝑚𝑎𝑙𝑒 𝐴𝑔𝑒 × 𝐹𝑒𝑚𝑎𝑙𝑒 𝑂𝑟𝑖𝑔𝑖𝑛 + 𝐴𝑑𝑗𝑢𝑠𝑡𝑚𝑒𝑛𝑡 𝑣𝑎𝑟𝑖𝑎𝑏𝑙𝑒𝑠 + 𝑟𝑎𝑛𝑑𝑜𝑚 𝑖𝑛𝑡𝑒𝑟𝑐𝑒𝑝𝑡𝑠
where Y was the dependent variable (the breeding parameter). Details of specific adjustment variables (mentioned in the section Adjustment variables) implemented in the starting model of each breeding parameter can be found in Appendix F.

Starting from this generic structure, models were simplified to develop all possible combinations of explanatory variables (including additive, polynomial terms and the interaction, see all model combinations in Appendix G). This procedure was applied to our six breeding dependent variables. For each breeding dependent variable, we ranked models using the Akaike Information Criterion corrected for small sample size (AICc). The best model (ΔAICc = 0; Burnham & Anderson 2002) was selected to estimate β coefficients and associated 95% confidence intervals for fixed-effect covariates. All analyses were conducted using R 3.3.3 (R Development Core Team, 2017) and the packages lme4 1.1-12 [START_REF] Bates | Fitting Linear Mixed-Effects Models Using lme4[END_REF], Ordinal 2015.6-28 [START_REF] Christensen | ordinal -Regression Models for Ordinal Data[END_REF], MuMIn 1.13.4 (Barton 2015), and nest survival (Herzog 2009).

Normality and homoscedasticity of the random residuals were checked graphically from the starting models, and no overdispersion [START_REF] Zuur | Mixed effects models and extensions in ecology with R, Statistics for Biology and Health[END_REF]) was observed for the Poisson distribution model.

Results

On average, females were monitored over 2.97 ± 1.50 years (3.58 ± 1.76 years for wild-born females and 2.70 ± 1.25 years for released females) and 86% of females were monitored over more than 1 year (88% for wild-born females and 85% for released females).

We analysed the effects of dependent covariates on our six breeding parameters based on a common statistical approach. To keep the presentation of the results clear and straightforward, we directly present the coefficients from the best selected models (ΔAICc = 0) in Table 1. AICc based model ranking selections, which compiled six tables, can be found in online appendix G.

According to the best models (Table 1), clutch size and egg volume declined along the breeding season. Daily nest survival was higher during the nest initiation stage than during the incubation stage, while daily brood survival increased with the age of the brood. For all breeding parameters concerned (clutch size, egg volume, daily nest survival and daily brood survival), the carry-over trapping effect had 95% confidence intervals overlapping 0, whereas direct trapping of females while rearing chicks had a very strong negative impact on the daily survival of the brood (Table 1). For all breeding parameters (with the exception of the daily nest survival rate), the Age effect was selected in all best models, indicating that one-year-old females had poorer breeding performances than older ones (Table 1).

The origin factor was selected in the best model sets for nesting effort, egg volume and daily brood survival rate. In Table 1, for the female origin factor, the results show the differences with respect to a reference level at the intercept. We defined wild-born females as the reference level at the intercept in Table 1, but all levels of the female origin factor were set as the reference level at the intercept in order to assess potential significant differences. Captive-bred females initiated fewer nests and produced smaller eggs than wild-born females (Figs 3A,B). However, confidence intervals of these effects did not overlap 0 only for captive-bred females released in spring relative to wild-born females (Table 1). Captive-bred females released in autumn had intermediate values of nesting effort and egg volume, but with confidence intervals difference to the two other groups overlapping 0 (Table 1). Finally, captive-bred females released in spring had lower daily brood survival than wild-born females and captive-bred females released in autumn (Fig. 3C, Table 1), and the confidence intervals of the difference between captive-bred females released in spring and wild-born females was close to negative (Table 1).

We were particularly interested in the interaction term of the age of nesting females and their origin, which quantifies whether any effect of origin depends on age (i.e., short-term) or not (i.e., long lasting). However, none of the best models (Table 1) included this interaction, providing no support for a more marked effect of the origin in one-year-old females as compared to older females.

Discussion

The assessment of the breeding performances of released individuals is of major importance in translocation programmes [START_REF] Bertolero | Conservation diagnosis of reintroducing Mediterranean pond turtles: what is wrong?[END_REF]. In this study, we investigated the differences of breeding parameters between wild and captive-bred released Houbara females and revealed how release strategies may influence the magnitude of this differences. This is a prerequisite to adapting and improving translocation strategies and thus to increase the chance of translocation success [START_REF] Armstrong | Directions in reintroduction biology[END_REF].

Most breeding parameters (except nest survival) consistently showed age-specific variation, suggesting that one-year-old females have reduced performances compared to older ones, in agreement with a previous study (Bacon et al. 2017a). Individual heterogeneity was taken into account in our statistical framework (see [START_REF] Nussey | The rate of senescence in maternal performance increases with early-life fecundity in red deer[END_REF][START_REF] Van De Pol | Age-Dependent Traits: A New Statistical Model to Separate Within-and Between-Individual Effects[END_REF]. This suggests that the observed effect of age is best explained by longitudinal improvement of individuals (increased breeding experience or an increase in breeding investment, [START_REF] Curio | Why young birds reproduce less well[END_REF], than by population level processes (e.g., selection of the best individuals over time). The lack of age effects on the daily nest survival may be explained by the weight of external factors, such as environmental conditions [START_REF] Morales | Annual productivity and individual female reproductive success in a Great Bustard Otis tarda population[END_REF] or predation, which is a prominent cause of nest failure in birds [START_REF] Martin | Avian Life History Evolution in Relation to Nest Sites, Nest Predation, and Food[END_REF], including Houbara (Bacon et al. 2017a[START_REF] Koshkin | Effects of habitat and livestock on nest productivity of the Asian houbara Chlamydotis macqueenii in Bukhara Province, Uzbekistan[END_REF], potentially overriding any age-dependent effects.

Regardless of age effects, our results indicate that, for some of the breeding parameters studied, wild-born individuals performed slightly better than individuals translocated in autumn and better than individuals translocated in spring. More precisely, wild-born females initiated more nests with larger eggs when compared to captive-bred females released in spring. In addition, the daily brood survival of captive-bred females released in spring was lower than for wildborn females. These differences have a marked size effect. It corresponds to a decline of ~47% in nesting effort and a decline of ~40% in brood survival (from hatching to fledging) when considering adult females released in spring relative to wild-born adult females (Figs 3A,C).

Conversely to our expectations, we did not observe any interaction between the age of the female and its origin. This result is consistent with the long lasting (i.e., permanent) negative post release effect on breeding parameters observed in other translocated species (see, e.g., [START_REF] Sarrazin | Breeding biology during establishment of a reintroduced Griffon Vulture Gyps fulvus population[END_REF]. It is also consistent with previous findings in Houbara that released captive-bred females exhibit age-dependent patterns of variation in breeding performance, similar to those observed in many vertebrate populations in the wild (i.e., a bell-shaped trend divided into youth, adulthood and senescence; see discussion in Bacon et al. 2017a).

As differences between wild-born and captive-bred individuals were long lasting, we cannot exclude the possibility that genetic causes related to captive-breeding issues partly explain the observed differences. However, strong and restrictive guidelines are conducted in the ECWP captive breeding programme [START_REF] Chargé | Does recognized genetic management in supportive breeding prevent genetic changes in life-history traits?[END_REF]. These guidelines aim to minimize the loss of genetic variation and to minimize the adaptation to captivity in the captive population in order to maximize the fitness of individuals released into the wild, while reducing the genetic consequences of the reinforcement programme on wild populations [START_REF] Chargé | Does recognized genetic management in supportive breeding prevent genetic changes in life-history traits?[END_REF]).

Nevertheless, the relaxation of sexual and natural selection and adaptation to captivity may partly explain the lower breeding performances of captive-bred individuals in the wild [START_REF] Araki | Genetic Effects of Captive Breeding Cause a Rapid, Cumulative Fitness Decline in the Wild[END_REF][START_REF] Frankham | Genetic adaptation to captivity in species conservation programs[END_REF]. Beside, owing to the substantial translocation effort (see online appendix B), a substantial proportion of wild-born individuals may be offspring from translocated birds (see general discussion in [START_REF] Barbanera | Genetic consequences of intensive management in game birds[END_REF], which may explain why the breeding performances of released and wild-born females are relatively similar. Therefore, it is likely that most differences between the two groups are related to non-genetic issues.

At the phenotypic level, previous research into captive breeding/rearing in birds has suggested that demographic differences between released and wild-born individuals can be explained by the lack of exposure of chicks hand-reared in captivity to adult behaviour [START_REF] Roche | Relative fitness of wild and captive-reared piping plovers: Does egg salvage contribute to recovery of the endangered Great Lakes population?[END_REF], permanently decreasing their fitness in the wild [START_REF] Buner | How do different releasing techniques affect the survival of reintroduced grey partridges Perdix perdix?[END_REF][START_REF] Buner | Experimental assessment of release methods for the re-establishment of a red-listed galliform, the grey partridge (Perdix perdix)[END_REF][START_REF] Pérez | Are parent-reared red-legged partridges (Alectoris rufa) better candidates for reestablishment purposes?[END_REF], Sánchez-Garcia et al. 2011). Although wild-born and released captive-bred Houbara females have a similar diet [START_REF] Bourass | Diet of released captive-bred North-African houbara bustards[END_REF] and exploit equivalent ecological niches [START_REF] Monnet | The realized niche of captive-hatched Houbara Bustards translocated in Morocco meets expectations from the wild[END_REF], released females may be less efficient at finding and exploiting resources [START_REF] Champagnon | Low survival after release into the wild: assessing "the burden of captivity" on Mallard physiology and behaviour[END_REF]. This could explain the lower investment of captive-bred females in reproduction (e.g., egg production), which may subsequently affect the survival of the chicks.

In species with precocial chicks (including nidifugous chicks), smaller eggs generally contain less yolk, giving rise to chicks that are in lower condition with lower nutrient reserves at hatching and that are more prone to early mortality [START_REF] Amundsen | Egg size and parental quality influence nestling growth in the shag[END_REF][START_REF] Blomqvist | Parental quality and egg size affect chick survival in a precocial bird, the lapwing Vanellus vanellus[END_REF][START_REF] Galbraith | Effects of egg size and composition on the size, quality and survival of lapwing Vanellus vanellus chicks[END_REF][START_REF] Williams | Intraspecific variation in egg size and egg composition in birds: effects on offspring fitness[END_REF]). The differences observed between wild-born and translocated females must be interpreted while taking into account that all individuals were fitted with transmitters. Markings devices can have significant impacts on some demographic, physiological and behavioural parameters [START_REF] Barron | Meta-analysis of transmitter effects on avian behaviour and ecology[END_REF], see discussion in Le [START_REF] Maho | Rovers minimize human disturbance in research on wild animals[END_REF]. However, because wild-born and released females are unlikely to exhibit differential responses to transmitter fitting, our main conclusions are unlikely to be affected by transmitter effects.

Our results indicate that the difference observed between wild-born and released females varies according to the period of release (autumn versus spring). Captive-bred females released in autumn had consistently higher nesting effort, egg volume and daily brood survival values than females released in spring (Table 1). Assuming that survival and reproductive parameters positively covary (e.g., "fixed heterogeneity", [START_REF] Cam | The Conundrum of Heterogeneities in Life History Studies[END_REF], such a pattern may be related to the death of the less performant females between their release and the breeding season among females released in autumn. This hypothesis is supported by previous findings showing that harsh environmental conditions in autumn and winter lower short-term survival rates compared to spring releases (Hardouin et al. 2014), representing a stronger environmental filter that primarily affects weak individuals.

It is important to note that not all breeding parameters differed between wild-born and captivebred females. The origin of females was not selected in the best models for clutch size, nest initiation date and daily nest survival analysis. Similar heterogeneous responses of multiple breeding parameters to female origin were observed by Hill & Roberston (1988) for a reinforced ring-necked pheasant (Phasianus colchicus) population, where wild-born and captive-bred released females had similar nesting effort and nesting success, yet chick survival to fledging and total productivity were higher for wild-born females.

Management implications

We observed consistent, permanent impairments in some breeding parameters of captive-bred released females relative to wild-born females. The magnitude of these differences highlights the importance of assessing demographic impacts of captivity and post release effects in translocated populations. In our study, the observed differences between wild-born and released females appears to be dependent on the translocation strategy (period of release), which implies that it could be buffered by adapting translocation strategies (favouring autumn releases).

Nonetheless, the current strategy (spring and autumn releases) seems suitable regarding the multiple objectives of the reinforcement programme, i.e., reinforcement of the wild population and supplementation for hunting (falconry). For autumn releases, environmental stochasticity may select individuals in good condition (Hardouin et al. 2014) that will later breed with relatively better parameters (as suggested by the present results). On the other hand, although outcomes from spring releases will be less certain in terms of breeding performances, higher survival rates in this season will help to maintain the population size at a level that can sustain hunting pressure (Hardouin et al. 2015a) and population growth rate. Finally, although differences were observed between wild-born and captive-bred females, our results confirm the capacity of captive-bred females to reproduce in the wild and to contribute to the dynamics of the population beyond their individual history. This partially supports the success of the reinforcement programme [START_REF] Bertolero | Conservation diagnosis of reintroducing Mediterranean pond turtles: what is wrong?[END_REF]) and provides a strong basis for more integrative and quantitative assessment of translocation success (Robert et al. 2015a).

  Prior to 2006, released Houbara were tagged on the left tarsus with an aluminium ring displaying a unique ID number. Beginning in 2006, all released individuals were subcutaneously tagged with a unique radio-frequency identification (RFID) microchip (TROVAN LID100 implantable transponder, DorsetID; see Hardouin et al. 2015b for details). In Eastern Morocco, 3 014 released birds (49% females) were equipped with transmitters (1 937 with VHF transmitters, 933 with satellite transmitters and 144 with GSM transmitters). The VHF transmitters used in this study were battery-powered necklace units with mortality signals (11 g model RI-2B-M, 20 g model RI-2D-M, Holohil System Ltd., Carp, Ontario, Canada) and solar-powered backpack (19 g model, Merlin Systems Inc., Boise, Idaho, USA). During the breeding season, females fitted with VHF transmitters were searched for at least once per week.

  et al. 2000 and personal communication from Stephen Dinsmore). Finally, as records of periodic nest visits were only standardized beginning in 2003, we discarded data collected prior to 2003.The same statistical approach was used to study brood survival, with the exception that the exposure period starts at hatching and ends 60 days later for successful broods (mean age at fledging, Hardouin et al. 2012, but see details inBacon et al. 2017a).

Acknowledgements

This study was funded by the Emirates Center for Wildlife Propagation (ECWP), a project of the International Fund for Houbara Conservation (IFHC). We are grateful to H.H. Sheikh Mohammed bin Zayed Al Nahyan, Crown Prince of Abu Dhabi and Chairman of the IFHC and H.E. Mohammed Al Bowardi, Deputy Chairman of IFHC, for their support. This study was conducted under the guidance of Reneco International Wildlife Consultants LLC., a consulting company managing ECWP. We are thankful to Dr. F. Lacroix, managing director, and G.

Leveque, project director, for their supervision. We sincerely thank all ECWP staff from the Ecology division who participated in data collection. We are thankful to Pascale Reding and Grégoire Liénart for producing our study area map. We are grateful to Doug Armstrong and two anonymous reviewers for their helpful and constructive comments on the manuscript, as well as to Dr Thomas Martin for improving the English text.

Hardouin, L.A., Nevoux, M., Robert, A., Gimenez, O., Lacroix, F., Hingrat, Y., 2012. Determinants and costs of natal dispersal in a lekking species. Oikos 121, 804-812. doi:10.1111/j.1600-0706.2012.20313.x Hardouin, L.A., Robert, A., Nevoux, M., Gimenez, O., Lacroix, F., Hingrat, Y., 2014. Meteorological conditions influence short-term survival and dispersal in a reinforced bird population. Journal Applied Ecology 51, 1494Ecology 51, -1503Ecology 51, . doi:10.1111Ecology 51, /1365Ecology 51, -2664.12302 .12302 Hardouin, L.A., Hingrat, Y., Nevoux, M., Lacroix, F., Robert, A., 2015a


Peter W Glynn

Marvin K Nakayama

Bruno Tuffin

Comparing Regenerative-Simulation-Based Estimators of the Distribution of the Hitting Time to a Rarely Visited Set

We consider the estimation of the distribution of the hitting time to a rarely visited set of states for a regenerative process. In a previous paper, we provided two estimators that exploited the weak convergence of the hitting time divided by its expectation to an exponential as the rare set becomes rarer. We now add three new estimators, based on a corrected exponential, a gamma, and a bootstrap approach, the last possibly providing less biased estimators when the rare set is only moderately rare. Numerical results illustrate that all of the estimators perform similarly. Although the paper focuses on estimating a distribution, the ideas may also be applied to estimate risk measures, such as a quantile or conditional tail expectation.

INTRODUCTION

The hitting time of a rarely visited set A of states is of interest in many areas, such as dependability and finance. While rare-event simulation is often used to study the hitting time (Rubino and Tuffin 2009b), most of the literature focuses on estimating its mean. But its quantiles or entire distribution can often provide more useful information. For example, if the hitting time represents the failure time of a product that a company plans to sell, a quantile of its distribution may be employed to define an appropriate warranty length. The time to ruin of an insurance company may be modeled as a hitting time, and we may be interested in its distribution; e.g., see Section 1.3.5 of [START_REF] Kalashnikov | Geometric Sums: Bounds for Rare Events with Applications[END_REF]. But the study of rare-event-simulation techniques for computing the distribution of hitting times has not received a lot of attention.

We now address this problem for a stochastic process having a regenerative structure; i.e., the process "probabilistically restarts" at an increasing sequence of regeneration times [START_REF] Kalashnikov | Topics on Regenerative Processes[END_REF]. We then can express the mean hitting time as a ratio, which may be utilized to design efficient simulation methods; see [START_REF] Goyal | A Unified Framework for Simulating Markovian Models of Highly Dependable Systems[END_REF], [START_REF] Nicola | Fast simulation of dependability models with general failure and repair processes[END_REF][START_REF] Glynn | On the Estimation of the Mean Time to Failure by Simulation[END_REF], and Rubino and Tuffin (2009b).

It is often the case [START_REF] Kalashnikov | Geometric Sums: Bounds for Rare Events with Applications[END_REF]) that when the probability p of hitting A before regeneration goes to zero, the distribution of the hitting time divided by its mean converges to an exponential with mean 1. We can exploit this property to approximate the distribution of the hitting time by an exponential, with the problem now being reduced to estimating the mean hitting time, which may be efficiently estimated, as noted above. [START_REF] Glynn | Using Simulation to Calibrate Exponential Approximations to Tail-Distribution Measures of Hitting Times to Rarely Visited Sets[END_REF] develop this idea to provide estimators of the distribution, quantile and conditional tail expectation (CTE), each requiring only an estimator of the mean hitting time. [START_REF] Glynn | Using Simulation to Calibrate Exponential Approximations to Tail-Distribution Measures of Hitting Times to Rarely Visited Sets[END_REF] further propose a second family of methods, resulting in so-called convolution estimators. This approach uses the fact that the scaled sum of lengths of the cycles (i.e., the process between regenerations) before the one hitting A converges weakly (as p shrinks) to an exponential, whose mean can be estimated. We then convolve the exponential with the distribution of the time to reach A within a cycle, which is estimated by importance sampling (IS; Asmussen and Glynn 2007, Chapter VI).

While the above estimators are both computationally efficient and can have small variance, a drawback is that they are also biased. The weak convergence on which they are based requires p → 0. But when this probability is not very small, bias could potentially surpass variance. For this reason, we propose in this paper three other estimators, which may be more accurate: i) a convolution estimator for which the exponential approximation for the sum of lengths of cycles before the one hitting A is replaced by a refined approximation of [START_REF] Blanchet | Uniform Renewal Theory with Applications to Expansions of Random Geometric Sums[END_REF]; ii) an estimator that substitutes the exponential approximation with one based on a gamma, giving one extra degree of freedom in the inference of the distribution, which may be helpful for moderate p; iii) a bootstrap estimator that does not use any distributional approximations but rather resamples simulated lengths of independent cycles. We focus on estimating the distribution of the hitting time, but quantile and CTE estimators can be derived as well. Our numerical experiments show that all of the distribution estimators are close in terms of performance.

The rest of the paper proceeds as follows. Section 2 sets the mathematical framework and reviews the exponential and convolution estimators of [START_REF] Glynn | Using Simulation to Calibrate Exponential Approximations to Tail-Distribution Measures of Hitting Times to Rarely Visited Sets[END_REF]. The next three sections develop the new estimators, with the hope that they can reduce the error when the probability p is not too small. Section 3 describes the modified convolution estimator based on a refinement of the exponential approximation by [START_REF] Blanchet | Uniform Renewal Theory with Applications to Expansions of Random Geometric Sums[END_REF]. Section 4 explains the gamma estimator, and we present the bootstrap estimator in Section 5. Section 6 compares all estimators through numerical experiments on two standard regenerative processes (a model of highly reliable systems used in dependability analysis, and an M/M/1 queue) for which the exact values are known, so we can analyze the error. Section 7 gives concluding remarks and directions for future work.

MODEL AND PREVIOUS RESULTS

Consider a stochastic process X = [X(t) : t ≥ 0] evolving on a state space S ⊆ ℜ d . The goal is to estimate the distribution/quantiles of the hitting time T = inf{t ≥ 0 : X(t) ∈ A } of a subset A ⊂ S of states.

Assume process X is regenerative, with 0 = Γ 0 < Γ 1 < Γ 2 < • • • as the sequence of regeneration times, so X "probabilistically restarts" at each Γ i . For i ≥ 1, define τ i = Γ i -Γ i-1 as the time between regenerations, and let τ be a generic copy of τ i . The process (X(Γ i-1 + s) : 0 ≤ s < τ i ) between regenerations i -1 and i is called the ith (regenerative) cycle of X. The pairs (τ i , (X(Γ i-1 + s) : 0 ≤ s < τ i )), i ≥ 1, are independent and identically distributed (i.i.d.) (Kalashnikov 1994, Section 1.3). Define also T i = inf{s > 0 : X(Γ i-1 +s) ∈ A } as the first time to reach A after the (i -1)th regeneration, and let M = sup{i > 0 : T i > τ i } be the number of cycles completed before the first in which A is hit. We then can express

T = S +V ≡ M ∑ i=1 τ i + T M+1 , (1)
where the geometric sum S is independent of V . Let F be the cumulative distribution function (cdf) of T , denoted as T ∼ F. Also, let G and H be the cdfs of S and V , respectively. Independent of M, each τ i in S has a cdf K that is the conditional cdf of τ given τ < T ; i.e., K(x) = P(τ ≤ x | τ < T ).

We also assume that the probability p = P(T < τ) to reach A before a regeneration is small. To study a limiting behavior, we parameterize (and index) the model by ε, and our goal is to study the cdf [START_REF] Glynn | Using Simulation to Calibrate Exponential Approximations to Tail-Distribution Measures of Hitting Times to Rarely Visited Sets[END_REF] note that in different contexts [START_REF] Kalashnikov | Geometric Sums: Bounds for Rare Events with Applications[END_REF],

F ≡ F ε of T ≡ T ε when p ≡ p ε → 0 as ε → 0. With µ ε = E ε [T ε ],
lim ε→0 P ε (T ε /µ ε ≤ t) = 1 -e -t , ∀t ≥ 0, (2)
so T ε /µ ε converges weakly to an exponential as ε → 0. An interesting consequence of (2) is that the computation of the cdf F ε of hitting times reduces asymptotically (as ε → 0) to computing its mean. We next describe two standard examples of asymptotic regimes in which p ε → 0 as ε → 0.

• For the process with X(t) representing the number of customers at time t in a GI/GI/1 queue with first-in-first-out discipline, we may be interested in computing the distribution of the hitting time

T = T ε of the set A ≡ A ε = { 1/ε , 1/ε + 1, .
. . } of states. Regenerations occur at the beginnings of each busy period [START_REF] Kalashnikov | Topics on Regenerative Processes[END_REF], Example 1.2.2), and it is clear that p ε → 0 as ε → 0. • For a highly reliable Markovian system (HRMS), studied among others in [START_REF] Shahabuddin | Importance Sampling for Highly Reliable Markovian Systems[END_REF], [START_REF] Nakayama | General Conditions for Bounded Relative Error in Simulations of Highly Reliable Markovian Systems[END_REF], L'Ecuyer and Tuffin (2012), Rubino and Tuffin (2009a), let c be the number of types of components, where each type l = 1, 2, . . . , c, has a given redundancy n l . Each component of each type is subjected to failures and repairs, with failure and repair times exponentially distributed.

We can then define X as a continuous-time Markov chain on a state space S , where a state in S includes information on the number of operational components of each type (along with any necessary queueing information about failed components waiting for repair). The set S of states is decomposed into a set of states in which the system is defined as working and the set A of failed states. Regenerations occur on returns to the fully operational state, and T represents the hitting time to A . Failure rates are small, of order O(ε), with respect to repair rates, of order Θ(1), with a repair possible from any non-fully operational state. (Recall that a function f

(ε) is O(g(ε)) if | f (ε)/g(ε)| remains bounded when ε → 0, and it is Θ(g(ε)) if | f (ε)/g(ε)
| is bounded and also bounded away from 0, when ε → 0.) Again, we have p = p ε → 0 as ε → 0 [START_REF] Shahabuddin | Importance Sampling for Highly Reliable Markovian Systems[END_REF].

We adopt the following notational convention in the rest of the paper. For an unknown quantity such as F, we use a tilde, as in F, to denote a (non-simulation) approximation to F, typically derived from a weak-convergence result, such as (2). A simulation-based estimator of F has a hat, as in F.

Exponential Approximation Estimator

By (2), the computation of the hitting time's cdf reduces asymptotically (as ε → 0, where we often omit the subscript ε to simplify notation) to estimating its mean µ = E[T ]: ∀t ≥ 0,

F(t) = P(T ≤ t) ≈ F exp (t) ≡ 1 -e -t/µ . (3)
For regenerative systems, we can express µ as a ratio [START_REF] Goyal | A Unified Framework for Simulating Markovian Models of Highly Dependable Systems[END_REF])

µ = E[T ∧ τ] E[I (T < τ)] ≡ ζ p (4)
with x ∧ y = min(x, y) and I (•) the indicator function. This allows estimating µ by measure-specific importance sampling (MSIS; [START_REF] Goyal | A Unified Framework for Simulating Markovian Models of Highly Dependable Systems[END_REF], which independently estimates ζ and p in (4). To do this, we specify a proportion 0 < γ < 1 used to allocate a total of n independent simulated cycles as follows.

• Generate n CS ≡ γn independent cycles by crude simulation (CS; i.e., no IS), giving i.i.d. observations 4) is a rare event [START_REF] Glynn | Using Simulation to Calibrate Exponential Approximations to Tail-Distribution Measures of Hitting Times to Rarely Visited Sets[END_REF], we employ IS to estimate p = E[I (T < τ)] = I (T < τ) dP. Specifically, rather than sampling using the original probabilistic dynamics of P, IS instead simulates under another probability measure P , and we can apply a change of measure to express p = I (T < τ)L dP = E [I (T < τ)L], with L = dP/dP as the likelihood ratio and E the expectation operator under P . This motivates estimating p by

T i ∧ τ i , i = 1, 2, . . . , n CS , used to estimate the numerator ζ in (4) by ζ n = (1/n CS ) ∑ n CS i=1 T i ∧ τ i . • Because {T < τ} in the denominator of (
p n = 1 n IS n IS ∑ i=1 I (T i < τ i )L i , (5)
where

(I (T i < τ i ), L i ), i = 1, 2, . . . , n IS ≡ (1 -γ)n, are i.i.d. copies of (I (T < τ), L) under P .
Applying MSIS leads to the following estimator [START_REF] Glynn | Using Simulation to Calibrate Exponential Approximations to Tail-Distribution Measures of Hitting Times to Rarely Visited Sets[END_REF].

Estimator 1 The exponential estimator of the cdf F(t) of T is

F exp,n (t) = 1 -e -t/ µ n , (6)
where

µ n = ζ n / p n is the MSIS estimator of µ.
As in [START_REF] Goyal | A Unified Framework for Simulating Markovian Models of Highly Dependable Systems[END_REF], we may select the proportion γ to minimize the variance per unit of computational budget (or work-normalized variance) of the estimator µ n . Also, [START_REF] Nakayama | Efficient Estimation of the Mean Hitting Time to a Set of a Regenerative System[END_REF] develop another regenerative estimator of µ based on the expression µ

= E[T ;T <τ]+E[τ;T >τ] E[τ]
, which can alternatively be used and could yield more accurate results, where E[X; A] = E[XI (A)] for an event A.

From the cdf estimator (6), [START_REF] Glynn | Using Simulation to Calibrate Exponential Approximations to Tail-Distribution Measures of Hitting Times to Rarely Visited Sets[END_REF] further estimate the q-quantile ξ = F -1 (q), 0 < q < 1, and the CTE γ

= E[T |T > ξ ] by ξ exp,n = F -1 exp,n (q) = -ln(1 -q) µ n and χ exp,n = (1 -ln(1 -q)) µ n , (7)
respectively, where the CTE estimator uses the fact that when T is exactly exponential, the memoryless property implies its CTE is µ + ξ . Moreover, given a confidence interval (CI) for µ based on µ n , we can easily obtain CIs for ξ and γ through (7), which are linear transformations of µ n [START_REF] Glynn | Using Simulation to Calibrate Exponential Approximations to Tail-Distribution Measures of Hitting Times to Rarely Visited Sets[END_REF]). In the following, we focus on estimating only the cdf F of T and for space reasons will not further consider estimators for quantiles and CTE; see [START_REF] Glynn | Using Simulation to Calibrate Exponential Approximations to Tail-Distribution Measures of Hitting Times to Rarely Visited Sets[END_REF] for numerical results on the latter two.

Convolution Estimator

Because T = S +V in (1) has the geometric sum S ∼ G independent of V ∼ H, the cdf F of T satisfies

F(t) = G H(t) = H(t -s) dG(s) (8)
where denotes the convolution operator. This suggests convolving estimators of G and H to obtain an estimator of F, which Glynn et al. ( 2018) develop as follows.

• For each t ≥ 0, we often have [START_REF] Kalashnikov | Geometric Sums: Bounds for Rare Events with Applications[END_REF]) that

G(t) ≈ G exp (t) = 1 -e -t/η for p ≈ 0, where η = E[S] = E[M] • E[τ | τ < T ].
As M is geometric with probability mass function (pmf) 4). We apply a change of measure to get

P(M = k) = (1 -p) k p for k ≥ 0, its mean is E[M] = (1 -p)/p, which we estimate by (1 -p n )/ p n , where p n is from (5). As E[τ | τ < T ] = E[τI (τ < T )]/(1 -p) is unknown, we estimate it via (1/((1 -p n )n CS )) ∑ n CS i=1 τ i I (τ i < T i ) using the n CS i.i.d. observations of (τ i , I (τ i < T i )) of (τ, I (τ < T )) from crude simulation. We then estimate η by η n = 1 p n n CS n CS ∑ i=1 τ i I (τ i < T i ), resulting in G exp,n (t) = 1 -e -t/ η n as an exponential estimator of G(t). • The cdf H of V satisfies H(x) = P(V ≤ x) = P(T ≤ x | T < τ) = P(T ≤ x, T < τ)/p by (
P(T ≤ x, T < τ) = E[I (T ∧ τ ≤ x, T < τ)] = E [I (T ∧ τ ≤ x, T < τ)L]. Thus, from i.i.d. observations (T i ∧ τ i , I (T i < τ i ), L i ), i = 1, 2, . . . , n IS , of (T ∧ τ, I (T < τ), L)
under IS measure P , we estimate H(x) by (using ( 5))

H n (x) = 1 p n n IS n IS ∑ i=1 I (T i ∧ τ i ≤ x, T i < τ i )L i = ∑ n IS i=1 I (T i ∧ τ i ≤ x, T i < τ i )L i ∑ n IS j=1 I (T j < τ j )L j . ( 9
)
Taking the convolution of G exp,n and H n , we get [START_REF] Glynn | Using Simulation to Calibrate Exponential Approximations to Tail-Distribution Measures of Hitting Times to Rarely Visited Sets[END_REF] Estimator 2 The convolution estimator of cdf F(t) of T is

F conv,n (t) = 1 - 1 p n • n IS n IS ∑ i=1 I (T i < τ i ) L i e -(t-(T i ∧τ i )) + / η n . ( 10
)
For a fair comparison with the exponential estimator F exp,n in (6), we defined the convolution estimator F conv,n in (10) using the same sample sizes n CS = γn and n IS = (1γ)n for CS and IS, respectively. But we could also choose the MSIS parameter γ differently for the convolution estimator.

CORRECTED CONVOLUTION ESTIMATOR

While Section 2.2 obtained F conv,n in (10) by convolving G exp,n and H n , we can construct additional estimators of F by instead taking the convolution of other estimators of G and H. In particular, we now propose to do this by exploiting an approximation, developed by [START_REF] Blanchet | Uniform Renewal Theory with Applications to Expansions of Random Geometric Sums[END_REF], to the cdf G of a (slightly different) geometric sum S . Specifically, let S = ∑ M i=1 Y i , where (Y i ) i≥1 are i.i.d. copies of a random variable Y and M = M + 1 for M as in ( 1), where M is independent of the Y i . Thus, M has the alternative definition of a geometric with pmf P(M = k) = (1p) k-1 p for k ≥ 1, in contrast to our geometric M in (1), whose support starts from k = 0. If Y has a strongly nonlattice distribution with moment-generating function Φ(v) ≡ E[e vY ] < ∞ for some v > 0, then Theorem 2 of [START_REF] Blanchet | Uniform Renewal Theory with Applications to Expansions of Random Geometric Sums[END_REF] provides rigorous justification for the "corrected" exponential approximation

P(S > t) ≈ c(p)e -θt (11)
as p → 0, where θ is the unique nonnegative solution of p) , with φ the derivative of Φ. While the right side of (11) approaches a true tail cdf for all t as p → 0, for a fixed p > 0, it may be greater than 1 for some (small) t. [START_REF] Blanchet | Uniform Renewal Theory with Applications to Expansions of Random Geometric Sums[END_REF] suggest two alternative options to estimate θ and c(p):

Φ(θ ) = (1 -p) -1 and c(p) = p (1-p) 2 θ φ (θ ) ≡ e r(
1. Use a numerical root-finding method (Press et al. 2007, Chapter 9) to get θ and then estimate c(p); 2. Develop expressions for θ and r(p) = ln(p)ln[(1p) 2 θ φ (θ )] in powers of p:

θ = ∞ ∑ k=1 χ k p k and r(p) = ∞ ∑ k=0 δ k p k ,
with the χ k and δ k computed via the implicit-function theorem, giving as a variation from [START_REF] Blanchet | Uniform Renewal Theory with Applications to Expansions of Random Geometric Sums[END_REF] -

χ 1 = 1/E[Y ], -χ 2 = 2(E[Y ]) 2 -E[Y 2 ] 2(E[Y ]) 3 -χ 3 = 3(E[Y 2 ]) 2 -6E[Y 2 ](E[Y ]) 2 +6(E[Y ]) 4 -E[Y ]E[Y 3 ] 6(E[Y ]) 5 -δ 1 = 2(E[Y ]) 2 -E[Y 2 ] 2(E[Y ]) 2 -δ 2 = 12(E[Y ]) 4 -12E[Y 2 ](E[Y ]) 2 +15(E[Y 2 ]) 2 -8E[Y 3 ]E[Y ] 24(E[Y ]) 4 -. . .
We now would like to use (11) to approximate the cdf G of S = ∑ M i=1 τ i from (1) (i.e., Y i is replaced with τ i given τ i < T i ) and convolve it with the cdf H of V to approximate the cdf F of T = S +V . However, we cannot directly apply (11) because our S sums M terms, whereas S in (11) has M = M + 1 summands. To account for this difference, note that (also see p. 12 of Kalashnikov 1997) for t ≥ 0,

P(S > t) = P(S > t | M ≥ 1)P(M ≥ 1) + P(S > t | M = 0)P(M = 0) = (1 -p)P(S > t | M ≥ 1) = (1 -p)P(S > t) ≈ (1 -p)c(p)e -θt .
To estimate θ and c(p), we use the above second option as the moments of Y are easily estimated. Specifically, the kth moment is

E[Y k ] = E[τ k |τ < T ] = E[τ k I (τ < T )]/(1 -p), which MSIS estimates by 1 (1 -p n ) n CS n CS ∑ i=1 τ k i I (τ i < T i ).
This yields estimators χ k,n and δ k,n of χ k and δ k , respectively, so we estimate θ and r(p) by

θ n = χ 1,n p n + χ 2,n p 2 n + χ 3,n p 3 n and r n ( p n ) = δ 1,n p n + δ 2,n p 2 n ,
respectively, leading to the corrected exponential estimator of G(t) as

G corr,n (t) = 1 -(1 -p n )e -θ n t+ r n ( p n ) .
We also estimate the cdf H of V by H n in (9). As (8) shows that F can be expressed as a convolution of G and H, we then estimate F by (with

A i = T i ∧ τ i ) F corr,n (t) = 1 p n n IS n IS ∑ i=1 I (T i < τ i )L i t 0 I (A i ≤ t -x) d G corr,n (x) = 1 p n n IS n IS ∑ i=1 I (T i < τ i )L i (t-A i ) + 0 d G corr,n (x) = 1 p n n IS n IS ∑ i=1 I (T i < τ i )L i G corr,n ((t -A i ) + ).
Finally, we get the following expression: Estimator 3 The corrected convolution estimator of cdf F(t) of T is

F corr,n (t) = 1 p n n IS n IS ∑ i=1 I (T i < τ i )L i (1 -(1 -p n )e -(t-A i ) + θ n + r n ( p n ) ).
Again, for sake of comparison and as the optimal MSIS parameter γ is unknown in this setting, our numerical studies (Section 6) take n CS = γn and n IS = (1γ)n with γ optimized for estimator µ n in (6).

GAMMA-BASED ESTIMATOR

We next propose to extend the exponential approximation (3) of cdf F of T by a gamma approximation since a gamma distribution generalizes an exponential. Instead of estimating only the rate of the exponential, we now have to estimate two parameters (shape and rate) for a gamma. Using this generalized family leads to two degrees of freedom instead of one and may lead to smaller bias when p is not very small.

Recall that the cdf F Γ of a gamma distribution with shape parameter α > 0 and rate β > 0 is

F Γ (t; α, β ) = t 0 β α Γ(α) x α-1 e -β x dx,
where Γ(z) = ∞ 0 x z-1 e -x dx is the gamma function. We will approximate the cdf F of T by F Γ , but as α and β in F Γ are unknown, we need to estimate them. Section 7.2 of [START_REF] Casella | Statistical Inference[END_REF] describes several methods for doing this, including maximum likelihood and Bayesian estimators, but we restrict ourselves to the method of moments (MoM). Indeed, the first two techniques require i.i.d. copies T i of T , which may not be possible to obtain in a reasonable amount of time in a rare-event context because each T i typically requires a long simulation. MoM can avoid this issue, as we will explain below.

As F Γ has v = 2 unknown parameters, MoM starts by expressing the first v true (central) moments of F Γ as functions of the distribution's v parameters, and then solves for the parameters in the resulting v simultaneous equations in terms of the moments. Finally, in the expressions for the parameters, replace the true moments with their estimators to obtain estimators of the parameters. Specifically, the mean and variance of T ∼ F Γ are µ = E[T ] = α/β and σ 2 = E[T 2 ]µ 2 = α/β 2 , and solving for α and β yields

α = µ 2 σ 2 and β = µ σ 2 . ( 12
)
(Note that if the scaled T converges in distribution to an exponential as p → 0, then α = µ 2 /σ 2 → 1 as expected (characterizing an exponential distribution).)

We next provide representations for µ and σ 2 that MSIS can efficiently estimate. As in (4), the mean satisfies µ = E[T ∧τ] P(T <τ) = ζ p , which Section 2.1 describes how to estimate by MSIS. For the variance

σ 2 = E[T 2 ] -µ 2 , the only term left to handle is E[T 2 ],
and the regenerative property implies

E[T 2 ] = E[T 2 ; T < τ] + E[(τ + T -τ) 2 ; T > τ] = E[T 2 ; T < τ] + E[τ 2 ; T > τ] + E[(T -τ) 2 ; T > τ] + 2E[τ(T -τ); T > τ] = E[(T ∧ τ) 2 ; T < τ] + E[(T ∧ τ) 2 ; T > τ] + E[(T -τ) 2 | T > τ] P(T > τ) + 2E[τ(T -τ) | T > τ] P(T > τ) = E[(T ∧ τ) 2 ] + E[T 2 ](1 -p) + 2E[τ | T > τ] E[T -τ | T > τ](1 -p) = E[(T ∧ τ) 2 ] + E[T 2 ](1 -p) + 2E[τ; T > τ] µ as E[(T -τ) k | T > τ] = E[T k
] by the regenerative structure. Solving for E[T 2 ] and using (4) then yield Proposition 1 If the process X is regenerative, then

E[T 2 ] = 1 p E[(T ∧ τ) 2 ] + 2E[τI (T > τ)] E[T ∧ τ] p . ( 13
)
From ( 4) and ( 13), we then can obtain expressions for µ and σ 2 as functions of cycle expectations, so we can apply MSIS to estimate them and substitute into (12) to obtain our MoM estimators of α and β . Specifically, for a fixed total number n of cycles to simulate,

• use n CS = γn cycles with CS to estimate E[T ∧ τ], E[(T ∧ τ) 2 ]
, and E[τI (T > τ)];

• use n IS = (1γ)n cycles with IS to estimate p by p n in (5).

Our numerical studies (Section 6) again will use the MSIS parameter γ optimized for estimator µ n in (6).

Having now estimated the moments µ and σ 2 , we estimate the shape and rate parameters of the gamma distribution through (12). We may then compute the resulting cdf (resp., quantile) estimator via efficient algorithms (Press et al. 2007, Chapter 6) for gamma cdf evaluation (resp., inversion). In summary:

Estimator 4 The gamma-based estimator of the cdf F(t) of T is F Γ,n (t) = F Γ (t; α n , β n ), with α n and β n as estimators of the shape and rate parameters, as described above.

BOOTSTRAP ESTIMATOR

All of the previous estimators are based on limit theorems that hold as p → 0. But in an actual system, we do not have p → 0 but rather a fixed p > 0, so the resulting estimators have bias, which does not vanish even as n → ∞. To address these issues, we next consider another approach that applies a type of bootstrap (Casella and Berger 2002, Section 10.1.4). The method is computationally more demanding but may be worthwhile for moderate values of p or when the model is time-consuming to simulate.

The algorithm estimates cdf F by resampling generated observations of T = S +V , where we recall from (1) that S = ∑ M i=1 τ i , M is geometric with parameter p (with support starting from 0), and S ∼ G and V ∼ H are independent. Also, M, τ 1 , τ 2 , . . . are mutually independent in the geometric sum S. Recall that each summand τ i in S has cdf K, where K(x) = P(τ < x | τ < T ). Our algorithm generates observations of T in two steps. The first step (setup) applies MSIS to estimate the geometric parameter p and the cdfs K and H. The second step (bootstrap) generates observations of T as follows: initially sample a geometric M * with the estimated parameter p, then sample M * i.i.d. observations from the empirical estimate of the cdf K, and finally add their sum to a sample from the empirical estimate of H. Next are the details. Estimator 5 The bootstrap estimator F boot,n,n of the cdf F of T is constructed as follows.

Step 1: Set-up

Apply MSIS with a total sample size n = n CS + n IS to estimate p, K, and H. -Use CS to generate n CS i.i.d. pairs (τ i , I (τ i < T i )), i = 1, 2, . . . , n CS , of (τ, I (τ < T )) under the original measure P, and estimate K(x) by

K n (x) = ∑ n CS i=1 I (τ i < x, τ i < T i )/[∑ n CS j=1 I (τ j < T j )].
-Generate n IS cycles by IS to estimate p by p n in ( 5) and the cdf H by H n in (9).

Step 2: Semi-Parametric Bootstrap

Repeat the following resampling procedure n independent times to obtain n conditionally i.i.d. observations T * l , l = 1, 2, . . . , n , of T , given the data from Step 1; i.e., in each iteration l = 1, 2, . . . , n , do the following:

-

Generate M * l ∼ geometric( p n ), i.e., M * l has pmf P(M * l = k) = (1 -p n ) k p n , for k ≥ 0. -Generate M * l i.i.d. observations τ * l,1 , τ * l,2 , . . . , τ * l,M * l
from the empirical distribution K n , and let

S * l = ∑ M * l j=1 τ * l, j be their sum. -Generate V * l ∼ H n , conditionally independently of S * l . -Return T * l = S * l +V * l . After completing all n iterations, construct F boot,n,n with F boot,n,n (t) = (1/n ) ∑ n l=1 I (T * l ≤ t).
Recall that K is the conditional cdf of τ given τ < T , but Step 1 generates (via CS) each (τ i , I (τ i < T i )) not conditional on τ i < T i . Hence, to estimate K from our data, we use only those τ i with τ i < T i , so the empirical cdf K n assigns probability

I (τ i < T i )/[∑ n CS j=1 I (τ j < T j )] to each τ i sampled in Step 1. Also, H n in (9) gives each IS-generated T i ∧ τ i a mass I (T i < τ i )L i /[∑ n IS j=1 I (T j < τ j )L j ]. Thus, although
Step 1 entails simulating the complete stochastic model via CS and IS, Step 2 resamples the generated values from

Step 1, providing a computational savings. Also, the bootstrap estimator is perhaps mainly appropriate when p is not too small; otherwise, the geometric M * l will typically take on huge values, so generating S * l in Step 2 incurs large cost. But the situation when p is not very small is exactly when the exponential approximations to F or G may not be so accurate. While the bootstrap estimator avoids the bias that arises from a weak convergence, as in (2), not holding for a fixed p > 0, the bootstrap incurs another type of bias. As ratio estimators, K n and H n are typically biased, although this bias vanishes as n → ∞, in contrast to exponential-type estimators, whose bias does not change as n grows. For our numerical experiments in Section 6, we will again use n CS = γn and n IS = (1γ)n, and take n = n in the second step.

NUMERICAL COMPARISONS

We ran numerical experiments to compare our estimators of the cdf F of the hitting time T of A for two stochastic models. For the exponential, convolution, corrected convolution, gamma, and bootstrap estimators, we simulated a total of n = 10 4 independent cycles for MSIS, where the allocation parameter γ (proportion of CS cycles) minimizes [START_REF] Goyal | A Unified Framework for Simulating Markovian Models of Highly Dependable Systems[END_REF]) the work-normalized variance of the estimator µ n of µ used in (6). For comparison, we also constructed an empirical estimator of F from directly generating 10 4 i.i.d. values of T . When possible, we also computed the exact cdf F numerically (i.e., no simulation). Our goal is to analyze the bias and error of the estimators. To do this, we generated m = 10 5 independent copies F j,n (t k ), j = 1, 2, . . . , m, of an estimator F n (t k ) at several points t k . Then we estimated the bias, E[ F n (t k )] -F(t k ), and the mean square error (MSE), E[(

F n (t k ) -F(t k )) 2 ], by 1 m m ∑ j=1 F j,n (t k ) -F(t k ) and 1 m m ∑ j=1 F j,n (t k ) -F(t k ) 2 , (14)
respectively, where we recall that MSE decomposes as the sum of variance and squared bias. In our experiments, except for the empirical and bootstrap estimators, the differences in the computation times for all simulation methods are negligible.

Highly Reliable Markovian Systems

We first consider an HRMS, as discussed in Section 2. The system has c = 3 types of components with n l = 3 components of type l (1 ≤ l ≤ c). Each component has an exponentially distributed time to failure with rate ε. Any failed component has an exponentially distributed repair time with rate 1. There are ∑ c l=1 n l repairmen, so failed components never queue for repair. The system is considered down whenever fewer than two components of any one type are operational, and T is the system's time to failure (TTF). We implemented IS using the so-called zero-variance approximation (L'Ecuyer and Tuffin 2012).

We first examine plots of the estimated cdf F of T obtained by the various simulation methods for n = 10 4 . Figure 1 shows the results for ε = 10 -4 on the left and for ε = 10 -2 on the right, while Figure 2 is for ε = 10 -1 . For ε = 10 -4 , we do not show the exact value because it requires too much time to compute; we also omit the empirical (it takes more than one hour to reach the first failure) and the bootstrap (which is potentially of interest only for moderate values of p ε ). For the gamma estimator, the estimated shape and rate parameters are α = 1 (resp., 0.99999 and 1.0108) and β = 6.0292e-15 (resp., 5.7437e-07 and 3.2142e-03) for ε = 10 -4 (resp., 10 -2 and 10 -1 ).

For ε = 10 -2 and ε = 10 -4 , the curves overlap, showing the good behavior of all algorithms. For ε = 0.1 we can see slight differences: the exponential estimator seems less accurate, but we have to check whether it is a general characteristic or if it is due to random noise.

To better analyze the bias and errors, Figure 3 plots the bias and MSE, estimated using ( 14) from m = 10 5 independent experiments, when ε = 0.1, where (L, U) is a 95% CI for the true value F(t). It took 40 minutes to get the curves. As the CIs largely overlap, the differences appear to be very limited. We do the same when ε = 10 -2 in Figure 4, which do not show the empirical and bootstrap estimators because they take too much time to compute. All the methods give equivalent results in terms of bias and MSE. The squared bias is much smaller than the MSE, meaning that, at least for n = 10 4 , bias may not be much of an issue for the CIs, which are based only on variance and do not directly account for the bias.

6.2 M/M/1 Queue For a second example, we consider the process where X(t) denotes the total number of customers at time t in an M/M/1 queue, a special case of the GI/GI/1 queue described in Section 2, with arrival rate λ = 0.5 and service rate µ = 1. We want to estimate the cdf of the hitting time of A ε = { 1/ε , 1/ε + 1, . . . }.

For IS, we swap the arrival and service rates [START_REF] Parekh | A Quick Simulation Method for Excessive Backlogs in Networks of Queues[END_REF]. which takes longer to compute, has much smaller bias but larger MSE due to a larger variance than the other methods. Compared to the HRMS, the M/M/1 queue has more variable cycle lengths. (Indeed, HRMS has very short cycles as returns to the regenerative state are very likely, which is less true for M/M/1.) It is also worth noting that again, even for this moderately rare case, the squared bias contributes little to the MSE. For N = 20, we do not display the results for the bootstrap estimator because it is too computationally demanding. All algorithms are almost identical in terms of bias and MSE.

CONCLUDING REMARKS

We have introduced in this paper new cdf estimators of the hitting time to a rarely visited set of states for a regenerative process. We have numerically observed that the estimators are very close in performance even for a moderately rare situation. It is also interesting to observe that bias is not a significant part of the mean squared error. All this emphasizes the interest of the exponential approximation estimator. As a next step, we plan to apply the same techniques to estimate risk measures, such as quantiles and CTE.

Figure 1 :Figure 2 :

 12 Figure 1: Plots of the estimated and exact cdf of the TTF for (left) ε = 10 -4 with MSIS allocation parameter γ = 0.75370, and for (right) ε = 10 -2 with γ = 0.26390.

Figure 3 :

 3 Figure 3: Bias (left) and MSE (right) at various points t k for the HRMS with ε = 10 -1 when m = 10 5

Figure 4 :

 4 Figure 4: Bias (left) and MSE (right) at various points t k for the HRMS with ε = 10 -2 when m = 10 5

Figure 5 :

 5 Figure 5: Bias (left) and MSE (right) at various points t k for the M/M/1 queue with N = 10 and m = 10 5 independent replications

Figure 6 :

 6 Figure 6: Bias (left) and MSE (right) at various points t k for the M/M/1 queue with N = 20 and m = 10 5 independent replications