Van Hoan Nguyen

Rivier Claire

Effect of the processes (extrusion cooking and/or amylase addition) used to increase energy density of gruels on energy and nutrient intakes of 7 to 10-month-old Vietnamese infants Short running title: Effect of processes on intakes of young children…

Keywords: Instant flour, ready-to-cook flour, α-amylase, energy density, energy intake

Effect of extrusion cooking and/or amylase addition used to increase energy density of gruels on energy and nutrient intakes of 7

INTRODUCTION

In the developing world, 32% (178 million) of children under five (WHO standards) still suffer from stunting. [START_REF] Black | Maternal and child undernutrition study group. Maternal and child undernutrition: global and regional exposures and health consequences[END_REF] In Vietnam remarkable progress has been made in reducing child undernutrition in the last two decades, and, as a consequence, the prevalence of stunting decreased from 59.7% in 1985 to 27.9% in 2006 (NCHS standard). However, malnutrition is still a major public health concern in Vietnam, and the government has implemented a national plan of action to accelerate the reduction of stunting. [START_REF] Khan | A national plan of action to accelerate stunting reduction in Vietnam[END_REF] One of the necessary conditions to protect young children from malnutrition is that households have access to affordable complementary foods of appropriate energy and nutrient densities in addition to breast milk. Studies conducted in the framework of the Fasevie programme demonstrated that it is possible to produce nutritious infant flours from local raw materials using different processes: extrusion cooking with a very low cost extruder (VLCE) designed in collaboration with a Vietnamese institute (Centre of Technological Transfer and Consultancy-CTC) or the addition of amylase. [START_REF] Mouquet | Ability of a "very lowcost extruder" to produce instant infant flours at a small scale in Vietnam[END_REF][START_REF] Bruyeron | Infant and young child feeding in Vietnam : 10 years of experience and lessons from the Fasevie project[END_REF] By enabling partial starch hydrolysis, these processing methods enable the preparation of high energy density gruels.

The effects on quantity and energy intake of feeding children with high energy density gruels have been described in many studies in different countries. [START_REF] Bennett Va Morales | Effects of dietary viscosity and energy density on total daily energy consumption by young Peruvian children[END_REF][START_REF] Moursi | Does the consumption of amylase-containing gruels impact on the energy intake and growth of Congolese infants?[END_REF][START_REF] Traoré | Effects of the duration of the habituation period on energy intakes from low and high energy density gruels by Burkinabe infants living in free conditions[END_REF][START_REF] Islam | Effects of varied energy density of complementary foods on breastmilk intakes and total energy consumption by healthy, breastfed Bangladeshi children[END_REF][START_REF] Islam | Effects of energy density and feeding frequency of complementary foods on total daily energy intakes and consumption of breast milk by healthy breastfed Bangladeshi children[END_REF] However, the effect of the type of processing methods used to produce the infant flours that enable the preparation of high energy dense gruels has not been investigated.

The objective of this study was to measure the intakes (in quantity, energy and nutrients) of gruels prepared from flours produced by three processing methods: extrusion cooking, amylase addition, or a combination of the two, compared to those of a minimallyprocessed gruel prepared from a flour made of crude or roasted raw materials without extrusion cooking or amylase addition.

SUBJECTS, MATERIALS AND METHODS

Recruitment of participants

