Bayesian Preference Learning Tutorial

Executive Summary: Bayesian preference learning uses probabilistic models to infer underlying
“utilities” or scores of items from observed pairwise comparisons. Core ideas include Bayes’ theorem
(posterior ∝ likelihood × prior), probabilistic utility models (e.g. Bradley–Terry logistic or Thurstone–
Mosteller probit), and Bayesian inference (MAP, MCMC, variational inference). We introduce the
necessary mathematics (probability, Bayes’ rule, Gaussian processes), derive key formulas, and explain
intuition simply. We work through numerical examples and provide clear Python code: first a minimal
from-scratch implementation of a Bradley–Terry model (with MAP and MCMC using PyMC), then a
preference-learning example using a library (a Gaussian-process approach). Each code block is fully
annotated and easy to run. We compare inference methods and libraries in tables, and illustrate results
with plots: e.g. posterior distributions of item utilities, preference probability curves, and (conceptual)
learning curves. Throughout we cite foundational sources (standard Bayesian texts and tutorials) and
explain every concept assuming only basic calculus and Python knowledge.

Background Mathematics
Probability and Bayes’ Theorem: A basic notion is that of conditional probability and Bayes’ rule. If θ
represents unknown parameters (e.g. item utilities) and D is observed data (pairwise comparisons),
Bayes’ theorem states

P (D ∣ θ) P (θ)
P (θ ∣ D) = ,
P (D)

i.e. the posterior P (θ∣D) is proportional to the likelihood P (D∣θ) times the prior P (θ) (with P (D)
ensuring normalization)【46†L309-L317】. In words, we update our prior beliefs about θ by how well
they explain the data (via the likelihood). The posterior encodes all our information about θ after seeing
the data. In practice we often omit the normalizing constant and write P (θ∣D) ∝ P (D∣θ) P (θ).

The likelihood P (D∣θ) is the probability of the observed comparisons under parameters θ. The prior
P (θ) expresses any pre-data beliefs (e.g. a Gaussian prior expressing that utilities are likely moderate).
The posterior then combines them. For example, in a simple case with a continuous parameter θ, the
posterior density is

f (θ) L(D ∣ θ)
f (θ ∣ D) = ,
∫ f (u) L(D ∣ u) du

where f (θ) is the prior density and L(D∣θ) = P (D∣θ) the likelihood【46†L309-L317】. We may often
work with log-probabilities to turn products into sums (the log-likelihood plus log-prior).

Gaussian Processes (GPs): Gaussian processes are widely used as nonparametric priors over functions.
Informally, a GP defines a distribution over functions f (x) so that any finite collection of outputs has a
multivariate normal distribution【40†L144-L153】. Equivalently, a GP is specified by a mean function
m(x) (often taken 0) and a covariance (kernel) function k(x, x′ ). For inputs x1 , … , xn , the vector
(f (x1 ), … , f (xn )) is jointly Gaussian: f (x) ∼ N (m(x), K) where Kij = k(xi , xj )【40†L144-
L153】. This makes GPs convenient for regression: given data we update to a posterior GP that

1
interpolates the data. In preference learning, GPs can serve as priors on latent utility functions, giving
flexible (nonlinear) models of item utilities【36†L306-L309】【33†L268-L276】.

<p align="center">【80†embed_image】 Figure: A logistic sigmoid curve (blue) mapping utility differences
to preference probability.</p>

Preference Models and Likelihoods

A common framework is the Random Utility Model (RUM): each item i has a latent utility or score ui .
When two items i, j are compared, the model assumes noisy utilities (e.g. ui + ϵi and uj + ϵj ) and that
i is preferred if its noisy utility is larger. The noise assumption leads to familiar probabilistic models
【36†L279-L288】.

• Bradley–Terry (BT) model (Logistic): Each item i has a positive strength pi . In the simplest BT
parameterization, the probability that i beats j is

pi
P (i ≻ j) = .
pi + pj

Equivalently, letting βi = log pi (so pi = eβi ), this becomes a logistic form:

