Mingxiang Gao
email: mingxiang.gao@epfl.ch

Denys Nikolayev

Marko Bosiljevac
email: marko.bosiljevac@fer.hr

Zvonimir Šipuš
email: zvonimir.sipus@fer.hr

Anja K Skrivervik
email: anja.skrivervik@epfl.ch

Rules of Thumb to Assess Losses of Implanted Antennas

Keywords: implanted antennas, radiation efficiency, losses, fundamental limits. I

Fundamental limits for implanted antennas should be able to evaluate the bounds on losses with the knowledge of physical limitations from antenna and tissue parameters. In this work, rules of thumb to assess losses of deeply implanted antennas are presented toward improving the radiation efficiency. By means of rigorous and approximate approaches, the radiation efficiency of an elementary antenna within a spherical phantom is calculated as a function of multiple variables. According to the rules summarized subsequently, losses from different contributions may play a key role in specific situations, such as changes in the antenna type, phantom size, implant location and operating frequency.

INTRODUCTION

Since the development of implantable pacemakers in the 1950s, implantable electronics are becoming more acceptable and available in recent decades and offer breakthrough capabilities for a wide variety of medical diagnostics, treatments and therapies [START_REF] Fitzpatrick | Implantable Electronic Medical Devices[END_REF]- [START_REF] Dinis | Extending the limits of wireless power transfer to miniaturized implantable electronic devices[END_REF]. Today, more cutting-edge technologies allow the realization of implantable electronics with higher complexity, further miniaturization, etc. In order to achieve long-term operation of implants, wireless communication and powering through the lossy phantom have become an indispensable and challenging task, where the implanted antenna has become a crucial implantable electronic component [START_REF] Skrivervik | Implantable antennas: The challenge of efficiency[END_REF]. Compared to increasing the transmitting power for communication or powering, a comprehensive understanding of implanted antennas is generally more economical, since the improvement of the implanted antenna in terms of radiation efficiency will do more with less power under the guidance of appropriate design rules.

In free space, physical bounds on radiation efficiency of arbitrary antennas have been extensively studied [START_REF] Wheeler | Fundamental limitations of small antennas[END_REF]- [START_REF] Gustafsson | Physical limitations on antennas of arbitrary shape[END_REF]. However, for implanted antennas, the main source of radiation losses is in general due to the dissipation and scattering in the lossy biological tissues and the level of miniaturization required. To analyze the basic radiation properties of implanted antennas, a simplified body model can be a good choice [START_REF] Kim | Implanted antennas inside a human body: Simulations, designs, and characterizations[END_REF]- [START_REF] Bosiljevac | Propagation in finite lossy media: An application to WBAN[END_REF]. Although only approximate results can be obtained, they can provide more physical insights for implanted antennas, which will be beneficial to put forward the design rules. For instance, based on the spherical dyadic Green's function expansions and finitedifference time-domain code, Kim and Rahmat-Samii [START_REF] Kim | Implanted antennas inside a human body: Simulations, designs, and characterizations[END_REF] analyzed the electromagnetic characteristics of dipole antennas and low-profile patch antennas implanted in the human head and body. Poon et al. [START_REF] Poon | Optimal frequency for wireless power transmission into dispersive tissue[END_REF] showed that the optimal frequency can be achieved at GHz-range for an inbody device operating in the midfield region. Using spherical wave expansion method, Merli et al. [START_REF] Merli | The effect of insulating layers on the performance of implanted antennas[END_REF] studied the effect of implantable device encapsulation on radiation efficiency, which improves with increasing permittivity and thickness of encapsulation.

This paper aims to investigate the radiation efficiency of the deeply implanted antennas using canonical models, as was first proposed in [START_REF] Bosiljevac | Propagation in finite lossy media: An application to WBAN[END_REF]. Following this approach and using a spherical wave expansion (SWE) to estimate the radiation efficiency of implanted canonical antennas (elementary sources) [START_REF] Bosiljevac | Propagation in finite lossy media: An application to WBAN[END_REF]- [START_REF] Skrivervik | Fundamental limits for implanted antennas: maximum power density reaching free space[END_REF], this paper presents the evolution of the radiation efficiency as a function of multiple variables: the frequency, the radius of the phantom and source, and the depth of the implant. Based on the understanding of the loss mechanisms in these canonical cases, rules of thumb to assess the losses of deeply implanted antennas are summarized qualitatively for different situations, e.g. changes in the antenna type, phantom size, implant location and operating frequency.

II. MODEL AND ANALYSIS METHODS

To analyze the losses of implanted antennas within the human body, we follow here the method developed in [START_REF] Skrivervik | Fundamental limits for implanted antennas: maximum power density reaching free space[END_REF]: The antenna, an elementary electric or magnetic dipole, is implanted in a spherical homogeneous or stratified body phantom. The model is depicted in Fig. 1 for the case of a homogeneous phantom. The source is surrounded by a small lossless sphere, usually air, of radius r impl (roughly represents the antenna size). It is located at a distance r feed from the center of the spherical phantom, the latter having a radius r body .

A. Spherical Wave Expansion (SWE)

Regarding the spherical body phantom, our first solution procedure makes use of the SWE proposed in [START_REF] Bosiljevac | Propagation in finite lossy media: An application to WBAN[END_REF]. Depending on the design of the model structure, the source can be placed at any point inside of the spherical phantom, and the phantom can be made stratified (representing additional fat/skin layers, for instance) with the help of a mode matching technique presented in [START_REF] Fuchs | Scattering of spherically and hemispherically stratified lenses fed by any real source[END_REF]. In the SWE method presented in [START_REF] Bosiljevac | Propagation in finite lossy media: An application to WBAN[END_REF], for the case where an off-centered source is implanted in a spherical phantom, the structure of interest consists of two spherical structures, i.e the body model and the antenna model, which are connected using addition theorems. In this paper, a simplified method is adopted: the source is directly placed in the spherical phantom without the surrounding lossless sphere, and the calculation results are normalized by the corresponding value at the outermost point of the lossless sphere boundary. Since the radius of the antenna (i.e. the radius of the surrounding lossless sphere r impl ) is always relatively small, the field distribution at the its boundary will contain only the fundamental spherical mode. Therefore, according to the equivalence principle, such a simplified method is reasonable for most implanted antennas with small radius (i.e. r impl ≪ λ).

B. Approximations of the Radiation Efficiency for Implanted Antennas

When the source is placed at the center of a homogeneous spherical phantom, approximations of the radiation efficiency for implanted antennas were proposed in [START_REF] Skrivervik | Fundamental limits for implanted antennas: maximum power density reaching free space[END_REF]. Here only the final expressions are repeated for the case with an electric dipole as excitation. The total radiated power reaching free-space could be expressed as total radiated power entering the body near field propagating field reflections

P P e e e     . (1)
In this way, we calculate the radiation efficiency as the ratio of P total radiated power to P entering the body .

The specific expressions for losses due to these different contributions are expressed as

  2 near field * 2 -1 2 1 3 impl impl impl Re Re 2 1 j k e k k r r k r k                                  , (2)
2Z T Z Z   , (4)
where is the wave impedance in the considered medium, k denotes the complex wave number in this medium, and α is the real part of the propagation constant. Z air and Z body are the mode impedances of the spherical wave in the corresponding medium.

Compared with SWE, this approximation helps to understand and separate the contribution of losses in different situations, which will guide the design of the antenna to work on the primary one in order to improve the radiation efficiency.

III. RESULTS

Our first example is similar to the implant described in [START_REF] Nikolayev | Optimal radiation of body-implanted capsules[END_REF]. However, we use here a Hertzian dipole as an antenna. We consider an electric dipole surrounded by a bubble having a radius r impl = 5 mm implanted at the center point of a spherical body phantom with the radius r body = 10 cm. The phantom is homogeneously composed of muscle tissue, and its permittivity is provided by the four-region Cole-Cole model defined in [START_REF] Gabriel | The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues[END_REF]. Fig. 2 shows the radiation efficiency as a function of the frequency (from 100 MHz to 4 GHz) computed using SWE and the approximations of the radiation efficiency for implanted antennas, respectively. The results using two different methods are almost superimposed, and these results correspond well to the ones published in [START_REF] Nikolayev | Optimal radiation of body-implanted capsules[END_REF] (see Fig. 2(d) in [START_REF] Nikolayev | Optimal radiation of body-implanted capsules[END_REF]), although the sources in the latter work are shaped as cylindrical capsule antennas with similar dimensions. Using the same methods, we investigated the radiation efficiency for the same canonical case as a function of the frequency and the radius of the phantom r body . The source is still placed at the center point of a spherical phantom made of muscle, and r impl is kept at 5 mm. Electric and magnetic dipoles are considered. The results are shown in Fig. 3 for the electric dipole and in Fig. 4 for the magnetic dipole. The results computed via the approximations of the radiation efficiency for implanted antennas [Fig. 3 (a) and Fig. 4(a)] are very similar to the results by the SWE [Fig. 3 (b) and Fig. 4 (b)]. Deviations can be found when the radius of the phantom rbody is small. This happens when near-field losses and the losses due to reflections at the phantom-air interface are predominant.

According to Fig. 3, the following insight can be gained:  The frequency and r body correlatedly affect the value of the radiation efficiency. This is mainly because the frequency and the properties of the phantom (i.e. its size and composition) jointly determine the field region of the source where the phantom-air interface is located. Another important factor is that the permittivity of the lossy tissues, like muscle, is always a function of the frequency (especially, the imaginary part of permittivity that represents losses due to both displacement and ionic currents).  For low frequencies (e.g. ≲ 500 MHz in the case of electric dipole), losses due to the reactive near-field absorption will significantly decrease radiation efficiency; when the frequency increases, losses due to propagating field absorption begin to play a major role in reducing the radiation efficiency.  When the radius of the phantom r body is small, the contribution of reflections losses can't be ignored; as r body increases, the losses due to propagating field absorption become more dominant.  For the case of electric dipole implanted at the center point of the spherical body phantom with r body of several centimeters, the frequency that optimizes the radiation efficiency is in the GHz-range or sub-GHz range. This being said, we see that, as r body increases, the satisfactory values of this frequency start to decrease, but its descending speed will gradually slow down, as in this case when r body is larger than 5 cm. Similarly, the above analyses are mostly applicable to the case excited by a magnetic dipole, as shown in Fig. 4 (a) and (b). Some more insight can be obtained:

 Compared with the results of an electric dipole, the effect of the near-field losses produced by a magnetic dipole is significantly smaller. Therefore, we see that the optimal frequency regarding the radiation efficiency becomes lower than in the case of the electric dipole case, especially as r body increases to more than 5 cm.

The radius of the implanted antenna r impl is another key parameter for implanted antennas. Fig. 5 (a) and (b) show the distributions of the radiation efficiency as a function of the frequency and r impl for an electric and a magnetic dipole placed at the center point, respectively, in which r body is fixed to 5 cm. Comparing these results, we could find that:

 Whether for the implanted electric dipole or magnetic dipole, the increase in r impl reduces the losses especially the reactive near-field losses. Since the near-field losses for the electric dipole are relatively large, the frequency that optimizes the radiation efficiency start to increase when r impl becomes small enough [i.e. r impl < 5 mm in Fig. 5 (a)]. However, for the magnetic dipole, the acceptable range of the optimal frequency does not change much with the decrease of r impl . Finally, by the SWE, a more general case has also been investigated: an electric dipole deeply implanted in a spherical phantom made of muscle (r body =10 cm) with the offset distance r feed increasing from 0 to 5 cm. We notice that:

 When r feed increases to a relatively small length (i.e. smaller than 4 cm in this case), the frequency that optimizes the radiation efficiency is basically maintained in a fixed range under the tolerance of -1 dB, as shown in Fig. 6. This means that the rules and conclusions obtained before are also applicable to most deeply implanted antennas, despite they may not be placed at the center point of the spherical phantom.