The study was conducted in two villages in the vicinity of Tam Ky town in Quang Nam province. Before launching the study, a meeting was organized by the coordinator to inform local authorities (village chiefs, health workers, and a representative of the women's union) about the objectives and activities of the study.

All infants aged from six to 10 months during the period planned for the study in the two villages were invited to participate in a general health examination at the municipal health center. A total of 82 infants were eligible, i.e. were still breastfed, were not severely malnourished (weight-for-height and height-for-age z-scores>-3) and their parents were willing to bring them to the canteen located nearest their home twice a day throughout the study period. A letter explaining the objectives and the organization of the study was sent to parents who had to give their written consent. The parents were then invited to come to a meeting during which the objectives and the activities of the study were presented. Written consent was obtained from the mother (or father) of each infant. Nine canteens were set up that could receive around eight infants each day. They were open every day from 7 to 10 a.m. and from 3 to 6 p.m.

During the study, of the 82 eligible infants, 15 (18.3 %) did not eat the gruel regularly because they were sick (5), or because the caregivers were too busy (7) or because the family had relocated far away from Tam Ky (3). Finally, data from 67 infants were used for analyses. For ethical reasons, children excluded for health reasons were allowed to continue to eat the gruels like the other infants.

The protocol of the study was approved by the ethical committee of the National Institute of Nutrition, in Vietnam.

Study protocol

Infants were fed four types of gruels -one type twice a day-for four consecutive periods of four days (Table 1). The gruels were prepared by the canteen staffs who were previously trained for the purpose.

Composition and preparation of the infant flours

The composition on a dry basis of the four infant flours before adding sugar and amylase was based on the same blend of rice (62.2%), sesame (5.9%), soybean (18.5%), milk powder (11.5%) and other ingredients, among which a mineral and vitamin premix (1.9%). The blend was designed to obtain gruels with nutrient densities that meet the latest international recommendations for young children (Table 2). [START_REF] Fao | Human vitamin and mineral requirements, second edition. Report of a joint FAO/WHO expert consultation[END_REF][START_REF] Lutter | Nutrient composition for fortified complementary foods: proposed nutrient composition for fortified complementary foods[END_REF] In order to obtain gruels of the same sweetness level, which is assumed to influence their acceptability [START_REF] Vieu | Effects of energy density and sweetness of gruels on Burkinabe infant energy intakes in free living conditions[END_REF] , their sugar content on a dry basis was calculated so that the sugar content on a wet basis was 3.23 g/100 g of gruel using the targeted dry matter (DM) content of the gruel. The four types of flour, A, B, C and D, were produced using same batches of rice, sesame and soybean. Flour A corresponded to the instant infant flour Favina (Fasevie programme), produced by extrusion cooking of a blend of rice and sesame with addition of roasted soybean and other ingredients. Gruel A was prepared using the instant procedure, i.e. by adding water at 75°C. Flour B was a ready-to-cook flour resulting from the extrusion cooking of the blend of rice, soybean and sesame, and addition of the other ingredients. An industrial powdered αamylase (BAN 800MG, Novo) was added at a rate of 0.001% to obtain the partial starch hydrolysis required for the preparation of gruel of high ED and appropriate semi-liquid consistency. Flours C and D were prepared by blending crude rice, sesame and roasted soybean. BAN amylase was also added at a rate of 0.006% to flour C to enable the preparation of gruel of high ED. Gruels B, C and D were prepared by bringing the mix of water plus flour to the boil over gentle heat and then boiling it at 100°C for 5 min in the case of flour B, and for 10 min in the case of flours C and D. The slow increase in temperature enabled the action of amylase on gruel starch and boiling ensured the inactivation of the enzyme. [START_REF] Trèche | Use of amylases in infant food[END_REF] Gruel D corresponded to the minimally-processed gruel that was used as control.

Extrusion cooking of the rice-sesame mix (flour A) or rice-sesame-soybean mix (flour B) was performed using a VLCE EX-800, which is a simple single-screw autogenous extruder developed by the Fasevie programme [START_REF] Mouquet | Ability of a "very lowcost extruder" to produce instant infant flours at a small scale in Vietnam[END_REF][START_REF] Bruyeron | Infant and young child feeding in Vietnam : 10 years of experience and lessons from the Fasevie project[END_REF] and manufactured by CTC in Vietnam. Soybean hulls were removed in a DH20 dehuller (CTC, Vietnam) and when necessary (flours A, C and D), the soybeans were roasted. Extrudates, roasted soybeans or blends of crude rice and sesame were ground and blended in adequate proportions (particle size < 0.5 mm). The resulting flours were then mixed with the other ingredients of the formula (i.e. sugar, milk powder, vitamin and mineral premix, CaCO3, milk aroma, and when required, BAN amylase). The cost of the different flours was estimated taking into account raw materials and production costs (Table 2). The cost of flours A and B produced by extrusion cooking was slightly higher than that of flours C and D, but this additional cost remained limited thanks to the use of the VLCE produced in Vietnam. One batch of each type of flour was produced at the beginning of the study and packed in individual sachets containing the quantity required to prepare a bowl of gruel for one child (300 g). The sachets were then stored until use.

The four types of flours -along with their appropriate preparation procedure-enabled the production of the gruels A, B, C and D, with the targeted dry matter content, energy density and consistency shown in Table 3.

Preparation and distribution of gruels

In each of the nine canteens, the appropriate gruel was prepared by using the number of sachets equal to the number of children present and the corresponding quantity of water. The canteen staff poured about 300 g of gruel into a previously coded and weighed bowl, and the mothers were asked to spoon-feed the gruel to their child ad libitum. Clean pre-weighed serviettes were given to the mothers to collect any gruel their infant did not consume. Trained personnel, selected among local volunteers, weighed the child, the bowl and the serviette before and after the consumption of the gruel by the child (0.1 g precision), and were ready to intervene in case a child choked on the gruel.

Samples of each prepared gruel were collected for the determination of DM content (by dehydration in an oven until constant weight) and Bostwick flow measurements at 45°C, according to the method described by Mouquet et al. [START_REF] Mouquet | Characterization of the consistency of gruels consumed by infants in developing countries: assessment of the Bostwick consistometer and comparison with viscosity measurements and sensory perception[END_REF]

Nutritional value of infant flours

Samples of each of the four infant flours were collected for biochemical analyses. Protein content (N x 6.25) was determined by the method of Kjeldahl (standard NF. V03-050). [START_REF] Afnor | Directives générales pour le dosage de l'azote avec minéralisation selon la méthode de Kjeldahl (produits agricoles alimentaires)[END_REF] Lipids were extracted with ether oil using the HT6 Soxtec system (Tecator, Höganäs, Sweden) according to the official AOAC method 2003.06. [START_REF] Aoac | Crude fat in feeds, cereal grains, and forages. Randall/Soxtec/hexane extractionsubmersion method[END_REF] Insoluble fibre (Acid Detergent Fibre, ADF) contents were determined by the method of Van Soest 17 using a Dosi-fiber (Selecta, Barcelona, Spain). Ash contents were determined by calcination in a furnace at 530°C. Available carbohydrates were calculated by difference (100 -proteins -lipids -ADFash). Energy contents of the infant flours were then calculated using Atwater coefficients: 4 kcal/g DM (16.7 kJ/g DM) for proteins and available carbohydrates and 9 kcal/g DM (37.6 kJ/g DM) for lipids. Mineral (calcium, iron and zinc) contents were determined by atomic absorption spectrophotometry (SpectrAA 200, Varian, Victoria, Australia). All analyses were made in duplicate and the two values averaged. The energy content and main nutrient composition of each study infant flours are listed in Table 2.

Data management and statistical analysis

Data management and statistical analyses were carried out using SAS software (version 9.1; SAS Institute, Cary, NC, USA). The association between the canteen (from 1 to 9) and the characteristics of the gruels (DM content and Bostwick flow) was assessed by analysis of variance. Concerning intakes, linear mixed models (the MIXED procedure of SAS) were used to evaluate the effect of the type of gruel (A, B, C or D), the day rank of the diet period (from day 1 to day 4), the age of children (in months), the number of the canteen (from 1 to 9) and the time of day the meal was eaten (morning or afternoon) -the last variable was included only for the outcomes expressed per meal. Each model was adjusted for the number of the period (1 to 4) and for all other factors studied. Interactions between the type of gruel and the day rank of diet, the period, the age, and the time of day of the meal were also included. The primary response variables were the amount of gruel consumed and the energy intake (EI) expressed in g or kcal per kg of body weight (BW) and per meal or per day. Different nutrient intakes were also calculated and analyzed. When appropriate, Tukey's test was performed to examine the pairwise levels of significance. In all comparisons, differences were considered significant when p < 0.05.

RESULTS

Characteristics of infants

Among the 67 study subjects, 32 were male and 35 female. The mean age of the study subjects was 7.84 ± 1.33 months and their mean body weight was 7.79 ± 1.22 kg. All children were still breastfed during the study as this was a criterion for inclusion.

Characteristics of study gruels prepared in the canteens

The mean DM content and corresponding ED of the study gruels as well as the mean Bostwick flow as an indicator of the consistency are listed in Table 3. All values were close to those obtained in the laboratory (targets). Thus, as expected, the dry matter content and the Bostwick flow of gruels B and C did not differ significantly. For dry matter content, the standard deviation was low, indicating a relatively good control of the methods of preparation of gruels by the canteen staff. Concerning Bostwick flow, its variability was remarkably low for gruel D, but higher for the three high ED gruels A, B and C. Indeed, the increase in ED is based on partial starch hydrolysis during processing, and the extent of starch hydrolysis obtained through extrusion cooking or amylase action depends greatly on the processing conditions. It may also have been difficult to thoroughly control the temperature of gruel during the measurements made at the canteens and it is known that the consistency of highly concentrated gruels is highly dependent on temperature. [START_REF] Mouquet | Characterization of the consistency of gruels consumed by infants in developing countries: assessment of the Bostwick consistometer and comparison with viscosity measurements and sensory perception[END_REF] Despite satisfactory control of the preparation of the different types of gruel, there was a significant effect of the canteen on both DM content and Bostwick flow for the three high ED gruels. In most cases, a higher Bostwick flow was linked to a lower DM content. The effect of the canteen can possibly be explained by the know-how of the canteen staff or the environmental conditions when water was being heated for the preparation of the instant flour A, or during cooking of gruels B and C, which may have had an effect on the quantity of water that evaporated.

Amounts of gruel consumed and energy intakes

Influence of the day of the period Neither the effect of the day (p = 0.46), nor the interaction of day with the type of gruel (p = 0.64) were significant. As a result, this variable was removed from the models.

Influence of the canteen

The effects of the canteen on the amount of gruel consumed and energy intake were both highly significant. This can be partially explained by differences in the characteristics of gruels due to the way the canteen staff prepared them, but also by the atmosphere of each canteen (from whether the canteen staff smiled or not, or the layout of the canteen facilities).

Influence of the age of the child

The effect of the age of the child on the amount of gruel consumed per meal was not significant (p = 0.67) (interaction non significant p = 0.57).

Influence of the time of day of the meal

Children received two meals of gruel per day, the first in the morning, and the second in the afternoon. The analysis of the effect of the rank of the meal in the day revealed significant differences (p = 0.018) between the two meals, the mean amount of gruel consumed being slightly higher during the second meal of the day. This difference was approximately the same for each type of gruel (interaction between rank of the meal and type of gruel was non significant p = 0.51).

Influence of the type of gruel

Mean amounts of gruel consumed and energy intakes expressed per meal or per day are shown in Table 4. The effect of the type of gruel was highly significant (p < 0.0001) for all outcomes. The largest amount consumed was obtained with the low ED gruel (D) and the amount was significantly lower in the case of the three high ED gruels (A, B and C). Despite this reduction, energy intakes were significantly higher in the case of the three high ED gruels compared to that of the low ED gruel (Table 4). Thus, the energy intake from the two meals of gruel during the day increased by 38%, 66% and 44% for gruels A, B and C respectively, compared to that of gruel D. There were also some differences between the three high ED gruels: the mean amounts of gruels prepared from flours produced by extrusion cooking consumed were significantly higher than that of gruel C prepared from flour with addition of amylase but without any extrusion cooking. And this was true even in the case of gruel B which had the same ED as gruel C. Consequently, the energy intake of gruel B was significantly higher than that of the other two high ED gruels, A and C.

The mean percentages of recommended energy intake were calculated for each type of gruel taking into account the latest figures for total energy requirements [START_REF] Butte | Energy requirements derived from total energy expenditure and energy deposition during the first 2 years of life[END_REF] and low and average energy intake from breast milk [START_REF] Who | Complementary feeding of young children in developing countries: a review of current scientific knowledge[END_REF] (Table 5). Considering average breast milk intake, the consumption of gruels A, B or C provided on average more than 100% and at least 72% of the total energy requirements, respectively, for infants aged 6-8 months and 9-11 months, but only 79% and 56% in the case of gruel D. These figures were much lower when low breast milk intake was taken into account, in this case suggesting the need to increase the frequency of meals. By looking at individuals, only a few children among the 67 in our study were able to cover 100% of the recommended energy intake from complementary foods when only consuming the gruel twice a day. Among the 6 to 8-month-old children, nearly half managed when consuming high ED gruels A, B, C and considering average breast milk intake, but very few children in the 9 to 11-month-old group did. When children get older, the contribution of breast milk decreases, and more than two meals of gruel per day are needed in addition to other foods such as snacks. The increase in DM content by using processing that partially hydrolyses starch in order to obtain high ED gruel has also the considerable advantage of increasing the intake of all the nutrients contained in the infant flour. Thus, in our study, as shown in table 4 for a few examples of nutrients, the nutrient intakes from the high ED gruels A, B and C were substantially higher than those from gruel D. With two meals of gruel D, only the protein requirements were fully satisfied whereas the mean calcium intake only just reached the minimum of the range of recommended daily intake, which means that about half the children did not have sufficient calcium intake. With the high ED gruels, the highest mean nutrient intakes were obtained with gruel B prepared from the flour produced by extrusion cooking and the addition of amylase, values for protein and calcium intakes being significantly higher than those of the three other gruels. Compared to the recommendations for daily nutrient intakes from complementary foods, on average, two meals of gruel A, B or C provide enough nutrients to cover 100% of lipid, protein, calcium and iron requirements.

DISCUSSION

The main objective of this study was to compare energy and nutrient intakes from gruels whose ED was increased using three processing methods (extrusion cooking, addition of amylase or both) and against a low ED gruel as control. Before we comment on the differences between the intakes of the four types of gruel, it should be noted that the mean amount of gruel consumed per meal ranged from 11.8 to 17.7 g/kg of BW and was consequently much lower than the assumed functional gastric capacity of children estimated at 30 g/kg of BW, which is the value most calculations are based on when recommendations for complementary feeding are established. [START_REF] Who | Complementary feeding of young children in developing countries: a review of current scientific knowledge[END_REF][START_REF] Dewey | Update on technical issues concerning complementary feeding of young children in developing countries and implications for intervention programs[END_REF] However, even lower values have been reported in other contexts. [START_REF] Traoré | Effects of the duration of the habituation period on energy intakes from low and high energy density gruels by Burkinabe infants living in free conditions[END_REF][START_REF] Vieu | Effects of energy density and sweetness of gruels on Burkinabe infant energy intakes in free living conditions[END_REF][START_REF] Trèche | Complementary foods in developing countries: importance, required characteristics, constraints and potential strategies for improvement[END_REF] In these studies as well as in ours, children were fed ad libitum and decided by themselves when to stop eating. This suggests that filling the gastric capacity may not be the only factor that determines satiety in children.

The reducing effect of the increase in gruel ED on the amount consumed observed in our study has often been reported and authors also all observed a consecutive increase in energy intakes. [START_REF] Moursi | Does the consumption of amylase-containing gruels impact on the energy intake and growth of Congolese infants?[END_REF][START_REF] Traoré | Effects of the duration of the habituation period on energy intakes from low and high energy density gruels by Burkinabe infants living in free conditions[END_REF][START_REF] Islam | Effects of energy density and feeding frequency of complementary foods on total daily energy intakes and consumption of breast milk by healthy breastfed Bangladeshi children[END_REF] Generally, the amount consumed is inversely linked to the concentration (and therefore to the ED) of the gruel, as was the case with gruels D, A and C in our study. However, the amount of gruel B consumed was higher than that of gruel C although there was no significant difference between the characteristics of the two gruels, either in DM content or in Bostwick flow. Consequently, the difference in the amounts consumed can be attributed to better acceptability or appetence of gruel B, produced by extrusion cooking. During the study, canteen staff also reported that gruel C was not accepted by the children or their mothers as well as the other types of gruel. Indeed, the many biochemical reactions that occur during the extrusion cooking step have been reported to improve the palatability of cereal and soybean based products. [START_REF] Molina | Some characteristics of whole corn:whole soybean (70:30) and rice:whole soybean (70:30) mixtures processed by simple extrusion cooking[END_REF]23 Finally, energy intakes from gruel B prepared by extrusion cooking and addition of amylase were significantly higher than that of gruel A prepared only by extrusion cooking, due the lower ED of the latter (94 vs. 122 kcal/100 g for gruel B), and also slightly higher than that of gruel C prepared only by addition of amylase, due to lower acceptability.

Despite the fact that the infant flour B is not instant and requires a short cooking time that may be perceived as a little more constraining by caregivers, the addition of amylase enables the preparation of a very high ED gruel, significantly higher than the gruel prepared with the instant flour (A). In addition, in contexts where the good hygienic quality of the water is not always ensured, which is the case in many developing countries, the short cooking time needed is a guarantee that the children will receive a safe gruel. Compared to gruel C with only amylase added, the required cooking time for gruel B was shorter, which led us to qualify this type of flour of "ready-to-cook". The additional production cost due to the use of extrusion cooking (+ 9.5 % compared to flour C) remains reasonable thanks to the use of a VLCE. But it could be higher if other types of single screw extruders are used.

In conclusion, the three different processing methods (extrusion cooking using a VLCE, addition of amylases, or both) used for the preparation of high energy density gruels from a fortified rice-sesame-soybean blend proved to be efficient in increasing energy and nutrient intakes by children, compared to a simpler processing method without an extrusion cooking step or the addition of amylase. Among the three processing methods tested, the one that combined extrusion cooking and the addition of amylase resulted in the highest intakes due to higher ED associated with better acceptability. Thus the preparation of ready-to-cook infant flour by combining extrusion cooking and amylase addition is a new way of using the VLCE that may extend its use for the production of nutritious infant flours beyond Vietnam to other developing countries. recommended energy intake from complementary foods (REI-CF) ¶ taking into account the mean intake for the four days of consumption

Mustakas

Table 1 .

 1 GC, Albrecht, WJ, Bookwalter GN, Mc Ghee JE and Griffin EL. Extruderprocessing to improve nutritional quality, flavor, and keeping quality of full-fat soy flour. Food Technol. 1970;24: 1290-1296. Study design

				Period 1			Period 2			Period 3			Period 4
		Day of consumption † 1	2	3	4	1	2	3	4	1	2	3	4	1	2	3	4
	Canteens 1, 2	19 infants	Gruel A			Gruel C			Gruel B			Gruel D
	Canteens 3, 4	13 infants	Gruel B			Gruel A			Gruel D			Gruel C
	Canteens 5, 6	13 infants	Gruel C			Gruel D			Gruel A			Gruel B
	Canteens 7, 8, 9	22 infants	Gruel D			Gruel B			Gruel C			Gruel A
	†	Two meals of gruel per day, one in the morning, one in the afternoon

Table 2 .

 2 Nutritional value and estimated cost of flours used in the studyafter adjusting sugar content

				Type of gruel
		A	B	C	D
		Instant	Extruded + BAN	Raw + BAN	Raw
	Nutritional value (per 100 g DM † )
	Energy value (kcal)	425	423	429	423
	Protein (g)	14.6	14.4	15.2	13.5
	Lipid (g)	8.5	8.1	9.3	7.7
	Calcium (mg)	259.3	282.6	276.2	242.3
	Iron (mg)	20.6	22.0	26.0	23.3
	Zinc (mg)	7.7	8.4	9.8	8.6
	Cost of flour including raw material and production cost (US $ per kg)	1.27	1.27	1.16	1.07
	† DM = dry matter

Table 3 .

 3 Characteristics of the gruels distributed to the infants

		Type of gruel
	A	B	C	D
	Instant	Extruded + BAN	Raw flour + BAN	Raw flour
	n = 72	n = 72	n = 72	n = 72

Table 4 :

 4 Mean amounts of gruel, energy and nutrient intakes as a function of the type of gruel † 11 † p < 0.0001 for all outcomes. All values are mean ± std error, adjusted for canteen, period, and meal (only in the case of data expressed /meal). ‡ Values in a row with different superscript letters are significantly different, p < 0.05 (Tukey's test). § For infants aged 6-11 months with average breast milk intakes (Lutter and Dewey[START_REF] Lutter | Nutrient composition for fortified complementary foods: proposed nutrient composition for fortified complementary foods[END_REF] ).

			Type of gruel		Recommended
		A	B	C	D	daily intake from
		Instant	Extruded + Ban	Raw flour + Ban	Raw flour	complementary foods §
	Amount of gruel ‡
	(g/meal)	118.4 ± 2.5 b	110.0 ± 2.6 b	93.4 ± 2.6 c	140.8 ± 2.5 a
	(g/kg BW/meal)	15.1 ± 0.3 b	13.8 ± 0.3 c	11.8 ± 0.3 d	17.7 ± 0.3 a
	(g/day)	238.1 ± 6.5 b	220.7 ± 6.6 b	186.7 ± 7.0 c	283.2 ± 6.5 a
	Energy intake ‡
	(kcal/meal)	111.8± 2.6 b	134.3 ± 2.6 a	116.7 ± 2.7 b	80.8 ± 2.6 c
	(kcal/kg BW/meal)	14.1 ± 0.3 b	16.8 ± 0.4 a	14.7 ± 0.4 b	10.1 ± 0.3 c
	(kcal/day)	224.8 ± 6.7 b	269.5 ± 6.8 a	233.5 ± 7.1 b	162.7 ± 6.6 c	176
	Lipid (g/day) §	4.5 ± 0.1 b	5.1 ± 0.1 a	5.1 ± 0.1 a	2.9 ± 0.1 c	4.8
	Protein (g/day)	7.7 ± 0.2 b	9.2 ± 0.2 a	8.3 ± 0.2 b	5.2 ± 0.2 c	3.0 -4.5
	Calcium (mg/day)	139.1 ± 4.2 b	178.9 ± 4.3 a	142.7± 4.6 b	100.7 ± 4.2 c	100 -200
	Iron (mg/day)	11.0 ± 0.4 b	14.6 ± 0.4 a	13.6 ± 0.4 a	9.0 ± 0.4 c

Table 5 .

 5 Energy requirements from complementary foods and mean percentage of recommended energy intake as a function of the type of gruel (taking into account 2 meals/day)

	Age group

†

Butte et al.

[START_REF] Butte | Energy requirements derived from total energy expenditure and energy deposition during the first 2 years of life[END_REF]

; ‡ WHO

19 §

ACKNOWLEDGEMENTS

This work was conducted in the framework of the Fasevie programme, realised by the GRET (Groupe de Recherche et d'Echanges Technologiques), IRD (Institut de Recherche pour le Développement) and the NIN (National Institute of Nutrition of Vietnam) with the financial support from Regional Council of Nord-Pas de Calais, General Council of Côtes d'Armor, the French Ministry of Foreign Affairs, the French Committee for Unicef, the European Union and Sight & Life Foundation.

The authors would like to thank the mothers and infants who participated in this study.

AUTHOR DISCLOSURES

Nguyen Van Hoan, Claire Mouquet-Rivier, Sabrina Eymard-Duvernay, Serge Trèche, no conflicts of interest.


Julie Oger

Emmanuel Lesigne

Philippe Leduc

A Random Field Model and Decision Support in Industrial Production

Keywords: kriging, Bayesian inference, Gaussian processes mixture prior, multivariate tdistribution, uncertainty analysis, manufacturing yield evaluation, decision support

We propose a new tool of decision support in front of a globally unknown phenomenon which is modeled by a random field representing simultaneously our knowledge and our lack of information. This tool is the distribution of a random variable called failure risk probability. Before giving the precise definition of this object, we describe an industrial context in which the decision problem occurs and we discuss Bayesian random field model constructions.

Introduction

Extended abstract

In competitive industries, a reliable yield forecasting is a prime factor to accurately determine the production costs and therefore ensure profitability. Indeed, quantifying the risks long before the effective manufacturing process enables fact-based decision-making. From the development stage, improvement efforts can be early identified and prioritized. In order to measure the impact of industrial process fluctuations on the product performances, the construction of a failure risk probability estimator is presented in this article. The complex relationship between the process technology and the product design (non linearities, multi-modal features...) is handled via random process regression (interpolation). A random field encodes, for each product configuration, the available information regarding the risk of non-compliance. After a brief presentation of the Gaussian model approach, we describe a Bayesian reasoning avoiding a priori choices of location and scale parameters. The Gaussian mixture prior, conditioned by measured (or calculated) data, yields a posterior characterized by a multivariate Student distribution.

In the second part of the article we describe the way from a random model to the failure risk probability. Our approach is to consider as random all unknown, inaccessible or fluctuating data. In order to propagate uncertainties, a fuzzy set approach provides an appropriate framework for the implementation of a Bayesian model mimicking expert elicitation. The underlying leitmotiv is to insert minimal a priori information in the failure risk model. The relevancy of this concept is illustrated with theoretical examples. Note that this article comes with supplementary material available on line (hal-00914192).

Settings and motivations

During the fabrication of an industrial product, the manufacturing process cannot be entirely controlled. There is an intrinsic variability (materials, machinery, environment ...) which can significantly deteriorate the characteristics of the manufactured pieces, inducing non-functional parts. Such variations are critical when dealing with complex technical systems, such as for example integrated circuits developed in the microelectronic industry. Indeed, if this fact is covered up, the resulting circuits are (in many cases) either consistently out of specifications or unnecessarily overdesigned.

To investigate the influence of fluctuating parameters, one can perform a corner analysis applying Design of Experiments principles. However, the conclusions drawn from such measurements campaigns are, in general, qualitative and consequently incomplete. Besides, time and money required by physical prototyping are often prohibitive. As an alternative, numerical (deterministic) models implemented in engineering simulation software offer a way to compute the relevant features (thermal, mechanical, electrical ...) of a device. Thus, nowadays, engineers can virtually explore various design configurations and get a deep insight of the final product performances. In this context, statistical studies can be conducted to evaluate the effect of process tolerances on the manufacturing yield. The most popular is the Monte Carlo (MC) method, which consists in sampling configurations of the product parameters according to their probability distributions and to count the failure events. The process is monitored "in line" to roughly assess distributions of environmental variables. This method is easy to implement but its efficiency depends directly on the complexity of the deterministic model. When the numerical simulator considered is time and/or memory demanding, we only get partial information. For instance, the duration of "Finite Elements" analyses, from several hours to several days, is not compatible with a brute Monte Carlo approach.

Data are sparse and it is therefore necessary to propagate the information in the factor space using an analytical representation (emulator). Monte Carlo sampling is then applied to this analytical representation, a fast surrogate of the computer code (simulator), see [START_REF] Pfingsten | Model-based design analysis and yield optimization[END_REF]. In their article, outputs are considered as a single outcome of a Gaussian random field. Our work focuses on a mixture of Gaussian fields, and this choice will be justified in the sequel.

Classical modeling approaches for failure rate (or yield) estimation only keep a small part of the model available information: mean, quantile... Besides, after the model has passed the validation tests, these outputs are taken at face values. We believe such a procedure is hazardous in the specific field of risk assessment. Indeed, it does not measure the impact of the uncertainty introduced by any modeling stage on the only quantity of interest, namely the risk of failure. For example, how can the decision-maker relate the model acceptance criteria to the accuracy of the failure risk estimate? Answering this question in particular is difficult and probably fruitless. Our work proposes a novel and general solution to address these issues. Once the random model has been determined, and given the probability distribution of the product parameters, we go beyond Pfingsten's approach to define the failure risk probability. The predictive uncertainty is not deduced from the posterior model, calculating for instance the conditional variance. The failure risk itself is probabilistic and randomness mirrors the model uncertainty.

Let us be more precise and introduce the basic objects used in the sequel. We consider that the product (meaning each individual manufactured piece) under study is characterized by a number D of numerical factors. Each factor can vary in a given interval what allows the definition of the factor space X ⊆ R D . Each set of factors x (∈ X) determines a numerical value y (x), and the specifications imposed on the product apply to the value of this response y. We can derive from these specifications the out of specifications space A ⊆ R: the product characterized by the set of factors x does not satisfy the specifications if and only if y (x) ∈ A. The factor space X and the out of specifications space A are considered as known, but regarding the deterministic function x → y (x), we have only very partial information.

Besides, a probability distribution P is given on the factor space X. This distribution reflects the factors variability and is considered as known.

As we said previously, the knowledge of y (x) is rarely available for all x. We only have access to a restraint number of data. Suppose we know n deterministic response values (y i = y (x i )) 1≤i≤n respectively for factor set values (x i ) 1≤i≤n . With these data, our goal is to define a random variable named failure risk probability. Its distribution should help for the robust estimation of the product manufacturability.

We begin with the construction of a random field model (Y x ) x∈X of the unknown response (y (x)) x∈X .

Model construction

Several methods are described in the literature to infer a model (Y x ) x∈X from a limited number of available deterministic data

Y x1 = y 1 , Y x2 = y 2 , . . . , Y xn = y n .
Linear regression analysis is the method of choice in the scientific community for estimating relationships among predictor variables. When the phenomenon studied is complex though, the order of the model (usually polynomial), necessary to correctly fit the data, leads to numerous unknowns. This drawback is magnified in a high-dimensional factor space. Indeed, for a p-variable polynomial of degree d, p + d d coefficients should be determined. If the number of data available is less than this limit, the model is singular. Thus, the player is required to make a "blind bet" to enter the game... In addition, the assumptions justifying the statistical model are rarely strictly respected in practice. For example, how should we understand random errors uncorrelated with zero mean when data are obtained through computer experiments, reproducible in essence? These difficulties are overcome by the Gaussian Process (GP) model. The original kriging method was formalized by [START_REF] Matheron | Estimating and choosing. An essay on probability in practice[END_REF] in the geo-sciences. First, a GP indexed by X is selected and then conditioned by the data. This approach has been widely used for different applications such as, geostatistics studies [START_REF] Berger | Objective bayesian analysis of spatially correlated data[END_REF]), optimization [START_REF] Emmerich | Single-and multiobjective evolutionary optimization assisted by gaussian random field metamodels[END_REF]), wind fields modeling [START_REF] Cornford | Modelling frontal discontinuities in wind fields[END_REF]) or design sensitivity analysis [START_REF] Pfingsten | Model-based design analysis and yield optimization[END_REF]), justified by better predictive performances than several other regression methods [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]). Indeed, GP modeling has several interesting properties. Belonging to the interpolation methods family, this model exactly reproduces the observed data set, there is no residual at the observation points. As a consequence, it is an appropriate tool for the analysis of computer experiments [START_REF] Sacks | Design and analysis of computer experiments[END_REF]). GP can be determined (in theory) even for data set of small size, a key property when information is missing which is a common situation when the number of factors (D) is large. Moreover, it is a very versatile model, able to describe non-continuous as well as non-differentiable surface responses. This noteworthy feature is particularly useful to handle response discontinuities, which may occur due to the numerical solving scheme of computer codes (especially when meshing algorithms are involved). Finally, the probabilistic nature of the predictions can be interpreted as a model for uncertainty, a confidence interval representing a degree of belief.

