Nonparametric estimation of utility functions
Mengyang Gu Debarun Bhattacharjya Dharmashankar Subramanian

input@path
 






















 




















 



definitionDefinition
theoremTheorem
lemmaLemma
propositionProposition
corollaryCorollary
remarkRemark
exampleExample
assumptionAssumption

#1equation (<ref>)
p.d.f.  c.d.f. 
p.d.f.s c.d.f.s
m.g.f.  m.g.f.s 
⊥⊥                        
 ( [ ) 
diag
loglog
tr
tantan





 


























 



































<cit.>
<cit.>
(<ref>)
7













equationsection
plain
thmTheorem[section]






Nonparametric estimation of utility functions














[addr1]Department of Applied Mathematics and Statistics, Johns Hopkins University, Baltimore, MD 
    e1 
    


[addr2]IBM T.J. Watson Research Center, Yorktown Heights, NY  
 e2,e3 
 













Inferring a decision maker's utility function typically involves an elicitation phase where the decision maker responds to a series of elicitation queries, followed by an estimation phase where the state-of-the-art is to either fit the response data to a parametric form (such as the exponential or power function) or perform linear interpolation. We introduce a Bayesian nonparametric method involving Gaussian stochastic processes for estimating a utility function. Advantages include the flexibility to fit a large class of functions, favorable theoretical properties, and a fully probabilistic view of the decision maker's preference properties including risk attitude. Using extensive simulation experiments as well as two real datasets from the literature, we demonstrate that the proposed approach yields estimates with lower mean squared errors. While our focus is primarily on single-attribute utility functions, one of the real datasets involves three attributes; the results indicate that nonparametric methods also seem promising for multi-attribute utility function estimation.
 
 
 
 
 
 
 
 




















Gaussian stochastic process
reference prior
robust estimation
utility function






§ INTRODUCTION



Making decisions under uncertainty using decision theory requires that beliefs about uncertainties be represented by probabilities and preferences over outcomes be summarized by utilities. The latter is often accomplished by assessing a decision maker's utility function over a single attribute (such as monetary units) or over multiple attributes, depending on the decision situation under consideration.

The vast literature on modeling utility functions describes schemes where the decision maker responds to elicitation queries; these responses are then subsequently used to infer the decision maker's utility function. We distinguish between these phases and refer to them as elicitation and estimation respectively. Posing elicitation questions in a specific form and then estimating the functional form that best represents the decision maker's preferences often go hand in hand.
In this paper, we contribute to the literature on estimation by introducing a Bayesian nonparametric approach for modeling utility functions. Specifically, we demonstrate the advantages of using a Gaussian stochastic process (henceforth GaSP) over the state-of-the-art, where parametric methods and the nonparametric method of linear interpolation are most popular. We show that linear interpolation corresponds to a subclass of our GaSP model using a specific covariance function in which the process turns to be a Wiener process (<cit.>); moreover, parametric methods can be incorporated into the GaSP model through the mean function.


While the GaSP model has been popular in domains such as nonlinear regression and classification in machine learning (<cit.>), spatial statistics (<cit.>]) and computer model emulation and calibration (<cit.>), it has not been applied to the problem of estimating utility functions, to the best of our knowledge.
There is some related work in the artificial intelligence literature on preference modeling with GaSP. For instance, <cit.> and <cit.> use GaSP for preference learning over a set of items. In contrast, our work is on estimating utility functions based on elicitation schemes from the operations research and management science literature (see e.g. <cit.>).


A fundamental challenge in utility elicitation is that the number of assessed queries is often very small, since elicitation burden increases significantly and responses become less reliable when more questions are posed.
The main contribution of this work is a robust estimation approach using GaSP models that reduces predictive errors when there are a small number of observations, in comparison with approaches such as parametric methods and linear interpolation. Furthermore, we provide a coherent model to accurately quantify the uncertainty associated with the utility elicitation process, as compared to the point predictions that are popular in the literature. Another benefit of such a coherent model is that we are able to assess risk attitudes of decision makers more accurately than current adopted methods.

Our proposed methods for utility elicitation are available online as an R package (<cit.>).





The remainder of the paper is organized as follows. In sections 2 and 3, we introduce basic notation and terminology, placing our contribution in the context of related literature in operations research and management science on utility elicitation and estimation. Here we motivate our research effort and also discuss some theoretical advantages of nonparametric methods. In section 4, we introduce our GaSP formulation for utility function estimation and also show how to compute derivatives for assessing risk attitude. In section 5, we make a practical case for our proposed method by comparing it with prevalent methods for estimating single-attribute utility functions using a few assessed points on the utility function. We conduct several simulation experiments and also study a real-world dataset (<cit.>).

In section 6, we delve into estimating multi-attribute utility functions, demonstrating that GaSP performs better than the estimation methods that were deployed for a three-attribute dataset (<cit.>). Finally, we conclude the paper in section 7.




§ UNCERTAINTY QUANTIFICATION FOR UTILITY ELICITATION



Consider outcomes 𝐱 in a separable metric space 𝒳 (e.g. ℝ^n).
<cit.> showed that preferences over  uncertain outcomes in such a space are complete, transitive and continuous in 𝒳 iff there exists a continuous utility function representation U: 𝒳→ℝ.

Behavioral literature indicates that people tend to construct their preferences during the elicitation process, often providing inconsistent responses to queries (<cit.>). The implication is that a decision maker's utility function U should perhaps be considered an approximate representation of their preferences; for this and other reasons we shall discuss shortly, we make the following distinction:


When the decision maker answers all the preference elicitation questions consistently with the same underlying (and typically unknown) utility function, the assessment is said to be noise-free, otherwise it is noisy.