IV. CONCLUSION

In this paper, we have investigated the radiation efficiency of an antenna deeply implanted in a lossy phantom using canonical models. According to the results for different situations, rules of thumb to assess losses of implanted antennas have been presented. We summarized the effects of frequency and phantom radius on the losses due to three contributions, i.e. reactive near-field absorption, propagating field absorption and reflections. We then analyzed the optimal frequency (the frequency yielding the highest overall radiation efficiency) and the according radiation efficiency as a function of the type of the source (electric or magnetic), the size of the phantom, the depth of the implant and the size of the implant. The established guidelines facilitate the choice of the antenna type and dimensions and provide simple rules to check the feasibility of a given antenna specification.

Fig. 1 .

 1 Fig. 1. View of the homogeneous spherical body phantom with the offcentered excitation.

Fig. 2 .

 2 Fig. 2. Spectrum of the radiation efficiency excited by an electric dipole placed at the center point of the spherical body phantom.

Fig. 3 .

 3 Fig. 3. Distributions of the radiation efficiency for an electric dipole as a function of the frequency and the radius of the phantom r body via (a) the SWE and (b) the approximations of the radiation efficiency for implanted antennas.

Fig. 4 .

 4 Fig. 4. Distributions of the radiation efficiency for a magnetic dipole as a function of the frequency and the radius of the phantom r body via (a) the SWE and (b) the approximations of the radiation efficiency for implanted antennas.

Fig. 5 .

 5 Fig. 5. Distributions of the radiation efficiency for (a) an electric dipole and (b) a magnetic dipole as a function of the frequency and the radius of the implanted antenna r impl , computed using the approximation of the radiation efficiency for implanted antennas.

Fig. 6 .

 6 Fig. 6. Distributions of the radiation efficiency for an electric dipole as a function of the frequency and the offset radius of the source.

ACKNOWLEDGMENT This work was supported in part by the China Scholarship

Council.


Jasper Goethals
email: jasper.goethals@ugent.be

Günter Vermeeren

Denys Nikolayev

Margot Deruyck

Luc Martens

Wout Joseph

Combined Antenna-Channel Characterization for Wireless Communication from Horse Hoof to Base Station

Keywords: channel characterization, modeling, path loss, link budget, horse, LoRA

This paper presents a design of an antenna in the complex environment of the horse hoof for sub-gigahertz (868 MHz) communication. The influence of the leg and the ground on the performance were examined by means of finitedifference time-domain simulations. Furthermore, an adaptation was presented to increase the efficiency of the antenna. The first device has a center frequency of 858 MHz and a bandwidth of 99.5 MHz. The total efficiency is 3.12%. Adding the ground to the model results in the ground absorbing most of the radiated power which leads to a total efficiency drop to 1.6%. The center frequency stays the same (at 862 MHz). The bandwidth decreases to 70 MHz. When the horseshoe is connected to the device, the antenna is not tuned anymore, yet the the total efficiency stays 3.08%. This shows that connecting the device to the horseshoe leads to better radiation efficiency. When using LoRa technology, this setup can reach 1631 m if the hoof is in the air. When the hoof is on the ground, only a range of 115 m is estimated.

I. INTRODUCTION

Owning a horse is a serious undertaking. The health of this valuable asset is of utmost importance. To monitor its health, one can put a mobile monitoring device on the animal. There already exist products with this function [START_REF]Smart Halter | Alerts Equine Distress, Foaling, Colic Signs of Horse[END_REF], [START_REF]EquiFone Foaling Alarm[END_REF], but the performance is not up to the expected standards. The overarching goal of this work is a part of a horse monitoring system. The device should monitor the horse 24/7 and automatically detect anomalies in its behavior while the horse is in its usual environment which include the pasture. This automatic detection of anomalies decreases the response time for a veterinarian intervention which leads to a reduce of the risk of complications and medical costs.

In this paper, the goal is to establish a wireless communication channel between the transmitting monitor device inside the hoof of the horse and off-body hoof antenna to a gateway, e.g. installed in the farm. To cover sufficiently long ranges during training and recreational activities and to be able to directly communicate from the horse to a gateway, radio technologies with suitable propagation characteristics need to be used. Lower frequency waves typically penetrate objects better at lower energy consumption costs, thereby improving the achievable range. However, lower radio frequencies typically require larger antenna sizes, and limit the maximum achievable throughput, thereby limiting device placement options and the amount of monitoring data, respectively [START_REF] Balanis | Antenna theory: analysis and design[END_REF]. Trade-offs between range, antenna size, and throughput result in using radio frequencies from the 433 MHz ISM (Industrial, Scientific and Medical) bands or the EU 868 MHz SRD (Short Range Device) bands (jointly referred to as sub-GHz bands).

For deciding upon optimal sub-GHz technology and antenna placement, an accurate propagation analysis is required, specifically related to the wireless body area networks (WBAN) propagation characteristics. The physical channel for off-body communication (from the sensor in the hoof to the private/cellular operated gateway) is required. Most research about off-body propagation is focused on humans [START_REF] Ambroziak | Impact of radio wave polarisation on off-body communications in indoor environments[END_REF], [START_REF] Van Roy | Dynamic Channel Modeling for Multi-Sensor Body Area Networks[END_REF]. However, [START_REF] Gabriel | The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues[END_REF], [START_REF] Kwong | Wireless Sensor Networks in Agriculture: Cattle Monitoring for Farming Industries[END_REF] concluded that radio performance will be significantly affected by an animal's body. Kwong et al. [START_REF] Kwong | Practical considerations for wireless sensor networks in cattle monitoring applications[END_REF] proposed a simple two-path model for off-body communication, but the model has not been validated through measurements. Benaissa et al. [START_REF] Benaissa | Experimental characterisation of the off-body wireless channel at 2.4 GHz for dairy cows in barns and pastures[END_REF], [START_REF] Benaissa | Internet of animals: characterisation of LoRa sub-GHz off-body wireless channel in dairy barns[END_REF] investigated off-body path loss and communication for dairy cows in stables at 2.4 GHz and 868 MHz. For on-body communications, the work of Reusens [START_REF] Reusens | Characterization of On-Body Communication Channel and Energy Efficient Topology Design for Wireless Body Area Networks[END_REF] and Roy [START_REF] Van Roy | Dynamic Channel Modeling for Multi-Sensor Body Area Networks[END_REF] are relevant for UWB (Ultra Wide Band) and 2.4 GHz frequencies.

With respect to the impact of metal in the channel (horsehoe near the antenna), Tanghe et al. [START_REF] Tanghe | Intra-, Inter-, and Extra-Container Path Loss for Shipping Container Monitoring Systems[END_REF] investigated the influence of metallic environments on the performance of antennas. Besides the metal environment, also the animal's body will affect the antenna design and its optimization [START_REF] Nikolayev | Optimal radiation of body-implanted capsules[END_REF] [START_REF] Nikolayev | Impact of tissue electromagnetic properties on radiation performance of in-body antennas[END_REF]. To obtain accurate link budgets, one needs to model the antennas and the channel jointly, which is not available for the considered use case. For this paper, we aim to use 868 MHz. This opens the possibility to use radio technologies like LoRaWAN, Sixfox and other low-power IoT technologies. For this, the channel of the configuration will be characterized and a link budget will be calculated based on simulations.

To the best of the authors' knowledge, this is the first time off-body horse hoof channel is been studied. The novelty of this paper is the antenna-channel design in the unique environment of animal tissue of the horse leg combined with the metallic horseshoe surrounding the antenna. Also the influence of the presence of the ground will be investigated on the propagation and the antenna. In this section, the general approach will be discussed. First an antenna operating in the hoof pad is designed while enclosed with the horse leg and the horseshoe. The parameters of interest are the return loss, bandwidth (BW), center frequency and total efficiency of the antenna. The BW is here defined as the range of frequencies which have a |S 11 | below -10 dB. In this step, the influence of the electrical properties of the horse leg on the antenna will be studied.

A. Antenna Design

The antenna design has to satisfy the spacial constrains of the device. In this work, the assumption is that the device will be installed in a large hoof pad with dimensions ca. 150 × 150 × 4 mm 3 . A border of 2 cm has to be kept clear. Here the iron will press the pad against the hoof.

1) Antenna: A planar inverted F-antenna is chosen as the antenna for the device. Fig. 1 shows the antenna lay-out. The target center frequency is 868 MHz (BW should be wide enough for LoRa, 250 kHz). The design is based on the works of Lizzi et al. [START_REF] Lizzi | Miniature Multiband Inverted-F Antenna over an Electrically Small Ground Plane for Compact IoT Terminals[END_REF]. The structure is easy to adapt. The placement of the feed can be used to match the antenna and the length of the arm determines the center frequency. With a well understood antenna, the examination of the channel can be performed. Another advantage of this antenna are the volumetric properties. They fulfill the current constrains of the device. The dimensions of the antenna are 43 × 25 × 1 mm 3 . The antenna itself is 15 mm long. The dielectric is FR-4 (ε r = 4, σ = 0.0023 S/m). The metal tracks and the ground are perfect electrical conductors (PEC). The antenna is embedded in the pad (ε r = 2.5, σ = 0 S/m) (Fig. 2a). This encapsulation will increase the robustness of the parameters of the antenna to the influence of the horse leg and the ground.

2) The configuration: Fig. 2b shows the complete configuration. The 20 cm long leg and the 7 mm thick horseshoe sandwiches the hoof pad wherein the device lays. The leg made out of bone and muscle (homogeneous both times) have been considered with parameters from human tissue [START_REF] Gabriel | The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues[END_REF]. Another simulation will be done with the horse leg out of muscle. The horseshoe is made out of PEC. This is the main configuration that is examined.

B. Connecting the Iron Shoe

In this configuration, the ground of the device is connected to the horseshoe with a wire through the nails. Fig. 2d shows the configuration. The effect of the addition of this connection will be examined.

C. Influence of the Soil

An interesting situation to examine is the performance of the antenna when the horse leg rests on the ground (Fig. 2c. Here, the underground is taken as dry sand (ε r = 2, σ = 10 S/m) with the dimensions of 500×500×50 mm 3 . The type of soil can very a lot, going from air to wet muddy dirt. When approaching pure water, an effective relative permittivity of 90 can be reached [START_REF] Risman | Retro-modelling of a dual resonant applicator and accurate dielectric properties of liquid water from -20 • C to +100 • C[END_REF], which could substantially influence the properties of the device. These extreme cases will not be examined here.

D. Path Loss Model and Link Budget

1) Path Loss model: To calculate the loss of the power as a function of the distance, a path loss model is proposed [10] [17]. The received power P r (d) after a distance d is given by

P r (d) = P t G t G r H 2 t H 2 r d 4 L (1)
with P t the power of the transmitter (Tx), G t/r the gain of the receiver (Rx)/Tx, H t/r the height w.r.t. the ground of Rx/Tx, d the distance between Rx and Tx and L the losses not accounted for in this path loss model. In this paper, these losses are shadowing losses (M s ) and multi-path fading (M f ).

In conclusion, the maximum range (in dB) is given by with S the sensitivity of the considered radio frequency (RF) technology.

d max = S -P t -G r -G t -M f -M s + 2H t + 2H r 4 (2)
2) RF Technologies: Here, the LoRaWAN technology [18] will be examined. The sensitivity of the chip is inversely correlated with the data throughput. LoRa link budget values are listed in Table I for a bandwidth of 125 MHz.