A major and often neglected problem of this method is the a priori choice of parameters (mean, variance and correlation function), left to the responsibility of the data analyst. In order to reduce the arbitrary nature of expert elicitation, we propose to replace the Gaussian field by a mixture of Gaussian fields, therefore avoiding the a priori choice of mean and variance. For this purpose, with the aim of introducing minimal a priori information into the model, mean and variance are considered as random variables with uniform probability distributions. We will show that, in doing so, any posterior random vector follows a multivariate t-distribution (or multivariate Student distribution). We note that a similar result has been obtained with a hierarchical Bayesian model approach described, for example, in [START_REF] Santner | The design and analysis of computer experiments[END_REF].

In the classical approach of yield evaluation based on a random model, the random variable

ω → P ({x ∈ X | Y x (ω) ∈ A})
is considered as the model of the unknown quantity P ({x ∈ X | y (x) ∈ A}) and estimated by Monte-Carlo methods. These Monte-Carlo methods can be optimized, as it is for example presented in the recent article by [START_REF] Auffray | Bounding rare event probabilities in computer experiments[END_REF]. Our approach is different, as it is explained below, and detailed in Section 3.

Failure risk probability viewed as a random variable

The random field (Y x ) x∈X represents the information as well as the uncertainty regarding the values of the response function y. Once this field is built, we can compute for each factor point x ∈ X, the failure probability P (Y x ∈ A).