eβi 1
P (i ≻ j) = = .
e i +e j
β β 1 + e i −βj )
−(β

P (i≻j)
Hence the log-odds log P (j≻i) = βi − βj . Thus BT is a logistic (sigmoid) model over utility
differences【22†L309-L317】. Intuitively, if βi is much larger than βj , then P (i ≻ j) will be near
1; if they are equal, the probability is 0.5. In practice we often set one β to 0 for identifiability
(only differences matter).

• Thurstone–Mosteller (Probit) model: Thurstone’s model assumes Gaussian noise on utilities.

Concretely, if the difference of true utilities is δ = ui − uj , and each has i.i.d. normal noise, the
probability P (i≻ j) = Φ(δ/σ 2), where Φ is the standard normal CDF. Thus P (i ≻ j) =
Φ(βi − βj ) if we absorb constants into β . Like BT, it depends only on the difference βi − βj , but
now through a probit (Gaussian) link. In fact, BT and Thurstone are special cases of Generalized
Linear Models for preferences: BT uses a logistic link (coming from Gumbel noise) whereas
Thurstone uses Gaussian noise (probit)【36†L279-L288】. Both produce similar S-shaped curves.

These models handle pairwise binary outcomes. For multiple comparisons, the likelihood of observed
data is the product over pairs: e.g. if item i beats j wij times out of nij contests, the likelihood
contribution is P (i ≻ j)wij (1 − P (i ≻ j))nij −wij . We will focus on the logistic BT form for concreteness
(easy interpretation and implementation).

<p align="center">【89†embed_image】 Figure: A typical sigmoid (“S”) curve (red) mapping differences in
utility to win probability (similar shape for probit).</p>

2
Bayesian Inference Methods
Having specified a likelihood P (D∣θ) (with parameters θ , e.g. all the βi ’s), and a prior P (θ), we seek the
posterior P (θ∣D). Exact calculation is often infeasible (no closed-form), so we use approximate
inference:

• Maximum a Posteriori (MAP) Estimation: MAP finds the mode (maximum) of the posterior
density. Equivalently, it maximizes log P (θ∣D) = log P (D∣θ) + log P (θ). If the prior is
uniform, MAP reduces to maximum likelihood estimation (MLE). In Bayesian stats the MAP is
simply the most probable parameter value under the posterior. As Wikipedia notes, MAP “is the
mode of the posterior density” and effectively regularizes MLE by the prior【49†L193-L200】.
Practically, one can use optimization (e.g. gradient ascent or [Link] ) on the log-
posterior to find MAP estimates. MAP is fast and gives a point estimate, but does not quantify
full uncertainty.

• Markov Chain Monte Carlo (MCMC): MCMC methods draw samples from the posterior
distribution. Classic algorithms include Metropolis–Hastings and Hamiltonian Monte Carlo (used
by Stan/PyMC). The sampler constructs a Markov chain whose equilibrium distribution is
P (θ∣D), so after “burn-in” we obtain correlated samples θ(1) , … , θ(N ) approximately
distributed from the true posterior. These samples allow us to estimate any posterior quantity
(means, intervals, predictive checks) and fully capture uncertainty. As Stan’s manual states, “Stan
uses MCMC techniques to generate draws from the posterior distribution for full Bayesian
inference.”【56†L111-L119】. MCMC can handle complex posteriors, but is computationally
intensive (especially as data or model complexity grows).

• Variational Inference (VI): VI turns inference into optimization. We posit a simple variational
family of approximate posteriors q(θ∣ϕ) (e.g. mean-field Gaussian) and choose parameters ϕ to
make q close to the true posterior (by minimizing KL divergence). The result is an analytic
approximation to P (θ∣D). Stan/PyMC implement algorithms (e.g. ADVI) that output “draws”
from the fitted variational posterior【56†L116-L124】. VI is usually faster than MCMC and scales
to larger data, but it can underestimate uncertainty and its accuracy depends on the choice of
variational family. For Gaussian-process preference models, sparse variational GPs are a popular
scalable approach.

Below is a concise comparison of inference methods:

Method Output Complexity Pros/Cons Use-Cases

Quick
Low; like Fast, simple; ignores
estimates,
MAP optimization, uncertainty
Point estimate large-scale
(optimization) often quadratic (overconfident)
settings,
costs in #params 【49†L193-L200】
initialization

High; each Exact (in limit),

Small/medium
sample ~cost of quantifies uncertainty
MCMC Samples from data, full
likelihood; can 【56†L111-L119】;
(sampling) posterior Bayesian
be parallel slow convergence,
inference
chains tuning needed

3
 Method Output Complexity Pros/Cons Use-Cases

Fast, scalable;
Large data,
Moderate; approximation error
Approximate models
VI optimization of (often
posterior amenable to
(approximation) variational underestimates
(parametric) mean-field
parameters variance)【56†L116-
structure
L124】

<div markdown="1"> Inference Example (Bradley–Terry): For a simple BT model with 3 items A,B,C,
suppose we observe wins A≫B, A≫C, B≫C. The likelihood factors into Bernoulli trials. In a Bayesian fit,
each βi might have a Gaussian prior. The log-posterior is

log P (β∣D) = ∑ [wij log σ(βi − βj ) + (nij − wij ) log(1 − σ(βi − βj ))] + ∑ log P (βi ),
pairs (i,j) i

where σ(x) = 1/(1 + e−x ). We could maximize this (MAP) or sample from it (MCMC). In practice one
often fixes one β to 0 (gauge) and estimates the others. The fitted β reveal a ranking: the highest β is
the most preferred item. </div>

Derivations and Worked Examples

Bradley–Terry Derivation: Starting from positive strengths pi , pj , define P (i ≻ j) = pi /(pi + pj ).
Equivalently let βi = log pi . Then

eβi 1
P (i ≻ j) = = ,
e i +e j
β β 1 + e i −βj )
−(β

a sigmoid function【22†L309-L317】. This makes clear why only differences βi − βj matter. One can
also derive the log-likelihood for one comparison: ℓ(β) = log P (i ≻ j) = βi − log(eβi + eβj ). Taking
gradient (for MAP/MLE) gives the famous update: if item i beats j , the gradient w.r.t. βi is 1 − P (i ≻ j)
and w.r.t. βj is −1 + P (i ≻ j)【76†L173-L180】. Intuitively, if i wins, we “nudge” βi up (and βj down)
by how surprising that win was.

Thurstone Derivation: Assume ui = u∗i + ϵi , uj = u∗j + ϵj , where ϵi , ϵj ∼ N (0, σ 2 ). Then P (i ≻

j) = P (ui > uj ) = P (ϵj − ϵi < u∗i − u∗j ). But ϵj − ϵi ∼ N (0, 2σ 2 ). So

u∗i − u∗j
P (i ≻ j) = Φ( ),
2σ
where Φ is the CDF of N (0, 1). Setting βi = u∗i /( 2σ) yields P (i ≻ j) = Φ(βi − βj ). This is the
probit link. Thus Thurstone’s model gives a similar S-curve to BT; indeed the logistic and probit curves
are very close (differ by <0.01 for many x)【25†L185-L192】.

Bayes Update (example): Suppose prior P (βi ) is N (0, τ 2 ). Given data D , the posterior density is

P (β∣D) ∝ [∏ P (D∣β)] ∏ exp(−βi2 /(2τ 2 )).

One could write out the unnormalized log-posterior and maximize it. For example with 2 items (only

4
β1 − β2 matters) and 10 comparisons where item 1 won 7 times, the log-likelihood contribution is
7 log σ(Δ) + 3 log(1 − σ(Δ)) with Δ = β1 − β2 . Adding −β12 /(2τ 2 ) − β22 /(2τ 2 ) from the Gaussian
prior completes the objective to maximize.

Numerical Example: Consider 3 items with “true” (unknown) utilities β ∗ = [2, 1, 0]. We generate 10
pairwise contests for each pair {A,B}, {A,C}, {B,C} based on the logistic model above, e.g.

− βB = 1, P (A ≻ B) = 0.73.
• A vs B: if βA
• A vs C: βA − βC = 2, P (A ≻ C) = 0.88.
• B vs C: βB − βC = 1, P (B ≻ C) = 0.73.
Randomly sampling gives, say, 8 out of 10 A>B wins, 9/10 A>C, 7/10 B>C.
Using this data and a flat prior, we could compute the posterior (via MCMC) or MAP estimate of
β . The MAP might be close to [2,1,0] (modulo shift).
Such computations are done in the implementation below.

Implementation in Python
Below we provide two implementations: (1) a minimal implementation of a pairwise-preference model
(Bradley–Terry) using basic Python and PyMC for inference; and (2) an example using a higher-level
library (Gaussian-process-based) for preference learning. Each code block is fully annotated, with simple
explanations and instructions.

(1) Minimal Bradley–Terry Model with MAP and MCMC

We implement the BT model for n items. The code first sets up example data, then builds a Bayesian
model in PyMC. We use a simple Gaussian prior on each β , then define the likelihood of observed wins
using the logistic function. We find the MAP estimate and draw posterior samples (via MCMC) to
illustrate both methods.

import numpy as np
import pymc as pm

# Example data: 3 items A,B,C (indices 0,1,2)

# Suppose we observe comparisons: A beats B 8 times, B beats A 2 times, etc.
wins = []
player1 = []
player2 = []
# A vs B: A wins 8, B wins 2
for _ in range(8):
[Link](0); [Link](1); [Link](1) # 1 means player1
wins
for _ in range(2):
[Link](0); [Link](1); [Link](0) # 0 means player1
loses (so B wins)
# A vs C: A wins 9, C wins 1
for _ in range(9):
[Link](0); [Link](2); [Link](1)
for _ in range(1):
[Link](0); [Link](2); [Link](0)
# B vs C: B wins 7, C wins 3

5
 for _ in range(7):
[Link](1); [Link](2); [Link](1)
for _ in range(3):
[Link](1); [Link](2); [Link](0)

player1 = [Link](player1, dtype=int)

player2 = [Link](player2, dtype=int)
wins = [Link](wins, dtype=int)

- We list each comparison as (player1, player2, result). player1[i] and player2[i] are indices of
the two items in the i-th trial; wins[i] =1 if player1 won, 0 if player2 won. In this example, item 0=A,
1=B, 2=C.

with [Link]() as bt_model:

# Prior on utilities (one beta per item)
beta = [Link]('beta', 0, sigma=5, shape=3)
# Logistic win probability for each observed comparison:
# sigmoid(beta[player1] - beta[player2])
p = [Link]('p', [Link](beta[player1] - beta[player2]))
# Likelihood: each comparison result is Bernoulli with prob p[i]
obs = [Link]('obs', p=p, observed=wins)
# Find MAP (optional)
map_est = pm.find_MAP()
print("MAP estimate of beta:", map_est['beta'])
# Run MCMC sampling
trace = [Link](draws=1000, tune=1000, chains=2, random_seed=1,
progressbar=False)

- We use a wide Normal(0,5) prior on each utility β. The probability that item player1[i] beats
player2[i] is the logistic function of the difference of their betas ( sigmoid(beta[p1] -
beta[p2]) ). We observe each outcome via [Link] . pm.find_MAP() performs
optimization to find the posterior mode (MAP); [Link]() draws from the posterior via MCMC.
- Dependencies: This requires PyMC (v4). Install with pip install pymc .
- Running: Save this code in a script or Jupyter cell and run. It will print the MAP estimate of β and
return a trace with posterior samples.

For example, the printed MAP might be something like [2.1, 1.0, 0.2] (depending on data and
random seed), reflecting A>B>C. The trace object can be summarized to get means/intervals.

# Posterior summary
beta_samples = [Link]['beta']
mean_beta = beta_samples.mean(dim=['chain','draw']).values
hpd_beta = [Link](beta_samples, hdi_prob=0.95).values
for i in range(3):
print(f"Item {i} beta: mean={mean_beta[i]:.2f}, 95% CI=({hpd_beta[0,i]:.
2f},{hpd_beta[1,i]:.2f})")