The spreading factor SF corresponds with the length in time of a symbol [START_REF] Lorawan | The Things Network[END_REF]. The sensitivity S is than calculated as

S = -174 + 10log 10 (BW ) + N F + SN R (3)
with BW the bandwidth in kHz, NF the noise factor which is dependent on the hardware, and SNR the signal to noise ratio.

The link budget also depends on the power of the transmitter. In Belgium, the maximum allowed power is 25 mW or 14 dBm in the 868 MHz ISM band [START_REF] Lorawan | The Things Network[END_REF].

3) Link budget: The link budget will be calculated for four configurations: once for the horse leg in the air, once for the leg on the soil, once with the connected device in the air, and once with the connected device on the ground. The range in both cases will be calculated for the maximum gain and once for a interesting direction (Fig. 3a). As interesting direction, a slice of the gain at 5 • from the plane earth is chosen. A trigonometry calculation shows that the angle will be around this value for longer distances (>100 m) from the base station. From the slice, the maximum gain will be taken. The orientation of the antenna w.r.t. the earth will be different for hoof on ground case and the hoof in air case (shown in Fig. 3b), as a horse turns its hoof 90 • when it lifts it in the air.

E. FDTD simulations

The calculations of the parameters of interest will be done by Sim4Life, a finite-difference time-domain (FDTD) simulator [START_REF]Sim4Life[END_REF]. For the calculation of the BW and center frequency, a Gaussian signal was excited over the edge source at the feed strip. The center of the Gaussian was 800 MHz and the BW was 250 MHz. To simulated an infinite flat earth, the soil was brought into the perfectly match layer for the hoof on the ground configuration. The simulation were then performed with a harmonic excitation of 868 MHz to calculate the gain pattern, efficiency and the |S 11 |. The number of cells of the grid were 440 912 kCells and 668 222 kCells for the BW simulations without or with ground, respectively. 472.004 kCells and 630.808 kCels for the gain pattern simulations without or with ground, respectively.

III. RESULTS

A. Configuration 1: The Hoof in the Air

When the leg is in the air, the antenna has a center frequency of 858 MHz and a bandwidth of 99.5 MHz (only 250 kHz is needed for LoRa), as shown in Fig. 4. At 868 MHz, |S 11 | = -19.57 dB. The maximum gain is -10.0 dBi. The gain pattern is shown in Fig. 5a. The gain in the interesting direction is -11.0 dBi. The total efficiency is 3.12%. The low efficiency is due to the absorption of the leg, hence a diminishing of the radiation efficiency, but is satisfying for WBAN's where often the efficiency is <1% are obtained. The antenna in free space without horseshoe and horse leg has a radiation efficiency of 50%.

Varying the electric properties of the leg (muscle instead of bones as above) changes the characteristics of the antenna. The center frequency with a muscle phantom changes to 784 MHz with a bandwidth of 99.8 MHz, as shown in Fig. 4. This does not satisfy the requirement of Section II anymore. This shows that a good knowledge of the horse leg is needed to design the device. The device should also be robust for different kind of horses and different horses in general. For this, a study of the electrical properties of the horse's hoof area and its variations will be conducted in the future.

B. Configuration 2: Connecting with the Iron

When the device is connected to the horseshoe via a wire and nails, the center frequency shift of 8 MHz to 886 MHz. The bandwidth decreases to 81MHz. This is shown in Fig. 4. At 868MHz, |S 11 | = -12 dB. The total efficiency, 3.08%, stays the same despite a lower matching. This shows that connecting to the horseshoe results in a slight gain in radiated efficiency. The maximum gain is -9.99 dBi. The maximum gain in the interesting direction of 5 • is -13.7 dBi. The gain pattern is shown in Fig. 5c. Connecting the shoe also changes the shape of the gain pattern. power is radiated in the upper half of the space. The soil absorbs the radiation, which explains the lower efficiency of 1.6%, but does not alter the antenna characteristics.

D. Configuration 4: Connected device on the ground

As final configuration, the device is now connected to the horseshoe and the hoof is put on the ground. The center frequency shifts 2 MHz to 884 MHz and the BW stays the same at 82 MHz. At 868 MHz, the |S 11 | = dB (Fig. 4). The total efficiency is 3.3%. The maximum gain is -6.3 dBi and -12.7 dBi at 5 • (Fig. 5g and5h). We conclude that it is interesting to connect the device with the horseshoe when the hoof is on the ground.

E. Range per configuration

We calculated the gain in sections III-A to III-D. The range of the device can be calculated with these results. As mentioned, the maximum allowed radiated power is 14 dBm. As base station antenna, a dipole is considered, this has a gain of 2.15 dBi. The maximum data rate is from Table I with a bandwidth of 125 kHz is 5469 bps. This results in a sensitivity of -120 dBm. Table II links the link budget calculations with the effective hoof channel. When the hoof is in the air, we estimate that the antenna will be 10 cm off the ground. Otherwise, it will be 7mm off the ground, the thickness of the horseshoe iron. The other path losses are taken from [START_REF] Benaissa | LoRa in-to-out-body link budget for animals[END_REF]. The configurations are hoof in the air, hoof in the air with the horseshoe connected, hoof on the ground with normal setup and hoof on the ground with connected setup respectively. The calculations show that it is better to send the data when the hoof is in the air (>1300 m when in the air compared to <200 m when on the ground). A difference of 41 dB between in the air and on the ground is due to the path loss model, this model punishes low-to-the-ground situated antennas. Yet, it is interesting to see that connecting the device with the horseshoe increases the range with more

IV. CONCLUSION

In this paper, the characterization of the horse hoof environment for sub-gigahertz communication with an inverted-F antenna is investigated using FDTD simulations. The dielectric properties of the hoof have to be understood to design a robust antenna. Connecting the horseshoe with the device can be beneficial to the performance, it can gain up to 20% in range. The best way for a efficient communication is sending data when the hoof is in the air. Future work will consist of improving models by constructing a live like horse model. Characterization of the horse hoof will aid the designer to make a robust device for every targeted horse. In later work, the antenna should be converted to dual-band to allow 2.45 GHz communication. Also the influence of different types of soil is part of future work. Lastly, manufacturing the antenna will validate the results. The preliminary experimental validation will be presented at the conference.

Fig. 1 .

 1 Fig. 1. The layout of the antenna

  Fig. 2. The total configuration and the alternative configuration

Fig. 3 .

 3 Fig. 3. Clearification link budget calculations

C. Configuration 3 :Fig. 4 .

 34 Fig. 4. S 11 of the different configurations: conf1 = horse leg without ground, conf1_m = conf1 with muscle properties, conf2 = the connected device with the horseshoe without ground, conf3 = conf1 with the ground, conf4 = conf2 with the ground, the blue vertical line notes the targeted frequency of 868 MHz

Fig. 5 .

 5 Fig. 5. Gain patterns

TABLE I DATARATE

 I TABLE LORA AT A BW OF 125 KHZ

	Spreading factor SNR limit (dB) Bitrate (bps)
	7	-7.5	5469
	8	-10	3125
	9	-12.5	1758
	10	-15	977
	11	-17.5	537
	12	-20	293

ACKNOWLEDGMENT

This work was executed within the imec.icon project Hoof-MATE, a research project bringing together academic researchers and industry partners. The Hoof-MATE project was co-financed by imec and received project support from Flanders Innovation & Entrepreneurship (project nr. HBC.2018.0536).


Yonatan-Carlos Carranza-Alarcon
email: yonatan-carlos.carranza-alarcon@hds.utc.fr

Sébastien Destercke
email: sebastien.destercke@hds.utc.fr

Imprecise Gaussian Discriminant Classification

Keywords: Discriminant Analysis, Robust Bayesian, Classification, Near-Ignorance

Gaussian discriminant analysis is a popular classification model, that in the precise case can produce unreliable predictions in case of high uncertainty (scarce or noisy data set). While imprecise probability theory offer a nice theoretical framework to solve such issues, it has not been yet applied to Gaussian discriminant analysis. This work remedies this, by proposing a new Gaussian discriminant analysis based on robust Bayesian analysis and near-ignorance priors. The model delivers cautious predictions, in form of set-valued class, in case of limited or imperfect available information. Experiments show that including an imprecise component in the Gaussian discriminant analysis produce reasonably cautious predictions, and that set-valued predictions correspond to instances for which the precise model performs poorly.

Introduction

In machine learning, the classification task consists in seeking to identify to which label (among a finite set K of such labels) a new unlabelled instance x ∈ X p belongs. The reliability of this precise prediction (or single decision) may depend heavily on prior beliefs (e.g. assumptions made by data analysts, such as asymptotically unbiased estimators) and the nature of training data sets (e.g. in small amounts [START_REF] Kitchin | Small data in the era of big data[END_REF]2] and/or with high degree of uncertainty 1 ), both will be referred as imperfect information 2 . A well-known precise generative classifier model used to perform the classification task is the Gaussian discriminant analysis (GDA) [6, §4.3] A classifier model is called precise when it performs pointwise predictions (or precise estimations) in the form of single class labels, even in extreme cases, regardless of the available information we have a about an instance. In these cases, it may be useful to provide set-valued, but more reliable predictions, especially for sensitive applications (e.g. medical diagnosis, control systems, etc.) where we cannot afford to make mistakes (see illustration in Figures 3(a) and 3(b)).

Imprecise probabilities (IP) [START_REF] Augustin | Introduction to imprecise probabilities[END_REF] can mitigate the impact of imperfect information, taking into account the lack of evidence by replacing precise estimates (or a single probability distribution) with imprecise estimates (or a set of probability distributions, most often in the form of convex set) in order to make cautious setvalued decisions. In this paper, we adopt such an approach to propose a novel method of (cautious 3 ) imprecise classification.

Cautious classification is a relatively new trend in machine learning which do not aim to do "better" than their precise counterparts, nor to implement a rejection option (i.e., not classifying at all) in case of ambiguity [START_REF] Herbei | Classification with reject option[END_REF], but to highlight those hard cases for which information is insufficient to isolate a single reliable precise prediction, and to propose a subset of possible predictions. We can find in the literature three "main" ways to access cautious classifier models: (1) using a classical precise classifier but deriving a set-valued predictions from them [START_REF] Cheng | Label ranking with partial abstention based on thresholded probabilistic models[END_REF] (e.g. partial reject [START_REF] Ha | The optimum class-selective rejection rule[END_REF], conformal prediction [START_REF] Shafer | A tutorial on conformal prediction[END_REF]) (2) making data imperfect (coarse or impartial observations) and then building a corresponding imperfect robust model, and finally (3) a cautious classifier under IP from which set-valued predictions follow naturally (such as robust frequentist inference [START_REF] Cattaneo | Statistical decisions based directly on the likelihood function[END_REF][START_REF] Cattaneo | Fuzzy probabilities based on the likelihood function[END_REF] or Bayesian inference [START_REF] Walley | Statistical reasoning with imprecise Probabilities[END_REF][START_REF] Walter | Generalized bayesian inference under prior-data conflict[END_REF][START_REF] Quaeghebeur | Imprecise probability models for inference in exponential families[END_REF]). We retain here the latter as it considers imprecision as parts of its basic axioms, rather than the other approaches where imprecision is added ex-post.