The reader sees here that we denote by P the probability measure associated to the random field (meaning that the random field is defined on an abstract probability space (Ω, P)). On the other hand, we recall that we denote by P the probability measure assigned to the factor space X.

We decide that there is a failure risk when the failure probability exceeds a fixed threshold α ∈ [0, 1]. So, the failure risk probability is defined by R (α) :

= P (P (Y . ∈ A) > α), that is, R (α) = P ({x ∈ X | P (Y x ∈ A) > α}).
Setting the correct value for the accident threshold α is a tricky problem. On the one hand, the decision will significantly impact the risk assessment of the product. On the other hand, this choice is eminently subjective since it depends heavily on the risk attitude of the individual. That is why we consider α as a random variable. As a consequence, R (α) is also a random variable and has the distribution

R := 1 0 δ R(α) η (dα)
where η is the probability distribution of α and δ t denotes the Dirac mass at point t. We will see hereafter that the uniform distribution is a choice for η which provides an interesting property of mean value preservation.

In practice, R (α) cannot be computed analytically, it is approximated via a Monte Carlo simulation. According to standard Bayesian methods for sampling study, R (α) follows a beta distribution, denoted here by β α . Finally, the distribution of the failure risk probability is defined, if η has been chosen uniform, by

R := 1 0 β α dα.
This reasoning will be described with more details in the sequel.

Contents

In Section 2, several constructions of random field models are discussed. After a reminder about Gaussian fields (2.1), we show that a random field prior with unknown mean (2.2) and variance leads to a conditioned random field following a multivariate t-distribution (2.3), and we discuss on the model implementation (2.4). In Section 3, the construction of the density of the failure risk probability is described, beginning with the elicitation model (3.1). We provide arguments in favor of a risk-neutral attitude, that is a uniform distribution for the accident threshold (3.2). We continue with practical considerations (3.3) to conclude with the description of our global strategy (3.4). Each section ends with illustrating examples.

From Gaussian to Student fields

The unknown function x → y (x) is modeled as a random function. In the absolute sense, if we disregard the variability that is sometimes introduced by numerical solving schemes, the function is deterministic. However, as the numerical solution to a hard problem (described, for instance, by partial differential equations), y does not have in general a closed-form expression. Consequently, the data analyst looking for a behavioral representation, can legitimately think of y as a "black-box" function: it is unknown for any particular configuration until the computer code is actually run. It hence makes sense to postulate a prior model for y, expressing our initial belief regarding y features. Bayesian updating then combines the evidences acquired, the data output (y i = y (x i )) 1≤i≤n , and the prior distribution to yield the posterior distribution.

Contrary to what is often implicitly stated in the literature, we believe that it is far from obvious to get a relevant prior representation for y. Standard priors reflect more of their mathematical tractability than a real understanding of the phenomenon under study. Incorporating doubtful prior information into the model may yield erroneous and overly confident forecasts, a dangerous cocktail in risk assessment. For that reason, we focus in the following on a prior based on quite few assumptions about how the model output relates to its inputs; it is deliberately weakly informative.

Gaussian field

In this first sub-section, we point out that the conditioning of a Gaussian field generates a new Gaussian field and we recall classical formulae giving the conditional densities.

A Gaussian field (Y x ) x∈X is characterized by the fact that all finite dimensional marginal distributions are Gaussian and by the following data:

• mean values: µ (x) := E (Y x ) • covariances: ρ (x, x ) := cov (Y x , Y x ).
Theorem 1. Let x 1 , x 2 , . . . , x n ∈ X and y 1 , y 2 , . . . , y n ∈ R. The conditional distribution of the Gaussian field (Y x ) x∈X given (Y xi = y i ) 1≤i≤n is still Gaussian with

• mean values: E Y x | (Y xi = y i ) 1≤i≤n = µ (x) + ρ (x, (x i ))Σ -1 (y -µ ((x i ))) T • covariances: cov Y x , Y x | (Y xi = y i ) 1≤i≤n = ρ (x, x ) -ρ (x, (x i ))Σ -1 ρ (x , (x i )) T where ρ (x, (x i )) := (ρ (x, x 1 ) , ρ (x, x 2 ) , . . . , ρ (x, x n )), Σ := (ρ (x i , x j )) 1≤i,j≤n is a positive- definite matrix, y := (y 1 , y 2 , . . . , y n ) and µ ((x i )) := (µ (x 1 ) , µ (x 2 ) , . . . , µ (x n )).

This is a classical result in Probability Theory.

A standard solution to our initial problem is to use this theorem with a given a priori Gaussian field (Y x ) x∈X where neither the mean value µ = µ (x) nor the variance σ 2 = ρ (x, x) depend on x. In this case, we will denote k (x, x ) = ρ (x, x ) /σ 2 , the correlation coefficient between Y x and Y x . If the covariance depends only on the difference between x and x , the Gaussian field is said stationary or homogeneous [START_REF] Abrahamsen | A review of Gaussian random fields and correlation functions -2nd edition[END_REF]). Moreover, if the covariance only depends on the Euclidean distance between x and x , the Gaussian field is said isotropic.

In practice, the mean µ, the variance σ 2 and the correlation function k are usually a priori defined by the expert or estimated by means of calibration methods. The most commonly used is the maximum likelihood method which will be discussed in Section 2.4.2.

Note that here and in the sequel, we avoid particular choices of the correlation k, vectors x 1 , x 2 , . . . , x n and numbers y 1 , y 2 , . . . , y n which could lead to degenerate distributions.

Gaussian mixture with random mean

In this sub-section, the prior model is a mixture of Gaussian processes, which differ by translation. This way, it is not required to set a priori the mean value. The posterior distribution of this parameter is deduced from the data via Bayesian inference. We will see that a (non informative) improper prior for the mean parameter can be defined as the limit of a (weakly informative) proper uniform distribution. The limiting step is totally justified here since, on the one hand, it leads to proper posteriors and, on the other hand, inferences have a positive miscalibration (see [START_REF] Gelman | Prior distributions for variance parameters in hierarchical models[END_REF] for details), that is an overestimate (on average) of the variance. The conditioned random field is Gaussian and we give formulae for the mean and variance. Doing so, we recognize expressions initially established by [START_REF] Sacks | Design and analysis of computer experiments[END_REF] and also derived by [START_REF] Santner | The design and analysis of computer experiments[END_REF].

