text:
A Theory-Based Evaluation of Nearest Neighbor Models Put Into Practice

In the $k$-nearest neighborhood model ($k$-NN), we are given a set of points $P$, and we shall answer queries $q$ by returning the $k$ nearest neighbors of $q$ in $P$ according to some metric. This concept is crucial in many areas of data analysis and data processing, e.g., computer vision, document retrieval and machine learning. Many $k$-NN algorithms have been published and implemented, but often the relation between parameters and accuracy of the computed $k$-NN is not explicit. We study property testing of $k$-NN graphs in theory and evaluate it empirically: given a point set $P \subset \mathbb{R}^\delta$ and a directed graph $G=(P,E)$, is $G$ a $k$-NN graph, i.e., every point $p \in P$ has outgoing edges to its $k$ nearest neighbors, or is it $\epsilon$-far from being a $k$-NN graph? Here, $\epsilon$-far means that one has to change more than an $\epsilon$-fraction of the edges in order to make $G$ a $k$-NN graph. We develop a randomized algorithm with one-sided error that decides this question, i.e., a property tester for the $k$-NN property, with complexity $O(\sqrt{n} k^2 / \epsilon^2)$ measured in terms of the number of vertices and edges it inspects, and we prove a lower bound of $\Omega(\sqrt{n / \epsilon k})$. We evaluate our tester empirically on the $k$-NN models computed by various algorithms and show that it can be used to detect $k$-NN models with bad accuracy in significantly less time than the building time of the $k$-NN model.

Introduction
The k-nearest neighborhood (k-NN) of a point q with respect to some set of points P is one of the most fundamental concepts used in data analysis tasks such as classification, regression and machine learning. In the past decades, many algorithms have been proposed in theory as well as in practice to efficiently answer k-NN queries [e.g., 1, 7-10, 17, 19, 22, 25, 26, 29, 35]. For example, one can construct a k-NN graph of a point set P , i.e., a directed graph G = (P, E) of size n = |P | such that E contains an edge (p, q) for every k-nearest neighbor q of p for every p ∈ P , in time O(n log n + kn) for constant dimension δ [8]. Due to restrictions on computational resources, approximations and heuristics are often used instead (see, e.g., [9,10] and the discussion therein for details). Given the output graph G of such a randomized approximation algorithm or heuristic, one might want to check whether G resembles a k-NN graph before using it, e.g., in a data processing pipeline. However, the time required for exact verification might cancel out the advantages gained by using an approximation algorithm or a heuristic. On the other hand, testing whether G is at least close to a k-NN graph will suffice for many purposes. Property testing is a framework for the theoretical analysis of decision and verification problems that are relaxed in favor of sublinear complexity. One motivation of property testing is to fathom the theoretical foundations of efficiently assessing approximation and heuristic algorithms' outputs.
Property testing [31], and in particular property testing of graphs [20], has been studied quite extensively since its founding. A one-sided error -tester for a property P of graphs with average degree bounded by d has to accept every graph G ∈ P and it has to reject every graph H that is -far from P with probability at least 2/3 (i.e., if graphs that are -far are relevant, it has precision 1 and recall 2/3). A graph H of size n is -far from some property P if more than dn edges have to be added or removed to transform it into a graph that is in P. A two-sided error -tester may also err with probability less than 1/3 if the graph has the property. The computational complexity of a property tester is the number of adjacency list entries it reads, denoted its queries. Many works in graph property testing focus on testing plain graphs that contain only the pure combinatorial information. However, most graphs that model real data contain some additional information that may, for example, indicate the type of an atom, the bandwidth of a data link or spatial information of an object that is represented by a vertex or an edge, respectively. In this work, we consider geometric graphs with bounded average degree. In particular, the graphs are embedded into R δ , i.e., every vertex has a coordinate x ∈ R δ . The coordinate of a vertex may be obtained by a query.

Main Results
Our first result is a property tester with one-sided error for the property that a given geometric graph G with bounded average degree is a k-nearest neighborhood graph of its underlying point set (i.e., it has precision 1 and recall 2/3 when taking -far graphs as relevant). Theorem 1. Given an input graph G = (V, E) of size n = |V | with bounded average degree d, there exists a one-sided error -tester that tests whether G is a k-nearest neighbourhood graph. It has query complexity c · √ n k 2 ψ δ / 2 , where ψ δ is the δ-dimensional kissing number and c > 0 is a universal constant.
We emphasize that it is not necessary to compute the ground truth (i.e., the k-NN of P ) in order to run the property tester. Furthermore, the tester can be easily adapted for graphs G = (P ∪ Q, E) such that P ∩ Q = ∅ and we only require that for every q ∈ Q, E contains an edge (q, p) for every k-nearest neighbor p of q in P . This is more natural when we think of P as a training set and Q as a test set or query domain. To complement this result, we prove a lower bound that holds even for two-sided error testers. Theorem 2. Testing whether a given input graph G = (V, E) of size n = |V | is a k-nearest neighbourhood graph with one-sided or two-sided error requires max( n/(8 k), kψ δ /6) queries.
Finally, we provide an experimental evaluation of our property tester on approximate nearest neighbor (ANN) indices computed by various ANN algorithms. Our results indicate that the tester requires significantly less time than the ANN algorithm to build the ANN index, most times just a 1 /10-fraction. Therefore, it can often detect badly chosen parameters of the ANN algorithm at almost no additional cost and before the ANN index is fed into the remaining data processing pipeline.

Related Work
We give an overview of sublinear algorithms for geometric graphs, which is the topic of research that is most relevant to our work. As mentioned above, the research on k-NN algorithms is very broad and diverse. See, e.g., [15,32] for surveys. Testing whether a geometric graph that is embedded into the plane is a Euclidean minimum spanning tree has been studied by Ben-Zwi et al. [3] and Czumaj and Sohler [11]. In [3], the authors show that any non-adaptive tester has to make Ω( √ n) queries, and that any adaptive tester has query complexity Ω(n 1/3 ). In [11], a one-sided eror tester with query complexityÕ( n/ ) is given. In a fashion similar to property testing, Czumaj et al. [14] estimate the weight of Euclidean Minimum Spanning Trees inÕ( √ n · poly( )) time, and Czumaj and Sohler [12] approximate the weight of Metric Minimum Spanning Trees inÕ(n · poly( )) time for constant dimension, respectively. Hellweg et al. [21] develop a tester for Euclidean (1 + δ)-spanners. Property testers for many other geometric problems can, for example, be found in [13,28].