Bayesian methods incorporate some prior beliefs in the form of probability distribution defined on unknown parameters of the model. Such beliefs typically comes from expert opinions of persons that are knowledgeable in the context of the problem. However, it is also well-known that the elicitation of prior beliefs can be absent or hard to obtain during the study of a problem, especially in when learning classifier. A classical way out of this problem is to use so-called non-informative prior, that allow one to obtain a posterior not including any prior knowledge [2]. Yet, the use of such prior is not without problem within the Bayesian theory, as they are not coherent in the sense of De Finetti. In addition, it has been argued and shown that using truly vacuous prior information while remaining coherent usually lead to vacuous posterior predictions [14, 17, §7.4, §5.6.2] (i.e. our model would not be able to learning from data). Moreover, it may seem strange that an absence of prior should lead to a fully precise, completely informed posterior. Walley have therefore proposed to use a set of non-informative prior distributions, called near-ignorance prior [14, §4.6.9], to solve this issue. These near-ignorance prior must respect certain properties [18, §2] so as not to obtain vacuous predictions. Hence, one of our motivations in this paper is to not to use a single prior distribution, but a set of prior distributions (or credal set [START_REF] Levi | The enterprise of knowledge: An essay on knowledge, credal probability, and chance[END_REF]) to reflect our lack of knowledge and obtain cautious predictions. Let X × K be the space of observations and possible labels, withX ∈ X = R p a random vector and Y ∈ K = {m 1 , ..., m K } the set of categories. The main goal of GDA is to estimate the theoretical conditional probability distribution (c.p.d) P Y =m k |X of the class Y = m k given the observation X via Bayes' theorem as follows

P Y =m k |X = P X|Y =m k P Y =m k m l ∈K P X|Y =m l P Y =m l . (1)
Thus, quantifying P Y =m k |X is equivalent to quantify P X|Y =m k and the marginal distribution P Y . In precise probabilistic approaches, this is typically done by using maximum likelihood estimation (MLE) and by making some parametric assumptions about the probability density P Y =m k |X (i.e. Gaussian probability distribution (g.p.d)) in order to find a plausible estimate (see Section 3.1). However, such precise estimates usually have trouble differentiating different kinds of uncertainties [START_REF] Senge | Reliable classification: Learning classifiers that distinguish aleatoric and epistemic uncertainty[END_REF], such as uncertainty due to ambiguity (mixed classes in some areas of the input space) and uncertainty due to lack of knowledge or information (limited training data set inducing biases in estimates [START_REF] Braga-Neto | Is cross-validation valid for small-sample microarray classification?[END_REF]). In both cases, it may be useful to provide setvalued, but more reliable predictions, especially for sensitive applications where we cannot afford to make mistakes (see illustrations in Figures 3(a) and 3(b) ). Section 2 recalls the basics of the precise classification setting, adopting the viewpoints of statistical decision theory [6, §2] and expected utility [20, §2.2]. It also introduces the corresponding extensions of these tools to the IP setting, in particular the maximality criterion [START_REF] Troffaes | Decision making under uncertainty using imprecise probabilities[END_REF] that extend the precise utility-based decision-making to the IP context.

In Section 3, we describe the estimation of the conditional distribution in the case of the precise GDA, using a frequentist inference approach. We then extend this precise parametric estimation to imprecise estimation in a robust Bayesian inference context, using the IP near-ignorance model proposed by Benavoli et al. [START_REF] Benavoli | Prior near ignorance for inferences in the k-parameter exponential family[END_REF] to do so and obtaining estimates in the form of a convex set . Coupling this imprecise estimation with the maximality criterion, we present our Imprecise GDA (IGDA) model and its different variants in Section 4.

In Section 5, we perform a set of experiments on different datasets using our imprecise model and compare it to its precise counterparts. We show that the cautious predictions are useful, in the sense that they concern instances for which the precise classifier often makes mistakes, often include the true class within the predicted set or these same instances, and are not overly imprecise. Furthermore, we briefly discuss (focusing on computational issues) in Section 6 the extension of our method to other settings, namely to the case where the class proportions P Y =m k are also imprecisely estimated, and where the criteria to minimise is not the raw number of errors (corresponding to a 0/1 loss function) but a generic loss function. In this paper, we will use mathematical notations of table 1.

Symbols

Description Data related notations

X p ⊆ R p Input space of dimension p. K = { m 1 , ..., m K } Output space of size K. (•) T Transpose Operator • Euclidian norm X Matrix n × p of all instances of dataset. y = (y 1 , . . . , y n ) T
Vector n × 1 of all label of dataset.

P (Y = m k |X = x) ∼ P Y |x Conditional probability of m k given x Φ = {ϕ |ϕ : X × Θ → K}
Learning space models

L 0/1 (y, ϕ(x)) = 1 y =ϕ(x) zero-one loss function Discriminant Analysis n k Number of observations of category m k (x i,k , y i,k ) n k i=1 = {(x 1,k , y 1,k ), . . . , (x n k ,k , y n k ,k )} Observations of category m k . x k = 1 n k n k i=1 x i,k Empirical mean of category m k σ j m k = 1 n k -1 n k i=1 (x j i,k -x j k ) 2 , ∀j ∈ {1, . . . , p} Empirical variance of category m k S m k = 1 N -n k n k i=1 (x i,k -x k )(x i,k -x k ) T Empirical covariance matrice of category m k S = 1 (N -K) K k=1 n k i=1 (x i,k -x k )(x i,k -x k ) T Empirical total covariance matrix π y = { π y=m k | π y=m k = n k /N, m k ∈K π y=m k = 1} Empirical marginal distribution P Y .

Preliminaries and basic reminders

In this section, we remind some notions of classical statistical learning and decision-making used to build a precise classification model, as well as basic notions needed to also deal with sets of probabilities.

Classification setting

Let D = {(x i , y i )|i = 1, . . . , N } be a training data set issued from X p × K, such that x i ∈ X p are 75 regressors or features (input space) and y ∈ K is the response variable or class (output space). We denote n k the number of observations that belong to the label m k , and so N = K k=1 n k . The goal of classification is to build a predictive model ϕ : X → K that predicts a label m k ∈ K given a new unlabelled instance (x, •) / ∈ D. A well-known approach for that is called inductive learning which involves two steps; (1) the learning phase that estimates or induce4 a model from observed data, and (2) the decision inferred from the learned model. Learning consists in determining the optimal model ϕ ∈ Φ among the chosen set Φ of models, using to do so a set of training data (x i , y i ) N i=1 generated from an unknown joint probability distribution P X,Y (see figure 1 for an illustration). After getting an "optimal" model ϕ, we must decide what is the label of a new unlabelled instance (x, •). This latter decision step can be handled by tools issued from decision theory. In machine learning, induction is often seen as the task of determining a decision function that will minimise the risk of getting misclassifications. The cost or risk of a misclassification is generally quantified through a function L : K × K → R penalising every bad decision, known as loss function, with L (y, ŷ) the loss incurred when predicting ŷ if y is the observed, true value. The optimal model is most often defined as the one minimizing the expected loss.

Definition 1 (Risk minimizing [6, §2.4]). Given a general loss function L (•, •), the optimal model is defined as the one minimizing the average loss of getting missclassification.

φ := argmin ϕ(X)∈Φ E X ×K [L (Y, ϕ(X))]
(2)

If loss function is defined instance-wise, then, Equation (2) can also be expressed as the minimization of conditional expectation [6, eq. 2.21]:

φ := arg min y∈K E Y |X [L(y, ϕ(X))] (3)
Classical accuracy corresponds to a zero-one loss function, where all missclassification are penalised identically, i.e. L 0/1 (y, ŷ) is equal to 1 if y and ŷ are different and 0 otherwise. Therefore, given L 0/1 , we can reformulate the risk minimization as the well-known Bayes classifier, which would choose the learning model maximizing the conditional probability (a.k.a. maximum a posterior (MAP) probability) given a new unlabeled instance x:

ϕ(x) = argmax m k ∈K P (Y = m k |X = x) (4)
Hence, in a precise probabilistic approach, the main task is to estimate the conditional distribution P Y |X , from which can be obtained the optimal decision. An alternative way of looking at this decision-making problem is to pose it as a problem of inferring preferences between the labels, as follows:

Definition 2 (Precise ordering [20, pp. 47]). Given a general loss function L(•, •) and a conditional probability distribution P Y |x , m a is preferred to m b , denoted by m a m b , if and only if:

E P Y |x [L(•, m a )|x] < E P Y |x [L(•, m b )|x] (5)
Definition 2 tells us that exchanging m b for m a would incur a positive expected loss, due to the fact that expectation loss of m b is greater than m a , therefore m a should be preferred to m b for a given new unlabelled instance x. In the particular case where we use the loss function L 0/1 , it is easy to prove that:

m a m b ⇐⇒ P (Y = m a |X = x) > P (Y = m b |X = x) (6)
where P (Y = m a |X = x) is the unknown conditional probability of label m a and a new unlabeled instance x. Therefore, given a set of labels K, we can then establish a complete preorder making pairwise comparisons (see figure 2) as follows:

m i K m i K-1 • • • m i1 ⇐⇒ P (Y = m i K |X = x) > • • • > P (Y = m i1 |X = x). (7)
We can then pick out one of the undominated labels, i.e., one with maximal probability . In the case where we find 1 < r ≤ K equal maximal conditional probabilities, which is unlikely in practice but not impossible, they can be considered as indifferent and chosen randomly. It should be noted that, whatever the quantity of data used to induce the model or the specific new instance x we face (that may come from a poorly populated region), we will always get (up to indifference) a unique undominated model. In contrast, the IP approach where we consider sets P Y |X may result in partial orders having multiple undominated and incomparable labels.

Classification with imprecise probabilities

Often, the decision maker can be faced with unreliable or hard situations where making a single decision may give rise to serious mistakes (e.g. cancer screening). The hardness of such situations can for instance be due to the lack of sufficient evidence or information (i.e. uncertainty in data). Such cases could be dealt more reliably via a cautious decision (i.e. a set of plausible choices, see e.g., Figure 3).

Group B Group A As shown by Equation 7and Example 1, usual statistics and probabilities usually model uncertainty in data with single distribution P, canonically ending up in a unique undominated label. While it is possible to implement decision rules providing set-valued predictions in such settings [START_REF] Ha | The optimum class-selective rejection rule[END_REF], several authors [START_REF] Walley | Statistical reasoning with imprecise Probabilities[END_REF][START_REF] Dempster | A generalization of bayesian inference[END_REF] have argued that a single distribution cannot always faithfully represent lack of information. Following them, we consider the framework of imprecise probabilities to account for such lack of information.