We propose to consider a prior random field Y x = U + W x where U is a real random variable following a uniform distribution on an interval [-m, m] and where (W x ) x∈X is a centered Gaussian field with constant variance. Moreover, we suppose that U and (W x ) x∈X are independent.

The parameters characterizing the Gaussian field (W x ) x∈X are the variance σ 2 and the correlation function k (recall that the mean is zero). So, for all

x, x ∈ X, cov (W x , W x ) = σ 2 k (x, x ) and, for all x ∈ X, k (x, x) = 1. Let x 1 , x 2 , . . . , x n ∈ X and y := (y 1 , y 2 , . . . , y n ) ∈ R n . Denote Σ := (k (x i , x j ))
1≤i,j≤n the positive-definite matrix of correlations and k (x) := (k (x, x j )) 1≤j≤n the correlation vector.

Theorem 2. The conditional distribution of the random field (Y x ) x∈X knowing that (Y xi = y i ) 1≤i≤n is given by explicit formulae for the densities of finite dimensional marginals. When the parameter m goes to infinity, this conditional distribution becomes Gaussian. In particular, when m → ∞, the univariate conditional distribution of the random variable Y x becomes Gaussian with mean

µ + k (x)Σ -1 (y -µ1) T with µ := yΣ -1 1 T 1Σ -1 1 T (1)
and variance

σ 2   1 -k (x)Σ -1 k (x) T + 1 -1Σ -1 k (x) T 2 1Σ -1 1 T    (2)
where

1 = (1, 1, . . . , 1) ∈ R n . Remark. The mean µ + k (x)Σ -1 (y -µ1)
T can also be written

1Σ -1 1Σ -1 1 T 1 -k (x)Σ -1 1 T + k (x)Σ -1 y T .
Note that an expression similar to (1) is obtained in Section 2.5.

Proof of Theorem 2. We look for the distribution of the random field

(Y x ) x∈X given (Y xi = y i ) 1≤i≤n .
Let r be a positive integer and (t 1 , . . . , t r ) ∈ X r . We have:

(Y t1 , . . . , Y tr , Y x1 , . . . , Y xn ) = (U, . . . , U ) + (W t1 , . . . , W tr , W x1 , . . . , W xn )
and (W t1 , . . . , W tr , W x1 , . . . , W xn ) follows the distribution N (0, ∆) with

∆ := σ 2 Σ 2 k (t) k (t) T Σ
, a positive-definite matrix

where k (t) := (k (t i , x j )) 1≤i≤r,1≤j≤n and Σ 2 := (k (t i , t j )) 1≤i≤r,1≤j≤r .

Denote by f the density of the random vector (Y t1 , . . . , Y tr , Y x1 , . . . , Y xn ):

f (ζ) = 1 2m m -m 1 √ 2π n+r |∆| exp - 1 2 (ζ -u) ∆ -1 (ζ -u) T du
where ζ := (y t1 , y t2 , . . . , y tr , y x1 , y x2 , . . . , y xn ) and u :

= (u, u, . . . , u) ∈ R n+r . The conditional density of (Y t1 , Y t2 , . . . , Y tr ) given (Y xi = y xi ) 1≤i≤n is z := (y t1 , y t2 , . . . , y tr ) → f (z, y) R r f (z , y) dz =: g (z) .
After simplification, we get:

g (z) = m -m exp -1 2 z y -u ∆ -1 z y -u T du R r m -m exp -1 2 z y -u ∆ -1 z y -u T dz du
.

By monotone convergence, some calculations show that:

g (z) -→ m→+∞ exp   -1 2   z y ∆ -1 z y T - 1∆ -1 z y T 2 1∆ -1 1 T     R r exp   -1 2   z y ∆ -1 z y T - 1∆ -1 z y T 2 1∆ -1 1 T     dz .
Within the exponential in the numerator of this expression we identify a non-negative second degree polynomial of the variable z. We recognize a Gaussian distribution (of dimension r).

In the limit situation m → ∞, we note that the distribution of the conditioned random field is welldefined and is Gaussian. Some calculations show that, at the point x ∈ X, the one-dimensional marginal Gaussian distribution of the field has mean

µ + k (x)Σ -1 (y -µ1) T with µ := yΣ -1 1 T 1Σ -1 1 T and variance σ 2   1 -k (x)Σ -1 k (x) T + 1 -1Σ -1 k (x) T 2 1Σ -1 1 T    .

Gaussian mixture with random mean and variance

We go one step further, avoiding the choice of the mean and variance. The prior model is a mixture of Gaussian processes, which differ by affine transformation. Once again, the conditioned random field is well-defined selecting (weakly informative) proper uniform priors for the unknown mean and standard deviation. When the distributions supports respectively tend to the entire and positive real lines, the multivariate Student distribution arises as the random field posterior distribution. We give explicit formulae for the location and scale parameters. [START_REF] Santner | The design and analysis of computer experiments[END_REF] and [START_REF] Kato | Improved prediction for a multivariate normal distribution with unknown mean and variance[END_REF] describe similar models with different priors. For example, with Jeffreys priors for the mean and standard deviation (scale-invariance property), Santner and al. come to the same posterior except for the number of degrees of freedom. We justify our choice arguing that a positive miscalibration is always preferred by the risk-averse individual.

Multivariate Student distribution

First, recall the definition of a multivariate Student distribution. We refer to the book of [START_REF] Kotz | Multivariate t-distributions and their applications[END_REF]. A p-variate Student distribution (or a multivariate t-distribution) has density:

t → 1 ( √ πν) p |Σ| Γ ν+p 2 Γ ν 2 1 + 1 ν (t -µ) Σ -1 (t -µ) T -ν+p 2 (t ∈ R p )
where the positive integer ν is the number of degrees of freedom, Σ is the p × p positive-definite matrix of scale parameters and µ is the 1 × p vector of location parameters. For ν > 1, the mean vector of the Student distribution is well-defined and equals to µ. For ν > 2, the covariance matrix of the Student distribution is well-defined and equals to ν ν-2 Σ. This is a multi-dimensional generalization of the Student distribution. When ν = 1, the distribution is a multivariate Cauchy distribution. When ν goes to infinity, the distribution tends to a multivariate Gaussian distribution.

Student field posterior

Now, we propose to consider a prior random field Y x = U +V W x , where U is a real random variable following a uniform distribution on an interval [-m, m], where V is a real and positive random variable following a uniform distribution on an interval [ , 1/ ] and where (W x ) x∈X is a centered normalized Gaussian field. Moreover, we suppose that U , V and (W x ) x∈X are independent.

The parameter characterizing the Gaussian field (W x ) x∈X is the correlation function k (recall that the mean is zero and the variance is 1). We suppose here that n ≥ 3. Let x 1 , x 2 , . . . , x n ∈ X and y := (y 1 , y 2 , . . . , y n ) ∈ R n . Denote Σ := (k (x i , x j )) 1≤i,j≤n the positive-definite matrix of correlations and k (x) := (k (x, x j )) 1≤j≤n the correlation vector.

Theorem 3. The conditional distribution of the random field (Y x ) x∈X knowing that (Y xi = y i ) 1≤i≤n is given by explicit formulae of densities of finite dimensional marginals. When the parameter m goes to infinity and goes to zero, for n > 2, this conditional distribution becomes a multivariate Student distribution. In particular, when m → ∞, → 0 and n > 2, the univariate conditional distribution of the random variable Y x becomes a Student distribution with n -2 degrees of freedom, with location parameter

µ + k (x)Σ -1 (y -µ1) T with µ := yΣ -1 1 T 1Σ -1 1 T and scale parameter 1 n -2 (y -µ1) Σ -1 y T   1 -k (x)Σ -1 k (x) T + 1 -1Σ -1 k (x) T 2 1Σ -1 1 T    where 1 = (1, 1, . . . , 1) ∈ R n .
The proof of Theorem 3 results from the same reasoning than the proof of Theorem 2 with more complex formal calculations. It is developed in the supplementary materials available on line (hal-00914192).

Remark. The expression (y -µ1) Σ -1 y T can also be written (y -µ1) Σ -1 (y -µ1)

T . Note that a similar expression is obtained in Section 2.5.

Model implementation

Discussion on correlation functions

The a priori choice of the correlation function is one of the major difficulties in random field modeling [START_REF] Sacks | Design and analysis of computer experiments[END_REF]) and the discussion is unfortunately often avoided. It is sometimes justified by expert knowledge [START_REF] Cornford | Modelling frontal discontinuities in wind fields[END_REF]) in particular applied cases. [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] give a complete list of correlation functions, classified according to practical considerations. The most popular correlation families are isotropic: Exponential, Matern and Rational Quadratic classes [START_REF] Abrahamsen | A review of Gaussian random fields and correlation functions -2nd edition[END_REF]). Indeed, an isotropic random field does not suffer the curse of dimensionality since the number of parameters in the model does not depend on the factor space dimension (D).

Determining appropriate values for the parameters of the correlation function is the purpose of calibration methods. The most popular of them is the maximum likelihood estimation (MLE) described in the next sub-section. [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], [START_REF] Stein | Interpolation of spatial data: some theory for kriging[END_REF] and [START_REF] Robert | The Bayesian choice: from decision-theoretic foundations to computational implementation -2nd edition[END_REF] describe this method and discuss about its capabilities and limitations such as numerical issues for too big samples, multiple optima for too small samples and over-fitting problems for which the sample is well learned but the unknown function values are poorly predicted everywhere else in the factor space.

In order to overtake these different limitations, a good way seems to try different choices of correlation function; each alternative is then evaluated by checking if known values of the responses are compatible with the confidence intervals predicted by the random field. The literature describes several types of cross-validation methods measuring the predictive capability of a model [START_REF] Currin | A bayesian approach to the design and analysis of computer experiments[END_REF]). It is of course necessary to use a set of data which is not involved in the construction of the a posteriori model. Because we cannot, when information is scarce, afford not to include all the gathered data in the final model, the partition into training and validation sets is only temporary, for the purpose of the cross-validation stage. In that case, the test is therefore performed on an "incomplete" posterior model version. [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] describe an alternative method which consists in choosing a parameterized family of correlation function and a prior distribution on its parameters, to construct a hierarchical model. This method involves analytical approximations of integrals. Markov chain Monte Carlo (MCMC) methods are popular solutions to make these computations, see [START_REF] Robert | The Bayesian choice: from decision-theoretic foundations to computational implementation -2nd edition[END_REF] for detailed description of these methods. A major drawback is here the cost of such calculations.

Maximum Likelihood Estimation

The classical maximum likelihood estimation (MLE) method defines an estimator of the unknown parameter vector θ of a probability distribution f θ . This estimator is the value θ max which maximizes the density distribution of a random sample calculated at the observed value of this sample.

