Emmanuelle Anceaume
email: emmanuelle.anceaume@irisa.fr

Aina Ravoaja
email: aina.ravoaja@irisa.fr

INCENTIVE-BASED ROBUST REPUTATION MECHANISM FOR P2P SERVICES

Keywords: peer-to-peer, reputation, ballot stu ng, bad mouthing, opinions, incentive peer-to-peer, r putation, collusion, opinions, incitation

In this paper, we address the problem of designing a robust reputation mechanism for peer-to-peer services. The mechanism we propose achieves high robustness against malicious peers from individual or collusive ones and provides incentive for participation. We show that the quality of the reputation value of trustworthy and participating peers is always better than the one of cheating and non participating ones. Finally we formally prove that, even when a high fraction of peers of the system exhibits a collusive behavior, a correct peer can still compute an accurate reputation mechanism towards a server, at the expense of a reasonable convergence time.

I R I S A 1 Introduction

With the emergence of e-commerce in open, large-scale distributed marketplaces, reputation systems are becoming attractive for encouraging trust among entities that usually do not know each other. A reputation system collects, distributes, and aggregates feedback about the past behavior of a given entity. The derived reputation score is used to help entities to decide whether a future interaction with that entity is conceivable or not. Without reputation systems, the temptation to act abusively for immediate gain can be stronger than the one of cooperating. In closed environments, reputation systems are controlled and managed by large centralized enforcement institutions. Designing reputation systems in P2P systems has to face the absence of such large and recognizable but costly organizations capable of assessing the trustworthiness of a service provider. The only viable alternative i s to rely on informal social mechanisms for encouraging trustworthy behavior 8 . Proposed mechanisms often adopt the principle that "you trust the people that you know best", just like in the word-of-mouth system, and build transitivity trust structures in which credible peers are selected 18, 1 9 , 2 0 . However such structures rely on the willingness of entities to propagate information. Facing free-riding and more generally under-participation is a well known problem experienced in most open infrastructures 2 . The e ciency and accuracy of a reputation system depends heavily on the amount of feedback it receives from participants.

According to a recognized principle in economics, providing rewards is an e ective w ay t o improve feedback. However rewarding participation may also increase the incentive for providing false information. Thus there is a trade-o between collecting a sizable set of information and facing unreliable feedback 6 . An additional problem that needs to be faced with P2P systems, is that peers attempt to collectively subvert the system. Peers may collude either to discredit the reputation of a provider to lately bene t from it bad mouthing, or to advertise the quality of service more than its real value to increase their reputation ballot stu ng. Lot of proposed mechanisms break down if raters collude [START_REF] Dellarocas | Reputation mechanisms[END_REF] .

In this paper we address the robust reputation problem. Essentially this problem aims at motivating peers to send su ciently honest feedback in P2P systems in which peers may free-ride or be dishonest. This work has been motivated by a previous one in which the proposed architecture is built on top of a supervising overlay made of trusted peers 3 . The mechanism we propose achieves high robustness to attacks from individual peers or from collusive ones, and provides incentive for participation. This is accomplished by a n aggregation technique in which a bounded number of peers randomly selected within the system report directly observed information to requesting peers. Observations are weighted by a credibility factor locally computed. Incentive for participation is implemented through a fair di erential service mechanism. It relies on peer's level of participation, a measure of peers' contribution over a xed period of time, and on the credibility factor, assessing the con dence one has in a peer.

Our results are promising: We prove that through su cient and honest cooperation, peers increase the quality of their reputation mechanism. We show that the reputation estimation e ciently lters out malicious behaviors in an adaptive w ay. Presence of a high fraction of malicious peers does not prevent a correct peer from computing an accurate PI n 1816

Anceaume & Ravoaja reputation value, at the expense of a reasonable convergence time. Furthermore, the tradeo between the sensitivity of the mechanism facing up malicious peers and the duration of the computation is tuned through a single input parameter. These properties, combined with the incentive scheme, makes our mechanism adapted to P2P networks. Finally, we provide a full theoretical evaluation of our solution.

The rest of the paper is organized as follows. In Section 2 related work is reviewed. Section 3 presents the model of the environment, and the speci cation of the robust reputation problem. Section 4 presents the incentive-based mechanism. Section 5 analyses its asymptotic behavior, its resistance to undesirable behavior and its convergence time.

Related Work

There is a rapidly growing literature on the theory and applications of reputation systems, and several surveys o er a large analyze of the state of art in reputation systems 14, 1 0 , 7 . According to the way ratings are propagated among entities and the extent of knowledge needed to perform the needed computations, reputation systems fall into two classes, namely centralized or distributed. An increasing number of online communities applications incorporating reputation mechanisms based on centralized databases has recently emerged. The eBay rating system used to nd traders allows partners to rate each other after completion of an auction. Despite its primitive reputation system, ebay is the largest person-to-person online auction with more than 4 millions auctions open at a time 15 . Regardless of this success, centralized approaches see for example, 20, 17 often pay little attention to misbehaving entities by assuming that entities give honest feedback to the requesting entity. More importantly, they rarely address non-participation and collusive behaviors.

Regarding decentralized p2p architecture, several research studies on reputation-based P2P systems have emerged. Among the rst ones, Aberer and Despotovic 1 propose a reputation mechanism in which trust information is stored in P-Grid, a distributed hash table-based DHT overlay. Their mechanism is made robust by guaranteeing that trust information is replicated at di erent peers, and thus can be accessed despite malicious entities. However, the e ciency of their approach relies on peers propensity to fully cooperate by forwarding requests to feed the P-Grid overlay. Additionally, as for most of the DHT-based approaches, peers have t o store data they are not concerned with. Thus, malicious peers may discard it to save private resources, leading to a loss of information. Other systems relying on the trust transitivity approach face false ratings by assuming the presence of speci c faithful and trustworthy peers e.g. 12 , or by w eighting second-hand ratings by senders' credibility 6 , 1 8 , 1 9 . Opposed to the aforementioned works, Havelaar reputation system 9 , exploits long-lived peers by propagating reports between sets of well de ned peers identi ed through hash functions. A report contains the observations made during the current round, the aggregated observations made by the predecessors during the previous round, and so on for the last r rounds. By relying on such an extensive aggregation, false reports hardly in uence the overall outcome. Furthermore by using hash functions collusion is mostly prevented. The e ciency of their approach mainly relies on the readiness of peers Irisa to store and propagate large amount of data, and to remain in the system for relatively long periods of time. To m o t i v ate peers to participate, Jurca and Faltings 11 propose an incentive-compatible mechanism by i n troducing payment for reputation. A set of brokers, the R-agents, buy and sell feedback information. An entity can receive a p a yment only if the next entity reports the same result. Weakness of such an approach is the centralization of the whole information at R-agents, and its robustness against malicious R-agents. Finally, Awerbuch et al. 4, 5 give l o wer bounds on the costs of the probes made by honest peers to nd good objects in eBay-like systems, and propose algorithms that nearly attain these bounds.

In contrast to these works, we propose a fully distributed mechanism based on local knowledge that provides malicious and non-participating entities an incentive for participation and honest behavior.

Model

A P2P Service Model

We consider a P2P service system where service p r oviders or servers repeatedly o er the same service to interested peers. We assume that the characteristics of a server capabilities, willingness to o er resources, etc are aggregated into a single parameter θ called type. This type in uences the e ort exerted by the server through a cost function c. The e ort determines the Quality o f the Service QoS provided by the server. We assume that the e ort exerted by a server is the same for all the peers that solicit him and takes its value within the interval [0, 1].

De nition 1 E ort. The e ort of a service p r ovider s is a value q * s that determines the quality of the service o ered t o t h e p eers that interact with s.

After each i n teraction with server s, each client or peer has an imperfect observation of the e ort exerted by s. Peers may h a ve di erent tastes about a server QOS. But basically these observations are closely distributed around s's e ort. Thus, we reasonably assume that an observed quality of service takes its value within the interval [0, 1] and follows a normal distribution of mean q * s and variance σ * s .

De nition 2 Observed Quality of Service Level. The Observed Quality of Service Level of a service provider s observed by peer p at time t is a value obs s p (t) which is drawn from a normal distribution over [0, 1] with mean q * s . The value 1 resp. 0 characterizes the maximal resp. minimal satisfaction of p.