•? P (ŷ * |X = x * ) ≈ 0.5 (a) Precise decision-making Group B Group A •? ŷ * ⊆ {A, B} (
In what follows, we will introduce some essential concepts about imprecise probabilities, as well as an extension of the precise decision approach given in Definition 2 relying on the convex set of distributions P y|X . As in this work we want to obtain more reliable decision by allowing partial predictions, i.e. cautious decisions, we will focus on extensions where we get a set-value label Y instead of a precise label y.

Basic notions about imprecise probabilities

Imprecise probability theory often (and will in our case) consists in representing our uncertainty by a convex set P X of probability distributions [START_REF] Walley | Statistical reasoning with imprecise Probabilities[END_REF][START_REF] Augustin | Introduction to imprecise probabilities[END_REF] (i.e. a credal set [START_REF] Levi | The enterprise of knowledge: An essay on knowledge, credal probability, and chance[END_REF]), defined over a space X rather than by a precise probability measure P X [START_REF] Taylor | Introduction to measure and integration[END_REF]. As they include precise distributions as special cases, such convex sets of distributions provides richer, more expressive models of uncertainty that allow us to better describe uncertainty originating from imperfect or scarce data.

Given such a set of distribution P X and any measurable event A ⊆ X , we can define the notions of lower and upper probabilities P X (A) and P X (A), respectively as:

P X (A) = inf P ∈P X P (A) and P X (A) = sup P ∈P X P (A) (8)
where P X (A) = P X (A) only when we have sufficient information about A.

Estimations of parameters in the context of imprecise probabilities is usually more complicated as we consider a set P X of distributions instead of a single distribution P X . In our case, such complications will be limited, as we will rely on previous works providing efficient generalized Bayesian inference methods for exponential families (which include Gaussian distributions), that we will present in Section 3.2 For theoretical developments of the next subsection about decision, we will assume that we know the set P Y |X of conditional distributions over the classes.

Decision making under imprecise probabilities

Within IP theories, we can find different methods extending the decision criterion given in Definition 2 (more details in [START_REF] Troffaes | Decision making under uncertainty using imprecise probabilities[END_REF]). For classifying a new instance x, we will make use of the maximality criterion [7, §8.6] that has strong theoretical justifications [14, §3.9.5] and often remains applicable in practice [START_REF] Zaffalon | The naive credal classifier[END_REF][START_REF] Yang | Cautious classification with nested dichotomies and imprecise probabilities[END_REF][START_REF] Benavoli | Classification with imprecise likelihoods: A comparison of tbm, random set and imprecise probability approach[END_REF]. This one extends Equation ( 5) in a robust way, requiring that a preference holds only if it holds for every model. More precisely, the criterion of maximality is defined as follows:

m a M m b ⇐⇒ inf P Y |x ∈P Y |x E P Y |x [L(•, m b ) -L(•, m a )|x] > 0 (9) if L(•, •) is 0/1 loss function, this gives: m a M m b ⇐⇒ inf P Y |x ∈P Y |x P (Y = m a |X = x) -P (Y = m b |X = x) > 0 (10)
Equation ( 10) amounts to asking that Equation ( 5) is true for all possible probability distributions in P Y |x . In practice, M can be a partial order with several maximal elements, in which case the prediction becomes imprecise due to high uncertainty in the model. Note that when N → ∞, imprecise and precise models will usually coincide. The prediction Y M resulting from M is defined as:

Y M = m a ∈ Y ∃m b ∈ Y : m b M m a ( 11
)
Example 2. Given the label set K = {m a , m b , m c }, we could have the following plausible partial ordering:

B = {m a M m b , m c M m b }
where Y M = {m a , m c } is the predicted set obtained from set B of comparisons by the criterion of maximality (figure 4). In the next section, we expose how the Gaussian discriminant model can be made imprecise to characterise the conditional distributions P Y |X by a set P Y |X . We will first recall the precise model, before making it imprecise.

Gaussian discriminant analysis model

As mentioned in the introduction, a classical way to estimate the distribution P Y |X is by using Bayes' theorem. Making use of Equation ( 1), we will discuss the precise and imprecise approach, respectively in Sections 3.1 and 3.2.

Statistical inference with precise probabilities

There are many ways to model P X|Y =m k , but in this work, we focus on parametric discriminant analysis which assume that P X|Y =m k follows a multivariate Gaussian distribution N (µ m k , Σ m k ) with unknown mean µ m k and covariance matrix Σ m k , i.e.:

G m k := P X|Y =m k ∼ N (µ m k , Σ m k ) (12)
whose probability density function is written

P (X = x|Y = m k ) = 1 (2π) p/2 |Σ m k | 1/2 e -1 2 (x-µm k ) T Σ -1 m k (x-µm k ) . ( 13
)
The marginal distribution is defined as a multinomial π y := P Y , where P (Y = m k ) = π m k . So, under a 0/1 loss function, the optimal prediction becomes:

ϕ(x|θ m k ) := arg max m k ∈K log π m k - 1 2 log |Σ m k | - 1 2 (x T -µ m k ) T Σ -1 m k (x T -µ m k ) (14)
where

Θ = {θ m k |θ m k = (π m k , Σ m k , µ m k ), ∀m k ∈ K}
is the parametric space from which comes our estimate. In Table 2, we remind different discriminant models arising from of the last equation, and corresponding to various constraints imposed to the covariance matrices of the conditional distributions. [START_REF] Marco | The euclidean distance classifier: an alternative to the linear discriminant function[END_REF] Unit-variance feature indep.: In frequentist inference, usual estimation of parameters of ( 14) is obtained by MLE using a subset Those estimates do not account for the quantity of data they are based on, which may be low to start with, and may also vary significantly across classes, especially in case of imbalanced data sets. To solve this issue, we propose in the next section an imprecise discriminant model, based on the use of imprecise probabilities and using results from Benavoli et al. [START_REF] Benavoli | Prior near ignorance for inferences in the k-parameter exponential family[END_REF].

Discriminant analysis model Assumptions (∀m

k ∈ K) Parametric space (∀m k ∈ K) Parametric Gaussian conditional distribution P y|X Linear Discriminant [6, §4.3] Homoscedasticity: Σ m k = Σ Θ = {θ m k |θ m k = (π m k , Σ, µ m k )} Quadratic Discriminant [6, §4.3] Heteroscedasticity: Σ m k = Σ k Θ = {θ m k |θ m k = (π m k , Σ k , µ m k )} Naive Discriminant [6, §6.63] Feature independence: Σ m k = σ T k I Θ = {θ m k |θ m k = (π m k , σ k , µ m k )} Euclidean Discriminant
Σ m k = I Θ = {θ m k |θ m k = (π m k , µ m k )}
D m k = {(x i,k , y i,k =m k )|i = 1, . . . , n k } ⊆ D of

Statistical inference with imprecise probabilities

To estimate P X|y and P Y in the form of convex sets of distributions, we will use robust Bayesian inference under prior near-ignorance models. Before describing our imprecise estimation, we make three general assumptions for our imprecise Gaussian discriminant model:

1. Normality of conditional probability distribution P X|Y =m k := G m k , as in the classical case.

A precise estimation of marginal distribution P

Y := π y . 3. A precise estimation of covariance matrix Σ k := Σ k = S k or S.
In Section 6, we will discuss the relaxation of assumption 2, considering a set of distributions P Y .

Robust Bayesian inference

The estimation of parameters in Bayesian inference relies mainly on two components; the likelihood function and the prior distribution, from which posterior inferences can then be made on unknown parameters of the model, in our case θ m k .

In the particular case of P X|Y =m k , the likelihood function is the product of conditional probabilities

n k i P x i,k |y i,k
,θm k and the prior distribution P θm k models our knowledge about θ m k = (Σ m k , µ m k ). In this paper, we focus on estimating imprecise mean parameters (i.e. θ m k = µ m k ), assuming a (precise) estimation of Σm k , for reasons of computational complexity that will be discussed in Section 7. Thus, the posterior on the mean is such that

P (µ m k | D m k ) ∝ n k i P (X = x i,k | µ m k , y i,k = m k )P (µ m k ). ( 15
)
To simplify, we will from now on remove the subscript m k , always bearing in mind that these estimations are related to a group of observations labelled m k .

Near-ignorance on Gaussian discriminant analysis

Near-ignorance models allow us to provide an "objective inference" approach, representing ignorance about unknown parameter and letting the data speak for themselves. In their work, Benavoli et al in [START_REF] Benavoli | Prior near ignorance for inferences in the k-parameter exponential family[END_REF] propose a new near-ignorance model based on a set of distribution M , which aims to reconcile two approaches, namely, the re-parametrization invariance and the Walley's near-ignorance prior. For that, Benavoli A. et al define four minimal properties, which must be satisfied whenever there is no prior information about the unknown parameter, on the set of distributions M (more details in [18, §2]).

(P1) Prior-invariance, that states that M should be invariant under some re-parametrization of the parameter space (e.g. translation, scale, permutation, symmetry, etc).

(P2) Prior-ignorance, that states that M should be sufficiently large for reflecting a complete absence of prior information w.r.t unknown parameter, but no too large to be incompatible with property (P3). (P3) Learning from data, that states that M should always provide non-vacuous posterior inferences, in other words, it should learn from the observations. (P4) Convergence, that states that the influence of M on the posterior inference vanishes when increasing number of observations, i.e. n → ∞, requiring consistency with the precise approach at limit. Benavoli et al [START_REF] Benavoli | Prior near ignorance for inferences in the k-parameter exponential family[END_REF] provide a set of conjugate priors M for regular multivariate exponential families [28, §3.3.4] (FExp) that satisfies the last four properties under quite weak assumptions. Borrowing from [START_REF] Benavoli | Prior near ignorance for inferences in the k-parameter exponential family[END_REF], we can define this set of prior distribution M as follows:

Definition 4 (Prior near-ignorance for k-parameter exponential families [18, §4, eq. 16]). Let L be a bounded closed convex subset of R k strictly including the origin ( [18, lem. 4.5]).

L = ∈ R k : i ∈ [-c i , c i ], c i > 0, i = {1, . . . , d} (16)
Let w ∈ W = R k be a real-valued parameter with a density having the following functional form that belongs to FExp:

p(w) = exp( T r) exp( T w)1 Wr (w) ( 17
)
where belongs to L and r ∈ R k is a real value. The set of prior distributions (c.f. [18, th. 4.6]) can be written as follows:

M w = w ∈ W | p(w) ∝ exp( T w), = [ 1 , . . . , k ] ∈ L (18)
Since our Gaussian probability distribution P X|y=m k given by Equation ( 12) belongs to FExp, we can use the set of prior distributions M µ of Equation ( 18) in order to get a set of posterior distributions M µ n having the same functional form (FExp) [17, §5.2]:

M µ n = µ x n , ∝ N T Σ + nx n n , 1 n Σ µ ∈ R n , ∈ L (19)
where

x n = 1 n n i=1
x i and ∈ L. We can then estimate the lower and upper values of the unknown µ parameters, so for every dimension i ∈ {1, . . . , d} [START_REF] Benavoli | Prior near ignorance for inferences in the k-parameter exponential family[END_REF]:

inf M µ n E[µ i |x n , ] = E[µ i |x n , ] = -c i + nx n n (20)
sup

M µ n E[µ i |x n , ] = E[µ i |x n , ] = c i + nx n n (21)
As a result, we will have for each label m k a convex space of plausible values for the mean µ m k which can be represented by the hyper-cube

G m k = µ m k ∈ R d µ i,m k ∈ -c i + n k x i,n k n k , c i + n k x i,n k n k , ∀i = {1, ..., d} . (22)
Remark 1. The convergence property (P4) ensures us that no matter the initial value of our convex space L, when the number of observations tends to infinity, n → ∞, their influence on the posterior inference

200 of μ will disappear, i.e G m k ----→ n→∞
x n , and will become the asymptotic estimator of the precise Gaussian distribution.

On the basis of the set G m k previously calculated, we can simply consider the following set of conditional probability distributions P X|y=m k (or set of predictive distributions) for every label m k on K :

P X|y=m k = P X|Y =m k P X|Y =m k ∼ N (µ m k , Σ m k ), µ m k ∈ G m k (23)
In what follows, we study how we can incorporate the sets of distributions P X|Y =m k in Gaussian discriminant analysis, using maximality (Definition 3) to get our (possibly) imprecise classification.

Imprecise Classification with L 0/1 loss function

Let us now present our approach to make cautious classification by using sets of conditional distribution given by Equation ( 23) and obtained from a near-ignorance model. Using the maximality criterion, to know whether m a M m b , we need to solve Equation ( 10) by applying Bayes' theorem: inf

P Y ∈P Y inf P X|ma ∈P X|ma P X|m b ∈P X|m b P (X = x|Y = m a )P (Y = m a ) -P (X = x|Y = m b )P (Y = m b ) > 0 ( 24
)
since the marginal P (X = x) = m l ∈Y P (X = x|Y = m l )P (Y = m l ), which is the same positive constant of normalisation for each probability, can be omitted.

As conditional distributions sets P X|Y =m k are independent of each others, we can rewrite Equation (24) as follows (cf. [24, eq. 4.3]): inf

P Y ∈P Y P (X = x|Y = m a )P (Y = m a ) -P (X = x|Y = m b )P (Y = m b ) > 0 ( 25
)
where P (P ) is the infimum (supremum) conditional probability. Also, applying Assumption 2 and the fact that every π y > 0, solving Equation ( 25) is reduced to finding the two values

P (X = x|Y = m a ) = inf P X|ma ∈P X|ma P (X = x|Y = m a ), (26)
P (X = x|Y = m b ) = sup P X|m b ∈P X|m b P (X = x|Y = m b ) (27)
As P X|y=m k is a set of Gaussian distributions, the solutions of Equations ( 26) and ( 27) are respectively obtained for the following values of the means

µ ma = arg inf µm a ∈Gm a - 1 2 (x -µ ma ) T Σ -1 m b (x -µ ma ), (28)
µ m b = arg sup µm b ∈Gm b - 1 2 (x -µ m b ) T Σ -1 m b (x -µ m b ), (29)
where Σ -1 m b is the inverse of the covariance matrix (Assumption 3). Depending on the internal structure of the precise covariance matrix Σk , solving for ( 28) and ( 29) may be more or less computationally challenging. We will consider two main different imprecise discriminant models: (1) with non-diagonal covariance matrix and (2) with diagonal covariance matrix.

Gaussian discriminant model with dependent features

Similarly to the distinction made in the precise case, we will consider two different variants of the nondiagonal case.

Case 1. Imprecise Quadratic discriminant analysis (IQDA): if we suppose that the covariance structures of all groups of observations are different, that is

Σ m k = S m k , ∀m k ∈ K .
Case 2. Imprecise linear discriminant analysis (ILDA): if we assume that all groups of observations have the same covariance structure, that is Σm k = Ŝ, ∀m k ∈ K .

In those cases where the covariance matrix contains collinear columns, Σm k will not be invertible, in which case we use the singular value decomposition (SVD) method for computing the pseudo-inverse of covariance matrix. Before studying the computational issues of IQDA and ILDA, i.e. Equations ( 28) and ( 29), we will illustrate the last case (ILDA) in Example 3.

Example 3. The interest of modelling an imprecise mean is to be able to detect areas where we should be cautious and predict sets of labels rather than a single one. For example, in Figure 5, we simulated two groups of observations x ma * and x m b * (i.e. binary case), each with two non-correlated regressors and different means: We also drew the (precise) mean of each group, i.e. µ ma and µ m b , as solid points, and a black dot (•) representing a new unlabelled instance x as well as positions of solutions of Equations ( 28) and (29). In Figure 5, we observe (in purple) an area of uncertainty generated by the imprecise mean and the maximality criterion. Let us now discuss the problem of solving Equations ( 28) and [START_REF] De Angelis | Quadratic programming with box constraints[END_REF]. Expressing G m b as constraints, the solution µ m b of (29) can be written as

xm a 1 xm a 2 ∼ N (( 0.25 0.5 ) , ( 1 0 0 1 )) xm b 1 xm b 2 ∼ N ( 0.5 -1.0 , ( 1 0 0 1 )) L = ∈ R 2 : i ∈ [-c i , c i ], c i = 2
µ m b = arg sup - 1 2 μT m b Σ-1 m b μm b + q T μm b s.t. -c j + n m b x j,n n m b ≤ μj,m b ≤ c i + n m b x j,nm b n m b , q T = -x T Σ-1 m b , ∀j = {1, ..., d} (BQP)
This optimisation problem is well-known as a box-constraint quadratic program (BQP) [START_REF] De Angelis | Quadratic programming with box constraints[END_REF], as (1) the constraint space G m k is a convex space, and (2) Σ -1 m k is a positive (semi)-definite matrix, pending the fact that the covariance matrix Σm k does not have multicollinearity problems [START_REF] Johnson | Positive definite matrices[END_REF]. Computing an optimal global solution of (BQP) in polynomial time is easy using modern optimisation libraries (e.g. using the CvxOpt python library [START_REF] Andersen | Cvxopt: A python package for convex optimization[END_REF]), as we have to maximize a concave function (or, equivalently, minimise a convex one).

Finding µ ma in equation ( 28) is much more difficult, as one seeks to solve the optimization problem

µ ma = arg inf μma ∈Gm a - 1 2 µ T ma Σ-1 ma μma + q T μma (NBQP)
That comes down this time to maximizing a convex function over box-constraints (G ma ), which is known to be NP-Hard [START_REF] Pardalos | Quadratic programming with one negative eigenvalue is np-hard[END_REF]. To solve it, we use a brand-and-bound (B&B) algorithm [START_REF] Burer | Globally solving box-constrained nonconvex quadratic programs with semidefinite-based finite branch-and-bound[END_REF][START_REF] Xia | Globally solving non-convex quadratic programs via linear integer programming techniques[END_REF], that employs a finite branching based on the first-order Karush-Kuhn-Tucker5 conditions and polyhedral semidefinite relaxation in each node of the B&B tree (more details in [START_REF] Burer | Globally solving box-constrained nonconvex quadratic programs with semidefinite-based finite branch-and-bound[END_REF]).

Gaussian discriminant model with independence features

When the number of features becomes high, and the associated optimisation problem quite timeconsuming to solve, it may be interesting to consider some additional assumptions which will significantly reduce the inference complexity. In what follows, we will assume that features x j i are independent conditional on the label m k .This translates in the fact that covariance matrices become diagonal matrices, i.e.

Σ m k = σ T m k I with σ T m k = (σ 1 m k , .
. . , σ p m k ) a p-dimensional vector containing the variance of each feature, which can be interepreted as weights of the features. Therefore, we can rewrite Equations ( 28) and ( 29) as follows:

µ ma = arg inf µm a ∈Gm a - 1 2 w ma x -µ ma 2 ( 30
)
µ m b = arg sup µm b ∈Gm b - 1 2 w m b x -µ m b 2 (31)
where

w m k = (w 1 m k , . . . , w p m k ) T such that w j m k = 1/σ j m k , ∀j ∈ {1, .
. . , p}, in this scenario, we will consider two new models.

Case 3. Imprecise naive discriminant analysis (INDA): this case is similar to the Naive Bayes classifier, as we simply consider the assumption Σ m k = σ T k I where σ k are the empirical variance estimator obtained from a group of observation belonging to the label m k .

Case 4. Imprecise Euclidian discriminant analysis (IEDA): this is a case more specific than INDA, where we assume that for every j ∈ {1, . . . , p} we have σ j k = 1, meaning that the measure used to evaluate the probability of a label given a new instance is proportional to the Euclidian distance between the instance and the corresponding mean. The Euclidean classifier is one of the simplest existing classifier, and is the supervised counterpart of the standard k-means method.

We show below that when the covariance matrix is diagonal, optimisation problems [START_REF] Johnson | Positive definite matrices[END_REF] and [START_REF] Andersen | Cvxopt: A python package for convex optimization[END_REF] become very easy (i.e. linear in p, O(p)) to solve. Proposition 1. For two vectors x, w ∈ R p , and a box-convex space on R p :

G = µ ∈ R p µ j ∈ [µ j , µ j ], ∀j ∈ {1, . . . , p}
-the infimum weighted distance subject to constraints G is:

inf µ∈G - 1 2 w T x -µ 2 = - 1 2 p j w j max j {(x j -µ j ) 2 , (x j -µ j ) 2 } ( 32
)
-and the supremum weighted distance subject to same constraints is:

sup µ∈G - 1 2 w T x -µ 2 = - 1 2 p j w j 0 if x j ∈ [µ j , µ j ] min j {(x j -µ j ) 2 , (x j -µ j ) 2 } otherwise (33)
Proof. Since each element of the sum is positive, we can interchange the infimum operator with summation, and calculate the supremum of each component as follows:

inf µ∈G - 1 2 p j w j (x j -µ j ) 2 ⇐⇒ - 1 2 p j w j sup µ j ∈[µ j ,µ j ] (x j -µ j ) 2 (34)
where the supremum can be calculated as follows:

sup

µ j ∈[µ j ,µ j ] (x j -µ j ) 2 = max j {(x j -µ j ) 2 , (x j -µ j ) 2 } ( 35
)
In the second case and for similar reasons, we can also put the supremum operator inside of summation and calculate of infimum value of each component:

sup µ∈G - 1 2 p j w j (x j -µ j ) 2 ⇐⇒ - 1 2 p j w j inf µ j ∈[µ j ,µ j ] (x j -µ j ) 2 (36)
where the infimum of squared subtraction of each element is: inf

µ j ∈[µ j ,µ j ] (x j -µ j ) 2 = 0 if x j ∈ [µ j , µ j ] min j {(x j -µ j ) 2 , (x j -µ j ) 2 } otherwise ( 37
)
The next section presents some experiences with different data sets and different precise and imprecise models.

Experiments setting

In this section, we provide experimental results evaluating the performance of our different imprecise Gaussian discriminant models (cf. Section 4).

How can we choose parameter c i ?

The choice of parameters c i determines the amount of imprecision in our posterior inference. It should be large enough to guarantee more reliable predictions when missing information, but small enough so as to provide informative predictions when possible. Therefore, in the absence of prior information and for symmetry reasons, we will consider a symmetric box around 0, as follows:

L = ∈ R k : i ∈ [-c, c], c > 0, i = {1, . . . , d} . (38)
In order to fix a value of c, there exists different approaches already mentioned in Section 4.3 of [START_REF] Benavoli | Prior near ignorance for inferences in the k-parameter exponential family[END_REF]. One can for example rely on the rate of convergence of the lower and upper posterior expectations [START_REF] Walley | Statistical reasoning with imprecise Probabilities[END_REF]:

∀i E[µ i |x n , ] -E[µ i |x n , ] = 2c n ----→ n→∞ 0 ( 39
)
meaning that for small values of c, we would reach a faster convergence of Equation ( 39) to a precise posterior inference (as precise models). A value of c ≤ 0.75 is recommended by [18, §4.3, §8], however since we are in a classification problem, we will select an optimal value of c through cross-validation on the training samples. More precisely, we restrict c to the interval [0.01, 5], discretised into [0.01, 0.02, ..., 5], with the optimal value decided by cross validation on the training samples. A typical empirical evolution of the accuracy measures used in the next sections is shown in Figure 8

Data sets and experimental setting

We perform experiments on 15 data sets issued from UCI machine repository [START_REF] Frank | UCI machine learning repository[END_REF](cf. Table 3), following a 10×10-fold cross-validation procedure. We aim to compare the performance of our imprecise Gaussian classifier model approach with the existing precise models (c.f. Table 2).

Owing to small amounts of samples in some groups of observations (belonging to a specific label m k ) of some data sets, the QDA model can suffer from a phenomenon known as ill-posed covariance matrix (i.e. n m k < p), and in such cases even calculating the pseudo-inverse of Σ mm estimated covariance matrice using SVD method cannot solve the problem. This affects the performance of our classifiers, getting highly significant drop (e.g. in Table 4, glass and yeast data sets). Therefore, in this case specific, we used a basic regularized method for estimated covariance matrice named Regularization QDA (or RQDA) [START_REF] Friedman | Regularized discriminant analysis[END_REF][START_REF] Friedman | The elements of statistical learning[END_REF]:

Σ m k (α) = α Σ m k + (1 -α)I, (40)
where Σ mm is the estimated covariance for a group of observations, I a identity matrix and α the regularization factor.

Comparing indeterminate predictions given in the form of a subset Y of plausible labels against just one determinate prediction y is a hard problem that mostly depends on the circumstances or the context in which a decision-marker may or may not accept partial predictions (or cautious decision) instead of a unique, risky decision. A good evaluation should reward cautiousness provided by Y when it allows to include the true observed label, but not so much as to systematically privilege imprecision over precision. In other words, we need an evaluation metric that seeks a compromise between cautiousness and informativeness. To do this, we adopt the evaluation metric proposed and theoretically justified in [START_REF] Zaffalon | Evaluating credal classifiers by utility-discounted predictive accuracy[END_REF], called utility-discounted accuracy, which makes it possible to reward the imprecision in a more or less strong way. It is written as follows:

u(y, Y ) = 0 si y / ∈ Y, α |Y | -1-α |Y | 2 autrement.
(41) [START_REF] Zaffalon | Evaluating credal classifiers by utility-discounted predictive accuracy[END_REF] shows that a value α = 1 amounts to not reward cautiousness and to confuse it with randomness, while α → ∞ does not penalize non-informativeness, as the vacuous prediction (i.e. Y = K ) would always get 275 a full, guaranteed reward. We will use the usual values u 65 with α = 1.6 and u 80 with α = 2.2 (as in [START_REF] Yang | Cautious classification with nested dichotomies and imprecise probabilities[END_REF]). To have an intuition about these measures, let us simply recall that the u 65 (u 80 ) measure rewards a binary correct prediction with 0.65 (0.80), while a purely random, non-cautious guesser picking one of the two possible label would reward it with 0.50. It therefore gives a "reward" of 0.15 (0.30) for rightful cautiousness.

Experimental results

The average results obtained according to u 65 and u 80 utilities, and the average execution time to predict the label of a new unlabeled instance are shown in Table 4 First, we can see that including some cautiousness can increase our accuracies on most data sets, by picking the right values of c. This increase is sometimes noticeable, for example in the vehicle (i), wine-quality (l), wall-following (o) and vowel data sets (j). All of this, keeping a time execution reasonable in view of the problems to be solved (e.g. a non-convex, NP-hard problem), and without an optimized implementation. As expected, assuming independence between the features (i.e., diagonal covariance matrices) significantly reduces the computational time, making it negligible, but overall reduces performances, as the assumptions are often violated in a stronger way.

In order to highlight the major role of cautiousness of an imprecise classifier model, we show in Figure 6(c) and 6(b) how, in the IRIS data set, our IQDA and ILDA models create different areas of decision boundaries (not to be confused with rejection area), where each area has a different combinations of subset of labels Y ⊆ K , in contrast to precise classifier model (LDA), in Figure 6(a), where it creates one area per label. We can clearly see that the two classifiers behave quite differently. In particular, ILDA will induces regions delimited by piece-wise linear functions, while IQDA will induces regions delimited by piece-wise quadratic functions. Also, in Figure 8, we show the evolution of utility-discounted accuracy (i.e. u 65 and u 80 of vowel dataset), with a standard deviation calculated by a 10-fold cross-validation on the training dataset, according to the imprecision of estimators µ. As expected we notice that when c reaches a too high value, the overall model performances decrease, as it becomes too imprecise with respect to our attitude towards cautiousness (modelled through utility (41)). The rest of experiments are in Appendix A. An imprecise classifier should abstain (i.e. by providing a set of plausible choices) on those hard instances, that is the instances where the precise classifier makes a unusual high amount of mistakes. In Figure 7, we verify that our imprecise classifiers follow this desirable behaviour on most data sets, for the u 80 measures (conclusions for the u 60 are similar, but not displayed to gain some space). Figure 7(a) displays the percentage of time the true label is in the prediction of ILDA, given that the prediction was imprecise, versus the accuracy of LDA on those same instances. The same graphs for the QLDA, IEDA and INDA methods are given by Figure 7(b), Figure 7(c) and Figure 7(d), respectively. We notice that on those hard instances where precise classifiers are wrong, our imprecise classifiers successfully overcome them, getting the ground-truth value into partial predictions (most often > 80%). A typical and quite remarkable example of this is the dermatology data set (h) for the linear case, where the accuracy on the imprecisely classified instances drop to 30% for the precise classifier (to be compared to an average of 96% on all instances), while the imprecise classifier always include the true class. Moreover, the fact that u 80 is higher indicates that the overall amount of imprecision remains acceptable. Our approach therefore seems to be able to well robustify the very simple, linear decision frontiers of the ILDA models. Before considering some generalisation of the presented methods, we would also like to mention that the imprecise probabilistic approach will in general induces decision frontiers that are different from classical rejection rule. Figure 5(b) illustrates this well: rejection regions in a binary setting are most often equivalent to require to predict {a, b} whenever P ({a}|x) ∈ [0.5 -, 0.5 + ] for some . This means that in the case of LDA, the rejection regions will be delimited by two parallel lines, corresponding to the iso-density points x for which P ({a}|x) = 0.5and P ({a}|x) = 0.5 + . In contrast, we can clearly see in Figure 5(b) that the boundaries are not linear, but piece-wise linear.

Imprecise prior marginal and generic loss functions

In this section, we will discuss about two new variants of IGDA model: (1) relaxing the Assumption 2, i.e. P Y := π, with the purpose of putting a set of probability distributions P Y instead, and (2) dealing with generic loss function instead classical L 0/1 loss function. We will evaluate the impact of this two new variants in our IGDA model in terms of added computational complexity.

Imprecise prior marginal

The first extension we will consider is to make imprecise the marginal distribution, considering a set P Y rather than a precise distribtuion, in the same vein as we have made the conditional distribution P X|Y imprecise. For the time being, we will still work with the L 0/1 loss function. Since the conditionals are still independent of each other, solving the maximality criterion amounts to solve Equation ( 25), that we recall here inf

P Y ∈P Y P (X = x|Y = m a )P (Y = m a ) -P (X = x|Y = m b )P (Y = m b ) (42)
with m a M m b if this is positive. This equation can be solved easily, as it is a linear form in P (Y = m a ), P (Y = m b ), meaning that we can either use linear programming over the constraints induced by P Y , or find the extreme point (e.g., by enumeration) of P Y for which the solution is obtained. Recall t The problem then amounts to estimate P Y . A quite popular choice to do so is to use an Imprecise Dirichlet Model (IDM) [START_REF] Walley | Inferences from multinomial data: learning about a bag of marbles[END_REF][START_REF] Bernard | An introduction to the imprecise dirichlet model for multinomial data[END_REF]. However, as Benavoli et al. has already mentioned in [18, §4.2], the set of prior distributions of IDM does not correctly satisfy (P1) Prior-invariance property and permutation invariance of near-ignorance model. So, to remain consistent with our previous estimates, we explore another solution proposed by Benavoli et al..

Let Y be a discrete random variable on a finite space of labels K with probability distribution P Y and let the parameters π m k , ∀m k ∈ K be the unknown non-negative chances, i.e. P (Y = m k ). The Corollary 4.10 in [START_REF] Benavoli | Prior near ignorance for inferences in the k-parameter exponential family[END_REF] propose adding some contraints in the space L in order not to favour some chances π m k over others. They then consider the following set of prior distributions:

P π = π 1-1 m1 π 2 -1 m2 . . . π -d-1 i=1 i-1 m K , 1 ≤ 2c, d-1 i=1 i ∈ [-c, c] . (43)
It is also shown [START_REF] Benavoli | Prior near ignorance for inferences in the k-parameter exponential family[END_REF]Eq. 24] that, after combining this set with the likelihood, the lower and upper expectations of the chances of observing a given subset A of categories result in