6
- This prints the posterior mean and 95% credible interval for each β. We expect item0’s β highest,
item2’s lowest.

# Optional: visualize posterior distributions (requires matplotlib)

import [Link] as plt
for i in range(3):
[Link](beta_samples[:,:,i].stack(draws=("chain","draw")), bins=30,
alpha=0.5, label=f'β_{i}')
[Link]('Utility β')
[Link]('Frequency')
[Link]()
[Link]('Posterior distributions of item utilities')
[Link]()

- Expected output: Three overlapping bell-shaped histograms, one for each β (as illustrated conceptually
above). The highest curve corresponds to the item most often winning.

Explanation: This code implements the BT model and uses Bayesian inference. The MAP step
( find_MAP ) solves for β maximizing the posterior, while [Link] uses MCMC to draw samples
(the console will show sampling progress). The final prints and plot summarize the posterior beliefs
about utilities.

(2) Library-based Preference Learning (Gaussian Process)

For a more advanced model, we can use a library. One example is GPro (a GP-based preference learning
library) or GPyTorch with custom likelihoods. Here we sketch using GPro’s ProbitPreferenceGP
(Gaussian-process with probit link) from its documentation【33†L268-L276】. (In practice you’d install
via pip install GPro and then:)

# Example using GPro library (pseudo-code, install GPro for real use)
import numpy as np
from [Link] import ProbitPreferenceGP