Estimation of the expected behavior of a server is based on its recent past behavior, that is, its recent i n teractions with the peers of the system. Such a restriction is motivated by game theoretic results and empirical studies on ebay that show that only recent ratings are meaningful 7 . Thus, in the following, only interactions that occur within a sliding window PI n 1816 of width D1 are considered. This approach is also consistent with Shapiro's work 16 i n which i t i s p r o ven that in an environment where peers can change their e ort over time the e ciency of a reputation mechanism is maximized by giving higher weights on recent ratings and discounting older ratings. Using a sliding time window is approximately equivalent t o this model. Every time a peer p desires to interact with a server s, p asks for feedback from peers that may h a ve directly interacted with s during the last D time units. We adopt the terminology witness to denote a peer solicited for providing its feedback. If P s p (t) represents the set of peers k whose feedback has been received by p e e r p by time t, then the reputation value of a server is de ned as follows:

De nition 3 Reputation Value. The reputation value r s p (t) of server s computed by peer p at time t is an estimation of the e ort q * s exerted b y s based on the feedbacks provided by the peers in P s p (t).

Speci cation of Undesirable Behaviors

In practice, peers may not always reveal their real ratings about other peers. They can either exaggerate their ratings by increasing or decreasing them, or they can simply reveal outright ratings to maximize their welfare. This behavior is usually called malicious, and can either be exhibited by a node independently from the behavior of other peers, or be emergent of the behavior of a whole group. By providing false ratings, malicious peers usually try to skew the reputation value of a server to a value which is di erent from its true e ort. Let q be this value, and d be the distance between the true e ort of the server and the false rating d = |q * s -q|. Then, we c haracterize the behavior of a peer by w s q = 1d α , with α a positive real value which represents the sensitivity o f w s q to the distance between the e ort and the expected observation given by q. A malicious peer tries to skew the reputation value to q by sending ratings that are distributed around q.

De nition 4 Malicious. A p eer p is called malicious if it lies on the reputation of peer s or creates s's reputation out of thin air. Formally :

E(obs s p (t)) = q = q * s , ∀t.
De nition 5 Collusive Group. A group of peers is called a collusive group if all the peers of this group behave maliciously towards a same goal. Formally, the set C i s a c olluding gro u p i f :

E(obs s k (t)) = q = q * s , ∀t, k ∈ C.
Another common behavior in P2P systems is peers non-participation. There are two main reasons why peers may not participate: either because they believe that their work is redundant with what the others in the group can do, and thus their participation can hardly in uence the group's outcome, or because they believe that by not contributing they maximize their own welfare note that information retention could be another pretence of not participating, however this is out of the scope of the paper. The latter behavior depicts what is typically called free-riding, while the rst one is described in the Collective E ort Model CEM as social loa ng 13 . Note that although e ects of both behaviors are similar, i.e., non-participation, their deep cause is di erent. Peers exhibiting one of these two behaviors are called in the following non-participating peers and are characterized as follows:

De nition 6 Non Participating. A peer is called non participating if it exerts less e ort on a collective task than it does on a comparable individual task or consumes more than its fair share o f c ommon resources.

A peer is called correct if during the time it is operational in the system it is neither malicious nor non-participating. Note that a malicious peer may not participate, on the other hand, a non participating one is not malicious.

Speci cation of the Robust Reputation Problem

Within this context, we address the problem of evaluating the reputation of a service provider in a dynamic environment in which peers are not necessary correct. This problem is referred in the sequel as the robust reputation problem. A solution to this problem should guarantee the following two properties. The rst one states that eventually correct peers should be able to estimate the reputation value of a target server with a good precision. The second one says that with high probability, correct peers have a better estimation of the reputation value of a target server than non correct ones. Formally: Property 1 Reputation Value -Accuracy. Eventually, the reputation of server s, evaluated by any correct peer re ects s's behavior with precision . That is, let β ∈]0, 1[ be some xed r eal, called in the sequel con dence level, then:

∃t s.t. ∀t ≥ t, Prob(|r s p (t) -q * s | ≤ ) ≥ 1 -β Let |E(r s p (t)) -q *
s | be the bias of the reputation value r s p (t) estimated by peer p. Suppose that two peers p and q interact with the same target servers at the same time, solicit the same witnesses, and get the same feedbacks at the same time. That is from the point of view of their interaction p and q are indistinguishable. However p is correct while q is not. Then, we h a ve: Property 2 Incentive-Compatibility. Eventually, the bias of the reputation value of server s estimated b y p is greater than or equal to the one estimated b y p eer q. That is, for a given level of con dence β, w e h a v e :

∃t s.t. ∀t ≥ t, Prob(|E(r s p (t)) -q * s | ≥ |E(r s q (t)) -q * s |) ≥ 1 -β

PI n 1816 4 The Reputation Mechanism

We propose a distributed reputation service which builds a social network among peers. Brie y, e v ery peer records the opinion about the late experiences it has had with a target server. Peers provide their information on request from peers willing to interact with that server. Providing a feedback based on direct observations also called rst-hand observations prevents the rumors phenomenon, i.e., the propagation of opinions about others, just because they have been heard from someone else 18 , however is better adapted to applications with relatively low c hurn. Upon receipt of "enough" feedback, the requesting peer aggregates them with its own observations if any to estimate the reputation of the target server, and provides this estimation to its application. Information is aggregated according to the trust the requesting peer has in the received feedback. Pseudo-code of the reputation mechanism is presented in Algorithm 1. The e ciency of the reputation mechanism fully depends on i the number of received feedbacks i.e., aggregating few feedbacks is not meaningful and thus not helpful, and ii the quality of each of them i.e., the trustworthiness of the feedback. The contribution of this work is the design of a reputation mechanism that enjoys both properties. The analysis presented in Section 5 shows the importance of each factor on the convergence time and accuracy of the reputation mechanism.

The solution we propose is a reputation mechanism, and therefore independent of the rewarding strategy used by the application built on top of this mechanism. That is, the willingness of a peer to interact with a server results from the application strategy, not from the peer's one. Clearly, the strategy of the application is greatly in uenced by the reputation value but other factors may also be taken into account.

Collecting Feedbacks

When a peer decides to evaluate the reputation value of a service provider, it asks rst-hand feedback from a set of witnesses in the network. Finding the right set of witnesses is a challenging problem since the reputation value depends on their feedback. Our approach for collecting feedbacks follows the spirit of the solution proposed by Y u et al 19 , in which feedbacks are collected by constructing chains of referrals through which peers help one another to nd witnesses. We adopt the walking principle. However, to minimize the ability of peers to collude, witnesses are randomly chosen within the system. We assume, in the following, that the network is regular. Speci cally, our approach is based on a random walk RW sampling technique. We use the random walk technique as shown in Algorithm 1. Function query is invoked by the requesting peer that wishes to solicit x witnesses through r random walks bounded by ttl steps. The requesting peer starts the random walks at a subset of its neighbors, and runs them for ttl steps. Each peer p involved in the walk is designated as witness, and as such sends back to the requesting peer its feedback.

When a peer q receives a request from p to rate server s, i t c hecks whether during the last sliding window of length D, it has ever interacted with s. In the a rmative, p sets its feedback t o F s p (t) = obs s p (t 0 ), t 0 ,. Because of non-participation volunteer or because of a crash, random walks may fail: it su ces that one of the peers in the walk refuses to participate or crashes to prevent the walk from successfully ending. If we assume that in the neighborhoodofany peer, a fraction µ of the peers do not participate, then among the r initial peers that start a random walk, the expected number of peer that may "fail" their random walk is µr. Then during the next step, µ(1µ)r walks may "fail", and so on until the TTL value is reached. In consequence, only x feedbacks may be received, with x = ttl t=1 (1µ) t r. By setting r to

x ttl t=1 (1-µ) t
the requesting peer is guaranteed to receive a t l e a s t x feedbacks see line 5 in Algorithm 1.

In addition to its feedback, each peer sends to the requesting peer p through the witness message, see lines 19 and 46 the identity of the next potential witness on the walk, i.e., the peer it has randomly chosen among its neighbors. Sending this piece of information allows p to know the identity of all potential witnesses. As will be shown in Section 4.4, this allows to detect non participants if any and then to motivate them to participate by applying a tit-for-tat strategy. As for feedbacks, non participation may prevent the requesting peer from receiving witness messages. A similar analysis to the preceding one shows that if r random walks are initiated then y = ttl-1 t=0 (1µ) t r witness messages will be received. Note that a requesting peer can adapt its collect policy according to its knowledge of the target server, or of its neighborhood. Speci cally, to get x witnesses, a peer can either increase ttl and restrict r, or increase r and lower ttl. Enlarging ttl would be more sensitive to colluding peers that bias the random walk. However, this technique would increase the set of crawled witnesses, and thus would a ord new peers the opportunity to be known by other peers and consequently to increase both their participation and their credibility. Conversely, enlarging r would crawl only peers in the neighborhood of the requesting peer. However, this technique would increase the chance to nd a path that does not contain colluding peers 2 .