E m k ∈A π m k n, ŷn = min 1, 1 n m k ∈A n k + c := P Y (A), (44)
E m k ∈A π m k n, ŷn = max 0, 1 n m k ∈A n k -c := P Y (A), ( 45
)
where n is the total number of observations in the data set, i.e n = |D| = N . We will then consider the probability set P Y = P P Y (A) ≤ P (A) ≤ P Y (A), ∀A ⊆ K (46) Such a model, which corresponds to take a neighbourhood around the empirical distribution using the total variation distance (i.e., L ∞ norm) has been recently investigated by Miranda et al. [START_REF] Miranda | A unifying frame for neighbourhood and distortion models[END_REF], showing for instance that it induced a 2-monotone lower probability, but was not a specific case of probability intervals, in contrast with the IDM model. Using this fact, we know that the result of Equation ( 25) will be obtained by the Choquet integral, which results in this particular case in In particular, this shows that there would be no differences if we considered only the projections of P Y over its singletons, which amounts to consider the bigger set

P Y = P (Y = m k ) = π m k π m k ∈ max 0, n k -c n , min 1, n k + c n ∀m k ∈ K . (47)

Generic loss function

The zero-one loss function is the default loss function used in classification problems (where we consider that penalty of being wrong is the same for every kind of error). However, in many practical problems different errors will have different impacts, and this is especially true for sensitive applications in which imprecise probabilistic approaches could be useful.