Under the same name, this method has been adapted to the identification of the parameters of an unknown random field (Y x ) x∈X when considering a family of values Y x1 , Y x2 , . . . , Y xn . It is important to have in mind that the observation set is a single outcome of the random field so that the validity of the MLE method, in this context, is not obvious. However, this topic is outside the scope of this article.

Considering now the probabilistic modeling of a deterministic phenomenon, the MLE provides a practical procedure to set the parameters of the a priori random field (Y x ) x∈X knowing Y xi = y i , 1 ≤ i ≤ n. This is a classical approach in Gaussian field modeling, see [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], [START_REF] Stein | Interpolation of spatial data: some theory for kriging[END_REF] and [START_REF] Robert | The Bayesian choice: from decision-theoretic foundations to computational implementation -2nd edition[END_REF].

Even if our justifications are incomplete, let us describe the MLE for the random model described in Section 2.3.

First, set a parametric family of correlation functions k θ depending on the parameter vector θ, which defines also the parametric correlation matrix Σ θ = (k θ (x i , x j )) 1≤i,j≤n . The MLE consists in choosing θ which maximizes the density of the random vector (Y x1 , Y x2 , . . . , Y xn ) according to m and :

f m, (θ) = 1 2m 1 1 - m -m 1 1 √ 2π n v n |Σ θ | exp - 1 2v 2 (y -u) Σ θ -1 (y -u) T du dv
with y := (y x1 , y x2 , . . . , y xn ) and u := (u, u, . . . , u) ∈ R n . The maximum likelihood estimator θ max (m, ) maximizes f m, (θ). It maximizes also:

fm, (θ) = m -m 1 1 v n |Σ θ | exp - 1 2v 2 (y -u) Σ θ -1 (y -u) T du dv.
We don't know any analytic expression of fm, (θ) and a fortiori of θ max (m, ) but since we are interested in large values of m and little values of , we consider the limit value f∞,0 (θ). So we propose to study:

f∞,0 (θ) = +∞ -∞ +∞ 0 1 v n |Σ θ | exp - 1 2v 2 (y -u) Σ θ -1 (y -u) T du dv.
A short calculation gives

f∞,0 (θ) = √ π2 n-3 2 nσ 2 θ 2-n 2 Γ n-2 2 |Σ θ | 1Σ θ -1 1 T with σ 2 θ := 1 n (y -µ θ 1) Σ θ -1 (y -µ θ 1) T and µ θ := 1Σ θ -1 y T 1Σ θ -1 1 T .
Maximizing f∞,0 (θ) is equivalent to minimize -ln f∞,0 (θ) and after simplifications, we get:

θ max (∞, 0) minimizes n ln σ 2 θ + ln (|Σ θ |) + 2 ln 1Σ θ -1 1 T σ 2 θ .
Observe the difference with the Gaussian case in which the MLE method proposes analytical expressions for the mean and the variance of the field and where the parameter vector θ of the correlation function is estimated by minimizing:

n ln σ 2 θ + ln (|Σ θ |) .

Comparison between Student and Gaussian fields

Let us come back for a while to the Gaussian field model described in Section 2.1, choosing constant mean and variance. If the mean µ and variance σ 2 are determined by the MLE method considering the data set (y (x i ) = y i ) 1≤i≤n , then the following formulae are obtained:

µ = yΣ -1 1 T 1Σ -1 1 T and σ 2 = 1 n (y -µ1) Σ -1 (y -µ1) T .
This is a classical result in this area. See for example [START_REF] Currin | A bayesian approach to the design and analysis of computer experiments[END_REF] or [START_REF] Stein | Interpolation of spatial data: some theory for kriging[END_REF].

Following Theorem 1, the a posteriori mean and variance of the Gaussian random variable Y x are respectively

µ + k (x)Σ -1 (y -µ1) T and σ 2 1 -k (x)Σ -1 k (x)
T .

It is interesting to compare these values with those given in Theorem 3. Note that, for the same correlation function, the mean remains the same and the variance is greater. We can illustrate this difference on a simple example in one dimension. Figure 1 shows, for each field after conditioning, the mean and the boundaries of an interval computed with a confidence of 0.90. As expected, the confidence intervals are both null at the 5 sampled points since the distributions at these points degenerate in Dirac distributions. The confidence interval of the Student posterior is larger than the Gaussian posterior equivalent, suggesting a lower degree of belief in the prediction. The Student field seems more reliable, since we take into account all the possible values of the mean µ and the variance σ 2 , that were arbitrarily fixed in the case of the Gaussian field.

Decision support

We present here a new approach in order to define a failure risk probability. As it is briefly explained in Introduction, it applies to any deterministic phenomena on which we have only partial information, when this information and our uncertainty are described by a random field model.

Failure risk probability

We use here notations introduced at the end of Sub-section 1.1. The real failure risk is the measure of the real failure set {x ∈ X | y (x) ∈ A} where X ⊆ R D is the factor space provided with a probability distribution P . In a lot of practical cases, the function x → y (x) is not available though; it can only be sampled at a few x locations. According to Section 2, a random field (Y x ) x∈X defined on a probability space (Ω, P) is derived from a set of observations in order to approximate the function y.

This modeling stage is not without consequences on the approach retained to calculate the failure risk probability. Since the knowledge on the function y is low, it seems reasonable to provide the failure risk prediction with a confidence measure. Thus a decision-maker could easily evaluate the model quality and rule on the manufactured product robustness more safely. This problem can be formulated as follows: how can we propagate the uncertainty inherent to any prediction based on a random field to the failure risk estimation?

If the phenomenon under study was really a random field (Y x ) x∈X , we could consider the random variable

Z := ω → P ({x ∈ X | Y x (ω) ∈ A}) on (Ω, P)
as the failure risk and its distribution as the failure risk probability. This is the definition usually retained in risk assessment studies. See for example the work of [START_REF] Auffray | Bounding rare event probabilities in computer experiments[END_REF].

Strictly speaking, the model (Y x ) x∈X cannot describe the reality, which is perfectly determined (but unknown). It correctly represents our knowledge at the observation points. Besides, we have to give an interpretation to the process randomness. If we retain the subjective conception of probability theory, randomness results from the incomplete knowledge of the quantity y (x). Indeed, several realizations of the conditioned random field meets the requirements of the modeling problem, that is interpolates the observations. Within this framework, the measure P quantifies a degree of belief regarding the model forecasts. Randomness becomes a way to assess the quality, i.e. the predictive capability of the model. The expression P (Y x ∈ A) = 0.5 does not mean that the point x has a one-in-two chance of belonging to the real failure set; it means that, regarding the available information, one random field realization out of two states that x belongs to the real failure set. The model is therefore inadequate at point x.

Even if the previous definition of the failure risk by the random variable Z is absolutely correct mathematically, we claim that it cannot be used, in this particular case, to model incomplete knowledge. Indeed, the distribution of Z can be a Dirac mass even in the case of a bad model. An example is necessary here to clarify our point. According to the Lemma demonstrated in Appendix A.1, if the measure P is continuous and if the random variables Y x are pairwise independent and satisfy P (Y x ∈ A) = , then the distribution of Z is a Dirac mass in . As a consequence, the particular case where is near 1/2, which is the sign of a poor model, is not incompatible with the fact that the distribution of the random variable Z is a Dirac mass1 . However, a Dirac mass distribution, or a distribution close to a Dirac mass, represents a phenomenon for which the information is sure, or close to sure. This does not reflect the degree of insecurity associated to an inadequate model. Consequently, we conclude that the distribution of the random variable Z = P (Y . (ω) ∈ A) does not provide relevant information regarding the model quality.

The random field alone is not a knowledge representation; it cannot encode the mental procedures involved in belief assessments. An additional analysis is hence necessary. We certainly do not have at our disposal a strict condition of membership to the real failure set. Nevertheless, we can easily compute the following failure probability, a random variable seen as a membership function M:

(X, P ) → [0, 1] x → M (x) := P (Y x ∈ A) .
The quantity M (x) is the proportion of the random field outcomes concluding that x belongs to the failure set. Such a function, only involves marginal distributions, and can be straightforwardly evaluated. It associates to each point x ∈ X a real number in the interval [0, 1] measuring a grade of membership of x to the failure set. Thus, the crisp failure set, unknown, can be only imprecisely characterized by the membership function M. Due to incomplete knowledge, we can only access to the fuzzy version of the failure set. According to [START_REF] Zadeh | Fuzzy sets[END_REF], "such a framework provides a natural way of dealing with problems in which the sources of imprecision is the absence of sharply defined criteria of class membership rather than the presence of random variables." By definition, the notion of "belonging" is not well-defined for a fuzzy set. For a given grade of membership M (x), a decision-maker may consider that x belongs to the failure set whereas another may not, depending on its own risk tolerance. In order to account for the subjective nature of the decision, we introduce the level α, 0 ≤ α ≤ 1, and agree to say that "x belong the failure set" if M (x) > α. The threshold α determines the status to give to uncertain predictions. These data are therefore aggregated in two classes: "good" and "poor". Let us define the α-level failure set:

{x ∈ X | M (x) > α} .
It is a crisp subset of X and can be consequently measured:

R (α) := P ({x ∈ X | M (x) > α}) .
It is the failure risk probability estimation by a decision-maker whose risk tolerance is α. We implicitly define here a causal model: the risk probability is conditional on the risk tolerance α [START_REF] Pearl | Do we need higher-order probabilities, and if so, what do they mean?[END_REF]). The decision-maker risk tolerance being unknown, we choose to consider the threshold α as random. It takes its values in [0, 1] and we denote by η its distribution. In this framework, the failure risk probability is finally defined as a random variable whose distribution is the image of the measure η by the function R. In other words, this distribution is

R := 1 0 δ R(α) η (dα) .

Uniform distribution for the tolerance threshold

In this sub-section, we discuss the choice of the distribution η of the threshold α.

A natural output in risk assessment is the average failure probability, E (P (Y . ∈ A)) in our case. It is an estimator of the P -measure of the real failure set, based on knowledge of the values (y i = y (x i )) 1≤i≤n . As we will see in the sequel, the causal model introduced in Section 3.1 can be interpreted as an operator acting on the space of probability measures on [0, 1], which links the distributions of the random variables R (failure risk probability) and M (failure probability). It seems clever to keep the first moment of the latter when the transform is applied. In that way, the mean value of the failure risk probability remains an estimator of the failure set measure. This remark leads us to the first proposition.

If the distribution of the failure probability is uniform, we cannot conclude about the model relevance, defined as its ability to discriminate the factor space (i.e. identify the failure set). In such situation, the distribution of the failure risk probability should also be uniform, meaning that the state of uncertainty has been preserved. Based on this desired invariance property, we state the second proposition.