Preliminaries
Let d, δ, , k ≥ 0, d ≥ k be fixed parameters. In this paper, we consider property testing on directed geometric graphs with bounded average degree d.
Definition 1 (geometric graph). A graph G = (V, E) with an associated function coord : V → R d is a geometric graph, where each vertex v is assigned a coordinate coord(v). Given v ∈ V , we denote its degree by d(v) and the set of adjacent vertices N(v) := {u | (v, u) ∈ E}.
The Euclidean distance between two points x, y is denoted by dist(x, y). For the sake of simplicity, we write dist(u, v) := dist(coord(u), coord(v)) for two vertices u, v ∈ V . When there is no ambiguity, we also refer to coord(v) by simply writing v. We denote the size of the graph G = (V, E) at hand by n = |V |.
Definition 2 (k-nearest neighborhood graph). A geometric graph G = (V, E) is a k-nearest neighbourhood (k-NN) graph if for every v ∈ V , the k points u 1 , . . . , u k ∈ V that lie nearest to v according to dist(·, ·) are neighbors of v in G, i.e., (v, u i ) ∈ E for all i ∈ [k] (breaking ties arbitrarily).
Let G = (V, E) be a geometric graph. We say that a graph G is -far from a geometric graph property P if at least dn edges of G have to be modified in order to convert it into a graph that satisfies the property P. We assume that the graph G is represented by a function , a degree function d G : V → N that outputs the degree of a vertex and a coordinate function coord G : V → R δ that outputs the coordinates of a vertex.
Definition 3 ( -tester). A one-sided (error) -tester for a property P with query complexity q is a randomized algorithm that makes q queries to f G , d G and coord G for a graph G. The algorithm accepts if G has the property P. If G is -far from P, then it rejects with probability at least 2/3. The motivation to consider query complexity is that the cost of accessing the graph, e.g., through an ANN index, is costly but cannot be influenced. Therefore, one should minimize access to the graph.
} define the subset of knn(v) that is not adjacent to v. If wit(v) = ∅ or d(v) < k, we call v incomplete, and we call elements of wit(v) the witnesses of v.
If G is -far from being a k-nearest neighborhood graph, an -fraction of its vertices are incomplete.
Lemma 5. If G is -far from being a k-nearest neighborhood graph, at least dn/(2k) vertices are incomplete.
Proof. Assume the contrary. For every incomplete vertex v, delete k edges such that the distance to the property does not increase and insert the missing edges from v to its k nearest neighbors. By the assumption, the total number of inserted or deleted edges is less than dn/(2k) · 2k. Therefore, G is -close to being a k-nearest neighborhood graph.
The main challenge for the property tester will be to find matching witnesses for a fixed set of incomplete vertices. The following result from coding theory for Euclidean codes bounds the maximum number of points q i that can have the same fixed point p as nearest neighbor.