The notion of a `true' underlying utility function can be viewed as a theoretical construct and one that has been discussed often in the literature. One way to model the uncertainty in preference elicitation responses is to include a random response error to a systematic component, either as additive error (Laskey and Fischer (1987)) or as random parameters of the utility function
(<cit.>, Fischer et al. (2000a)).
Another approach is to treat the utility function as inherently stochastic (Becker et al. 1963).
In practice, decision analysts are well aware of these issues and have devised measures to counter inconsistencies (see for example <cit.>).



Making a noise-related distinction enables us to express different sources of uncertainty in the elicitation process. One source of uncertainty is prediction uncertainty, representing the system's uncertainty about the decision maker's utility at an unassessed 𝐱. This uncertainty is present regardless of whether decision makers are consistent with their answers. The second source of uncertainty depends on the decision maker. When they consistently answer questions with the same underlying utility function, the elicited utility u(𝐱) is identical to U(𝐱) at each assessed 𝐱. We refer to an estimation method that agrees with the elicited utility at each assessed 𝐱 as an interpolator. It is ideal for a method to be an interpolator for a consistent decision maker. Note that the second source of uncertainty only appears when decision makers do not answer questions consistently, leading to noisy assessments.








Preference elicitation for decision making under risk is typically conducted using gambles. Consider gamble (𝐱_a, p; 𝐱_b) that results in outcome 𝐱_a with probability p and outcome 𝐱_b with probability 1-p, where  𝐱_a, 𝐱_b ∈𝒳.
In an elicitation query, the decision maker must evaluate two or more gambles presented to them, e.g. they may be asked to compare gamble (𝐱_a, p; 𝐱_b) with the degenerate gamble (𝐱_c). If the assessment is noise-free and if the evaluation is done by expected utility theory (EUT), they prefer the first gamble iff V_EUT(𝐱_a, p; 𝐱_b) = p U(𝐱_a)+(1-p)U(𝐱_b) ≥ U(𝐱_c), for some underlying utility function U (<cit.>).
A subsequent estimation task must be performed to infer U from the responses.




There is significant empirical evidence from the descriptive literature on prospect theory demonstrating that people tend to overweight low probabilities and underweight high probabilities, thus under prospect theory, the gamble (𝐱_a,p;𝐱_b) evaluates to
V_PT(𝐱_a, p; 𝐱_b)= ω (p) U(𝐱_a) +ω (1-p)  U(𝐱_b), where ω (.) is a probability weighting function and the reference point is assumed to be 0 (<cit.>). Prospect theory also explains observed behavior such as loss aversion and diminishing sensitivity relative to the reference point.


<cit.> and others argue that descriptive violations of expected utility bias utility elicitations, therefore the design of elicitation schemes and estimation of utility functions should be conducted based on prospect theory rather than expected utility theory, resulting in more accurate estimates of U. We are sympathetic to this view but as we explain in the next subsection, our research goal is purely that of estimating U, provided elicitation responses in the following form:

 Assessed tuples are assessments of points on the utility function, (𝐱_i, u(𝐱_i) ) for i=1,…,n, 𝐱_i ∈𝒳.


We propose a GaSP approach to estimating a utility function given assessed tuples as input.
Our method can thus be applied to different data sets and is agnostic to the underlying theory. In fact, we demonstrate GaSP estimation using elicitation schemes based on both expected utility theory and prospect theory. Furthermore, it is also agnostic to the generative structure of randomness in the case of probabilistic preferences, i.e. whether the data is generated by an additive noise model, or a utility function with random parameters, or indeed any stochastic utility model.
As mentioned earlier, we consider both noise-free and noisy assessments. Note that for noise-free assessments, assessed tuples (𝐱_i, u(𝐱_i) ) = (𝐱_i, U(𝐱_i) )

for underlying utility function U.

We refer the reader to <cit.> for a comprehensive review of preference elicitation schemes involving single attributes; <cit.> also describe elicitation schemes for multi-attribute problems.


 





§ UTILITY ESTIMATION


The goal of any estimation task associated with preference elicitation is to use responses to the elicitation queries to infer the decision maker's utility function U(𝐱). Here we discuss parametric estimation and linear interpolation, highlighting potential limitations.



 §.§ Parametric Utility Estimation


Parametric estimation is a popular approach for estimating utility functions involving a single attribute x (<cit.>).
The most common parametric forms are the exponential and power functions. The exponential utility function takes the form a - b sgn(ρ) exp ( -x/ρ ), where a and b>0 are constants and sgn(ρ) is the sign of the risk tolerance parameter ρ≠∞.
The power utility function is of the form a + b sgn(α) sgn(x) |x|^α, with constants a and b>0, where sgn(α) and sgn(x) are the signs of α≠ 0 and x.
Note that the power function exhibits discontinuous behavior around 0 for α < 0; see <cit.> for details.
The limiting cases for the exponential and power functions are the linear and logarithmic functions; these two families of functions are the only ones that satisfy constant risk aversion and constant relative risk aversion respectively (<cit.>).


For multi-attribute utility functions, the most popular approach is to check for utility independence assumptions that decompose the multi-dimensional function into an additive or multiplicative function of one-dimensional marginal utility functions (<cit.>). These marginal functions typically take the form of the afore-mentioned standard parametric models.

There has also been literature on utility functions that do not enforce utility independence assumptions.
To account for utility dependence, for instance, parametric forms have also been used for conditional utility functions, i.e. the utility of a subset of attributes given that the other attributes are set at fixed consequences (see for example <cit.>).
A more recent approach is to use a copula (<cit.>), which has the benefit that the level of dependency between marginal utility functions can be represented in an easily interpreted parametric way.
In our GaSP model, utility dependence between attributes is modeled through a product covariance function, as discussed in section <ref>.

It is not hard to see that a parametric method will not be an interpolator unless the decision maker's underlying utility function follows the parametric class being used. The lack of this desirable property could result in large elicitation errors when the parametric class is misspecified, as we show in Section <ref>.
















































 §.§ Linear Interpolation


An alternate approach that is popular in the empirical literature on estimating single-attribute utility functions is that of piece-wise linear interpolation across assessed tuples (<cit.>). Such an approach is essentially equivalent to the predictive mean of the extended Wiener process, defined as a stochastic process W_t with independent, normally distributed increments W_t-W_s ∼𝒩(0,t-s) for t≥ s≥ 0 with continuous sample paths (<cit.>). A Wiener process is typically defined to have initial value W_0=0 but we relax this assumption. The predictive distribution formalized in the following lemma. (Please see Appendix for all proofs.)

 
Assume W_t, t ∈𝒯 follows a Wiener process. Assume we have observations W_t_1,...,W_t_n with 0< t_1< t_2<...<t_n. For any t_i≤ t_*≤ t_i+1, for any 1≤ i<n and i ∈ℕ , the predictive distribution of W_t_* given  W_t_1,...,W_t_n is
 
    W_t_*| W_t_1,...,W_t_n∼𝒩(μ_*, V_* ),

where μ_*=(t_i+1-t_*)W_t_i+ (t_*-t_i)W_t_i+1/t_i+1-t_i and V_*=(t_i+1-t_*)(t_*-t_i)/t_i+1-t_i.






 


Lemma <ref> states that if the utility function is modeled as a Wiener process for a single attribute t,  the posterior mean μ_*=[W_t_*| W_t_1,...,W_t_n] for assessing the utility at a point t_* (that has not been assessed) is equivalent to linear interpolation between two nearby assessed tuples.




 

Indeed, a Wiener process is a special case of GaSP (which we will formally define in the next section) with initial value W_0=0, mean zero and covariance function (W_s,W_t)=min(s,t) at any s,t ≥ 0, and continuous sample path. However, a Wiener process is not differentiable everywhere and consequently using the posterior mean (i.e. linear interpolation) poses problems for estimating the coefficient of risk aversion.




[Linear interpolation and overconfidence]
Consider Figure <ref>, which displays the function y= 3sin(5π t)+cos(7π t) (treated as unknown) whose values are assessed at 12 equally spaced points in [0,1].
The left panel shows the prediction (blue curves) by the extended Wiener process (for which we do not assume W_0=0). Not only does the prediction show discrepancy at places where the derivatives of the function change, it is also clearly overconfident as the 95% confidence interval covers regions that are a lot smaller than the nominal 95%. In comparison, the right panel is the prediction by the method we propose with the same 12 assessed points, using the default setting in the RobustGaSP R Package (<cit.>). The improvement arises due to the use of other covariance functions with objective Bayesian inference, discussed in the next section in detail. While this example illustrates a shortcoming of linear interpolation using a generic function that is measured at a finite number of values in its domain, the above weakness is relevant for our purposes of estimating a utility function using experimentally assessed tuples.







Another limitation of a Wiener process estimation approach is that it is only defined in a one-dimensional domain and thus is limited to single-attribute utility function estimation. Although there is some literature on interpolation in multi-attribute problems (see for example <cit.>), it is a challenging task due to the curse of dimensionality. A more general GaSP approach with suitable covariance functions built on the space of multiple attributes may be able to better handle the estimation task.



 §.§ Quantile-Parameterized Distributions


Quantile-parameterized distributions (QPD) have recently been introduced for modeling uncertainties in decision analysis (<cit.>). This approach characterizes a continuous probability distribution based on a number of assessed quantile/probability pairs.
Although QPDs have not been discussed in the context of utility elicitation, it is straightforward to apply the approach conceptually as assessed tuples are analogous to quantile/probability pairs.
Denoting these assessed tuples as (x_i, u_i) for i=1,..,n, where u_i is the assessed utility scaled from [0,1], the inverse CDF of a QPD takes the following form:

    F^-1(u)={ L_0      u=0,
    ∑^q_i=1  a_ig_i(u)     0<u<1, 
     L_1     u=1,
    .

with left handed limit L_0=lim_y→ 0^+F^-1(u), right handed limit L_1=lim_y→ 1^-F^-1(u), and g_i(·) referring to basis functions for i=1,...,q, q≤ n. When the domain of x is a real line, a simple choice is the Q-normal distribution, where the basis functions are g_1(u)=1, g_2(u)=Φ^-1(u), g_3(u)=yΦ^-1(u), g_4(u)=u, with Φ being the normal CDF (<cit.>). Note that the results depend entirely on the choice of basis functions; further research is required to explore suitable models for utility elicitation.
Here we consider QPDs solely as another benchmark for the GaSP model since they provide some flexibility to model a curve ranging from [0,1].







§ ESTIMATION WITH GASP


We introduce a Bayesian nonparametric method that takes assessed tuples as training data input, regardless of the underlying theory and assumptions used to derive them, and provide an estimated utility function û ( 𝐱 ), where 𝐱 could either be a single attribute or multiple attributes. Bayesian inference explicitly provides a formal way to quantify the uncertainty introduced by the estimation through the likelihood and prior, as discussed in detail in this section.




 §.§ Model Formulation



To set notation, let 𝐱=(x_1,...,x_p)^T be a vector of p different attributes and let  u(𝐱) be the utility evaluated at 𝐱. Let us consider a random utility function modeled in a general regression way with the following form,
   
    u(𝐱)= m(𝐱)+z(𝐱),

where m(𝐱) is the mean function, modeled as
    
    E [u( 𝐱)]= m(𝐱)=𝐡( 𝐱)θ =∑^q_j=1  h_j(𝐱)θ_j,

where 𝐡(𝐱) is assumed to be a q dimensional domain dependent basis function for any 𝐱∈𝒳, with unknown regression parameters θ_j for each basis function h_j(𝐱). h_j(𝐱) could be chosen, e.g., as a particular parametric form as specified in Section <ref> or as a polynomial function in 𝐱. For the additive residual term, instead of taking z(𝐱) as independent measurement errors as in <cit.>, we model z(·) as a stationary GaSP

    z(·)  ∼  GaSP(0,  σ^2c(·,  ·) )  ,

with variance  σ^2 and the pair-wise correlation function c(·, ·). In return, the joint distribution of any n inputs {𝐱_1,…,𝐱_n}∈𝒳, follows a multivariate normal distribution,

    ((z(𝐱_1),…,z(𝐱_n))^T| (σ^2,  𝐂) ) ∼𝒩 (0, σ^2   𝐂  ) ,

i.e. a normal distribution that is conditional on the unknown variance σ^2 and the Gram matrix 𝐂 (<cit.>) whose (i,j) element is c   (𝐱_i, 𝐱_j).
By definition, the covariance of the utility is
  
    (u( 𝐱_a), u( 𝐱_b))=σ^2c( 𝐱_a, 𝐱_b ),

for any 𝐱_a, 𝐱_b ∈𝒳. If p=1 (i.e. single attribute) and c(x_a,x_b)= min(x_a, x_b), GaSP becomes a Wiener Process when the process is defined on x ∈ [0,+∞), with initial value 0. To extend the definition to the case when p>1, the  isotropic assumption is sometimes made for modeling a spatial process (<cit.>), meaning that the correlation function c(𝐱_a, 𝐱_b ) is only a function of ||𝐱_a- 𝐱_b||  where  ||·|| is the Euclidean distance. However, the domain of attributes typically varies on completely different scales (e.g. between the price and comfort of a car), so the effect of the attributes on the correlations will be highly variable. Consequently, the assumption of isotropy may not be reasonable. Instead, the product correlation function is often assumed:


    c(𝐱_a,  𝐱_b)=∏^p_l=1c_l( x_al,  x_bl),

where c_l(·, ·) is a one-dimensional correlation function for the l^th attribute.  We list several frequently used correlation functions in Table <ref>. The difference between the above product correlation function and the isotropic assumption is that for the former, there are parameter(s) in each correlation  c_l(·, ·) that can control the smoothness of the utility function on this attribute (which could be learned from the data), while  the latter assumes the Euclidean distance.







The power exponential covariance and the Matérn covariance have been used in many applications. When ν_l=(2k+1)/2 where k∈ℕ,  Matérn correlation has a closed form. For example, when the roughness parameter  ν_l=5/2, the Matérn correlation is


    c_Mat(x_al,x_bl)=(1+√(5)d_l/γ_l+5d_l^2/3γ_l^2)exp(-√(5)d_l/γ_l),

where d_l=|x_al- x_bl|.  GaSP with the Matérn correlation defined in Equation (<ref>) means it is  twice mean square differentiable (<cit.>), allowing one to infer the risk attitude using derivatives of the GaSP as discussed later in Section <ref>. We use the Matérn correlation defined in Equation (<ref>) for the purpose of demonstration but we do not preclude the use of other correlation functions with suitable differentiable results in future applications. In practice, the roughness parameter (ν_l) is fixed at a chosen value (e.g. 5/2 for the Matérn), and the correlation function c_l(.,.) is essentially parameterized by the range parameters (γ_l).



















































 §.§ Posterior Predictive Distribution



Denote assessed tuples as ( 𝐱^𝒟, 𝐮(𝐱^𝒟) ), where 𝐱^𝒟 = {𝐱^𝒟_1, 𝐱^𝒟_2,…, 𝐱^𝒟_n } are n points in the domain of the multi-attributes. As we are empirically limited by a small number of assessed utility values, we seek the posterior predictive distribution of the utility function u^*(x^*) at any input x^* ∈𝒳 based on the assessed tuples ( 𝐱^𝒟, 𝐮(𝐱^𝒟) ). For simplicity, denote 𝐮^𝒟= (u(𝐱_1^𝒟),u(𝐱_2^𝒟),...,u(𝐱_n^𝒟))^T as the assessed utility points in the design. As per the chosen GaSP model, the likelihood is a multivariate normal likelihood:
  
    ℒ(𝐮^𝒟| θ,σ^2,γ )=(2πσ^2)^-n/2|𝐂|^-1/2exp{- ( 𝐮^𝒟- 𝐡(𝐱^𝒟) θ )^T (σ^2 𝐂)^-1 ( 𝐮^𝒟- 𝐡(𝐱^𝒟) θ )/2},

where 𝐡(𝐱^𝒟 ) is the n × q  basis design matrix with (i,j) element h_j(𝐱^𝒟_i).  The model parameters for posterior estimation are the mean parameter θ=(θ_1,θ_2,...,θ_q)^T, variance parameter σ^2 and range parameters γ=(γ_1,γ_2,...,γ_p)^T in the correlation function in Equation (<ref>).


 We complete the model by specifying reference priors, π, for the model parameters (<cit.>)

    π(θ ,σ^2,γ ) ∝|𝐈^*(γ )|^1/2/σ^2 ,


where 𝐈^*(·) is the expected Fisher information matrix, given as

    𝐈^*(γ ) = ( [      n - q    tr(𝐖_1)    tr(𝐖_2)        ...    tr(𝐖_p);             tr(𝐖_1^2) tr(𝐖_1𝐖_2)        ... tr(𝐖_1𝐖_p);                        tr(𝐖_2^2)        ... tr(𝐖_2𝐖_p);                                           ⋱          ⋮;                                              tr(𝐖_p^2) ])_(p+1) × (p+1),

where 𝐖_l = 𝐂̇_l𝐐,   1≤  l ≤ p with 𝐐 =  𝐂^ - 1 ( 𝐈_n - 𝐡(𝐱^𝒟 ){𝐡^T(𝐱^𝒟  ) 𝐂^ - 1𝐡(𝐱^𝒟  )}^ - 1𝐡^T(𝐱^𝒟  ) 𝐂^ - 1) and 𝐂̇_l is the derivative of 𝐂 with regard to γ_l.








 After integrating (θ, σ^2) using the above reference prior,  the marginal posterior of γ follows

    p(γ|𝐮^𝒟)  ∝ℒ(𝐮^𝒟|γ)|𝐈^*(γ )|^1/2,

where the marginal likelihood is

    ℒ(𝐮^𝒟 |γ ) ∝  |𝐂|^- 1/2 (|𝐡^T(𝐱^𝒟  ) 𝐂^ - 1𝐡(𝐱^𝒟   )|)^- 1/2( S^2)^-(n - q/2),

 with S^2=    ( 𝐮^𝒟 )^T 𝐐𝐮^𝒟.


The predictive distribution from a Bayesian approach is obtained by integrating the uncertainty from the parameters
  
    p(u(𝐱^*)| 𝐮^𝒟) =∫ℒ(u(𝐱^*)|  𝐮^𝒟, θ,σ^2,γ )p(θ,σ^2,γ|  𝐮^𝒟  ) dθ dσ^2 dγ.


Note that γ cannot be marginalized out explicitly  - there is no prior on γ that leads to a closed form marginal likelihood after integrating out γ. Numerical integration algorithms are both inefficient (as the computation requires the inverse of the correlation 𝐑, which has O(n^3) operations) and less stable. Instead, γ is often estimated by the  maximum marginal posterior mode
 
    (γ̂_1, …γ̂_p)=  arg max _γ_1,…, γ_p{ p(γ|𝐮^𝒟) }.

As discussed in <cit.>, this leads to a robust estimation of γ, yielding much better results than a non-robust method such as the maximum likelihood estimator (MLE). Note that the derivative of the reference prior is computationally intensive when the sample size is large. One can use some other objective priors that both have the robust estimation properties and the closed form derivatives (see e.g. <cit.>).  






With the above setup,  the predictive distribution of u(𝐱^*) at a new point 𝐱^* ∈𝒳, given the assessed tuples and the estimated range parameter γ̂, is a t-distribution (<cit.>)

    u(𝐱^*) |𝐮^𝒟  ,γ̂∼ t( û(𝐱^*),σ̂^2c^**,n - q).

with n-q degrees of freedom, where   û (𝐱^*), σ̂^2 and c^** are  given respectively as,

    û (𝐱^*) =   𝐡(𝐱^*)θ̂+𝐜^T(𝐱^*)𝐂^-1(𝐮^𝒟  -𝐡(𝐱^𝒟  )θ̂)  ,
    σ̂^2 =    (n-q)^-1( 𝐮^𝒟  -𝐡(𝐱^𝒟  )θ̂)^T𝐂^-1(𝐮^𝒟 -𝐡(𝐱^𝒟  )θ̂), 
    
     c^** =     c(𝐱^*, 𝐱^*)-𝐜^T(𝐱^*)𝐂^-1𝐜(𝐱^*) + (𝐡( 𝐱^*)-𝐡^T(𝐱^𝒟  )𝐂^-1𝐜(𝐱^*) ) ^T 
       ×(𝐡^T(𝐱^𝒟 )𝐂^-1𝐡(𝐱^𝒟  ))^-1(𝐡(𝐱^*)-𝐡^T(𝐱^𝒟  )𝐂^-1𝐜(𝐱^*) )  ,

with θ̂=(𝐡^T(𝐱^𝒟  )𝐂^-1 𝐡(𝐱^𝒟  ))^-1𝐡^T(𝐱^𝒟  )𝐂^-1𝐮^𝒟 as the generalized least squares estimator for θ and    𝐜(𝐱^*) = (c(𝐱^*,𝐱^𝒟_1 ), …, c(𝐱^*,𝐱^𝒟_n ))^T being the correlation of the GaSP at 𝐱^* and 𝐮^𝒟.


The predictive mean, û(𝐱^*) is typically used for estimation of the utility function. The following lemma specifies that such an estimation is an interpolator.




The predictive mean û(𝐱^𝒟_i) = U(𝐱^𝒟_i) for any i=1,...,n for noise-free assessments with utility function U(·).























For noisy assessments, we do not expect the prediction of GaSP to be exact at the assessed points.
To model such situations, an independent noise can be added in Equation (<ref>) by defining z̃(𝐱)= z(𝐱) +ε, where ε is independent white noise. Objective Bayesian inference for such a GaSP model is similar to the method discussed above (<cit.>). The afore-mentioned methods have been developed as an R package (<cit.>).






 §.§ Derivatives of GaSP


Differentiability is an important property of utility functions, e.g., the Arrow-Pratt measure of local risk aversion is given as λ(x) = -u”(x)/u'(x). Our proposal to use GaSP is helpful in this regard since the derivative processes are also GaSP when the covariance function is mean square differentiable (<cit.>).

Closed form expressions of derivatives of GaSP with Gaussian correlation (i.e. power exponential correlation with ν_l=2 in Table <ref>) are shown in <cit.> (chapter 3). Here we present the derivatives of the Matérn class correlation with roughness parameter equal to 2.5. For simplicity, we assume there is only a single attribute in this subsection, but the following result can be easily extended to directional derivative processes with regard to each attribute, i.e. ∂ u(𝐱)/∂ x_l.

Denote 𝐝^*=(|x^𝒟_1-x^*|,…, |x^𝒟_n-x^*|)^T and 𝐠=(g_1,…, g_n)^T with g_i={
1    x^*≥ x^𝒟_i

-1    x^*<x^𝒟_i
.. Let ∘ denote the element-wise product of each entry between two matrices. The following two lemmas present the joint distribution of the GaSP and its derivatives, and they can be derived through directly differentiating the covariance function (see e.g. <cit.>).







The joint distribution of GaSP 𝐮^𝒟 and GaSP derivative u' ( x^*) on any x^* ∈𝒳 is:
	
    ( 𝐮^𝒟
     u' ( x^*) 
    ) |θ ,σ ^2,γ,   ∼ N(( [     𝐡(𝐱^𝒟)θ; 𝐡'( x^*  )θ;             ]) ,σ ^2( [                 𝐂    𝐜^01(𝐱^𝒟, x^*);   𝐜^10(  x^*,𝐱^𝒟) c^11( x^* ,  x^*);                   ])).

 where 𝐡'(x^*)=∂𝐡( x^*)/∂ x^*, 𝐜^01(𝐱^𝒟, x^*)=-5 𝐠/3γ^2∘(𝐝^*+√(5)𝐝^* ∘𝐝^* /γ)∘ exp(-√(5)𝐝^*/γ),   and c^11( x^* ,  x^*)= 5/3 γ^2, and 𝐜^10(  x^*,𝐱^𝒟) =𝐜^01(𝐱^𝒟, x^*)^T.
















 The joint distribution of GaSP 𝐮^𝒟 and GaSP derivative u” ( x^*) on any x ∈𝒳 is:
	
    ( 𝐮^𝒟
     u” ( x^*) 
    ) |θ ,σ ^2,γ,   ∼ N(( [ 𝐡( 𝐱^𝒟  )θ;   𝐡”(x^*)θ;            ]) ,σ ^2( [                 𝐂    𝐜^02(𝐱^𝒟, x^*);    𝐜^20( x^*,𝐱^𝒟) c^22( x^* ,  x^*);                   ])).

 where h”( x^*)=∂^2 h( x^*)/∂ x^*∂ x^*, 𝐜^02(𝐱^𝒟, x^*)=5/3γ^2{5𝐝^* ∘𝐝^*/γ^2-1-√(5)𝐝^*/γ}∘ exp(-√(5)𝐝^*/γ), and c^22( x^* ,  x^*)=25/γ^4 and 𝐜^20( x^*,𝐱^𝒟) = 𝐜^02(𝐱^𝒟, x^*)^T.
 





Note that the above joint distribution is used to calculate the predictive distribution of the derivative. Indeed, by marginalizing out the mean and variance parameter, the predictive distribution of u”(x^*) given 𝐮 and  γ̂=(γ̂_1, …, γ̂_p), is a t-distribution


    u”( x^*) |𝐮^𝒟 ,γ̂∼ t( û”( x^*),σ̂^2c”,n - q),

with n-q degrees of freedom, where

    û” ( x^*)    =   𝐡”( x^*)θ̂+𝐜^02(𝐱^𝒟, x^*)𝐂^-1(𝐮^𝒟 -𝐡(𝐱^𝒟  )θ̂)  ,
    σ̂^2    =   (n-q)^-1(𝐮^𝒟 -𝐡(𝐱^𝒟  )θ̂)^T𝐂^-1(𝐮^𝒟 -𝐡(𝐱^𝒟  )θ̂), 
    
     c”   =     c^22( x^*,  x^*)-𝐜^02(𝐱^𝒟, x^*) 𝐂^-1𝐜^20(𝐱^𝒟, x^*)  + (𝐡”( x^*)-𝐡^T(𝐱^𝒟  )𝐂^-1𝐜^20(𝐱^𝒟, x^*)  ) ^T 
          ×(𝐡^T(𝐱^𝒟 )𝐂^-1𝐡(𝐱^𝒟  ))^-1(𝐡”( x^*)-𝐡^T(𝐱^𝒟  )𝐂^-1𝐜^20(𝐱^𝒟, x^*)  )  .


The predictive distribution of u'(x^*) is similar to Equation (<ref>). The point here is to demonstrate that the risk attitude can be obtained using the derivative processes with full assessment of the uncertainty associated with the analysis.




§ SINGLE ATTRIBUTE UTILITY FUNCTION ESTIMATION


So far we have primarily seen the theoretical benefits of estimating utility functions with GaSP. In this section, we explore practical ramifications on estimating single-attribute utility functions by first performing simulation experiments and then analyzing a real dataset.



 §.§ Simulation Experiments





  §.§.§ Comparing utility function estimates


In preference elicitation, only a limited number of questions can be posed, thus the number of assessed tuples n is typically small:
n=4, 7 and 10 are considered for these experiments. Out of sample mean squared error (MSE),
MSE=∑_i=1^n^*{û(x^*_i)-u(x^*_i)  }^2/n^*,

is utilized for comparison, where x^*_i ∈𝒳 is the i^th equally spaced held-out point and n^*=1,001 is used for testing throughout this section. We assume U(x_min)=0 and U(x_max)=1, where x_min=0 and x_max=10^5 are lower and upper bounds of the domain 𝒳 in this subsection. We compare the exponential function (Exp), power function (Pow), linear interpolation (LI) and quantile-parameterized distribution (QPD) method with the GaSP model.

For the parametric methods and QPD method, we estimated the parameters with the minimum least squares error. For the QPD, we choose the basis function to be g_1(u)=1, g_2(u)=Φ^-1_*(u), g_3(u)=uΦ^-1_*(u), g_4(u)=u. As the domain of x is [0,10^5], the usual Q-normal distribution is not a sensible choice. Instead, we let Φ^-1_*(u) be a normal distribution truncated at 0 and 10^5 centered at 0 with standard deviation 5× 10^4. This seems to perform the best among all basis functions considered.




In the first scenario, we assume that assessments are noise-free. Tables <ref> and <ref> display the out of sample MSE when the underlying utility function U is a power function and an exponential function respectively. Because the  underlying utility function is assumed to be the power utility function in Table <ref> and the power utility function is a subclass of models contained within the GaSP framework with mean basis h(x)=x^α and variance σ=0, we choose the mean function of the GaSP to be misspecified by selecting an inconsistent mean basis (x^0.5), demonstrating that GaSP performs reasonably well even in this scenario.













Tables <ref> and <ref> both indicate that the GaSP method outperforms the other methods for most of the cases. The discrepancy indicated by MSE of testing points by GaSP is usually several orders of magnitude less than that for parametric fitting, and it is usually ten to hundred times better than LI. Note that we only use n=4 and 7 (excluding two boundary points x_min and x_max), which is reasonably small for elicitation purposes, yet the predictive power achieved by GaSP is already realized. The GaSP method with 4 assessed points even performs better than its competitors with 7 assessed points, meaning that the cost for eliciting the utility function can drop significantly using GaSP method.   Moreover, when the sample size increases, the MSE of GaSP decreases, while the MSE for parametric fitting does not necessarily decrease as it is not consistent when the class of parametric utility is misspecified.

The power and exponential families can sometimes be too restrictive, while GaSP is better for these methods because it allows a more flexible structure. The LI and QPD methods also induce a flexible structure. Indeed, they outperform GaSP for two cases when the true utility function is exponential and when the sample size is small (n=4), as the GaSP method uses a misspecified mean function. However, as the sample size increases to n=7 for these two cases, GaSP outperforms the other methods by a lot, meaning that GaSP model converges faster than the LI and QPD methods. A more reasonable mean function would, of course, increase the prediction of the GaSP model. However, as for all cases shown in Tables <ref> and <ref>, the mean function of the GaSP model is misspecified, but the GaSP model still shows the best prediction results for most cases. This feature makes it suitable for a default method for estimation of utility functions.




The fitted curves for the different methods and for 4 cases are plotted in Figure <ref>. In almost all cases, the GaSP-fitted functions (blue curves) are close to the true functions (black curves). The GaSP model has very tight 95% confidence intervals in most regions, meaning that it is confident about the prediction. The confidence interval of the GaSP is large when x∈ [0,4× 10^4] and x∈ [0,2× 10^4] in right panels in Figure <ref>, implying that the uncertainty of the GaSP model is high, and the prediction by the GaSP model has comparatively large discrepancy to the real utility function in these two regions. It is an advantage of the GaSP model that it has an internal assessment of the uncertainty, therefore one can determine whether and where to obtain more assessments given a required level of precision.

In the second scenario, we assume noisy assessments, specifically that utilities are random with additive noise,
 
    u(x)   =U(x) +ϵ
    ϵ|σ^2_ϵ   ∼𝒩(0,σ_ϵ^2)

where U(x) is the underlying utility function and σ_ϵ^2 is the variance of the Gaussian noise. Since the assessed points are noisy, we simulated N=200 experiements to average out the design effect, and calculate the average MSE: AvgMSE=1/N∑_j=1^NMSE_j. Here the total held-out testing points for each case is thus thus n^*× N=200,200.




Compared to the results of noise-free assessments, all methods have comparatively large MSE for noisy assessments as shown in Table <ref>. Indeed, if the elicited utilities are dominated by  noise, none of the methods work as well as for the noise-free cases. Compared to parametric fitting to the rest of the methods, GaSP still has the smallest MSE in almost all cases. 




  §.§.§ Comparing derivatives


In <cit.>, the curvature of utility functions is classified as either concave, convex or of mixed type, using LI or parametric fitting. Such classification is global and characterizes the dominant risk attitude implied by the assessed utility function throughout the  entire domain. In the LI method, empirical derivatives of assessed utility points are utilized for estimation of curvature. Let P_i be an observed utility middle point between two neighboring elicited utility points P_i-1 and P_i+1. Denote S^-(P_i) and S^+(P_i) as slopes of the straight line between  P_i to P_i-1 and  P_i to P_i+1 respectively. Δ S(P_i) = S^+(P_i) -  S^-(P_i) is used as the estimate of convexity at point P_i. The utility function is typically estimated to be concave (convex) if more than ≈ 2/3  of elicited utility points are estimated to be  concave (convex), otherwise it is denoted as a mixed type.

 
 

In the GaSP method, one can compute the posterior probability of the second derivatives for any point x^* in the domain, p(u”(x^*) |u(x_1),…,u(x_n)), as shown in Section <ref>. When p(u”(x^*) ≤ 0|u(x_1),…,u(x_n))>0.5, the utility function is predicted to be concave at x^*, otherwise convex. In the following simulations, we consider n^*=10,000 equally spaced points and average them out to predict curvature in the GaSP model using the proportion of points that are predicted to be concave/convex.  The results for estimating global concavity by LI and GaSP are shown in Figure <ref>.

 

 

In the first row of Figure <ref>, the underlying utility functions are all concave. Due to the effect of noisy assessment, the proportion of concave points by LI is between 1/3 to 2/3 in most of the experiments, meaning that LI fails to identify the concavity of the functions, classifying them of mixed type instead. In the second row of Figure <ref>, when the underlying utility function is convex, LI fails to identify the convexity of the utility functions and again classifies a majority of points as the mixed type.

Compared to LI, GaSP predicts the concavity of the functions a lot more accurately. In the upper row, the proportions of concavity points are almost all close to 1 across  N=500 experiments.  In the second row, most of the points are predicted to be convex for a majority of experiments.  Using 2/3 as the threshold would correctly classify most of the utility functions using GaSP.

There are two main reasons why GaSP performs better. First, the GaSP prediction of concavity of point x^* utilizes information from all assessed tuples rather than just the two neighboring points as in LI. Second, the prediction of concavity by GaSP is averaged over many points (n^*=10,000 chosen here) rather than the limited number of n assessed tuples in LI.

GaSP estimation is therefore better at analyzing preference properties like risk attitude.


















 §.§ Real dataset from <cit.>


Let us compare the parametric and nonparametric methods using a real dataset from <cit.>, collected from a prospect theory based scheme. k=48 people answered a series of questions about comparisons between risky gambles. Due to the effect of loss aversion, the range of the loss domain (negative outcomes) is assumed to be [-1,0] and the range of the gain domain (positive outcomes) is constructed to be [0,0.25], with 11 and 7 assessed tuples in each domain respectively. The design and elicited scheme are described in  <cit.>.




In Figure <ref>, the assessed tuples and fits for the power function and GaSP along with its 95 % posterior credible interval are shown for 2 people. Although the power function is widely used in prospect theory (<cit.>), note that the least squares fit deviates from the assessed curve systematically over a wide range of the domain e.g.  x ∈ [-6× 10^4, -2× 10^4] for the first person.  For the second person, the power function only seems to fit around x=-3 × 10^4. The discrepancy is more likely to be caused by the assumption of the parametric form than the least squares estimation.

In comparison, GaSP seems more flexible for modeling utility functions and fits better in this real dataset.



To test the predictive performance of different methods, we randomly sample n^*_loss=4 and n^*_gain=3 assessed tuples in the loss and the gain domain for each person respectively and save them as the test set, while the remaining assessed tuples are used as the training set. We repeat this experiment N=500 times, and fit each method at each experiment.


We see in Table <ref> that the average out of sample MSE for GaSP is similar to that of LI but much smaller than those for the power and exponential functions.  The parametric fit, albeit straightforward in interpretation, induces extra assumptions on the shape of utility functions and may lead to a poor estimate.







§ MULTI ATTRIBUTE UTILITY FUNCTION ESTIMATION





To demonstrate the performance of the GaSP model for utility functions with multiple attributes, we study a real dataset with three attributes provided in <cit.>. This experiment was conducted among 22 Duke students about decisions pertaining to course selection involving three attributes. The attributes, denoted as 𝐱=(x_1,...,x_3), include the degree of interest, expected teaching quality and the average grade, each of which has 5 levels. The output, denoted as u(𝐱), is the utility rating of a course given a set of attributes.
Each volunteer is asked to rate the same 20 courses along the three attributes.



The RandMAU model is proposed in  <cit.> as a sampling model to capture the stimulus properties of attribute conflict and attribute extremity.
The model is intended to characterize potentially inconsistent responses to  preference elicitation questions through a stochastic model of the assessed utility. It is defined by:

    u(x_1,...,x_p)   =   ∑_i=1^pω_iu_i(x_i),  +ω∑_i=1^p ∑_j>i^p ω_iu_i(x_i)ω_ju_j(x_j)+...+ω^p-1∏_i=1^pω_iu_i(x_i),  
    
          1+ω   =   ∏_i=1^p(1+ωω_i),  
    ω_i    ∼_ind    Beta(r_i, u_i),

where p=3 is the number of attributes and u_i(x_i)=[ x_i-x_i0/x_i^*-x_i0]^α_i with x_i0 and x_i^* as the lower and upper limits for attribute  x_i. This model has 7 parameters (ω_1, ω_2, ω_3,ω, α_1, α_2, α_3) but when (ω_1, ω_2, ω_3) are known, ω can be uniquely solved (as in <cit.>) which leaves the model with 6 degrees of freedom.

<cit.> specify two ways of estimating parameters, namely the corner point and nonlinear least squares estimation methods. The first approach uses only 7 out of the 20 data points (for each participant) to fit the model while the second minimizes the square error using all assessed tuples.



We perform experiments to compare the afore-mentioned two approaches with our proposed GaSP approach for  modeling the average ratings (shown in Table 1 in <cit.>) based on the attributes. We compare the out of sample Means Square Error (MSE) and R^2 of our proposed GaSP estimation with the RandMAU-based methods of corner point and nonlinear least squares estimation, using the 13 data points that are not used in the corner points approach. Since the sample size is very small, each time we only leave a data point out and compute the MSE and R^2 for this point. The average out of sample MSE and R^2 are shown in Table <ref> and Figure <ref>.


In Table <ref>, we find that the prediction by the GaSP model is several times better than the previous inference methods because it is flexible enough to learn the structure of the utility function through data, while RandMAU is restrictive since each u_i(x_i) is assumed to be the power utility in Equation (<ref>) and thus cannot capture other shapes.  For the RandMAU methods, we find that the nonlinear least square estimates are better than the corner point estimates in terms of MSE and R^2. This is because  only 7 observations are used for prediction in corner point estimation, which is statistically inefficient as it simply throws away the rest of the data for fitting the model.


The plot in Figure <ref> shows the prediction of these 13 average ratings by the different methods vs. the real average ratings. Both corner point estimation and nonlinear least square estimation in RandMAU indicate imbalanced and biased predictions. Small values are all over predicted and larger values are under predicted, perhaps because the form in RandMAU is unable to capture the curvature and latent pattern of the data; the right most figure corresponding to GaSP seems to solve these difficulties.
























§ CONCLUSIONS


We have presented a nonparametric Bayesian approach involving GaSP for inferring a decision maker's utility function using assessed tuples as training data, regardless of the choice of elicitation protocol and underlying theory. We have argued theoretical benefits over parametric approaches around the desired property of interpolation for noise-free assessment. Unlike simpler approaches like linear interpolation, the proposed approach provides for differentiability whenever a differentiable correlation function is chosen.

The Bayesian aspect of our approach provides a probabilistic view of the decision maker's preference properties, such as the utility function and its risk aversion, while the nonparametric aspect provides the flexibility to fit a large class of functions.

Our numerical experiments with simulated, noise-free assessment data show that the GaSP model has lower MSE when compared with parametric fitting from a misspecified function class as well as the quantile parametrized distribution method and linear interpolation.

When the simulated data is corrupted with additive Gaussian noise, GaSP continues to enjoy the least MSE, though all approaches lead to comparatively higher MSE than the noise-free case.

We also demonstrated promising results with real data sets from the literature. These results provide evidence of the effectiveness of using Gaussian processes with respect to out of sample MSE, as well as their flexibility to model empirical utility assessment data.
While the proposed approach is only marginally better than linear interpolation for the single-attribute empirical data set, it provides a principled approach to estimate risk aversion coefficients using the posterior predictive view of the derivative processes. Analysis with the multi-attribute dataset indicates that the proposed estimation method is promising, but further experiments are needed to demonstrate effective estimation in higher dimensions as compared to existing approaches.



There are some unique issues that arise when applying Gaussian processes for estimating utility functions. For instance, a utility function is typically assumed to be monotonic in its arguments.

A potential avenue for future work is around constrained GaSPs, which would extend the proposed approach to respect monotonicity. Even without such considerations, GaSP estimation appears to out-perform the most popular state-of-the-art techniques.




§ APPENDIX


Proof of Lemma <ref>

The mean and the covariance of  the Wiener process are [W_t]=0 for any t≥ 0 and (W_s, W_t)=min(s,t) for any s, t≥ 0.  For any 0<t_i<t_*<t_i+1<0, the joint distribution of ( W_t_i ,  W_t_*,  W_t_i+1  )^T is
   
    ( [   W_t_i;   W_t_*; W_t_i+1;         ])   ∼𝒩( ([ 0; 0; 0;   ]),  ( [   t_i   t_i   t_i;   t_i   t_*   t_*;   t_i   t_* t_i+1;       ]) ),


By  properties of the multivariate distribution, the conditional distribution of W_t_*| W_t_i W_t_i+1 is

    W_t_*| W_t_i W_t_i+1∼𝒩(μ_*, V_*  ),

where μ_*=(t_i+1-t_*)W_t_i+ (t_*-t_i)W_t_i+1/t_i+1-t_i and V_*=(t_i+1-t_*)(t_*-t_i)/t_i+1-t_i. By the Markov property of the process,

    p(W_t_*| W_t_1,W_t_2,...,W_t_n)   =p(W_t_*| W_t_i,W_t_i+1,...,W_t_n) 
       = p(W_t_*, W_t_i,W_t_i+1,...,W_t_n)/p(W_t_i,W_t_i+1,...,W_t_n) 
       = p(W_t_*, W_t_i,W_t_i+1,...,W_t_n)/p(W_t_i+2,W_t_i+3,...,W_t_n| W_t_i, W_t_i+1   ) p(W_t_i,W_t_i+1)
       = p(W_t_*, W_t_i,W_t_i+1,...,W_t_n)/p(W_t_i+2,W_t_i+3,...,W_t_n|W_t_i,W_t_*, W_t_i+1) p(W_t_i,W_t_i+1)
       =p(W_t_i,W_t_*, W_t_i+1)/p(W_t_i,W_t_i+1)
       =p(W_t_*|W_t_i, W_t_i+1).

The result follows.


Proof of Lemma <ref>
For any i, note that when the prediction is at the design points 𝐱_i

    𝐜^T(𝐱^𝒟_i) 𝐂^-1 = 𝐞^T_i,

where  𝐞_i is an n dimensional vector with the i^th row as 1 and the others as 0. Thus

    û (𝐱^𝒟_i )=   𝐡(𝐱^𝒟_i )θ̂+𝐜^T(𝐱^*)𝐂^-1(𝐮^𝒟  -𝐡(𝐱^𝒟  )θ̂) 
    
    =   𝐡(𝐱^𝒟_i )θ̂ +  𝐮(𝐱^𝒟_i)- 𝐡(𝐱^𝒟_i  )θ̂
    
    =   𝐮(𝐱^𝒟_i).

































































































ba ]