First of all, we introduce some definitions and notations, useful for the subsequent demonstrations. As above, for each α ∈ [0, 1], we set: R (α) := P (M (x) > α) .

The function R is the complementary cumulative distribution function (or tail distribution) of the random variable M (x), defined on the probability space (X, P ). Now, we can also look at R as a random variable defined on the probability space ([0, 1] , η), and we define, for all t ∈ [0, 1]:

G (t) := η ({α ∈ [0, 1] | R (α) > t}) = η (R > t) .
So G is the tail distribution of the random variable R.

Moreover, we denote by K the cumulative distribution function of η: K (u) = η ([0, u]) for all u ∈ [0, 1]. Finally, let us recall the definition of the generalized inverse of the decreasing function R.

Definition (Generalized inverse). Let S a decreasing right continuous function on the interval [0, 1] with S (0) ≤ 1 and S (1) = 0. For all t ∈ [0, 1], we set:

S -1 (t) = sup {α ∈ [0, 1] : S (α) > t} .
At this stage, it is not difficult to check that G = K • R -1 . Let us define in a general setting the operator which maps G onto R. The probability measure η is assumed fixed. To any probability measure m on the interval [0, 1], we can associate another probability measure m on [0, 1] in the following way:

• firstly denote by F m the function F m (α) = m ((α, 1]), • secondly define F m = K • F -1
m . Denote by L η the operator, acting on the space of probability measures on [0, 1], which associates m to m.

Proposition 1.

1. The operator L η preserves the first moment of the probability (i.e. 2. If the probability η is uniform, then the operator L η is its own inverse, meaning that (m ) = m.

Proof of Proposition 1. The first moment of the probability measure m associated with F m is equal to:

1 0 F m (α) dα.
The following neat identity can be proved:

1 0 F m (α) dα = 1 0 F -1 m (α) dα.
Thus, the first moment of m is equal to the first moment of m if and only if:

1 0 K F -1 m (α) dα = 1 0 F -1 m (α) dα.
It can be easily verified that the identity map is the unique function K which satisfies this equality for all choice of m. This means that η has to be a uniform distribution. It is not difficult to check point 2.

Proposition 2. The uniform distribution η is the unique distribution for which the operator L η applied to a uniform distribution gives a uniform distribution.

The proof does not present any difficulty.

In the sequel, η is chosen as the uniform distribution. It results from proposition 1 that E η (R) = E (P (Y . ∈ A)), where E (P (Y . ∈ A)) is an estimator of the real failure risk P (y ∈ A).

At this stage of our presentation, the distribution of the failure risk probability is defined as follow. The relevance of this formula is illustrated in Appendix A.2, on the toy example introduced in Section 3.1.

R := 1 0 δ R(α) dα.

Monte Carlo approach

We describe in this sub-section why we have to add a new layer of Bayesian randomness, due to our practical approach of the calculation of the risk probability R(α).

In practice, we use a MC method to get a numerical approximation of R (α). More precisely, an importance sampling is performed: independent samples (x m ) 1≤m≤M , in the factor space X, are generated from the distribution P . For more details on the related preferential MC method, see 4.5 p. 19 of [START_REF] Caflish | Monte carlo and quasi-monte carlo methods[END_REF]. For each point x m , the failure probability M (x m ) = P(Y xm ∈ A) can be computed numerically since the Student distribution is tabulated. Denote by n (α) the number of MC draws

x m such that M (x m ) > α. The MC estimator of R (α) is the quotient n (α) /M .
Note that a classical Bayesian reasoning is used to model the uncertainty due to finite sampling. In this approach, the MC sampling is considered as a binomial experiment. This assumption is thoroughly justified for a good quality pseudo random number generator. As a consequence, R (α) is a random variable whose distribution is the beta distribution β α with shape parameters (n (α) + 1, M -n (α) + 1), assuming a uniform prior. Adding this additional probabilistic stratum in the model, we may appear a little bit pernickety. Nevertheless, it guards us against the basic (but widespread) fallacy which consists in presuming that a system is perfectly safe as long as no accident has been observed.

In conclusion, taking into account the numerical aspects of the problem, we define the distribution of the failure risk probability as

R := 1 0 β α dα.

Global strategy

In this sub-section, we gather together Sections 2 and 3, and we describe the global strategy that has been implemented in the practical application of our present work. From a set of virtual experiments, the computer code emulator (Student field) is derived. A stochastic simulation is then performed (uncertainty analysis) to extract the distribution of the failure risk probability. Standard statistical quantities such as the mean, the standard deviation and a confidence interval are available to the decision-maker for risk assessment. The confidence interval indicates the reliability of the estimate, the sources of uncertainty being the partial knowledge of the original function y and the finite MC sampling.

If the prediction quality is too low, new data points from the space X have to be added to the observations set {x 1 , x 2 , . . . , x n }. Let us describe succinctly a process of choice for the new data points. Note first that the MC simulation only requires model evaluations. Consequently, we have full scope to complete the observations set; it is in no way linked to the factors distribution P . We could sample data uniformly in X but, as each datum "costs the earth", it is strongly recommended to structure the data collection process. This issue is addressed by experimental design techniques, which intend to optimize the information gathering. We retain the differential entropy h (Y x ) := E (-log f (Y x )), where f is the density function of the random variable Y x , as a measure of the lack of information [START_REF] Cover | Elements of information theory[END_REF]). The information provided to the probabilistic model is increased if observations are made at locations in X where the entropy is maximum. We face a multimodal optimization problem : the objective is to find a set of local maxima for h (Y x ). It can be solved by conventional methods such as gradient descent or metaheuristics as evolutionary algorithms. We refer to [START_REF] Talbi | Metaheuristics: from design to implementation[END_REF] for a complete presentation of metaheuristics. Once the new observation points are identified, we are back to step one... The iterative procedure is repeated until the precision requirements regarding the failure risk probability are met or the due date for the risk analysis is reached.

Example and results

In order to illustrate the interest of the risk assessment scheme presented in this article, three theoretical examples described in Appendix B are considered. Using test functions with closedform expressions, the failure probability can be calculated exactly. Therefore, estimates can be compared to the true value.

The Gaussian mixture prior defined in Section 2.3 is kitted out with the anisotropic γ-exponential correlation function (see Appendix C) and trained using the MLE method described in Section 2.4.2. The resulting model has D + 1 parameters, a reasonable level of complexity. The factor distribution P is assumed uniform.

As an alternative to our model-based Monte Carlo (MMC) method, we consider the brute-force Monte Carlo (BMC) method plus the Bayesian model for finite sampling, recalled in Section 3.3. Let k the number of defectives samples, i.e. such that y (x) ∈ A, from M trials. The inferred failure risk probability follows the beta distribution with shape parameters (k + 1, M -k + 1). Note that this trivial model implicitly assumes a sampling of the factor space according to the distribution P (uniform). We used a quasi-random rather a pseudo-random source: (x m ) 1≤m≤M are chosen as elements of the Sobol low discrepancy sequence. Such sampling covers the factor space more evenly, a desirable property to obtain a relevant statistical population. For the MMC and the BMC methods, the mean value and the boundaries of the smallest confidence interval, stated at the 90% confidence level, are plotted as a function of the number M of evaluations, 10 ≤ M ≤ 1000.

Before going further, it is important to have in mind that a confidence interval only provides a statistical estimation of the error on the result. It obviously does not imply a strict condition of membership to the interval. Indeed, whatever the sampling method used, this one may not catch the essential features of the function y. As a result, uncertainty can be underestimated, a dangerous situation in risk assessment. This is particularly true for BMC, which does not take into account the spatial distribution of the data (geometrical structure of the space X). The third example well illustrates this point. This is the key difference with the MMC method we propose. Besides, and this is the cornerstone of this article, we aim at taking into account the uncertainty introduced by the interpolation stage. As a consequence, the comparison of the convergence properties of MMC and BMC should not only focus on confidence intervals. The uncertainty model is more relevant in the second case.

Quadric example

The characteristic of the MMC method is to warn the user if the emulator is poor. It is justified in practice, though, if it can outperform the BMC method at least in some cases.

So, let us start with the quadric example, considering a space X of dimension D = 5 and a true failure risk r t = 1.008 × 10 -1 (the theoretical expression is available in the supplementary materials). See Figure 2 for a comparison of the convergence rates of the model-based and bruteforce Monte Carlo. It turns out that, in this case, MMC definitely outperforms BMC since it estimates the real failure risk as accurate as BMC on 600 samples, requiring a third of the function evaluations.

Although MMC has an extra source of uncertainty, its confidence interval is smaller. At least 60 samples are necessary to get a confidence interval size inferior to 0.05. This lower bound soars to 400 data points using BMC. Besides, the MMC confidence interval always includes the failure risk true value. All this suggests that the prior assumptions of the model are appropriate and significantly impact the prediction accuracy, which is not surprising, the function under study being quite simple.

Sine example

The sine example (D = 3) is a trickier problem since the function to be modeled is oscillating. However, we expect the BMC method to perform correctly since the true failure risk is set to r t = 2.048 × 10 -1 , a relatively high value.

Both methods give quite similar results, plotted in Figure 3. The relative error on the mean failure risk predicted by MMC is less than 5% beyond 300 samples, to be compared to 400 samples for BMC. The MMC and BMC confidence intervals size is less than 0.05 respectively beyond 600 and 800 samples.

We note that, below 200 samples, the MMC confidence interval is large, which indicates that the conditioning of the prior model does not restrict much the set of interpolating functions. The propagation of information is therefore weak and we conclude that the γ-exponential correlation function does not really suit. In addition, a spurious pinch of the confidence interval shows up for 30 samples, due to strong variations of the correlation function parameters. The maximization of the marginal likelihood yields an overly confident estimate of the correlation function parameters posterior that would be obtained by classical Bayesian analysis. Thus, MLE reaches its limits when information is scarce, firstly because it becomes very sensitive to the training set and secondly because the Bayesian posterior can be hardly approximated by its mode.

Bell-shaped example

Consider the bell-shaped example (D = 5). The benchmark function, as a mixture of 10 multidimensional Gaussian functions, is hilly. The failure set is disconnected: it is the union of 10 small hyperspheres of different radii. The real failure risk is r t = 4.341 × 10 -3 . Below 700 draws, the statistical mean decreases as O n -1 , which is typical of an undersampling of the failure set. Indeed, the first defective part has been observed beyond 800 draws. At the end of the sampling process (M = 1000), the relative error on the mean failure risk is still 23% for MMC and 31% for BMC. BMC is here obviously overconfident: the upper boundary of the 90% confidence interval comes close to touching the risk true value. MMC prediction is much more robust: the confidence interval safely flanks r t and the ratio of the mean value to the confidence interval size unambiguously indicates that additional data are required.