# Suppose X is an array of item features (here just 1D positions for

illustration),
# and prefs is list of (i,j) indices where item i was preferred over j.
X = [Link]([[0.0],[1.0],[2.0]]) # 3 items on the real line
prefs = [(0,1),(0,2),(1,2)] # A>B, A>C, B>C (for simplicity)

model = ProbitPreferenceGP()
[Link](X, prefs) # fit the GP utility model
mu, var = [Link](X) # posterior means/variances of utilities

print("Posterior means of utility:", mu)

- This code uses GPro (Gaussian Process for preferences) as an example. It assumes item features X
(could be any descriptors; we used 1D positions). We input observed preferences prefs as index
pairs. fit learns the GP model, and predict gives the posterior mean utility for each item.

7
Dependencies: Install with pip install GPro . It relies on SciPy and NumPy.
Explanation: This approach infers a nonparametric utility function via a GP prior. The
ProbitPreferenceGP assumes a probit (Gaussian noise) likelihood for comparisons【33†L268-
L276】. It automatically computes a posterior over utilities. GPro’s documentation shows similar
examples (with plots of GP predictive curves)【33†L323-L332】.

# Plot preference probability curves (for one latent dimension):

import [Link] as plt
from [Link] import norm
delta = [Link](-3,3,200)
[Link](delta, 1/(1+[Link](-delta)), label='Logistic')
[Link](delta, [Link](delta), label='Probit (Gaussian)')
[Link]('Utility difference (β_i - β_j)')
[Link]('P(i preferred to j)')
[Link]()
[Link]('Preference probability vs utility difference')
[Link]()