Upper Bound
The idea of the tester is as follows (see Algorithm 1). Two samples are drawn uniformly at random: S , which shall contain many incomplete vertices if G is -far from being a k-nearest neighborhood graph and T , which shall contain at least one witness of an incomplete vertex in S . For every v ∈ S , the algorithm should query its degree, its coordinate as well as every adjacent vertex and their coordinates and calculate the distance to them. If d(k) < k or if one of the vertices in T is a witness of v, the algorithm found an incomplete vertex, and hence rejects. Otherwise, it accepts.
However, we have to deal with the case that some vertices in S have non-constant degree, say, Ω(1/ ), such that querying all their adjacent vertices would require too many queries. To this end, we prove that one can prune these vertices to obtain a subset S ⊆ S of low degree vertices that still contains many incomplete vertices with sufficient probability. T ← sample ln(10) · k · ψ δ · √ n vertices from V u.a.r. with replacement; Proof of Theorem 1 We prove that Algorithm 1 is an -tester as claimed by Theorem 1. Since Algorithm 1 does never reject a k-nearest neighbourhood graph, assume without loss of generality that G = (V, E) is -far from being a k-nearest neighborhood graph. Algorithm 1 only queries the neighbors of S, and therefore its query complexity is at most |S| · 100k/ = O( √ nk 2 / 2 ). It remains to prove the correctness.
In the following, let L := {v ∈ V | d(v) ≤ 100k/ } denote the set of all vertices in G that have low degree, let I denote the set of incomplete vertices in L, and let I S ⊆ S denote the set of incomplete vertices in S. By an averaging argument, |V \L| ≤ dn/(100k). It follows from Lemma 5 that L contains at least dn/(4k) incomplete vertices, and therefore we focus on finding incomplete vertices that have low degree. The following random variable identifies witnesses of vertices incomplete vertices in S ⊂ L.
Definition 7. Given u ∈ T , let W S (u) be a random variable that is 1 if u is a witness of an incomplete vertex v ∈ I S and 0 otherwise.
The proof of Theorem 1 follows from the following three claims. First, note that S is a uniform sample without replacement from L whose size |S| is random. However, |S| is sufficiently large with constant probability.
Proof. The expected cardinality of S \S is d √ n. Therefore, the probability that |S| is less than 20 √ n/ is at most 1/10 by Markov's inequality.
In the subsequent sections, we prove the following two claims. Given that S is sufficiently large, it will contain at least √ n incomplete vertices with constant probability.
Finally, we show that if S contains at least √ n incomplete vertices, then T will contain at least one witness of such an incomplete vertex with constant probability.
The correctness follows by a union bound over these three bad events.
Analysis of the Sample S: Proof of Claim 9 We bound the cardinality of S such that S contains at least √ n incomplete vertices.
Proof. Since S was sampled without replacement, the random variable |I S | follows the hypergeometric distribution. Let X be a random variable that denotes the number of draws that are needed to obtain √ n incomplete vertices in S, which therefore follows the negative hypergeometric distribution. By Lemma We apply Markov's inequality to obtain Analysis of the Sample T: Proof of Claim 10 We prove the following lower bound on the number of witnesses in G, which will imply a bound on |T | by k-reducing it to the case k = 1.
Proposition 12. Given a point set P ⊂ R δ , p ∈ P and k ∈ N, the maximum number of points q i ∈ P that can have p as k-nearest neighbor is bounded by k · ψ δ .
We note that this bound is tight, as shown in Lemma 15.
( * ) Pick a point q ∈ Q that lies furthest from w and let Proof of Proposition 12. We apply Definition 13 to p and prove that the size of Q at the beginning of the process is at most k · ψ δ , which proves the claim. At first we show that every vertex that is picked by ( * ) stays in Q: Let q 1 , q 2 be arbitrary points that are picked by ( * ) in the process of k-reducing, with q 1 being picked in an earlier iteration than q 2 . The latter implies dist(q 2 , p) < dist(q 1 , p). Assume that q 1 ∈ Q q 2 at the time q 2 is selected, and therefore q 1 is removed from Q. Since q 1 is deleted by (#), it holds that dist(q 1 , p) < dist(q 2 , p), which is a contradiction as q 1 has been selected before q 2 .
We continue to bound the maximum number of vertices that share their k-nearest neighbor: Because p is the nearest point for the remaining q ∈ Q, we apply Lemma 6 and conclude that at most ψ δ vertices are remaining in Q after k-reducing. Since every iteration of step (#) removed at most k − 1 points from Q, the cardinality of Q at the beginning of the process was at most Since at most k · ψ δ vertices can share a witness by Proposition 12, there are at least |I S | k·ψ δ distinct witnesses of vertices in S. We employ this bound to calculate the size of the sample T such that it contains at least one witness of an incomplete vertex in S with constant probability. Lemma 14. If |S| ≥ 10 √ n and |T | ≥ ln(10) · k · ψ δ · √ n, then Pr w∈T W S (w) = 0 ≤ 1 10 .
Proof. Since every vertex is sampled uniformly at random with replacement, the event that one vertex is a witness is a Bernoulli trial with probability Pr w∈V [W S (w) = 1] ≥ |I S | k·ψ δ · 1 n = |I S | k·ψ δ ·n . Therefore Pr w∈V [W S (w) = 0] ≤ 1 − |I S | k·ψ δ ·n . We have By Claim 9, Eq. (1) holds for |S| as chosen in Algorithm 1. In Eq. (2) we use the fact that 1 − x ≤ e −x and in Eq. (5) we use that all events W S (u) = 0 for u ∈ T are independent Bernoulli trials.
Finally, we observe that the factor k that is introduced in Proposition 12 is tight.
Lemma 15. For every δ ≥ 3, k ≥ 2, there exists a point set P ⊂ R δ such that there is a set of kψ δ points q i ∈ P that have the same k-nearest neighbor.

Lower Bound
We prove the first lower bound by constructing two (distributions of) graphs that are composed of multiple copies of the same building block. All graphs in one distribution are k-nearest neighborhood graphs, and all graphs in the other distribution are -far from the property. It suffices to show that no deterministic algorithm that makes o( √ n) queries can distinguish these two distributions with sufficiently high probability. Our building block is defined as follows.
Note that a line gadget is a k-nearest neighborhood graph itself. In the following, let k := k + 1. The graphs in the first distribution D 1 are composed of n/k line gadgets with sufficiently large pair-wise distances that maintain the k-nearest neighborhood property. The construction of the distribution of -far graphs D 2 is a bit more complicated. Basically, we want to move n/k line gadgets to the exact position of n/k other line gadgets such that in the resulting graph, n/k pairs of line gadgets share the same coordinates. However, we have to make sure that the algorithm is oblivious of this relocation with sufficiently high probability. Lemma 17. Testing whether a graph is a k-nearest neighborhood graph with two-sided error requires n/(8 k ) queries.
Proof. For the sake of simplicity, let n be a multiple of k . Let G = (V, E) be a graph that is composed of n/k line gadgets L 3k i for i ∈ [n/k ], and let D 1 be the uniform distribution over all vertex labellings of G. All graphs in D 1 are k-nearest neighborhood graphs. For every graph G ∈ D 1 , we define a random graph G as follows. Let S(G) = (s i ) i∈[2 n/k ] be a sequence of random numbers, drawn without replacement from the uniform distribution over [n]. For every i ∈ n/k , move the line gadget L 3k s i to the coordinates of L 3k s i+ n/k such that G contains no L 3k s i but two L 3k s i+ n/k afterwards. Note that G is -far from being a k-nearest neighborhood graph. Let D 2 be the uniform distribution over G∈D 1 G .
We may assume that if the tester queries for (a neighbor of) a vertex v, then the oracle returns the whole line gadget that v belongs to. This is only beneficial for the query complexity of the algorithm. We consider the knowledge graph of the algorithm, which is defined as the subgraph of the input graph that consists of the vertices, edges and non-edges that are revealed by the answers to the queries. Without loss of generality, we may assume that the algorithm only asks queries whose answers cannot be deduced from the knowledge graph. To prove the theorem, it is sufficient to show that for any deterministic algorithm that makes o( √ n) queries, the distributions of knowledge graphs that are obtained from distributions D 1 and D 2 , respectively, have small statistical distance.
First, note that every query of the algorithm to a graph from D 1 reveals a line gadget that is not in the knowledge graph yet. Let R be the set of revealed vertices after the i-th query q i . Then, the probability that q i+1 reveals a line gadget L 3k j is k /(n − k i) if L 3k j / ∈ R and 0 otherwise. We turn to D 2 now. Consider the i-th query q i . We claim that conditioned on the event that q 1 , . . . , q i+1 do not reveal two instances of a line gadget that is contained twice in the graph, then the line gadget that is revealed by q i+1 is distributed uniformly. Given G ∈ D 2 , let S j (G) be the j-th half of S(G) for j ∈ {1, 2} such that S(G) is the concatenation of S 1 (G) and S 2 (G), and let S 2 (G) denote the subset of S 2 (G) that has already been revealed. Let D 2 be the restriction of D 2 to the graphs that are compatible with the current knowledge graph, and the let be the support size of D 2 . We denote the event that a line gadget L 3k j is revealed by E j . Fix some arbitrary j. We have Therefore, the knowledge graph distribution of D 2 is exactly the same as the distribution of D 1 , i.e., uniform over all graphs that are compatible with current knowledge graph, as long as the algorithm does not reveal two line gadgets with the same coordinates. The probability that two line gadgets with the same coordinates are revealed by the first b = n/(8 k ) queries is upper bounded by the probability that for any pair of queries (i, j) ∈ [b] 2 , query i hits one of the k vertices in one of the n/k line gadgets L 3k s + n/k such that ∈ n/k , times the probability that query j hits the same gadget. In particular, by the union bound this probability is at most n k · k n · k n ≤ b 2 k n ≤ 1/4. Therefore, the total variation distance between the knowledge graph distributions is at most 1/4. In property testing, it is common to fix problem specific parameters such as the dimension and analyze the asymptotic behavior with respect to n and . However, it may be interesting that the computational complexity of a tester for k-nearest neighborhood graphs is at least linear in ψ δ . Lemma 18. Testing whether a graph of size n is a k-nearest neighborhood graph with two-sided error requires at least kψ δ /6 queries.
Proof sketch. Let P (q) be a point set as constructed in the proof of Lemma 15 with the common k-nearest neighbor located at q. Without loss of generality, assume that n = (c + 1) · |P | for some c ∈ N. Let q i = (i, 0, . . . , 0) ∈ R δ for i ∈ [c]. By scaling P accordingly, we can construct a set R = P (q 1 ) ∪ . . . ∪ P (q c ) such that for every i ∈ [c] and every p ∈ R ∩ P (q i ) \ {q i }, q i ∈ Q is the k-nearest neighbor of p. Let G be the k-nearest neighborhood graph of R.
We construct two graphs, G and H, and prove that one requires at least kψ δ /6 queries to distinguish between uniform distributions over all labellings of G and H, respectively. For i ∈ [ 2 c ] and some sufficiently small δ, we move q i to q i = (i − δ, 0, . . .) and insert a new point q i = (i + δ, 0, . . . , 0). Now, at least half of the vertices from P (q i ) have one of these points as k-nearest neighbor. Without loss of generality, let q i be this point. We obtain G by applying this modification to G . We obtain H from G by moving q i for every i ∈ [ 2 c ] to (−1, 0, . . . , 0) and recalculating the k-nearest neighbors for all these points. It follows that H is a k-nearest neighborhood graph, while G is -far from being a k-nearest neighborhood graph.
By the union bound, the probability to sample a point from Applying the union bound once again, it follows that the first kψ δ /6 queries to G or H will not contain any point from Q with probability greater than 1/3.

Experiments
As discussed above, property testing aims at distinguishing perfect objects and objects that have many flaws at very small cost. Given the output of an approximate nearest neighbor (ANN) algorithm, a natural use case for a property tester is to decide whether the nearest neighbor index computed by the ANN algorithm is accurate or resolves many queries incorrectly.
Although Algorithm 1 already gives values for the sizes of |S| and |T |, one would probably want to minimize the running time of the tester beyond worst-case analysis in practice. When used as a tool to assess an ANN index before actually putting it to work, it is also important that the tester actually reduces the total computation time compared to observing poor results at the end of the data processing pipeline (e.g., bad classification results) and starting over. Therefore, we seek to answer the following questions: Setup We implemented our property tester in C++ and integrated it into the Python framework ANN-Benchmarks [2,4]. The key-idea of ANN-Benchmarks is to compare the quality of the indices built by ANN implementations, with respect to their running times and query-answer times. To evaluate our property tester, we chose three algorithms with the best performance observed in [2]: KGraph [1] and hnsw and SW-graph from the Non-Metric Space Library [7,27]. All of the ANN algorithms are implemented in C / C++ and build upon nearest neighbor / proximity graphs. We computed the ground truth, i.e., a k-NN graph of the input data, for the Euclidean datasets MNIST (size 60 000, dimension 960, [24]), Fashion-MNIST (size 60 000, dimension 960, [30]) and SIFT (size 1 000 000, dimension 128, [18]) to evaluate the answers of the tester.
We ran our benchmarks on identical machines with 60 GB of free RAM guaranteed and an Intel Xeon E5-2640 v4 CPU running at 2.40 GHz (capable of running 20 concurrent threads) and measured CPU time. To minimize interference between different processes, a single instance of an ANN algorithm was run exclusively on one machine at a time.
The C++ source code of the property tester that was used for the experiments is available here [16]. The modified version of ANN-Benchmarks is available here [6].
Q1: Parameterization of the Property Tester We analyze how different choices for |S | and |T | in Algorithm 1 affect which quality of ANN indices the tester is likely to reject. All ANN algorithms were run ten times for each choice of parameters built into ann-benchmarks (as listed in [5]) and every dataset. Then the tester was run once for each output and for every choice from {0.001, 0.01, 0.1} × {0.05, 0.5, 5} for (c 1 , c 2 ) in |S | = c 1 · 8k √ n and |T | = c 2 · k √ n log(10), with oracle access to the resulting ANN index. We chose to evaluate the tester for k = 10 because indices that are very close to 10-NN graphs -which is the hard case for the property tester to detect -can be computed by the ANN algorithms in reasonable time, and we support this decision by an additional experiment for k = 50. The ground truth, i.e., a k-NN graph of each dataset, and the -distance (see Section 2) of each ANN index to a k-NN graph was computed offline.
We evaluate the recall of the property tester by distance of a tested ANN index to ground truth, where graphs that are no k-NN graphs are relevant (note that the tester always provides a witness when it rejects, so its precision is 1). Since the quality of an ANN index varies depending on the ANN algorithm's parameters and internal randomness, we group the computed ANN indices into buckets according to their distance to ground truth and depict the resulting recall on these classes in Fig. 1 for all datasets combined and for each dataset individually. As the oracle access that is provided to the property tester is oblivious of the underlying ANN algorithm, the figures show the combined results for all algorithms.
We observe that for distances and parameters that result in a reasonable overall recall, say, at least greater than 0.75, the property tester behaves comparable on all datasets. Since the property tester is guaranteed to have precision 1, even parameterizations with low recall on a small distance can be amplified by running the tester multiple times, possibly for different values of c 1 , c 2 . In summary, after choosing a target distance that the property tester should detect, the tested parameters seem suitable for data with dimensions up to roughly 800. For higher dimensions, it is likely advisable to apply dimensionality reduction techniques first before computing and using nearest neighbors in Euclidean space.
To get an indication of how the tester behaves for larger k, we conducted an additional experiment where we ran the tester on KGraph indices with k = 50. As one might expect, less indices are close to being a 50-NN graph than a 10-NN graph for the same sets of KGraph parameters (although the distance is normalized by k and therefore it allows more errors), but the results indicate that it is also easier for the property tester to spot errors. This suggests that, at least for KGraph, errors are spread quite uniformly in the index rather than they are concentrated on some vertices.
Q2: Performance Consider the following scenario: an algorithm that processes data employs an ANN index. The quality of the algorithm's result (e.g., the classification rate) depends on the quality of the ANN index. However, the best parameters for the ANN algorithm are not known, and conclusions about the quality of the ANN index can only be drawn by looking at the algorithm's final result, which may be a long costly way to go. Does it pay out to run the property tester on the ANN index and recompute the index using different parameters if the tester rejects? We address this question by measuring the tester's performance. However, whether to use the tester or not depends heavily on the cost incurred otherwise. Therefore, we compare the property tester against the minimum cost that every algorithm that uses an ANN index must invest before it can employ it or even just draw conclusions about its quality: the build time of the index. Figure 3 shows the time required by the property tester normalized (divided) by the time required to build the ANN index for each ANN algorithm and each dataset. There are two plots: one for all graphs that are between 0.005 and 0.01-close to a 10-NN graph, and one for all graphs that are between 0.01 and 0.02-close to a 10-NN graph.
In general, the running time of the property tester is always smaller than the build time for hnsw and SW-graph and at most five times the build time for KGraph. Mostly, it is even smaller than 1 /10 of the build time, and therefore running the property tester comes at almost no additional cost. For the runs of the tester on KGraph indices with k = 50, the testing time is also upper bounded by five times the build time and the tester time vs. build time ratio is 0.1 for k = 10 (restricted to MNIST and Fashion-MNIST) and k = 50.

Conclusion
We have studied the task of efficiently identifying NN models with low accuracy by exploring possibilities within the theoretical framework of sublinear algorithms and evaluated our approach by moving to experiments. In particular, we have proved that there is a one-sided error property tester with complexity O( √ nk 2 / 2 ), i.e., a sublinear (randomized) algorithm that decides whether an input graph G is a k-NN graph or requires many edge modifications to become a k-NN graph (i.e., precision 1 and recall 2/3 when taking -far graphs as relevant). We also proved that even a two-sided error property tester requires complexity Ω( n/( k)). Our experiments of the property tester on ANN indices computed by various algorithms indicate that testing comes at almost no additional cost, i.e., the testing time is significantly smaller than the building time of the ANN index that is tested.
From the perspective of applications, it would be desirable to analyze the tester for a more context sensitive notion of edit distance. For example, an edge to the (k + 1)-nearest neighbor of a point instead of an edge to its k-nearest neighbor might be a defect that is much less severe than an edge to the k 2 -nearest neighbor. It would be interesting to investigate what results can be obtained under established oracle access models, which are oblivious of the graph's structure, and whether other useful models can be devised.


extfieldsofstudy:
Computer Science
Mathematics

text:
Differentially Private Bayesian Inference for Exponential Families

The study of private inference has been sparked by growing concern regarding the analysis of data when it stems from sensitive sources. We present the first method for private Bayesian inference in exponential families that properly accounts for noise introduced by the privacy mechanism. It is efficient because it works only with sufficient statistics and not individual data. Unlike other methods, it gives properly calibrated posterior beliefs in the non-asymptotic data regime.

Introduction
Differential privacy is the dominant standard for privacy [1]. A randomized algorithm that satisfies differential privacy offers protection to individuals by guaranteeing that its output is insensitive to changes caused by the data of any single individual entering or leaving the data set. An algorithm can be made differentially private by applying one of several general-purpose mechanisms to randomize the computation in an appropriate way, for example, by adding noise calibrated to the sensitivity of the quantity being computed, where sensitivity captures how much the quantity depends on any individual's data [1]. Due to the obvious importance of protecting individual privacy while drawing population level inferences from data, differentially private algorithms have been developed for a broad range of machine learning tasks [2][3][4][5][6][7][8][9].
There is a growing interest in private methods for Bayesian inference [10][11][12][13][14]. In Bayesian inference, a modeler selects a prior distribution p(θ) over some parameter, observes data x that depends probabilistically on θ through a model p(x | θ), and then reasons about θ through the posterior distribution p(θ | x), which quantifies updated beliefs and uncertainty about θ after observing x. Bayesian inference is a core machine learning task and there is an obvious need to be able to conduct it in a way that protects privacy when x is sensitive. Additionally, recent work has identified surprising connections between sampling from posterior distributions and differential privacy-for example, a single perfect sample from p(θ | x) satisfies differential privacy for some setting of the privacy parameter [10][11][12][13].
An "obvious" way to conduct private Bayesian inference is to privatize the computation of the posterior, that is, to design a differentially private algorithm A that outputs y = A(x) with the goal that y ≈ p(θ | x) is a privatized representation of the posterior. However, using y directly as "the 32nd Conference on Neural Information Processing Systems (NIPS 2018), Montréal, Canada. posterior" will not correctly quantify beliefs, because the Bayesian modeler never observes x, she observes y; her posterior beliefs are now quantified by p(θ | y). This paper will take a different approach to private Bayesian inference by designing a pair of algorithms: The release mechanism A computes a private statistic y = A(x) of the input data; the inference algorithm P computes p(θ | y). These algorithms should satisfy the following criteria: • Privacy. The release mechanism A is differentially private. By the post-processing property of differential privacy [15], all further computations are also private. • Calibration. The inference algorithm P can efficiently compute or approximate the correct posterior, p(θ | y) (see Section 4 for our process to measure calibration). • Utility. Informally, the statistic y should capture "as much information as possible" about x so that p(θ | y) is "close" to p(θ | x).
Importantly, the release mechanism A is public, so the distribution p(y | x) is known. Williams and McSherry first suggested conducting inference on the output of a differentially private algorithm and showed how to do this for the factored exponential mechanism [16]; see also [17][18][19][20].
Our work focuses specifically on Bayesian inference when the private data X = x 1:n is an iid sample of (publicly known) size n from an exponential family model p(x i | θ). Exponential families include many of the most familiar parametric probability models. We will adopt a straightforward release mechanism where the Laplace mechanism [1] is used to release noisy sufficient statistics y [12,19], which are a finite-dimensional quantity that capture all the information about θ [21].
The technical challenge is then to develop an efficient general-purpose inference algorithm P. One challenge is computational efficiency. The exact posterior p(θ | y) ∝ p(θ)p(x 1:n | θ)p(y|x 1:n )dx 1:n integrates over all possible data sets [16], which is intractable to do directly for large n. We integrate instead over the sufficient statistics s, which have fixed dimension and completely characterize the posterior; furthermore, since they are a sum over individuals, p(s | θ) is asymptotically normal. We develop an efficient Gibbs sampler that uses a normal approximation for s together with variable augmentation to model the Laplace noise in a way that yields simple updates [22].
A second challenge is that the sufficient statistics may be unbounded, which makes their release incompatible with the Laplace mechanism. We address this by imposing truncation bounds and only computing statistics from data that fall within the bounds. We show how to use automatic differentiation and a "random sum" central limit theorem to compute the parameters of the normal approximation p(s | θ) for a truncated exponential family when the number of individuals that fall within the truncation bounds is unknown.
Our overall contribution is the pairing of an existing simple release mechanism A with a novel, efficient, and general-purpose Gibbs sampler P that meets the criteria outlined above for private Bayesian inference in any univariate exponential family or multivariate exponential family with bounded sufficient statistics. 1 We show empirically that when compared with competing methods, ours is the only one that provides properly calibrated beliefs about θ in the non-asymptotic regime, and that it provides good utility compared with other private Bayesian inference approaches.

Differential Privacy
Differential privacy requires that an individual's data has a limited effect on the algorithm's behavior. In our setting, a data set X = x 1:n := (x 1 , . . . , x n ) consists of records from n individuals, where x i ∈ R d is the data of the ith individual. We will assume n is known. Differential privacy reasons about the hypothesis that one individual chooses not to remove their data from the data set, and their record is replaced by another one. 2 Let nbrs(X) denote the set of data sets that differ from X by exactly one record-i.e., if X ∈ nbrs(X), then X = (x 1:i , x i , x i+1:n ) for some i. 1 There are remaining technical challenges for multivariate models with unbounded sufficient statistics that we leave for future work. 2 This variant assumes n remains fixed, which is sometimes called bounded differential privacy [23].
Definition 1 (Differential Privacy; Dwork et al. [1]). A randomized algorithm A satisfies -differential privacy if for any input X, any X ∈ nbrs(X) and any subset of outputs O ⊆ Range(A), We achieve differential privacy by injecting noise into statistics that are computed on the data. Let f be any function that maps datasets to R d . The amount of noise depends on the sensitivity of f .
We drop the subscript f when it is clear from context. Our approach achieves differential privacy through the application of the Laplace mechanism. Definition 3 (Laplace Mechanism; Dwork et al. [1]). Given a function f that maps data sets to R m , the Laplace mechanism outputs the random variable L(X) ∼ Lap f (X), ∆ f / from the Laplace distribution, which has density Lap(z; u, b) = (2b) −m exp (− z − u 1 /b). This corresponds to adding zero-mean independent noise u i ∼ Lap(0, ∆ f / ) to each component of f (X).
A final important property of differential privacy is post-processing [15]; if an algorithm A isdifferentially private, then any algorithm that takes as input only the output of A, and does not use the original data set X, is also -differentially private.

Private Bayesian Inference in Exponential Families
We consider the canonical setting of Bayesian inference in an exponential family. The modeler posits a prior distribution p(θ), assumes the data x 1:n is an iid sample from an exponential family model p(x | θ), and wishes to compute the posterior p(θ | x 1:n ). An exponential family in natural parameterization has density where η are the natural parameters, t(x) is the sufficient statistic, A(η) = h(x) exp η T t(x) dx is the log-partition function, and h(x) is the base measure. The density of the full data is where h(x 1:n ) = n i=1 h(x i ) and t(x 1:n ) = n i=1 t(x i ). Notice that once normalizing constants are dropped, this density is dependent on the data only directly through the sufficient statistics, s = t(x 1:n ).
We will write exponential families more generally as p(x | θ) to indicate the case when the natural parameters η = η(θ) depend on a different parameter vector θ.
Every exponential family distribution has a conjugate prior distribution p(θ; λ) [24] with hyperparameters λ. A conjugate prior has the property that, if it is used as the prior, then the posterior belongs to the same family, i.e., p(θ | x 1:n ; λ) = p(θ; λ ) for some λ that depends only on λ, n, and the sufficient statistics s. We write this function as λ = Conjugate-Update(λ, s, n); our methods are not tied to the specific choice of conjugate prior, only that the posterior parameters can be calculated in this form. See supplementary material for a general form of Conjugate-Update.

Release Algorithm: Noisy Sufficient Statistics
If privacy were not a concern, the Bayesian modeler would simply compute the sufficient statistics s = t(x 1:n ) and use them to update the posterior beliefs. However, to maintain privacy, the modeler must access the sensitive data only through a randomized release mechanism A. As a result, in order to obtain proper posterior beliefs the modeler must account for the randomization of the release mechanism by performing inference.
We take the simple approach of releasing noisy sufficient statistics via the Laplace mechanism, as in [12,13,19]. Sufficient statistics are a natural quantity to release. They are an "information bottleneck"-a finite-dimensional quantity that captures all the relevant information about θ. The released value is y = A( When t(·) is unbounded this quantity becomes infinite; we will modify the release mechanism so the sensitivity is finite (Sec. 3.3).

Basic Inference Approach: Bounded Sufficient Statistics
The goal of the inference algorithm P is to compute p(θ | y). We first develop the basic approach for the simpler case when t(x) is bounded, and then extend both A and P to handle the unbounded case. The full joint distribution of the probability model can be expressed as: where p(θ) = p(θ; λ) is a conjugate prior and the goal is to compute a representation of p(θ | y) ∝ s p(θ, s, y)ds by integrating over the sufficient statistics. We will develop a Gibbs sampler to sample from this distribution. There are two main challenges. First, the distribution p(s | θ) is obtained by marginalizing over the data sample x 1:n , and is usually not known in closed form. We will address this with an asymptotically correct normal approximation. Second, when resampling s within the Gibbs algorithm, we require the full conditional distribution of s given the other variables, which is proportional to p(s|θ)p(y | s). Care must be taken to make it easy to sample from this conditional distribution. We address this via variable augmentation. We discuss our approach to both challenges in detail below.
Normal approximation of p(s | θ). The exact form of the sufficient statistic distribution p(s | θ) is obtained by marginalizing over the data: In general, the exact form of this distribution is not available. In some cases, it is-for example if x ∼ Bernoulli(θ) then s ∼ Binomial(n, θ)-but even then it may not lead to a tractable full conditional for s.
Properties of exponential families pave the way toward a general approach that always leads to a tractable full conditional. By the central limit theorem (CLT), because s = i t(x i ) is a sum of iid random variables, it is asymptotically normal. It can be approximated as are the mean and variance of the sufficient statistic of a single individual. This approximation is asymptotically correct [25]. The quantities µ and Σ can be computed using well-known properties of exponential families [25]: where η = η(θ) is the natural parameter.
Note that we will not use this approximation for Gibbs updates of θ. Instead, we will compute the conditional p(θ | s) using standard conjugacy formulas. In this sense, we maintain two views of the joint distribution p(θ, s)-when updating θ, it is the standard exponential family model, which leads to conjugate updates; when updating s, it is approximated as p(θ)N (s; nµ, sΣ), which will lead to simple updates when combined with a variable augmentation technique.
Variable augmentation for p(y | s). We seek a tractable form for the full conditional of s under the normal approximation, which is the product of a normal density and a Laplace density: A similar situation arises in the Bayesian Lasso [22], and we will employ the same variable augmentation trick. A Laplace random variable z ∼ Lap(u, b) can be written as a scale mixture of normals by introducing a latent variable σ 2 ∼ Exp(1/(2b 2 )), i.e., the distribution with density 1/(2b 2 ) exp −σ 2 /(2b 2 ) and letting z ∼ N (u, σ 2 ). We apply this separately to each dimension of the vector y so that: s ∼ NormProduct nµ, nΣ, y, diag(σ 2 ) 6: The Gibbs Sampler. After the normal approximation and variable augmentation, the generative process is as shown to the right. The final Gibbs sampling algorithm is shown in Algorithm 1. Note that the update for θ is based on conjugacy in the exact distribution p(θ, s), while the update for s uses the density of the generative process to the right, so that p(s | θ, σ 2 , y) ∝ p(s | θ) p(y | σ 2 , s), which is a product of two normal densities where µ s and Σ s are are defined in Algorithm 1 [26]. The update for σ 2 follows Park and Casella [22]; the inverse Gaussian density is InverseGaussian( . Full derivations are given in the supplement.

Unbounded Sufficient Statistics and Truncated Exponential Families
The Laplace mechanism does not apply when the sufficient statistics are unbounded, because ∆ s = max x,y t(x) − t(y) 1 = ∞. Thus, we need a new release mechanism A and inference algorithm P.
We present a solution for the case when x is univariate. All elements of the solution can generalize to higher dimensions, except that one step will have running time that is exponential in d; we leave improvement of this to future work and focus on the simpler univariate case.
Release mechanism. Our solution is to truncate the support of the (now univariate) p(x | θ) to x ∈ [a, b], where a and b are finite bounds provided by the modeler. If the modeler cannot select bounds a priori, they may be selected privately as a preliminary step using a variant of the exponential mechanism (see PrivateQuantile in Smith [27]). 3 Then, given truncation bounds, the data owner redacts individuals where x i / ∈ [a, b] and reports the truncated sufficient statisticŝ is the indicator function of the set S. The sensitivity ofŝ is now ∆ŝ = max x,y∈R t (x) −t(y) 1 wheret(x) = 1 [a,b] (x) t(x). An easy upper bound for this quantity (see supplement) is: where t j (x) is the jth component of the sufficient statistics. The bounds [a, b] will be selected so this quantity is bounded. The released value is y ∼ Lap(ŝ, ∆ŝ/ ).
Inference: Truncated Exponential Family. Several new challenges arise for inference. The quantityŝ is no longer a sufficient statistic for the model p(x | θ), and we will need new insights to understand p(ŝ | θ) and p(θ |ŝ). Sinceŝ is a sum over individuals where x i ∈ [a, b], it will be useful to examine the probability of the event x ∈ [a, b] as well as the conditional distribution of x given this event. To facilitate a general development, assume a generic truncation interval [v, w], not necessarily equal to [a, b]. Let F (x; θ) = x −∞ p(x | θ)dx be the CDF of the original (univariate) exponential family model. It is clear that Pr(x ∈ [v, w]) = F (w; θ) − F (v; θ). The conditional distribution of x given x ∈ [v, w] is a truncated exponential family, which, in its natural parameterization is: Note that this is still an exponential family model (with a modified base measure), and all of the standard results apply, such as the existence of a conjugate prior and the formulas in Eq. (1) for the mean and variance of t(x) under the truncated distribution.
Random sum CLT for p(ŝ | θ). We would like to again apply an asympotic normal approximation forŝ, but we do not know how many individuals fall within the truncation bounds. The "random sum CLT" of Robbins [28] applies to the setting where the number of terms in the sum is itself a random variable. The sum can be rewritten asŝ = Proof. Each term of the sum has meanμ and varianceΣ, and the number of terms is N ∼ Binom(n, q). The result follows from Robbins [28].
Computingμ andΣ by automatic differentiation (autodiff). To use the normal approximation we need to computeμ andΣ. Lemma 1. Let p(x | θ) be a univariate exponential family model and letp(x | θ) be the corresponding exponential family model truncated to generic interval [v, w]. Then Proof. It is straightforward to derive from Eq. (2) thatÂ(η) = A(η) + log F (w; η) − F (v; η) . The result follows from applying Eq. (1) to this expression forÂ(η). See the supplement for derivation of A(η) and proof of this lemma.
We will use Equations (3) and (4)  When x is multivariate, analogous expressions can be derived forμ andΣ. The adjustment factors will include multivariate CDFs, with a number of terms that grow exponentially in d. This is currently the main limitation in applying our methods to multivariate models with unbounded sufficient statistics.
Conjugate updates for p(θ |ŝ). The final issue is the distribution p(θ |ŝ), which is no longer characterized by conjugacy becauseŝ are not the full sufficient statistics. We again turn to variable Conditioned on all other variables, each component is multivariate normal, so the sum s is also multivariate normal. We can therefore sample s and then sample from p(θ | s) using conjugacy. We will also need to drawŝ c separately to be used to update σ 2 .
The Gibbs Sampler. The (approximate) generative process in the unbounded case is: The Gibbs sampler to sample from this distribution is given in Algorithm 2. Note that in Line 8 we employ rejection sampling in which sufficient statistics are sampled until the values drawn are valid for the given data model, e.g., s must be positive for the binomial distribution. The RS-CLT algorithm to compute parameters of the random sum CLT is shown in Algorithm 3.

Experiments
We design experiments to measure the calibration and utility of our method for posterior inference.
We conduct experiments for the binomial model with beta prior, the multinomial model with Dirichlet prior, and the exponential model with gamma prior. The last model is unbounded and requires truncation; we set the bounds to keep the middle 95% of individuals, which is reasonable to assume known a priori for some cases, such as modeling human height.
Methods. We run our Gibbs sampler for 5000 iterations after 2000 burnin iterations (see supplementary material for convergence results), which we compare to two baselines. The first method uses the same release mechanism as our Gibbs sampler and performs conjugate updates using the noisy sufficient statistics [12,13]. This method converges to the true posterior as n → ∞ because the Laplace noise will eventually become negligible compared to sampling variability [12]. However, the noise is not negligible for moderate n; we refer to this method as "naive". For truncated models we allow the naive method to "cheat" by accessing the noisy untruncated sufficient statistics s. Thus the method is not private, and receives strictly more information than our Gibbs sampler, but with the same magnitude noise. This allows us to demonstrate miscalibration without highly technical modifications to the baseline method to be able to deal with truncated sufficient statistics.
The second baseline is a version of the one-posterior sampling (OPS) mechanism [11][12][13], which employs the exponential mechanism [29] to release samples from a privatized posterior. We release 100 samples using the method of [12], each with ops = /100, such that the entire algorithm achieves -differential privacy. Private MCMC sampling [11] is a more sophisticated method to release multiple samples from a privatized posterior and could potentially make better use of the privacy budget; however, private MCMC will also necessarily be miscalibrated, and only achieves the weaker privacy guarantee of ( , δ)-differential privacy for δ > 0, so would not be direct comparable to our method. OPS serves as a suitable baseline that achieves -differential privacy. We include OPS only for experiments on the binomial model, for which it requires the support of θ to be truncated to [a 0 , 1 − a 0 ] where a 0 > 0. We set a 0 = 0.1.
We also include a non-private posterior for comparison, which performs conjugate updates using the non-noisy sufficient statistics.
Evaluation. We evaluate both the calibration and utility of the posterior. For calibration we adapt a method of Cook et al. [30]: the idea is to draw iid samples (θ i , x i ) from the joint model p(θ)p(x | θ), and conduct posterior inference in each trial. Let F i (θ) be the CDF of the true posterior p(θ | x i ) in trial i. Then we know that In other words, the actual parameter θ i is equally likely to land at any quantile of the posterior. To test the posterior inference procedure, we instead compute U i as the quantile at which θ i lands within a set of samples from the approximate posterior. After M trials of the whole procedure we test for uniformity of U 1:M using the Kolmogorov-Smirnov goodnessof-fit test [31], which measures the maximum distance between the empirical CDF of U 1:M and the uniform CDF; lower values are better and zero corresponds to perfect uniformity. We also visualize the empirical CDFs to assess calibration qualitatively.
Higher utility of a private posterior is indicated by closeness to the non-private posterior, which we measure with maximum mean discrepancy (MMD), a kernel-based statistical test to determine if two sets of samples are drawn from different distributions [32]. Given m i.i.d. samples (p, q) ∼ P × Q, an unbiased estimate of the MMD is where k is a continuous kernel function; we use a standard normal kernel. The higher the value the more likely the two samples are drawn from different distributions.
Results. Figure 1a shows the results for three models and varying n and . Our method (Gibbs) achieves the same calibration level as non-private posterior inference for all settings. The naive method ignores noise and is too confident about parameter values implied by treating the noisy sufficient statistics as true ones; it is only well-calibrated with increasing n and when noise becomes negligible relative to population size. OPS is not calibrated because it samples from an over-dispersed version of p(θ | x). Figure 1b shows the empirical CDF plots for n = 1000 and = 0.01. Our method and the non-private method are both perfectly calibrated. The naive method's over-confidence in the wrong sufficient statistics causes its posterior to usually be too tight at the wrong value; thus the true parameter always lies in a tail of the approximate posterior, so too much mass is placed near 0 and 1. OPS shows the opposite behavior: its posterior is always too diffuse, so the true parameter lies close to the middle. For multinomial we show measures only for the parameter of the first category, but results hold for all categories. Figure 1c shows the MMD test statistic between each method and the non-private posterior, used as a measure of utility. Our method consistently achieves utility at least as good as the naive method for binomial and multinomial models. We omit OPS, which is never calibrated. For the exponential model (not shown) we did not obtain conclusive utility comparisons due to the lack of a naive baseline that properly handles truncation; the "cheating" naive method from our calibration experiments sometimes attains higher utility than our method, and sometimes lower, but this comparison is not meaningful because it receives strictly more information.

Conclusion
We presented a Gibbs sampling approach for private posterior inference in exponential family models. Rather than trying to approximate the posterior of p(θ | x 1:n ), we divide our procedure into a private release mechanism y = A(x 1:n ) and an inference algorithm P that computes p(θ | y). The release mechanism is designed to facilitate inference. We develop a general-purpose Gibbs sampler that applies to any exponential family model that has bounded sufficient statistics; a truncated version applies to univariate models with unbounded sufficient statistics. The Gibbs sampler uses general properties of exponential families to approximate the distribution of the sufficient statistics, and therefore avoids the need to reason about individuals. Promising lines of future work are to develop efficient methods for multivariate exponential families with unbounded sufficient statistics, and to develop methods for conditional models based on exponential families, such as generalized linear models.

A Properties of Exponential Families
A.1 Form of Conjugate-Update(λ, x 1:n ) Following Diaconis and Ylvisaker [24], the prior is where the parameters are λ = [λ 1 , λ 2 ] and sufficient statistics are [η, −A(η)] The posterior after observing x 1:n is Proof: The proof for Varp[t(x)] is similar.

B Derivation of σ 2 Gibbs update
We fully derive the Gibbs update for the noise variance σ 2 of the augmented model as stated in Park and Casella [22]. We represent the Laplace distribution with scale b = ∆ s / as a scale mixture of normals, i.e. a zero-mean normal with an exponential prior on the variance:

C Sensitivity of Sufficient Statistics in Truncated Model
Recall thatt(x) = 1 [v,w]  |t j (x)|, max x,y∈ [v,w] |t j (x) − t j (y)| D Proof of uniformity of CDF transform used by Cook et al. [30] Claim: Let X be a random variable with CDF F . The random variable U = F (X) is uniformly distributed.

Proof:
Pr(U ≤ u) = Pr(F (X) ≤ u) = Pr F −1 (F (X)) ≤ F −1 (u) = Pr(X ≤ F −1 (u)) = F F −1 (u) = u E Convergence of Gibbs Sampler Figure 2 shows the progress of sampled model parameters over the course of 500 iterations for both binomial and exponential models. For both models the samples quickly converge to the vicinity of the true parameter.