Reputation of a Server

Estimation of the reputation value of a target server is based on the QoS directly observed at the server if any and on the feedbacks received during the collect phase. The accuracy of the estimation depends on the way these informations are aggregated. The aggregation function we propose answers the following qualitative and quantitative preoccupations: First, to minimize the negative in uence of unreliable information, feedbacks are weighted by the credibility of their senders. Brie y, credibility i s e v aluated according to the past behavior of peers and re ects the con dence a peer has in the received feedback. Credibility computation is presented in the next subsection. Second, to prevent malicious nodes from ooding p with 2 Remark that selecting peers according to their credibility should be more e cient in the sense that only highly credible peers would be selected, however, newcomers may be penalized by this ltering. Furthermore, the resilience of the crawling technique to collusion highly relies on the way the graph of witnesses is constructed. Studying these issues is part of our future work.

PI n 1816 [START_REF] Josang | A survey of trust and reputation systems for online service provision[END_REF] Anceaume & Ravoaja fake feedback and thus from largely impacting the accuracy of its estimation, p keeps only a subset of each received feedback. More precisely, among the set of observations sent b y each witness over the last D time units, only the last f ones are kept, with f t h e s i z e o f t h e smallest non-empty set of non-default feedbacks received by p i.e., f = min k∈P s p (t) (

r s p (t) = 1 k∈P s p (t) c s p,k (t) k∈P s p (t) c s p,k (t).ρ k (t) if f = ∅ obs max otherwise 1 with, ρ s k (t) = 1 f (obs s k (t ),t )∈F s k (t) obs s k (t )

Trust in Witnesses

In this section, we tackle the issue of malicious peers. As remarked in the Introduction, malicious peers may alter the e ciency of the reputation mechanism by sending feedbacks that over-estimate or sub-estimate the observed QoS of a server to in ate or tarnish its reputation. This is all the more true in case of collusion. We t a c kle this issue by e v aluating peers credibility. Credibility is a 0,1 -valued function which represents the con dence formed by p e e r p about the truthfulness of q's ratings. This function is local and is evaluated on the recent past behavior of both p and q peers. It is locally used to prevent a false credibility from being propagated within the network. Speci cally, peer p estimates at time t how credible q is regarding server s as a decreasing function of the distance between q's feedbacks on s's e ort and p's direct observations on s's QoS. As for the reputation value computation, the distance is computed on the last f observations made by both p and q during the last D time units. Note that in case p has not recently observed s's QoS, then credibility of all its witnesses are set to a default value c 0 . Indeed, p cannot evaluate the distance between its own observations and those observed by witnesses. Determining c 0 value needs to solve the following trade-o : by a ecting a high value to the default credibility one increases the vulnerability of the system to the whitewashing phenomenon, that is, the fact that peers change their identity in order to reset their credibility to the default value. However, by Irisa setting this variable to a low v alue the mechanism tends to lter out new witnesses and thus, loses the bene t of the potential information a new peer can a ord, which clearly decreases the usefulness of the reputation mechanism. In order to cope with that, we set c 0 to the value of a decreasing function of φ, with φ an estimation of the number of whitewashers in the network. By adopting the notations of Equation 1, c s p,q (t) represents the credibility formed by p at time t about q regarding the target server s, and is given by:

c s p,q (t) = 1 -|ρ s q (t) -ρ s p (t)| α if f = ∅ c 0 otherwise 2
where |ρ s q (t)ρ s p (t)| α represents the distance between q and p's observations. Note that α is the variable introduced in Section 3. Then we h a ve the following lemma: Lemma 1. Credibility Accuracy Eventually, credibility of a peer q evaluated by any correct peer p re ects q's behavior with a precision . That is, let β ∈]0, 1[ be s o m e x e d r eal, there exists t such that, for all t ≥ t,

Prob(|c s p,q (t) -w s q (t)| ≤ ) ≥ 1 -β.
Proof. sketch Suppose rst that witness q is correct. Then by de nition w q = 1. As q truthfully reports the QoS of s it has observed during the last interval of length D, then obs s q (t)

follows a normal law o f m e a n q * s . By the law of large numbers,ρ s

q (t) = 1 |θ s q (t)| t ∈θ s q (t) obs s q (t )
eventually converges to q * s . By assumption, p is correct. Thus eventually, ρ s p (t) converges to q * s . By Equation 2, c s p,q (t) converges to 1 which concludes this case. Suppose now that peer q is malicious. By de nition see Section 3, its observations are distributed around q. Peer q's behavior is characterized by w q = 1-|q * s -q| α . Then, c s p,q (t) = 1-|ρ s q (t)ρ s p (t)| α converges to w q . Thus, for all p, there exists t such that, for all t ≥ t, Prob(|c s p,q (t)-w s q (t)| ≤ ) ≥ 1-β.

Incentive for Participation

Non participation may jeopardize the e ciency of the reputation mechanism. A certain amount of participation is required before reputation can induce a signi cant l e v el of cooperation. Facing non-participation in the reputation problem is challenging and has deserved few attention 19 . To m o t i v ate peers to send their feedback w e adopt a tit-for-tat strategy. We i n troduce the level of participation notion as the propensity of a peer for replying to a rating request. It is described by function l s p,q such t h a t l s p,q (t) represents the percentage of times q provided its feedback t o p's queries regarding server s's QoS over the last D time units, with l s p,q (t = 0) = l 0 = 1. Its computation is performed after p's collect phase see line 10 of the algorithm. Note that factor µ prevents correct peers from being penalized by walking breaks.. We apply the tit-for-tat strategy during the collect phase. When a peer p receives a rating request for s from peer q, then with probability l s p,q (t) p provides its feedback t o q, otherwise it sends a default feedback ⊥, ⊥ t o p r e v ent p from being tagged as non-participant. By PI n 1816 providing this default feedback, p lets q knows that its recent non-participation has been detected. Consequently, b y not participating, requesting peers drive correct witnesses providing them worthless feedback, which clearly makes their reputation mechanism useless. Hence there is a clear incentive for non participating peers to change their behavior. The following lemma proves that participation decreases the bias of the reputation value. As previously, let us consider two peers p and q such that both peers are indistinguishable from the point of view of the servers with which t h e y i n teract, that is both p and q observe t h e same QoS from these servers at the very same time, solicit and are solicited by the same set of peers at the same time; However, p is correct while q is non participating. Then we claim that: Lemma 2. Participation decreases the bias of the reputation value. That is,

|E(r s p (t)) -q * s | ≤ |E(r s q (t)) -q * s |
Proof. sketch Let us consider some correct peer k, s u c h t h a t k solicits both p and q during the same collect phase. By construction, l s k,p (t) = 1 since p is correct. On the other hand, l s k,q (t) < 1 by assumption. Suppose now that both p and q solicit peer k. Then by construction, k provides its rating to p with probability 1 while it provides it to q with probability 1l s k,q (t). By assumption both p and q have the same history in terms of solicitations, thus q cannot receive more ratings than p. From above, p receives k's rating with probability 1, and k is correct. Thus the quality of the feedbacks received by p is at l e a s t a s g o o d a s q's one. Which completes the proof.

Incentive for Truthful Feedbacks

We n o w address the problem of motivating peers to send truthful feedbacks. So far we h a ve presented strategies aiming at improving the quality of the reputation value estimation by aggregating more feedbacks and by w eighting feedback according to the credibility of their sender. We have shown that by using both strategies, utility of correct peers increases. However, none of these solutions have an impact on the e ort devoted by a witness to send a truthful feedback. To tackle this issue we use the credibility as a way to di erentiate honest peers from malicious ones. As for non-participating peers, when peer p receives a request to rate server s from a requesting peer q then p satis es q's request with probability c s p,q (t). By doing so, p satis es q's request if it estimates that q is trustworthy, otherwise it noti es q of its recent faulty behavior by sending it the (⊥, ⊥) feedback. As previously, b y cheating, a malicious peer penalizes itself by pushing correct witnesses to send meaningless feedbacks to it, leading to its e ective isolation. We claim that this social exclusion-based strategy motivates q to reliably cooperate. Lemma 3. High credibility decreases the bias of the reputation value. That is,

|E(r s p (t)) -q * s | ≤ |E(r s q (t)) -q * s |
Proof. Similar to the one of Lemma 2, by replacing l p,q (t) with c s p,q (t).

Irisa

Incentive-Based R obust Reputation Mechanism for P2P Services 13 Finally, to elicit su cient and honest participation, both strategies are combined, i.e., upon receipt of a rating request from peer q, with probability min(c s p,q (t), l s p,q (t)) p provides its feedback, otherwise it sends the default feedback (⊥, ⊥) see line 31 in Algorithm 1.

Theorem 1. The reputation mechanism described in Algorithm 1 is Incentive-Compatible in the sense of Property 2. Proof. sketch Given a non-participating peer, the result follows directly from Lemma 2. Let us consider a malicious peer p and a correct one q. Then, from Lemma 1, there exists a time t such that with high probability p's credibility i s -far from 1d α while q's credibility is -far from 1. Then by applying Lemma 3, we conclude the proof.

Analysis

Computing the reputation of a peer reduces to estimating, in the statistical sense, its effort. Our algorithm falls into the category of robust estimation algorithms. Indeed, robust estimation techniques consider populations where a non-negligible subset of data deliberately pollute the system. This analysis describes the asymptotic behavior of the reputation mechanism and its convergence time according to undesirable behaviors. In the following, we assume that a fraction γ of witnesses are malicious.

Asymptotic behavior

In this section we determine the accuracy of the reputation mechanism with respect to parameters α and . Thus for the purpose of this analysis, we assume that the number of aggregated feedbacks f is in nite. Recall that w q denote the characterization of q's behavior, with w q = 1 if q is correct, and w q = 1d α otherwise. By Lemma 1, c s p,q (t) = w q , when t → ∞. Moreover, the expected number of correct witnesses is (1γ)n + 1 while the expected number of malicious ones is γn, with n the expected number of witnesses Figures are plotted for n = 10. Thus, by replacing c s p,q (t) with their asymptotic values in Equation 1, the expected reputation value of a server s estimated by a correct peer p when t → ∞ is given by Equation 3:

r s p (t) t→∞ = 1 1 -γd α + 1 n 1 -γ + 1 n q * s + γ (1 -d α ) q .
3

The bias of the reputation value, when t → ∞ is given by the following equation: as shown by the crossed curves in Figure 1, our algorithm bounds the power of colluders whatever their percentage dotted curves. Indeed, witnesses' ratings are weighted by a decreasing function of d which lters out false ratings. Figure 2 shows the impact of α on the bias of the reputation value. As can be observed, the bias decreases with decreasing values of α, re ecting the sensitivity of the reputation value to the distance between direct observations and received feedbacks. Thus, decreasing values of α makes the reputation mechanism very sensitive to false feedbacks.

|r s p (t) t→∞ -q * s | = γd 1 -d α 1 -γd α + 1 n . Figures 1 a n
Theorem 2. The reputation value is -accurate, with > γd 1-d α 1-γd α + 1 n . Proof. The bias of the reputation estimation converges to γd 1-d α 1-γd α + 1 n see Equation 4. Then, given β and in ]0, 1[, there exists t such that, Prob(||r

s p (t) t→∞ -q * s |-γd 1-d α 1-γd α + 1 n | ≤ ) ≥ 1 -β. Thus, Prob(|r s p (t) -q * s | ≤ γd 1-d α 1-γd α + 1 n + ) ≥ 1 -β.
We conclude by setting

= γd 1-d α 1-γd α + 1 n + .
From the above, assuming an upper bound on γ, b y setting the maximal bias to and solving the corresponding equation one can derive a n upper bound on α under which the reputation value converges to the true e ort with an accuracy level of . Hence, for α with α ≤ α, the reputation value converges to the true e ort with an accuracy level of , with α given by:

α = ln n+1- √ (1-γ)n 2 +(2-γ)n+1 γn ln( ) 5
To conclude, one can always nd a value of α such that eventually the reputation value is accurate. This parameter, however, signi cantly in uences the convergence time of the algorithm. The next Section addresses this issue.

Convergence

In this section, we study the convergence time of the reputation mechanism. To d o s o , w e assume that the ratings of a malicious peer are drawn from a normal distribution with mean Irisa q and variance σ over [0, 1]. This assumption includes a wide range of possible behaviors. Indeed, a small value of σ depicts peers that try to rapidly skew the reputation value to q by giving reports tightly distributed around q. In contrast, a high value of σ depicts peers that try to hide their mischievious behavior to other peers by giving sparse reports. While the rst behavior is easily detected, the second one hardly skew the reputation value to q.

Recall that the reputation value is estimated by aggregating the last f interactions witnesses have had with the target server during a sliding time window of length D. Finding the optimal value of D is important. Indeed, it determines the resilience of the mechanism to e ort changes and the con dence level in the estimation. The optimal value of D is the one for which the estimation is at most -far from the true e ort with a g i v en con dence threshold β. To determine such a v alue, let us rst assume that f is known, and determine a lower bound on Prob(|r s p (t)q * s | ≤ ) ∀t ∈ D. Suppose that the credibility c s k (t) is -far from w s k for all the witnesses k. Then, because of Bayes' Theorem, we know that

Prob(|r s p (t) -q * s | ≤ ) ≥ Prob(|r s p (t) -q * s | ≤ | |c s k (t)-w s k (t)|≤ ,∀k∈P s p (t) ) • Prob(|c s k (t) -w s k (t)| ≤ , ∀k ∈ P s p (t)).
Remark that, assuming that the witnesses' credibility are -far from w s k , the probability that r s p (t) is -far from q * s is maximal under the following condition C: credibility of correct witnesses is minimal, i.e., equal to 1 -, and the one of malicious ones is maximal, i.e., equal to . Then, we h a ve:

Prob(|r s p (t) -q * s | ≤ ) ≥ Prob(|r s p (t) -q * s | ≤ | (C) ) • Prob(c s k (t) ≥ 1 -) (1-γ)n • Prob(c s k (t) ≤ ) γn
6 By Lemma 1, the probability that the witness's credibility is at most -far from w s k converges to 1 when t, and thus f , increase. Thus, this bound approaches Prob(|r s p (t)q * s | ≤ ) when f increases. Knowing the probability distribution of the reports, deriving a closed form of the lower bound can be done. Then, given a desired con dence threshold β for α ≥ α, solution of Equation 7provides two threshold values of d d1 and d2 on Figure 2 beyond which the false reports are eliminated:

= γd 1 -d α 1 -γd α + 1 n within [0, 1]. 7
Using Equation 6, we can see the e ect of α on the convergence time of the reputation mechanism to reach the exact e ort exerted by the server see Figure 3. Decreasing values of α signi cantly increases the convergence time while it decreases the bias of the reputation mechanism as shown in Figure 2. Thus a trade-o exists between the robustness of the mechanism and its convergence time. By setting the value of D, the application designer may d e r i v e the corresponding minimal numberofinteractions f and nally tune the value of α such that the desired con dence level is achieved within f steps through Equation 6. The resulting reputation estimation is less sensitive to false reports, but still eliminates peers that try to skew the reputation of the server to a value that is far from the true e ort, by ltering out extreme values of d see Figure 2. Finally, Figure 4 shows the impact of malicious peers on the convergence time. As can be seen, a relatively small percentage of malicious peers has a minor impact on the convergence time since the number of correct feedbacks is hardly in uenced by false ones. On the other hand, whenever a requesting peer has to face a large proportion of malicious peers, it can only rely on its own feedback to estimate the e ort exerted by the target server which clearly takes longer than when helped by correct witnesses. The same result applies for non-participating peers.

Conclusions

In this paper we have proposed a reputation mechanism that achieves high robustness to attacks and provides incentive for participation. This is achieved by an aggregation function in which a subset of the information provided by randomly chosen peers is kept and weighted by a con dence factor locally computed. We have proposed a simple and local incentive mechanism that guarantees a better quality of the reputation value estimation. Lessons learned from simulations are twofold: rst, decreasing values of α guarantees a greater sensibility of the mechanism to false ratings. It however increases the number of required feedbacks as well, and thus the time to get an accurate estimation of the e ort exerted by the target server. Second, the presence of a large number of malicious and non-participating peers does not prevent the mechanism from being accurate, however has an impact on its convergence time.

  A T O I R E S P U B L I C A T I O N I N T E R N E N o

Figure 1 :Figure 2 :

 12 Figures 1 a n d 2 show the bias of the reputation value with respect to d for increasing values of γ resp. increasing values of α. Recall that d = |q * s -q| re ects colluders' behavior. Unlike mean-based reputation value estimation in which the bias linearly increases with d

Figure 3 :

 3 Figure 3: Con dence level for γ = 5, n = 10, = .1, d = 1, σ * s = .2, and σ = .2.

Figure 4 :

 4 Figure 4: Con dence level for α = 3, n = 10, = .1, d = 1, σ * s = .2, a n d σ = .2.

  . ., obs s p (t l ), t l with obs s p (t i ) the QoS of s observed at time t i , where t i ∈ [max(0, t -D), t]. In case p has not recently interacted with s, p sends Irisa back t o q a default feedback F s p (t) = {(obs max , ⊥)}. As will be shown later, this feedback prevents p from being tagged non participant by q.

  Finally, if among all the witnesses including p none has recently directly interacted with s i.e., f = 0, then p a ects a maximal value obs max to s's reputation value. A ecting a maximal value re ects the key concept of the Dempster-Shafer theory of evidence which argues that "there is no causal relationship betwe e n a h ypothesis and its negation, so lack o f belief does not imply disbelief". In our context, applying this principle amounts in xing an an priori high reputation to unknown servers, and then updating the judgment according to subsequent i n teractions and observations 19 .We can now i n tegrate these principles within the aggregation function we propose. Let us rst introduce some notations: Let F s k (t) be the union of the last f non-default feedbacks received from k during the last D time units t ∈ [max(0, t -D), t]; P s p (t) be the set of witnesses k for which F s k (t) is non empty; ρ s k (t) represent the mean value of the observations drawn from F s k (t); and c s p,k (t) the credibility formed by p at time t about k regarding s. Then, at time t, p estimates s's reputation value as follows:

	with t ∈ [max(0, t -D), t].	|F s k (t)|)

  : p: requesting peer; s: target server; x: expected number of feedbacks output: r s p : estimation of s reputation value ttl ← default value; t ← getTime; let D be the the time interval maxt-D,0,t ; (t) messages are r eceived f r om x peers and witness(s, w, t) messages are r eceived New ← pick a random subset of r peers from p's neighbors; has interacted w i t h s at time t 0 , . . . , t l in the last D time units then pick o n e o f q's neighbor with probability 1 d , d = q's neighbors ; 43 send a rw p,s,ttl -1,t message to next; Estimation of the reputation value of server s by peer

	1 2	r ←	x (1-µ) l ; y ← ttl-1 l=0 (1 -µ) l r; l=1 ttl
	queryp,s,ttl,r,t; wait until F s q 5 5 3 from y peers; P s (t) ← q such that a feedback message is received from q ; 6 F s (t) ← F s q (t) for all q ∈ P (t) ; 7 W s (t) ← w for all witness(s, w, t) that have been received ; 8
	10 10	foreach k ∈ W (t) do
	11	l s p,k (t) ← min( t∈D	(|P |k (t)|
	12 13	end return r s p (t) to application;
	14	queryp,s,ttl,r,t begin
	22
	23		else
	24		F s p (t) ← {(obsmax, ⊥), p};
	25 26		end if send F s p (t) to p;
	27	end
	28	end
	upon receipt of a rw p,s,ttl,t message at peer q do with probability min(l s p,q (t), c s p,q (t) do if 33 33 29 31 31
	34		F s q (t) ← {(obs s
	35		else
	36		F s q (t) ← {(obsmax, ⊥), q};
	37		end if
	38	otherwise
	39		F s q (t) ← {(⊥, ⊥), q};
	40 41 42	end do send F s q (t) to p; if ttl = 0 then
	next ← 44 send a witness s,next,t t o p; 46 46
	47	end if
	48	end do Algorithm 1:

input t∈D |W |k (t)| (1µ) + µ.l 0 , 1);

15 forall next ∈ New do 16 send a rw p, s, ttl -1,t message to next; 17 send a witness s,next,t message to p; 19 19 if p has interacted with s at time t 0 , . . . , t l in the last D time units then 21 21

F s p (t) ← {(obs s p (t 0 ), t 0 ), .

. . , (obs s p (t l ), t l ), p}; q q (t 0 ), t 0 ), . . . , (obs s q (t l ), t l ), q}; p PI n 1816

D can have a n y pre-de ned length of time, i.e., a day, a w eek or a month. In the sequel, we suppose that D is insensitive to clock drift.Irisa

Anceaume & Ravoaja

Aina Ravoaja is partially supported by a grant from Brittany Region


Nicolas Desassis

Didier Renard

Hélène Beucher

Sylvain Petiteau

X Freulon

A pairwise likelihood

Keywords: Plurigaussian models, empirical variography, semiparametric estimation, pairwise likelihood (PL), underlying Gaussian Random Fields (GRFs) 1

Introduction

Regionalized categorical variables often appear in several scientific domains. For instance, in the earth sciences, some continuous soil properties (e.g the permeability, the grade of an element, ...) can be better described by first categorizing the rock types into lithofacies (or facies) which present a certain homogeneity with respect to the studied variable. Then, the continuous variables are studied separately in each category. In the scope of conditional simulations, the lithofacies are first simulated conditionally to the observed lithofacies, then the continuous variables are simulated inside each simulated category according to their associated spatial distribution. To model and simulate a categorical random field, the plurigaussian model is particularly appealing. Based upon the simple truncation of one [START_REF] Matheron | Conditional simulation of the geometry of fluvio deltaic reservoirs[END_REF] or several Gaussian Random Fields (Le Loc'h and [START_REF] Loc'h | Truncated plurigaussian method: theoretical and practical point of view[END_REF][START_REF] Loc'h | Improvement in the truncated gaussian method: combining several gaussian functions[END_REF], it allows to reproduce a wide range of patterns. Applications of the plurigaussian model can be found for mineral resources evaluation [START_REF] Carrasco | Application of the truncated gaussian simulation method to a porphyry copper deposit[END_REF][START_REF] Riquelme | Truncated gaussian and plurigaussian simulations of lithological units in mansa mina deposit[END_REF][START_REF] Talebi | Stochastic rock type modeling in a porphyry copper deposit and its application to copper grade evaluation[END_REF], in hydrology [START_REF] Mariethoz | Truncated plurigaussian simulations to characterize aquifer heterogeneity[END_REF]. In petroleum, some authors use the plurigaussian models in link with history matching [START_REF] Hu | Gradual deformation and iterative calibration of gaussian-related stochastic models[END_REF][START_REF] Liu | Automatic history matching of geologic facies[END_REF][START_REF] Romary | History matching of approximated lithofacies models under uncertainty[END_REF]. In this paper, we suppose that we directly observe the categorical variable.

When the underlying Gaussian Random Fields (GRFs) are supposed to be stationary, two ingredients are necessary to fully specify the plurigaussian model: the coding function (or truncation rule) which defines the sets associated to each category and which can vary along the space and the multivariate covariance function.

Concerning the coding function, some authors use a simple parametric form, for instance a cartesian product of intervals and they allow the threshold values to vary along the space. It is often the case vertically through the vertical proportion curves [START_REF] Felletti | Statistical modelling and validation of correlation in turbidites: an example from the tertiary piedmont basin (castagnola fm., north italy)[END_REF] but also laterally, for instance when using an auxiliary information as seismic data in the model. Other authors concentrate on the estimation of more complex coding functions constant in space [START_REF] Astrakova | Truncation map estimation for the truncated bigaussian model based on univariate unit-lag probabilities[END_REF]. Finally, [START_REF] Allard | Nonparametric diagrams for pluri-gaussian simulations of lithologies[END_REF] estimate complex coding function varying in space by using auxiliary information.

As mentionned by [START_REF] Mariethoz | Truncated plurigaussian simulations to characterize aquifer heterogeneity[END_REF], one of the main difficulty arising from the use of the plurigaussian model is the inference of the variogram models of the underlying GRFs. Indeed, the available empirical variograms are the variograms of the indicator functions of the categories (one simple variogram per category and the crossvariograms for all the bivariate combinations) while the variograms required by the model are the variograms of the underlying GRFs whose realizations are hidden by the truncation.

Until now, most of the methods to estimate the variogram of the underlying GRFs rely on the indicator variograms. For instance [START_REF] Mariethoz | Truncated plurigaussian simulations to characterize aquifer heterogeneity[END_REF] determine the variogram model of the underlying GRFs by using simulations. More precisely, they choose a parametric model for the underlying GRFs and they compute the parameters value such as the indicator simple variograms of the simulations are the closest to the data indicator variograms. The optimization is performed with simulated annealing. [START_REF] Armstrong | Plurigaussian simulation in geosciences[END_REF] exploit the mathematical relationships between the underlying GRFs variogram models, the coding function and the indicator simple and cross-variograms. Some industrial softwares (as Isatis R , 2014) also use these relations and the users must choose the parameters of the variogram models of the underlying GRFs by visual inspection of the resulting indicator variograms through a trial-and-error procedure [START_REF] Galli | The pros and cons of the truncated gaussian method[END_REF]. [START_REF] Emery | Simulation of geological domains using the plurigaussian model: New developments and computer programs[END_REF] performs the numerical integration of the Gaussian density by using its expansion into the normalized Hermite polynomials. All these methods are rather tedious as they have a high computational cost or require a lot of trials. [START_REF] Dowd | Plurigau: a computer program for simulating spatial facies using the truncated plurigaussian method[END_REF] and [START_REF] Xu | A flexible true plurigaussian code for spatial facies simulations[END_REF] propose to find the range parameters automatically by minimizing a squared differences with a grid-search but the choice of the covariance models of the underlying GRFs remains arbitrary and limited.

In this paper, we will supppose that the coding function is known and we will concentrate on the estimation of the variograms of the underlying GRFs. We propose an original methodology based on the pairwise likelihood (PL) maximization principle to directly compute the empirical variograms of the underlying and hidden GRFs. More precisely, we perform a semiparametric estimation by considering the variogram at a given lag (distance in the omnidirectional case or vector otherwise) as a parameter of the model. Then, we maximize the PL by grouping the pairs of points approximately separated by this lag in the same way as for the classical empirical variogram estimation [START_REF] Matheron | Traité de Géostatistique appliquée[END_REF]. We iterate this calculation on all distances (respectively vectors). Thereby, we obtain an empirical variogram which helps the user to choose a suitable valid model that can then be fitted by least squares or estimated with a likelihood based approach. Then, the simple and cross-variograms in the indicator scale can be derived and compared to the empirical variograms of the indicators to check the quality of the resulting models.

In the first part, we give the main notations of the paper and we recall the definition of the plurigaussian model. In section 2.2, the relationships between variograms of GRFs and variograms of indicators are recalled for comparison purposes. Then we present our method in section 3. First, we describe the general principle which should make the estimation of a complex multivariate spatial model possible. Then we describe with more details the implementation in the case where the underlying GRFs are considered as independent. To assess the efficiency of the method and to evaluate the uncertainty associated to the variogram estimation, a Monte-Carlo study is performed and its results are summarized in section 4. We finish with some perspectives offered by the method to simplify the inference of the hierarchical spatial models.

The data model 2.1 General formulation of the plurigaussian model

Let F = {f 1 , . . . , f K } be a finite set with K categories. For a set of n sites {x i } 1≤i≤n of a domain D ⊂ R d , we observe f = (f (x 1 ), . . . , f (x n )) a F-valued vector. We suppose that for a given location x ∈ D, the value f (x) is the realization of a F-valued random variable F (x).

To characterize the spatial distribution of F (.), we use the plurigaussian model as described in [START_REF] Armstrong | Plurigaussian simulation in geosciences[END_REF]. Let Y(.) = {Y(x), x ∈ D} be a q-variate centered and standardized GRF on D: for all x ∈ D, Y(x) = (Y 1 (x), . . . , Y q (x)) is a random vector with q components and for all N ∈ N and for all (x 1 , . . . , x N ) ∈ D N , the N × q-vector

(Y(x 1 ), . . . , Y(x N )) = (Y 1 (x 1 ), . . . , Y 1 (x N ), . . . , Y q (x 1 ), . . . , Y q (x N )) is a standard Gaussian vector with E[Y r (x i )] = 0 and Var[Y r (x i )] = 1 for all r ∈ [[1, q]] and i ∈ [[1, N ]].
In this paper, we supose that Y(.) is a second-order stationary multivariate function, i.e there exists a matricial cross-covariance function [START_REF] Wackernagel | Multivariate Geostatistics -An Introduction with Application, 3rd Edition[END_REF], for an introduction on multivariate spatial random fields).

C such as Cov(Y r (x), Y s (x )) = C rs (x -x) for (r, s) ∈ [[1, q]] 2 (see
Let C be a coding function on D such as, for all

x ∈ D, C(x) = (C 1 (x), . . . , C K (x)) where, for k ∈ [[1, K]], the subsets C k (x) form a (measurable) partition of R q .
The model is defined by the following equivalence

F (x) = f k if and only if Y(x) ∈ C k (x).
(1)

Note that the formulation given by ( 1) provides a quite general class of models. Indeed, it also contains the models defined by

F (x) = f k if and only if ϕ(Y(x)) ∈ Ck (x)
for any surjective function ϕ from R q to any set E where the sets Ck (x) for [[1, K]] form a partition of E. The subsets ϕ -1 ( Ck (x)) have to be some measurable sets of R q . This remark aims to highlight the fact that the marginal gaussianity of the random variables Y r (x) is arbitrary. Nevertheless, the Gaussian assumption is a convenient way to describe the spatial multivariate relationships of the underlying random fields. It also provides a multivariate random field easy to simulate (see e.g [START_REF] Lantuejoul | Geostatistical Simulation: models and algorithms[END_REF]).

We will note c(x) the set defined as:

c(x) = C k (x)
where k ∈ [[1, K]] is the index of the category at location x. In other words, f (x) = f k . In all the sequel, we will suppose that the classes C k (x) are known.

Indicators cross-variograms based methods

As mentionned in the introduction, most of the methods to choose the covariance model of the underlying GRFs rely on the indicator simple and cross-variograms or covariances. [START_REF] Armstrong | Plurigaussian simulation in geosciences[END_REF] or Isatis R (2014) use the mathematical relationships between the simple and cross-covariances (or variograms) of the GRFs and the simple and cross-covariances (or variograms) of the indicators of each category. In the current paper, we only use these relationships to check the quality of the results given by our proposed method. We recall these relationships below. For that purpose, we will note the random indicator function of the category f k ∈ F as follows:

1 1 f k (x) = 1 if F (x) = f k 0 otherwise
and 1 f k (x) the associated true value.

Variogram between two points

For (k, l) ∈ [[1, K]] 2
, one can define the variograms between indicators of facies k and l, between two locations x and x of D by:

γ kl (x, x ) = 1 2 E[(1 1 f k (x ) -1 1 f k (x))(1 1 f l (x ) -1 1 f l (x))]
When k = l, we have the simple variogram:

γ kk (x, x ) = E[1 1 f k (x)] + E[1 1 f k (x )] 2 -E[1 1 f k (x )1 1 f k (x)] (2)
When k = l, we have the cross-variogram:

γ kl (x, x ) = - E[1 1 f k (x )1 1 f l (x)] + E[1 1 f l (x )1 1 f k (x)] 2 (3)
We note Σ x and Σ x,x the respective correlation matrices of the vectors Y(x) and (Y(x), Y(x )). Furthermore, g Σ (u) stands for the centered and standardized Gaussian density with correlation matrix Σ computed for the vector u.

With these notations, we can establish the link between γ kl (x, x ) and the correlations between the underlying GRFs. Indeed, the expectation of the indicator of facies k (which corresponds to its proportion at location x) is equal to:

E[1 1 f k (x)] = C k (x)
g Σx (u)du ( 4)

and

E[1 1 f k (x)1 1 f l (x )] = C k (x) C l (x ) g Σ x,x ((u, v))dudv (5)
where each integration symbol represents an integration over a qdimensional space. These integrations and all the others mentionned in the current paper are integral of the Gaussian probability density function. They can be computed numerically with the efficient algorithm proposed by [START_REF] Genz | Numerical computation of multivariate normal probabilities[END_REF].

Note that it is sometimes useful to work with the non-centered covariances E[1 1 f k (x)1 1 f l (x )] which can be computed in the same way. Indeed, it has the advantage to capture asymetry in the model.

When the GRFs are stationary and the coding function C is constant over D, the variograms of all the involved random fields only depend on the lag between the points. Therefore, we can deduce the simple and cross-variograms of the indicator for a given lag from the variograms value of the underlying GRFs by using formulas (2), ( 3), ( 4) and ( 5). However, when the coding function varies over D, the theoretical simple and cross-variograms of the indicators for a given lag do not exist anymore. Nevertheless, it is still possible to compute the associated empirical variograms and compare them with an averaged version of the variograms between two points computed in the indicators domain as described below.

Variogram for a specific lag

For a given vector h ∈ R d , we will note x j -x i h when the pair (x i , x j ) is used to compute the empirical variogram at lag h (see e.g Chilès and Delfiner, 2012, for details). N (h) stands for the number of such pairs. γkl (h) = 1 2N (h)

x j -x i h (1 f k (x j ) -1 f k (x i ))(1 f k (x j ) -1 f k (x i )).
This is the quantity that many authors suggest to fit with the image of the variogram model of the GRFs in the indicators scale which is defined for (k, l) ∈ [[1, K]] 2 by:

γ kl (h) = 1 N (h) x j -x i h γ kl (x i , x j ). ( 6
)
This quantity depends on the spatial characteristics of Y(.) (defined through its multivariate cross-covariance function in the stationary case) and from the the set functions C k through the values at the observation locations {x 1 , . . . , x n } 3 Estimating the variogram by using pairwise likelihood

In this part, we describe a new methodology to perform the multivariate empirical variography of the underlying GRFs from the category observations. This methodology is based on the pairwise likelihood (PL) maximization. We first recall the principle of the more general composite likelihood based approach. Then we show how to apply it for the plurigaussian model. Finally, we describe more precisely the algorithm in two particular cases: the monogaussian case (q = 1) and the plurigaussian case in which the q GRFs are independent and the sets c(x) are cartesian products of real subsets.

Composite likelihood maximization

The PL approach belongs to the family of the composite likelihood methods (see e.g. [START_REF] Varin | An overview of composite likelihood methods[END_REF], for a comprehensive review). It is generally used to estimate a parameters vector θ of a statistical model, for instance when the usual maximization of the full likelihood is computationally cumbersome. In these cases, the full likelihood is replaced by a weighted product of marginal or conditional likelihoods. [START_REF] Lindsay | Composite likelihood methods[END_REF] defines the composite likelihood as follows: if W a random vector with multivariate density f (w; θ) and {A β , β ∈ I} is a set of marginal or conditional events with associated likelihoods L β (θ; w) ∝ f (w ∈ A β ; θ) for a finite set I, the composite likelihood is the weighted product

L C (θ; w) = β∈I L β (θ; w) λ β
where λ β are nonnegative weights to be chosen.

One of the advantages of the composite likelihood based approaches is that they enable to estimate only some components of the parameters vector θ. We use this idea below to derive a semiparametric estimator of the underlying variograms. Before that, we present an introductive example to show that, under suitable assumptions, the usual empirical variogram can be considered as a maximum of a composite likelihood.

The empirical variogram as a solution of an optimization problem

Suppose that we observe (y(x 1 ), . . . , y(x n )) derived from an intrinsec GRF with variogram γ.

Curriero and Lele (1999) model γ with a parametric form γ θ and estimate θ by maximizing the marginal likelihood based on pairwise differences:

L C (θ; w) = n-1 i=1 n j=i+1 f (v ij ; θ) λ ij (7)
where f (v ij ; θ) stands for the density of the increment V ij = Y (x j )-Y (x i ) and λ ij are some weights to choose.

Here we adopt a semiparametric approach, (as [START_REF] Im | Semiparametric estimation of spectral density with irregular observations[END_REF], with the spectral density). We suppose that θ contains the variogram values γ(h α ) for a set of n l lags h α , α ∈ [[1, n l ]]. Then we group the pairs of points (x i , x j ) according to their distance as follows: when their exists α such as x j -x i h α , we set λ ij to 1 and we consider that γ(x j -x i ) = γ(h α ); otherwise, if it is not possible to associate a pair (x i , x j ) to a lag, its weight is set to 0. Then, it is straightforward to show that the quantity which maximizes the associated marginal likelihood based on pairwise differences (7) is nothing but the traditional empirical variogram estimator of [START_REF] Matheron | Traité de Géostatistique appliquée[END_REF]. We use the same idea to estimate the underlying variograms in the plurigaussian model.

Adaptation to the plurigaussian case

We would like to estimate Σ(h α ), the cross-covariance matrices of the vectors

(Y 1 (x), Y 1 (x + h α ), . . . , Y q (x), Y q (x + h α )) for a set of n l separation vectors h α , α ∈ [[1, n l ]].
The particular composite likelihood which is adapted to this problem is the pairwise likelihood

L C (θ; w) = n-1 i=1 n j=i+1 f (w i , w j ; θ) λ ij
where f (w i , w j ; θ) are the bivariate densities of (W i , W j ) for all (i, j

) ∈ [[1, m]] 2 .
To estimate the covariance matrices Σ(h α ), we consider that θ contains all the unknown elements of the matrices Σ(h α ) for α ∈ [[1, n l ]]. Then, we group the pairs of sites according to their separation vector in the same way as the empirical variogram computation and we write the log PL as follows:

l C (θ; f) = n l α=1 x j -x i hα log p ij (Σ(h α )) (8)
where

p ij (Σ) = c(x i ) c(x j ) g Σ ((u, v))dudv
is the probability that F (x i ) = f (x i ) and F (x j ) = f (x j ) when the cross-covariance matrix of the vector (Y 1 (x), Y 1 (x+h α ), . . . , Y q (x), Y q (x+ h α ) is Σ. Again, the weights λ ij attached to a pair (i, j) have been set to 1 if there exists α ∈ [[1, n l ]] such as x j -x i h α and to 0 otherwise.

Then, the PL estimator is obtained by maximizing l C with respect to all the matrices Σ(h α ). Note that to satisfy the stationarity of the resulting model, the condition

Cov(Y r (x i ), Y s (x i )) = Cov(Y r (x j ), Y s (x j ))
is required for all locations x i and x j and all variable indices (r, s) ∈ [[1, q]] 2 . It implies that the 2q × 2q-matrices Σ(h α ) belong to the set noted S 2q defined by

B ∈ S 2q ⇔ b 2r,2s = b 2r-1,2s-1
where b r,s stands for the (r, s) th element of the matrix B. Furthermore, it is important to remark that all the matrices Σ(h α ) share some common terms to estimate, the ones corresponding to C rs (0). These two constraints on the global solution make the problem numerically difficult to solve.

In this paper, we focus on the simplified cases where the q GRFs are independent which is a quite ordinary assumption for the practitionners of the plurigaussian models. Then, C rs (h) = 0 for all (r, s) ∈ [[1, q]] 2 and all h ∈ R d . Therefore, the matrices Σ(h α ) do not share any common term to estimate simultaneously. Furthermore, the positive definiteness of the matrices Σ(h α ) is satisfied as soon as its non-zero off-diagonal elements belong to ] -1, 1[. It results that the maximization of the log PL can be achieved by solving n l simpler maximization problems: for each α ∈ [[1, n l ]], we can estimate the "parameters" q-vector θ α with r th element ρ r (h

α ) = Cor(Y r (x), Y r (x + h α )) by θα = arg max ρ∈]-1,1[ q x j -x i hα log c(x i ) c(x j ) q r=1
g ρr (u r , v r )du r dv r (9)

where ρ = (ρ 1 , . . . , ρ q ). In the monogaussian case (q = 1), the only quantity to estimate for a given lag h α is the spatial correlation of the underlying univariate GRF,

ρ(h α ) = Cor(Y (x), Y (x + h α )) or equivalently γ(h α ) = 1 - ρ(h α ).
Hence, the PL maximization problem ( 9) is reduced to a one dimensional optimization problem over a bounded interval; it can easily be solved, for instance with the golden section search algorithm [START_REF] Press | Numerical Recipes: The Art of Scientific Computing[END_REF].

For q > 1, the generalization is straightforward if we assume that all the sets C k (x) are cartesian products of subsets of R:

C k (x) = q r=1 T r k (x) with T r k (x) ⊂ R, By denoting t r (x) = T r k (x)
where k is such that f (x) = f k is the actual category at site x, we have:

p ij (Σ(h α )) = q r=1 tr(x i ) tr(x j ) g ρr(hα) (u, v)dudv.
In other word, each ρ r (h α ) is estimated by:

ρ r (h α ) = argmax ρ∈]-1,1[ x j -x i hα log tr(x i ) tr(x j ) g ρ (u, v)dudv.
which is equivalent to solve q problems similar to the monogaussian case.

C 1 (h) = e -h/20 C 2 (h) = e -(h/40) 2
4.1 The monogaussian case, q = 1 On a 1-dimensional regular grid with mesh size 1 and 2000 nodes, for i = 1, 2, 1000 realizations of standardized GRF Y i (.) with covariance C i (h) have been drawn. For each realization, y(.), one category among the set F = {black, red, green} is assigned to each node x of the grid according to the following rule:

f (x) =    orange if y(x) ∈ C 1 (x) = (-∞, s 1 (x)) black if y(x) ∈ C 2 (x) = (s 1 (x), s 2 (x)) green if y(x) ∈ C 3 (x) = (s 2 (x), +∞)
We consider two coding function cases:

• the constant case: s 1 (x) = -s 2 (x) are chosen such as the probability that P (Y (x) ∈ C i (x)) = 1 3 for all i = 1, 2, 3 and all x ∈ R d , • the varying case where s 1 (x) and s 2 (x) have been simulated once for all simulations.

Hence, 4 different categorical random fields have been considered by crossing all the possibilities (C 1 or C 2 vs. constant or varying coding function). The empirical variogram of the underlying GRF is computed by pairwise likelihood from the categories as described in section 3.3, for 150 distances ranging regularly from 1 to 150. For comparison purpose, in each case, the traditional empirical variogram has been computed directly on the realizations y(.) for the same set of distances.

The results are summarized on figure 3. In each case, the average over all the simulations display a negligible bias. As expected, the variability of the estimator increases with the distance. The variogram seems to be better estimated when computed from categories by PL than with the original Gaussian values despite the loss of information due to the truncation. The reason is that we provide additional information by fixing the sill to 1 in the computation by PL. If we compare the results with respect to the covariance model, the simulations with model C 1 always display more important statistical fluctuations than with C 2 . The comparison between the constant coding function case and the varying one shows that the statistical fluctuations around the mean are greater for the latest. It is probably due to the fact that with the constant proportions 1/3 of each category, a lot of transitions occur, bringing more information on the spatial correlation structure of the GRF than in the varying case. Indeed, for this latter case, the values of the probability of a given category are very strong in some areas, leading to few transitions and therefore less information on the hidden field.

4.2 The independent bigaussian case, q = 2 with C a constant cartesian product of intervals

In this example, we consider the same categories as previously. They are generated by assuming that Y 1 (.) and Y 2 (.) are independent with respective covariance function C 1 and C 2 . The categories are assigned to a point x according to the following rule:

f (x) =    black if y(x) ∈ C 1 (x) = (-∞, s 1 ) × R orange if y(x) ∈ C 2 (x) = (s 1 , +∞) × (-∞, t 1 ) green if y(x) ∈ C 3 (x) = (s 1 , +∞) × (t 1 , +∞)
where s 1 = t 1 = 0 such that P ((Y 1 (x), Y 2 (x) ∈ C 1 (x)) = 1 2 and P ((Y 1 (x), Y 2 (x) ∈ C i (x)) = 1 4 for i = 2, 3. A scheme of this coding function is displayed fig. 1 (b). Then 1000 simulations are performed on 800 locations chosen uniformely on the square [0, 1] × [0, 1] one time for all the simulations. A realization is displayed fig. 5.

For each simulation, an empirical omnidirectional variogram is computed for 30 distances ranging from 0 and 150 with a tolerance factor on the distance. The results are summarized fig. 5 and again, they are rather good compared to the empirical variograms computed directly from the Gaussian data. Note that for each simulation, the computation of the empirical variogram of the second Gaussian from y(.) has been computed by using only the subset of locations for which the first Gaussian is greater than 0.

Discussion

In this paper, we propose to use the pairwise likelihood principle to estimate empirical variograms of the underlying GRFs in the plurigaussian model. The use of a composite likelihood based approach to provide an exploratory data analysis tool is apparently original, even if, as above-mentionned, the usual empirical variogram can be seen as the solution of an optimization involving a composite likelihood.

Once the empirical variograms of the underlying GRFs have been computed, we can use them to choose a valid multivariate variogram model. They can be fitted by least squares, for instance by using the algorithm proposed in [START_REF] Desassis | Automatic variogram modeling by iterative least squares: Univariate and multivariate cases[END_REF] or their parameters can be estimated by a likelihood based method. The likelihood will probably remain intractable since it involves at least an integral on R n where n is the number of samples. Again, a composite likelihood based approach should be used instead. Again, the PL seems well suited.

To conclude, note that the method presented in this paper is implemented in the R-package RGeostats [START_REF] Renard | RGeostats: The geostatistical package[END_REF] in the function named vario.pgs. Some demonstration scripts are provided through a tutorial on the dedicated website.

Further researches will concentrate on the generalization of the approach presented in this paper to the case where no independence assumption is made between the underlying GRF. Indeed, more complex transitions between categories can be investigated, with more general multivariate spatial models (see [START_REF] Galli | An application of the truncated plurigaussian method for modeling geology[END_REF]. In that case, all the elements (except the diagonal) of the correlation matrices Σ(h α ) must be estimated with the constraints mentionned in section 3.3. This is computationally much more challenging. Another natural extension could be to adapt the local variogram kernel estimator proposed in [START_REF] Fouedjio | A generalized convolution model and estimation for non-stationary random functions[END_REF] to the plurigaussian context by using PL. Then, it would be possible to deduce a non-stationary model for the underlying GRFs, for instance with varying anisotropies.

Finally, the PL likelihood approach to compute empirical variograms seems to be a promising idea which could be applied to other similar context of hidden variable, or variable known after a transformation.

To cite some of them:

• compute the empirical variogram of the underlying GRF in the hierarchical geostatistical models (see e.g [START_REF] Diggle | Model based geostatistics[END_REF].

Some authors have already proposed a way to compute empirical variogram of underlying random fields in hierarchical models: [START_REF] Oliver | Estimating the risk of childhood cancer[END_REF] treats the binomial case, [START_REF] Monestiez | Geostatistical modelling of spatial distribution of balaenoptera physalus in the northwestern merditerranean sea from sparse count data and heterogeneous observation efforts[END_REF] the poisson case. However, these estimators are based on the method of moments and the distribution of the underlying random field is not specified. Thus, the underlying intensity can only be predicted by kriging but they can not be simulated. An approach based on the PL in a distribution based framework could be a good alternative;

• perform the multivariate empirical variography of the underlying GRFs when one has to deal with a continuous variable vs. discrete variable [START_REF] Emery | Conditional co-simulations of continuous and categorical variables for geostatistical applications[END_REF], or even two discrete variables [START_REF] Chagneau | A hierarchical bayesian model for spatial prediction of multivariate nongaussian random fields[END_REF][START_REF] Renard | Heterotopic bicategorical variables in plurigaussian simulation[END_REF];

• compute the empirical variogram of a variable at punctual level when the observations are some regularizations with different supports. ), average of the empirical variograms over all the simulations ( ), 25th and 75th percentiles ( ), 5th and 95th percentiles (

).

Figure 1 :Figure 2 :

 12 Figure 1: Example of coding functions for (a) q = 1 and (b) q = 2.

Figure 3 :Figure 4 :

 34 Figure 3: Simulation results in the monogaussian cases with covariance C 1 (left) and C 2 (right), in the constant coding function case (up) and in the varying case (middle). The down figures stand for the Gaussian case. Actual model ( ), empirical variogram of ten arbitrary simulations ( ), average of the empirical variograms over all the simulations (), 25th and 75th percentiles (), 5th and 95th percentiles ( ).

Illustration on simulationsIn this section, we present two simulation studies to assess the efficiency of the proposed method. We work with two covariance models:

Acknowledgments: We are very grateful to Christian Lantuejoul for the idea to directly compute the empirical variogram of the underlying GRF. We would also like to thank Geovariances and the School of Earth Science of the University of Queensland to have partly funded this research during the visit of the first author in Brisbane in 2012/2013.