- Expected output: An S-shaped (sigmoid) curve for both logistic and probit links (yellow vs red curves like
above). This illustrates that as the utility difference grows, preference probability smoothly goes from 0
to 1.

Comparison of Inference Methods and Libraries

Inference /
Type Use-Case Complexity Pros Cons
Library

Fast, simple;
Quick Low (convex gives single No uncertainty
MAP
Optimization point problem for best-fit values quantification;
(analytical)
estimate logit) 【49†L193- can overfit
L200】

Exact (samples)
posterior;
MCMC Full Computationally
High (esp. measures
(PyMC, Sampling Bayesian intensive;
large data) uncertainty
Stan) inference requires tuning
【56†L111-
L119】

Scalable, fast;
Large- Approximation
Moderate approximate
VI (PyMC, Approx. scale error; depends
(optimize posterior
GPyTorch) sampling Bayesian on variational
VI) 【56†L116-
models family
L124】

8
 Inference /
Type Use-Case Complexity Pros Cons
Library

Easy model
spec; many
Flexible likelihoods; Can be slow for
PyMC Moderate/
General Bayes built-in large GPs;
(Python) High
modeling samplers overhead
【56†L111-
L119】

Fast HMC Very fast MCMC Less mature

NumPyro Moderate/
General (JAX) on GPU/ via JAX; supports ecosystem; JAX
(Python) High
TPU MCMC/VI learning curve

Fast
Not specialized
GP with approximate
GPyTorch for preference
GP modeling custom High (O(n³)) GPs; GPU;
(PyTorch) (need custom
kernels flexible GP
link)
constructs

Direct
implementation
Limited to GP
GPro (GPro Preference Moderate of GP
GP probit; less
Preference with GP & (SciPy- preference
preference community
GP) probit based) models
support
【33†L268-
L276】

No direct
Scikit- Baseline Very fast;
Ranking (via preference API;
Learn ranking Low familiar
features) needs feature
(Logistic) model interface
design

Notes: MAP and MCMC are methods, not libraries, while PyMC, NumPyro, GPyTorch, etc. are tools
implementing them. GPs have cubic complexity in data size, so for large datasets one uses sparse or
approximate methods (e.g. inducing points with VI). PyMC/NumPyro allow easy MCMC/VI for custom
models (including BT or probit). GPro provides a high-level GP preference model (probit link) out of the
box【33†L268-L276】.

Visualizations
We already saw illustrative plots above:

• Posterior distributions of utilities: The histogram figure (in code example) shows each item’s
posterior. Items with higher latent utility have distributions shifted right.
• Preference probability curves: The sigmoid plots compare logistic vs Gaussian (probit) links.
Both smoothly map utility difference to preference probability (as shown by the yellow and red
curves above).
• Learning curves (conceptual): In preference learning, one could plot predictive accuracy or
posterior certainty vs number of comparisons. In general, more data shrinks posterior variance
and improves ranking accuracy. For example, with few comparisons the blue logistic curve would

9
 be “flatter” (more uncertain), and as data accumulate it steepens around true differences (not
shown here).

Further Reading
• Bradley–Terry Model: See the original papers (Bradley & Terry 1952) or Marden’s text【36†L279-
L288】. A clear introduction is Moran (2026) “Learning from pairwise preferences”【76†】.
• Gaussian Processes: The Stan manual provides a good exposition【40†L144-L153】. The tutorial
by Benavoli & Azzimonti (2024) comprehensively covers GP preference models【36†L306-L309】.
• Bayesian Inference: Gelman et al., Bayesian Data Analysis (3rd ed.) for fundamentals; and Bishop’s
Pattern Recognition for MAP/MCMC【46†L309-L317】【49†L193-L200】.
• Software: PyMC ([Link] and NumPyro ([Link] for general Bayesian
modeling; GPyTorch ([Link] for GPs; GPro ([Link] for GP
preference models.

10