The raw analysis of the prediction error of y gives no evidence to conclude that information is lacking. Indeed, the average relative error, introduced by the approximation E (Y x ) of y, is about 1.7%, which is not exactly the signature of a poor model. In contrast with global sensitivity analysis, discussed in [START_REF] Oackley | Probabilistic sensitivity analysis of complex models: a bayesian approach[END_REF], goodness-of-fit requirements do not only depend upon the visiting probability distribution P : good accuracy is also required in this case at the failure set boundaries. Thus, the decision regarding the model adequacy should not rest on a global criterion. In addition, focusing on the model performances to validate the failure risk prediction may lead to over-quality: the accident set A should also be taken into account. The probabilistic risk assessment scheme we propose naturally overcome these difficulties.

Conclusion

We introduce in this article a data sparing scheme in order to measure the robustness of a design to fabrication fluctuations. It was evolved in order to offset the lack of information, either because of the complexity of the data acquisition process or due to the dimensionality of the state space. The method is therefore particularly adapted to virtual experimentation, in which physically realistic computer experiments replace test runs to foresee product performances. Indeed, simulations duration are often prohibitive; the numerical simulator is considered as a "black-box" function and approximated by an analytical emulator. In this work, the regression is based on a Student process derived from low-informative priors (uniform) on the location and scale parameters of a Gaussian process.

Risk evaluation can be viewed as a binary classification problem, the product specifications defining the class membership. Thus, we use the behavioral model as a probabilistic classifier, by interpreting predictive probabilities as degrees of belief. In this scope, the failure risk probability naturally comes up as a random variable, randomness modeling the attitude of the decision-maker while exposed to uncertainty. As a result, the distribution of the failure risk probability is a straightforward measure of the impact of the manufacturing process variability on the product performances, in a way that reflects the various sources of uncertainty: incomplete knowledge of the "black-box" function, values of the model variables, numerical approximations... The relevance of the proposed methodology has been demonstrated on theoretical examples. This study shows that, for risk assessment, model-based Monte Carlo provides a more reliable estimator of the failure rate than brute-force Monte Carlo. Although it is not a crucial point in our argument, we noticed that introducing an intermediary analytical representation (i.e. incorporating prior knowledge) may come with a higher convergence rate of the MC sampling, thus reducing the computational effort for a given accuracy.

This work opens up new perspectives for the application of random field regression to engineering risk analysis and we can already give further lines of research:

• A natural extension of this work is the generalization of the approach in order to address the multi-responses problem. A solution is to build a random field over a mixed state space including the factors as well as the responses. The main drawback here is that the space dimension may be increased drastically. In addition, it may be difficult to find an appropriate a priori correlation function. A possibly clever alternative is to model each response independently and to compute the upper bound of the failure probability of the union event.

The "worst case" risk probability is then easily derived, applying the procedure prescribed in this article.

• We can study a more sophisticated regression model of the type k j=1 U j f j (x) + V W (where f j are deterministic known functions). It is the Student extension of some classical Gaussian models extensively described in the literature (see [START_REF] Santner | The design and analysis of computer experiments[END_REF]).

• The following ternary classification problem can be investigated. If the MC samples are separated into three groups:"good", "bad" and "undecidable", it is possible to define a deficiency risk. Such random variable is a quantitative measure of the random model inability to detect defective parts. It is intended to notify explicitly the decision-maker that the model is not informed enough to draw reliable conclusions.

• In practice, we also have to control carefully numerical calculations and approximations. A particularly resistant problem is the definite-positive matrix inversion when the matrix is closed to singular.

1. Why the distribution of the random variable Z := ω → P ({x ∈ X | Y x (ω) ∈ A}) is not the object that contains the information we need ?

In order to show that the desirable information is not contained in the distribution of Z, let us give an example. Let be any number between 0 and 1. We want to describe two extreme situations (somewhat caricatural, but it is possible to approximate them by realistic situations). In these two radically different situations, the distribution of Z will be a Dirac mass at .

(i) First is the case where the random field (Y x ) gives an exact (or excellent) description of the function x → y(x), and where P (y ∈ A) = . In this case, we have a perfect (excellent) model and, for all x ∈ X, P(Y x ∈ A) = 0 or 1.

(ii) Second is an example of bad model. Consider the case where the probability P on the set X is continuous, the random variables Y x are pairwise independent and satisfy P (Y x ∈ A) = for all x. Under these conditions, we have P (Y . (ω) ∈ A) = for almost all ω. (In Appendix A.1, we give the proof of this claim.)

2. What is the information given by the random variable called failure risk probability and denoted by R ?

A necessary condition for a good model is the fact that the distribution of R is concentrated. Coming back to our example above, we see that in case (i) the distribution of R is a Dirac mass at , whereas in case (ii) the distribution of R is (1 -)δ 0 + δ 1 . (We denote here by δ a the Dirac mass at the point a).

B Examples

Discussions in the Section 3.5 are based on three arbitrary chosen benchmark functions that we describe briefly here and more precisely in supplementary material (hal-00914192) available on line. In each of these examples, the out of specification space is of the form A := [m, +∞) for some m ∈ R.

B.1 Quadric example

Fix m > 0, a 1 , a 2 , . . . , a D > m 2 , X := [-1, 1] D and for x := (x 1 , x 2 , . . . , x D ) ∈ X:

y (x) := D i=1 a i x 2 i .

B.2 Sine example

B.3 Bell-shaped example

Let R be a positive integer, X = [0, 1] D and for x ∈ X :

y (x) := max 1≤i≤R f i (x, µ i , σ i )
where for all 1 ≤ i ≤ R, µ i ∈ R D , σ i ∈ R, σ i > 0 and

f i (x, µ i , σ i ) = 1 √ 2πσ i D exp - 1 2σ 2 i x -µ i 2
such that for the fixed threshold m, the failure areas associated to each of the f i do not overlap each other.

C Correlation Function

The correlation function chosen in the Section 3.5 belongs to γ-exponential family:

ρ (x, x ) := exp   - D i=1 x i -x i l i 2 γ 2   for x, x ∈ X.
It is stationary and anisotropic. Roughly speaking, the exponent γ ∈ (0, 2] controls the smoothness of the random process (it is mean square differentiable only when γ = 2) and the positive real numbers l 1 , l 2 , . . . , l D are characteristic length-scales, defining the influence hyper-ellipsoid of an observation point.

  Let X = [-5, 5]. Let a sample of 5 points (x 1 = -4, x 2 = -3, x 3 = -1, x 4 = 0, x 5 = 2) with the associated response values (y 1 = -2, y 2 = 0, y 3 = 1, y 4 = 2, y 5 = -1). We build for an arbitrary fixed correlation function, ρ (x, x ) := exp -100 |x -x | 2 , a Gaussian field and a Student field given the knowledge of the previous sample.

Figure 1 :

 1 Figure 1: Univariate illustration. Gaussian and Student field posteriors derived from five evaluations of the unknown function y (x).

  (dt) for all m) if and only if the probability η is uniform on [0, 1].

Figure 2 :

 2 Figure 2: Quadric example: convergence plots for MMC method (left) and BMC method (right)

Figure 3 :

 3 Figure 3: Sine example: convergence plots for MMC method (left) and BMC method (right)

Fix a 1

 1 , a 2 , . . . , a D non zero integers, X := [0, 1] D and for x := (x 1 , x 2 , . . . , x D )

See Appendix A.2 for a detailed description and discussion of this toy example.

Acknowledgments

This research has been developed through a PhD thesis, supported by STMicroelectronics and the Association Nationale de la Recherche et de la Technologie (http://www.anrt.asso.fr/) via a CIFRE funding.

The authors thank the anonymous referee of this article for stimulating questions which enable them to clarify some points of the presentation.

A Lemma

A.1 Lemma Lemma. If the probability P on the set X is continuous, the random variables Y x are pairwise independent and satisfy P (Y x ∈ A) = for all x, then P (Y . (ω) ∈ A) = for almost all ω.

Proof of the lemma.

From the fact that P is continuous, we deduce that for P ⊗ P almost all (x, x ), we have x = x , hence Y x and Y x are independent random variables. We obtain

which gives

As a consequence the random variable P (Y . ∈ A) is deterministic, its distribution is a Dirac mass. Since we have

we conclude that the distribution of the random variable P (Y . ∈ A) is the Dirac mass at .

A.2 A persuasive test case

The toy example, introduced in the section 3.1, is now detailed and discussed in order to answer these two frequently asked questions :


Pierre-Eric Allier

Ludovic Chamoin

Pierre Ladevèze

Pierre Ladevèze Pgd

PGD driven by the Constitutive Relation Error -Minimal CRE/PGD

teaching and research institutions in France or abroad, or from public or private research centers.

○ simulate to design and decide ○ resolution of PDEs with • a huge number of degrees of freedom (dofs)

• a huge number of parameters • even nonlinearities … ○ direct approaches impossible due to large computation times and large data sets ! ❶ -Problematic ▪ Principle ○ capture the main features of the behavior ○ separation of variables: approximation of the solution on a reduced basis ○ without prior knowledge of this reduced basis ○ alternative methods : POD / RB

• drawback: require partial knowledge of the reduced basis ❷ -PGD: model reduction

PDE: L (u(M, t)) = 0 u(M, t) ' m X i=1 i (M) i (t)
▪ A new PGD strategy with ○ improved PGD solution ○ reliable error reduction estimation ○ "m-1, #m-1 known, computation of a m th mode (i.e. $m, %m, &m)

▪ Fixed point algorithm ○ computes alternately %m and ($m, &m)

❺ -Feasibility study

⇢ c @u @t @ @x = (x vt) u(0, t) = u(L, t) = u(x, 0) = 0 vt & (M, t) ' m X i=1 @ i (t) @t Z i (M) Z ⌦ Z i @w @x + c i w dx = 0, 8w 8i

▪

  Mix of two tools ○ the PGD model reduction ○ the Constitutive Relation Error (CRE) estimation ▪ Idea: minimization of the error indicator ❹ -Cachan, CNRS, Université Paris Saclay -61 avenue du Président Wilson, F-94235 Cachan, France ▪ Numerical models for engineering

▪

  Outlook: extension to nonlinear problems ❼ -Conclusion ▪ Different PGD strategies ○ [Ladevèze 99, Ammar 10, Nouy 10] ▪ Versus a reference reduced solution ○ SVD decomposition of a FEM solution ▪ Measure of the deviation ○ highly refined space-time meshes: discretization error ignored ❸ -Steering the progressive algorithm with an immediate and accurate error estimation, in order to improve the PGD solution Towards a more robust and stable estimator PGD computational methods on a benchmark problem: introducing a new strategy based on CRE. AMSES (2015)