Equation Equation ( 9) can be written that uses the fact that the conditional probabilities P (X = x|Y = m k ) are all independent. As Equation (50) remains a linear form of the probabilities P (Y = m k ), it can be solved as previously, i.e., thourgh the use of linear programming or the identification of the extreme point for which the bound is reached.

If we now consider the credal set given by Equation ( 46) and induced by the constraints (44)-(45), we still have that this induces a 2-monotone lower probability, meaning that we can estimate (50) by using the Choquet integral.

All these remarks show that making the marginal probabilities imprecise or considering generic loss functions does not make the model more complex to use, as the computational complexity is not increased by much, especially when P Y has mathematical properties making computations easier (which is luckily the case for most IP models over multinomial distributions).

Conclusion

In this paper, we have generalized classical Gaussian discriminant models to the imprecise setting, mainly by allowing the estimated means of the conditional Gaussian distributions to become imprecise. This was achieved by a robust Bayesian procedure using sets of prior satisfying near-ignorance properties.

We have explored the computational issues associated to the predictions of such models, essentially showing that considering general covariance matrices ended up in practically manageable but theoretically difficult to solve problems, while considering diagonal covariance matrices essentially made the problem much easier to solve.

Experiments on various data sets shows that the method is providing quite satisfactory results, in the sense that the induced imprecision in the predictions is reasonable and mostly concerns instances that were wrongly classified by the precise methods. We have also discussed some possible extensions of our approaches, showing that such extensions would not add a prohibitive computational cost.

A natural next step would be to also make the covariance matrix estimate imprecise, possibly leaving the mean estimate precise in a first step. Computationally, this would be attractive, as the objective functions are a mainly linear of the covariances if the mean is left imprecise. The main problem would then be to derive a principled approach (i.e., using near-ignorance prior) that would deliver an easy-to-deal convex set of inverse covariance matrices.

Figure 1 :

 1 Figure 1: Learning model steps. Figure (a) show the initial training data, from which are induced the boundaries defining the decision function (b), then used to perform the predictions (c).

Example 1 . 1 P

 11 Given a set of labels K = {m a , m b , m c }, a new unlabeled instance x, and the probability estimates of the conditional distribution P Y |X : P (Y = m a |X = x) = 0.3 P (Y = m b |X = x) = 0.(Y = m c |X = x) = 0.6 the complete preorder between labels w.r.t estimated probabilities is m c m a m b .{m c } being the maximal label dominating other ones (Figure2), it is the predicted one.

m c m a m b Figure 2 :

 b2 Figure 2: Graph of complete preorder on labels

Figure 3 :

 3 Figure 3: Cautious vs precise decision-making.Figure (a) shows a precise model, where there are no regions where the model will output set-valued predictions, in contrast with (b) where such a region exists (in red).

  Figure 3: Cautious vs precise decision-making.Figure (a) shows a precise model, where there are no regions where the model will output set-valued predictions, in contrast with (b) where such a region exists (in red).

Definition 3 (

 3 Partial Ordering by Maximality Criterion [21, §3.2]). Let L(•, •) be a general loss function and P Y |X a set of probability distributions, then under the maximality criterion, m a is preferred to m b iff the cost of exchanging m a with m b have a positive lower expectation:

Figure 4 :

 4 Figure 4: Graph of partial order of set B .

  observations of training data. We have πm k = n k /N (frequency of m k ) and μm k = x k (sample mean of D m k ). Depending on whether we assume the model to have (1) dependent features, we will have an hetero-or homo-scedastic assumption, with respectively Σm k = S m k (sample covariance matrix of D m k ) or Σm k = S (within-class covariance matrix D), or to have (2) independent features, we will have features weighted proportionally to their inverse variance Σm k = σT k I or unweighed with all weights equal to 1, i.e. Σm k = I.

Figure 5 (

 5 Figure 5(a) illustrates this example and pictures the following things: groups of observations x ma * and x m b * with the symbols and , respectively, and the posterior convex estimates G (solid) of the means after injecting the information contained in the training data.We also drew the (precise) mean of each group, i.e. µ ma and µ m b , as solid points, and a black dot (•) representing a new unlabelled instance x as well as positions of solutions of Equations (28) and (29). In Figure5, we observe (in purple) an area of uncertainty generated by the imprecise mean and the maximality

Figure 5 :

 5 Figure 5: Imprecise boundary area and estimation.Figure 5(a) shows an example of imprecise estimation of means µ * , and Figure 5(b) shows an imprecise decision area of purple colour where the subset Y ={ma, m b } of labels is the imprecise decision, that is in this region ma and m b are incomparable.

Figure 6 :

 6 Figure 6: Figure 6(a) shows how a precise model divides the instance space in three single different zones by label (i.e {a}, {b}, {c}), the Figure 6(b) shows how an ILDA model divides the instance space in different zones as much as different combinations of a subset of labels (i.e {a}, {b}, {c}, {a, b}, {b, c}, and so on), and the Figure 6(c) shows how IQDA model can also divide in different zones with smooth curves instead.

Figure 7 :

 7 Figure7: Correctness of the different methods in the case of abstention versus accuracy of their precise counterparts, only on those instances for which an indeterminate prediction was given. Graphs are given for the u 80 accuracies.

Figure 8 :

 8 Figure 8: Figures shows performance evolution with a standard deviation region of three principales methods, (1)Figure 8(d) for ILDA model, (2) Figure 8(c) for IQDA model, and (3) Figure 8(b) for INDA model, w.r.t. utility-discount accuracy u 65 , u 80 and c tuning parameter on vowel dataset

infP

  Y ∈P Y P (X = x|Y = m a )P (Y = m a ) -P (X = x|Y = m b )P (Y = m b ) = P (X = x|Y = m a )P Y ({m a }) -P (X = x|Y = m b )P Y (Y = m b ) = P (X = x|Y = m a ) max 0, n ac n -P (X = x|Y = m b ) min 1, n b + c n

m

  a M m b ⇐⇒ infP X|m * ∈P X|m * P Y ∈P Y m k ∈K (L(m k , m b ) -L(m k , m a )) P (Y = m k |X = x) > 0,(48)which, if we denote by c b-am k := L(m k , m b ) -L(m k , m a ), gives ⇐⇒ inf P X|m * ∈P X|m * P Y ∈P Y m k ∈K c b-a m k P (X = x|Y = m k )P (Y = m k ) > 0 (49) ⇐⇒ inf P Y ∈P Y {k|c b-a m k >0} c b-a m k P (X = x|Y = m k )P (Y = m k ) + {k|c b-a m k ≤0} c b-a m k P (X = x|Y = m k )P (Y = m k ) > 0 (50)

  x = (x 1 , . . . , x p ) T

		New unlabeled instance to predict.
	ŷ	Precise output prediction.
	Ŷ	Set-valued output predictions.
	N	Number of training instances.
		Decision theory
	Θ	Parameter space
	P = {P

θ , θ ∈ Θ} Set of probability distributions P Y |X , P X|Y , P Y (Conditional) Probability distributions

Table 1 :

 1 Mathematical notations used in this paper

	70

Table 2 :

 2 Gaussian discriminant analysis models

Table 3 :

 3 Data sets used in the experiments

	#	name	# instances # features # labels
	a	iris	150	4	3
	b	wine	178	13	3
	c	forest	198	27	4
	d	seeds	210	7	3
	e	glass	214	9	6
	f	ecoli	336	7	8
	g	libras	360	91	15
	h	dermatology	385	34	6
	i	vehicle	846	18	4
	j	vowel	990	10	11
	k	yeast	1484	8	12
	l	wine quality	1599	11	6
	m	optdigits	1797	64	10
	n	segment	2300	19	7
	o wall-following	5456	24	4

for the four IGDA methods. It clearly shows that performances first increase in average with imprecision, but then degrades as imprecision becomes too large.

Table 4 :

 4 . Average utility-discounted accuracies (%) and time to predict in seconds.

		(a) LDA versus ILDA						(b) QDA versus IQDA
		LDA	ILDA		Avg.		QDA	RQDA	IQDA		Avg.
	#	acc.	u80	u65	Time	#	acc.	acc.	u80	u65	Time
	a b c d e f g h i j k l m n o avg.	97.96 ± 0.05 98.85 ± 0.36 94.61 ± 0.60 96.35 ± 0.25 62.15 ± 0.76 87.14 ± 0.37 64.45 ± 0.57 96.58 ± 0.35 77.96 ± 0.48 60.10 ± 0.68 58.92 ± 0.17 59.25 ± 0.27 95.40 ± 0.09 91.60 ± 0.09 67.96 ± 0.07	98.38± 98.99 ± 1.17 98.95 ± 1.26 97.16± 94.56 ± 1.08 94.05 ± 1.02 96.59 ± 0.23 96.51 ± 0.23 66.78 ± 0.73 58.87 ± 0.77 88.27 ± 1.43 87.72 ± 1.42 97.06 ± 0.62 96.94 ± 0.61 81.98 ± 0.91 79.59 ± 0.82 67.45 ± 0.48 62.41 ± 0.40 65.83± 60.31± 90.76 ± 0.35 89.70 ± 0.32 71.34± 66.65±	0.56 1.49 12.14 1.50 12.40 19.24 3.10 4.95 34.85 10.77	a b c d e f g h i j k l m n o avg.	97.29 ± 0.44 99.03 ± 0.45 89.43 ± 1.34 97.47 ± 3.37 96.66 ± 4.47 98.89 ± 2.22 94.64 ± 0.47 94.29 ± 2.86 7.15 ± 2.39 51.40 ± 9.79 46.19 ± 2.97 88.25 ± 5.97 34.04 ± 2.14 72.22 ± 6.21 82.47 ± 0.42 96.92 ± 0.88 85.07 ± 0.86 85.11 ± 2.63 87.83 ± 0.49 87.07 ± 3.49 13.18 ± 2.37 56.27 ± 2.29 55.62 ± 0.47 55.79 ± 5.35 87.18 ± 1.31 98.94 ± 0.85 64.69 ± 1.82 91.17 ± 1.61 65.87 ± 0.17 70.56 ± 2.63	98.08 ± 0.41 97.13 ± 0.42 99.39 ± 0.14 99.09 ± 0.13 91.77 ± 1.38 88.90 ± 1.32 95.20 ± 0.26 94.72 ± 0.24 64.38 ± 1.36 58.36 ± 1.30 87.34 ± 0.90 84.79 ± 0.87 84.24 ± 0.87 84.05 ± 0.88 87.96 ± 0.34 86.13 ± 0.27 89.96 ± 0.67 88.40 ± 0.70 65.85± 60.36± 71.79 ± 0.12 69.75 ± 0.12
			(c) NDA versus INDA					(d) EDA versus IEDA
		NDA	INDA		Avg.			EDA		IEDA		Avg.
	#	acc.	u80	u65	Time		#	acc.	u80	u65	Time
	a b c d e f g h i j k l m n o avg.	95.07 ± 0.44 95.73 ± 5.78 97.70 ± 0.58 98.39 ± 3.24 93.60 ± 13.93 95.53 ± 5.94 95.26 ± 0.33 89.95 ± 15.28 99.95 ± 18.59 90.38 ± 0.19 91.14 ± 5.88 88.21 ± 11.32 43.92 ± 1.36 51.25 ± 10.84 50.21 ± 10.41 82.39 ± 1.22 56.54 ± 18.69 56.54 ± 18.69 63.17 ± 1.43 65.50 ± 4.76 65.22 ± 5.30 85.52 ± 0.98 90.71 ± 4.57 90.45 ± 4.64 45.63 ± 0.89 46.33 ± 7.44 45.70 ± 7.14 67.26 ± 0.39 70.35 ± 7.44 71.70 ± 10.23 43.36 ± 0.51 49.51 ± 6.78 50.71 ± 5.81 54.83 ± 0.34 57.24 ± 8.27 57.33 ± 3.31 89.69 ± 0.19 85.17 ± 7.75 89.39 ± 7.01 79.83 ± 0.11 77.11 ± 9.23 78.75 ± 10.31 52.55 ± 0.12 51.00 ± 6.64 54.32 ± 5.09	0.46 × 10 -3 1.72 × 10 -3 4.09 × 10 -3 1.12 × 10 -3 4.98 × 10 -3 3.41 × 10 -3 10.35 × 10 -3 3.10 × 10 -3 1.19 × 10 -3 3.14 × 10 -3 2.50 × 10 -3 1.51 × 10 -3 4.44 × 10 -3 2.69 × 10 -3 1.29 × 10 -3		a b c d e f g h i j k l m n o avg.	91.60 ± 0.61 94.80 ± 4.24 46.65 ± 0.85 61.78 ± 3.93 81.09 ± 0.39 82.38 ± 7.42 90.38 ± 0.36 92.86 ± 8.05 46.26 ± 1.68 56.16 ± 12.06 49.08 ± 8.85 93.13 ± 4.75 52.04 ± 4.08 81.84 ± 9.30 87.29 ± 8.94 42.59 ± 0.04 43.34 ± 9.40 40.91 ± 7.08 49.36 ± 1.53 54.62 ± 7.21 49.27 ± 10.53 12.25 × 10 -3 0.29 × 10 -3 0.41 × 10 -3 1.31 × 10 -3 0.46 × 10 -3 0.67 × 10 -3 0.87 × 10 -3 51.22 ± 0.92 54.78 ± 10.44 52.87 ± 10.99 1.39 × 10 -3 28.03 ± 0.19 45.13 ± 6.03 38.33 ± 4.06 0.67 × 10 -3 58.08 ± 0.90 63.94 ± 5.28 58.94 ± 4.60 2.48 × 10 -3 31.27 ± 0.13 31.56 ± 3.10 31.39 ± 2.88 2.18 × 10 -3 19.72 ± 0.19 22.11 ± 5.64 23.40 ± 8.44 1.03 × 10 -3 88.17 ± 0.15 89.11 ± 1.95 88.83 ± 1.82 5.09 × 10 -3 22.10 ± 0.11 39.30 ± 5.87 30.65 ± 1.10 2.15 × 10 -3 57.90 ± 0.11 56.66 ± 7.01 58.47 ± 9.07 1.06 × 10 -3

In the sense that it goes from singular observations to a generic model.

Also known as KKT, which allows to solve problems of optimisation subject to non-linear constraints in the form of inequalities.

Acknowledgements

This work was carried out in the framework of the Labex MS2T, funded by the French Government, through the National Agency for Research (Reference ANR-11-IDEX-0004-02).

Appendix A. Experimental results

Complementary experimental results are shown in the