{'fraction_non_alphanumeric': 0.04381390005743825, 'fraction_numerical': 0.01707064905226881, 'mean_word_length': 4.459174714661984, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Valid distribution-free inferential models for prediction

A fundamental problem in statistics and machine learning is that of using observed data to predict future observations. This is particularly challenging for model-based approaches because often the goal is to carry out this prediction with no or minimal model assumptions. For example, the inferential model (IM) approach is attractive because it has certain validity guarantees, but requires specification of a parametric model. Here we show that a new perspective on a recently developed generalized IM approach can be applied to construct an IM for prediction that satisfies the desirable validity guarantees without specification of a model. One important special case of this approach corresponds to the powerful conformal prediction framework and, consequently, the desirable properties of conformal prediction follow immediately from the general IM validity theory. Several numerical examples are presented to illustrate the theory and highlight the method’s performance and flexibility.


Introduction
An important function played by a statistical model is that it identifies the unknown parameters about which inference is desired; for example, if the model is normal, then answers to relevant questions generally require inference on the mean and variance parameters. In the context of prediction, however, the future response is a well-defined quantity with or without a statistical model, so it is often expected that prediction be carried out without making specific model assumptions, such as normality. This modelfree perspective on prediction poses a challenge for those model-based approaches, such as Bayesian, fiducial, etc. For the most part, it is not possible to employ a model-based approach without a model, so the only option is to make the model wider and more flexible, which requires a high-or infinite-dimensional parameter. For example, in the now welldeveloped literature on Bayesian nonparametrics (e.g., Ghosal and van der Vaart 2017;Ghosh and Ramamoorthi 2003;Hjort et al. 2010), the likelihood is expressed in terms of an infinite-dimensional parameter, such as a density function, and in the currentlybeing-developed nonparametric fiducial literature (e.g., Cui and Hannig 2019), the datagenerating equation is expressed in terms of an infinite-dimensional parameter, such as a distribution function; see, also, Martin (2019a). An advantage of classical/frequentist approaches is that they can circumvent this challenging modeling step and construct methods for inference or prediction by taking advantage of other structures in the data. For example, in scalar prediction problems, it is straightforward to construct valid prediction intervals based on order statistics (e.g., Wilks 1941). A key observation is that these classical methods are based on certain distributional results-a "quasi-model"-derived from other weak, high-level assumptions about the data and structure of the problem. The question to be considered here is if those model-based approaches can be generalized to the extent that they can be applied at this quasi-model level.
Here we focus on a particular model-based approach, namely, inferential models (IMs) for inference and prediction; see, e.g., Liu (2013, 2015b) and Martin and Lingham (2016). This approach has close connections to fiducial inference (Fisher 1935;Taraldsen and Lindqvist 2013), generalized fiducial inference (Hannig et al. 2016), structural inference (Fraser 1968), Dempster-Shafer theory (Dempster 2008(Dempster , 1967(Dempster , 1968(Dempster , 2014Shafer 1976), and others, as will become clear as we proceed, but is related to Bayesian inference as well (Martin and Liu 2015a, Remark 4). The starting point of the IM construction is a so-called association, generically expressed as Y = a(θ, U ), that appropriately links the data Y , unknown parameter θ, and unobservable auxiliary variable U ; this is similar to the fiducial data-generating equation in (Hannig et al. 2016). In virtually every IM (or fiducial) application, this association is effectively an algorithm for simulating data from a posited model, which makes the approach "model-based." However, Martin (2015Martin ( , 2018 proposed a generalized association, motivated by the observation that the association need not fully characterize the data-generating process. Breaking the link between the association and the data-generating process creates an opportunity to construct an inferential model that, despite its name, is not model-based. In this paper, we carry out this construction in the context of prediction. We begin, in Section 2, with a brief introduction to IMs in general and a more detailed survey of the first IM-based thoughts on prediction in the context of a parametric statistical model. Next, in Section 3, we present our key insight, namely, that an IM can be constructed without requiring a statistical model or, alternatively, without requiring an explicit formula that describes the data generating process. That is, the aforementioned association, a critical ingredient to the IM construction, can be based on a suitable transformation of the observable and to-be-predicted data with some known distributional features. The particular implementation of this idea that we focus on here has some close connections to the powerful conformal prediction method of Vovk et al. (2005), Shafer and Vovk (2007), and Balasubramanian et al. (2014), which are now widely used in statistics and machine learning. Conformal prediction is mysterious in the sense that it does not obviously belong under the classical frequentist or Bayesian umbrellas, but our analysis here will shed light on what conformal prediction actually is, namely, a part of the IM framework. More specifically, we show • in Section 3.2 that the output of the conformal algorithm, often referred to as a "p-value" (e.g., Section 4.2 of Shafer and Vovk (2007)) is, in fact, in the language pioneered by Shafer (1976) and later adopted by Martin and Liu, just a plausibility function resulting from the general IM construction; • and in Section 3.3 that the prediction coverage probability guarantees, a distin-guishing feature of the conformal algorithm, is a direct consequence of the general IM validity theory.
After some technical remarks in Section 4, concerning optimality and an important special case of our developments, we present several examples in Section 5 to highlight the method's versatility: prediction of one-and two-dimensional observations, where the latter makes use of the notion of data depth (Liu et al. 1999;Tukey 1975), and then prediction of a response variable with covariates. Finally, we conclude in Section 6 with a brief discussion of some open problems.

Basic inferential models
Suppose the observable data Y ∈ Y is modeled by a distribution P Y |θ indexed by an unknown parameter θ ∈ Θ. Following Liu (2013, 2015b), an inferential model is a map from observed data Y = y, statistical model, etc. to a data-dependent function b y (·), taking values in [0, 1], where, for any hypothesis A ⊆ Θ, b y (A) represents the data analyst's degree of belief about A based on y. This definition encompasses many different approaches to statistical inference, e.g., Bayes, fiducial, etc. To avoid making "systematically misleading conclusions" (Reid and Cox 2015), it is essential that b y (A) is valid in the sense that b Y (A), as a function of Y ∼ P Y |θ , does not tend to be too large when the hypothesis A is false. More precisely, b y (A) should satisfy According to the false confidence theorem in Balch et al. (2019), the validity property in (1) cannot be achieved when b y is additive, like a probability distribution. Therefore, to achieve validity, non-additive functions b y must be considered. To our knowledge, the only general construction that achieves (1) is the one first presented in Martin and Liu (2013) and developed further in Martin and Liu (2015b). Since the focus in this paper is prediction, which is a little different than inference, we present the relevant IM details only for the prediction context; see Section 2.2 below.

IMs for prediction: first thoughts
Suppose the goal is to predict a future observation Y n+1 based on observed data Y n = (Y 1 , . . . , Y n ). Assume (Y n , Y n+1 ) are iid from P Y |θ , so that the old and new data are linked through a common parameter θ. Martin and Lingham (2016) showed that prediction can effectively be treated as a marginal inference problem, one where θ itself is a nuisance parameter. This suggests employing the two dimension-reduction strategies outlined in Martin and Liu (2015a,c), and the three-step construction of an inferential model for prediction is summarized below; for details, see Martin and Lingham (2016).
A-step. Write the joint association for Y n and Y n+1 as where a T and a are known functions, and P V,U n+1 is a known distribution. The left hand side equation above represents a dimension-reduced association between data, parameter and auxiliary variable V . For example, T (Y n ) could be a minimal sufficient statistic for θ. The key assumption is that the equation T (Y n ) = a T (V, θ) can be solved for θ; denote the solution by θ = θ(T (Y n ), V ). Plugging this solution into the second equation gives as the marginal association for Y n+1 . For fixed Y n , the right-hand side is a function of random variables (V, U n+1 ); write G Y n for its distribution function. Then (3) can be rewritten as P-step. Specify a suitable random set S ∼ P S that targets the unobservable auxiliary variable W in (4). One example is a symmetric random set Further discussion of this choice is given in Martin and Lingham (2016) and, in the context of our new developments here, in Section 3.2.2 below.
C-step. The C-step proceeds by combining the A-and the P-steps as follows. For the observed Y n = y n , write G −1 y n (S) for the data-dependent random set obtained as the image of S under G −1 y n . Then the plausibility function to be used for prediction of Y n+1 is a summary of the distribution of G −1 y n (S). That is, if A is some assertion about Y n+1 , then the plausibility of A, based on data y n , is This can be evaluated for any A, so we obtain a sort of distribution-like summary of our uncertainty about Y n+1 based on observations y n . For certain classes of assertions, it may be possible to summarize the plausibility function graphically. One such class of assertions are singletons, A = {ỹ}, for varyingỹ. We will refer to the plausibility function evaluated at singleton assertions as the plausibility contour and write pl y n (ỹ) ≡ pl y n ({ỹ}) = P S {G −1 y n (S) ỹ}.
Plots of the plausibility contour are shown in Figures 1-2 below. Also, when the random sets are nested, the plausibility function A → pl y n (A) is completely determined by its plausibility contour.
As discussed in Section 1, validity of inference and prediction is essential to the logic of statistical inference and, therefore, it is necessary that this validity property can be established for our proposed methods. In the present context, the desired result is That is, the plausibility contour value, pl Y n (Y n+1 ), as a function of the data Y n+1 is stochastically no smaller than Unif(0, 1). Intuitively, this is desirable because we would be prone to errors in prediction if the actual Y n+1 was determined to be relatively implausible based on data Y n . More specifically, one particularly relevant summary of our datadependent plausibility function/contour is a (prediction) plausibility region: P α (y n ) = {ỹ : pl y n (ỹ) > α}, for any fixed α ∈ (0, 1).
Then (6) above ensures that the frequentist prediction coverage probability of this plausibility region is at least the nominal level, uniformly over θ. Martin and Lingham (2016) demonstrate that (6) holds, at least approximately, for a wide range of problems.
3 Nonparametric IMs for prediction

Setup
The prediction approach in Martin and Lingham (2016), described above, is powerful but makes very specific model assumptions. These assumptions are undesirable because, in prediction settings, the goal is often to make as few model assumptions as possible. In particular, we might only be willing to assume that there is an exchangeable stochastic process Y 1 , Y 2 , . . ., where the common marginal distribution of the (possibly multivariate) Y i 's is continuous. Under these minimal model assumptions, the goal is to construct a valid IM for predicting Y n+1 based on observable data Y n = (Y 1 , . . . , Y n ).
The first obstacle is that the existing work on IMs has focused primarily on cases with a fully specified parametric model. To overcome this obstacle, we make use of the key insight in Martin (2015Martin ( , 2018, namely, that the association, as formulated in the A-step of the IM construction, need not fully characterize the data-generating process. Instead, the association can be built between suitable functions of the observable data, unknown quantities, and unobservable auxiliary variables.

A-step
Exchangeability of Y 1 , Y 2 , . . . is our only model assumption, so this will be central to our IM construction. But according to Martin (2015Martin ( , 2018, the association between the observable data, quantity of interest, and unobservable auxiliary variables can be generalized-it need not be formulated directly in terms of the observable data and an assumed data-generating process. That is, certain non-invertible transformations can be considered and, as long as those transformations inherit the exchangeability in the observable data, then we can construct a valid IM for predicting Y n+1 based on Y n .
Specifically, consider a transformation Y n+1 → T n+1 defined by the rule A is a specified real-valued function, symmetric in the elements of its first vector argument, and I k = {1, 2, . . . , k} for integer k ≥ 1. A simple example is which compares y i to the average of the entries in y n+1 −i . Shafer and Vovk (2007) refer to the function A as a non-conformity measure, and the interpretation is that A(y n+1 −i ; y i ) is small if and only if y i agrees with a prediction derived based on data y n+1 −i . Other examples of A are possible, e.g., Vovk et al. (2005), Lei et al. (2018), Hong and Martin (2019), some fairly sophisticated, especially in regression-like situations where the Y 's are response variables paired with corresponding covariates; see Section 5.3. The key point is that the symmetry in A guarantees that exchangeability is preserved when Y n+1 gets mapped to T n+1 . Where necessary in what follows, we will highlight T i 's dependence on the data Y n+1 by writing it as T i (Y n , Y n+1 ). The joint distribution of T n+1 is complicated, but we can write where the U i 's are exchangeable, marginally Unif(0, 1), with dependence structure that is too complicated to describe precisely, and the common marginal distribution function, G, is an infinite-dimensional nuisance parameter. Our goals are, first, to avoid directly dealing with G and, second, since we aim to predict the single data point Y n+1 , to reduce the dimension so that specifying a random set for the complicated and relatively highdimensional U n+1 is unnecessary.
The key observation is that, since G is monotone, the ranks of the T i 's are the same as those of the U i 's and, moreover, according to the following lemma, the latter ranks are marginally discrete uniform on I n+1 , denoted by Unif(I n+1 ).
Proof. The case where U 1 , . . . , U n , U n+1 are iid is straightforward: using iterated expectation, conditioning on the values of U (v−1) and U (v) , independence, and the standard result that uniform order statistics are equally spaced, on average, we get Then the more general exchangeable case follows from this, by appealing to de Finetti's theorem (e.g., Hewitt and Savage 1955), which says that exchangeability is equivalent to a certain form of conditionally iid. That is, by writing where F is the random distribution under which the U i 's are (conditionally) iid, and then repeating the same argument for the iid case above inside the square brackets, we get that the right-hand side above is equal to (n + 1) −1 , which proves the claim.
The quantity T n+1 is the one that holds the to-be-predicted value, Y n+1 , in special status, and this suggests a dimension-reduced (generalized) association, where r(·) is the ranking operator, that depends implicitly on T n+1 and, hence, on Y n+1 .
Here we assign rank 1 to the smallest value, rank 2 to the second smallest, and so on, because small values of the non-conformity measure A are "better" in a certain sense.
That the values are ranked in ascending versus descending order will be important in the P-step below. For now, we have completed the A-step: (11) defines generalized association-in the sense of Martin (2015Martin ( , 2018-that links the data Y n and the to-bepredicted value Y n+1 to an auxiliary variable V with known distribution.

P-step
Towards valid prediction of Y n+1 , here our intermediate goal is to specify a random set targeting the unobserved realization of the auxiliary variable V introduced above. The existing literature on random sets for this purpose has focused exclusively on cases where the auxiliary variable being targeted has a continuous distribution. In the present context, however, the auxiliary variable is discrete, so some special considerations are needed. Below we give only a brief summary of what specifically is needed for our present application, and refer the reader to Cella and Martin (2020) for the general details. The auxiliary variable V in our case takes values in I n+1 so we consider a random set S supported on a subset S of the power set, 2 I n+1 . Write P S for the distribution of S which, in this case, is rather simple: it is just a mass function defined on S. A summary of P S that is particularly relevant for our purposes here is its contour function, This is called the covering function of the random set S in Nguyen (2006)p. 40. The subscript "S" in γ S indicates which random set/distribution it pertains to; it is shorthand notation for the more appropriate but cumbersome γ P S . Following the developments in Martin and Liu (2015b), we are interested in the distribution of γ S (V ) as a function of a random variable V which, as in Lemma 1, has a Unif(I n+1 ) distribution. In particular, we seek a random set S such that the distribution of γ S (V ), as a function of V ∼ Unif(I n+1 ), is as close to (discrete) uniform as possible. For general problems, especially in the context of vector-valued auxiliary variables with a continuous distribution, this choice of random set S is a difficult one (e.g., Cella and Martin 2020) . Fortunately, thanks to our dimension-reduction steps and the structure imposed by our choice of A, the choice of S is straightforward here. Define the random set which is the push-forward of the Unif(I n+1 ) distribution through the set-valued mapping v → I v . This random set makes intuitive sense because, since the ranking operator is relative to ascending order, so that "rank equals 1" corresponds to a prediction consistent with the observed data, the random set S should almost surely include value 1. Furthermore, it can be shown that S in (12) is valid in the following sense.
Definition 1. A random set S is valid for predicting an auxiliary variable V ∼ Unif(I n+1 ) if γ S (V ), as a function of V , is stochastically no smaller than Unif{(n + 1) −1 I n+1 }, i.e., where a denotes the greatest integer less than or equal to a.
Validity of the random set S is an important property because the desired prediction calibration is an immediate consequence; see Theorem 1. Moreover, since the user is free to choose S, this validity property is not an assumption, it is completely under the control of the data analyst.
Lemma 2. The random set S defined in (12) is valid in the sense of Definition 1.
Proof. By direct calculation, we have , the claim follows immediately.
The proof of Lemma 2 shows that γ S (V ) is exactly Unif{(n + 1) −1 I n+1 } distributed, not just stochastically no smaller. This is an important characteristic because it implies that the set S in (12) is efficient in a certain sense. To see this, note first that the constant random set S ≡ I n+1 is valid but, because it is the full V -space, it cannot provide any valuable information. So, at least intuitively, we seek the "smallest" random set that is valid. If we measure the "size" of a random set S via its covering function, γ S , and if S is another random set that is smaller than S in the sense that then it is easy to check that S cannot be valid in the sense of Definition 1. Therefore, while validity guarantees certain prediction calibration properties (see Theorem 1), a random set that has exactly a uniform distribution, like that in (12), is the smallest such random set and so we call it efficient.
We conclude the P-step by introducing the random set S in (12) targeting the unobserved auxiliary variable V in the A-step from Section 3.2.1.

C-step
Next we combine the results of the A-and P-steps to construct a plausibility function on, in this case, the support of Y n+1 . If we express the result of the A-step with the v-indexed collection of sets Y y n (v) = {y n+1 : r(T n+1 (y n , y n+1 )) = v}, v ∈ I n+1 , and then combine this with S above, then the C-step yields the new random set on the Y n+1 -space, given by Y y n (S) = {y n+1 : r(T n+1 (y n , y n+1 )) ≤Ṽ },Ṽ ∼ Unif(I n+1 ).
Following the general construction of plausibility functions from these data-dependent random sets, it is easy to see that the corresponding prediction plausibility contour for Y n+1 is given by We immediately recognize the right-hand side above as the "p-value" formula in the conformal prediction approach; see, e.g., Section 4.2 of Shafer and Vovk (2007). Consequently, we have theoretical justification to treat the output from the conformal prediction algorithm as a genuine plausibility function arising from a valid prediction IM.

Prediction validity property
As a consequence of the above construction and the general theory of inferential models, we get the following prediction validity property.

Proof.
For a realization of the random set S ∈ 2 I n+1 , write S and S for the lower and upper bounds, respectively. By definition of the plausibility contour, we have If t n (α) ∈ I n+1 , then t n (α) = α and P V {γ S (V ) ≤ α} = α by Lemma 2; otherwise, t n (α) < α and P V {γ S (V ) ≤ t n (α)} < α, so the claim with "≤" holds generally.
This shows that the well-known validity of conformal prediction (e.g., Vovk et al. 2005) is a consequence of the general theory presented in, e.g., Martin and Liu (2015b). Indeed, it follows immediately from Theorem 1, that the 100(1−α)% prediction plausibility region P α (y n ) in (7) where the infimum is over all distributions P for the data process. This implies that prediction regions derived from the IM construction above will have frequentist coverage guarantees, for any sample size n and with no assumptions on the data-generating process beyond exchangeability.

Computation, visualization, and interpretation
Computation of the plausibility contour pl y n (y) in (13) is not trivial-it requires reprocessing the data for various candidate values of y in a grid. However, specification of the non-conformity measure is essentially the only task required of the data analyst. Once done, the computation of pl y n (y) proceeds as described in Algorithm 1.
A plot of y → pl y n (y), like in Figure 1(b), is a nice visual way of summarizing the information available in the data y n regarding the future observation Y n+1 , similar to the Bayesian posterior predictive density function. However, unlike in the Bayesian framework where integration is required to derive prediction regions, obtaining 100(1 − α)% prediction plausibility regions is straightforward: one just has to compare pl y n (y) to a threshold α to obtain P α (y n ). Moreover, the prediction region P α (y n ) has a simple and natural interpretation as the set of all y values such that the assertion "Y n+1 = y" is sufficiently plausible, given the observed data, y n . No other approaches offer prediction regions that assign such an interpretation to the individual elements they contain.

Optimality
Martin (2017) established a connection between valid IMs and confidence regions. He shows that, given a confidence region for some feature ϕ = ϕ(θ) of the full parameter, there exists a valid inferential model for θ whose corresponding marginal plausibility region for ϕ is at least as efficient as the given confidence region. Given the flexibility of the nonparametric IMs for prediction developed in Section 3, one may wonder if a similar result is possible for prediction intervals based solely on an exchangeability assumption. Vovk et al. (2005) and Shafer and Vovk (2007) present an optimality result in the context of conformal prediction. Specifically, for any valid and nested prediction interval, there exists a non-conformity measure A such that the corresponding conformal prediction interval derived from it is at least as efficient. Now, recall from Section 3 that the conformal prediction output is just a plausibility function arising from a valid prediction IM. Therefore, a result analogous to that in Martin (2017) for prediction is a consequence of the optimality argument in, e.g., Shafer and Vovk (2007, Sec. 4.4).

An important special case
An important special case arises when the Y i 's are scalar and, in the A-step of the nonparametric IM construction of Section 3.2.1, one considers T i = Y i in (8), i.e., In other words, no transformation Y n+1 → T n+1 is considered. The consequence is that the association in (10) becomes where the U i 's are exchangeable, marginally Unif(0, 1) and F is the continuous distribution function that characterizes the common marginal distributions of the Y i 's. As the dimension-reduced association in (11) is now in terms of r(Y n+1 ), the rank of Y n+1 , it can be rewritten as where Y (0) and Y (n+1) are defined as the infimum and supremum of the support of the marginal distribution of Y 1 , respectively.
Recall that the choice of the random set (12) in the P-step of the nonparametric IM construction (Section 3.2.2) was justified by the structure of the function A. There, we considered mappings A such that smaller values of them are "better" in the sense that y i agrees with the prediction derived based on data y n+1 −i . However, this is not necessarily the case for A in (15), as we are merely returning y i when comparing it to y n+1 −i ; see, also, Section 5.2. How to motivate a particular choice of random set S? Martin and Lingham (2016) consider different classes of assertions/hypotheses about Y n+1 , and the structure of these assertions can inform the choice the random set. For one-sided assertions, e.g., {Y n+1 ≤ỹ} or {Y n+1 ≥ỹ}, for varyingỹ, the natural choice of random set is similarly one-sided; see Cella and Martin (2020). Here, however, we will focus on singleton assertions, {Y n+1 =ỹ}, and the corresponding IM will produce two-sided prediction intervals. Specifically, for the P-step, we recommend the random set which consists of a random number of points closest to the midpoint n+2 2 . This is the discrete version of what Martin and Liu (2015b), call the "default" random set. Moreover, it is also easy to check that the random set S above is valid in the sense of Definition 1.
For the C-step, we combine S with the association (16), now represented as a collection to obtain the following new random set on the Y n+1 space:

The resulting plausibility function is
Using the same argument in the proof of Theorem 1, it can be verified that this plausibility function is valid in the sense of (14) and, consequently, prediction intervals derived from it have the nominal frequentist coverage probability. The prediction interval (7) based on the plausibility function in (19) will be an equal-tailed prediction bound for Y n+1 , always including the median of the observed data y n . See Section 5.1, Figure 1(b).

Connections to existing literature
Here we make a few remarks to connect the above developments to other related work in the literature. First, the reader is sure to recognize that the prediction intervals derived in Section 4.2 are not new, they agree with the classical intervals based on order statistics where, for pre-specified integers r and s such that 1 ≤ r < s ≤ n, the interval is a 100 s−r n+1 % prediction interval (Wilks 1941). It is satisfying to see that the IM framework, which has so far been focused on situations with a parametric statistical model, can be extended to cases without such a model and provide what would be considered a classical solution. Our analysis also sheds new light on that classical solution, in particular, it reveals that the latter has a "probabilistic" interpretation, beyond its more familiar interpretation as a frequentist procedure for constructing prediction intervals with error rate control. Indeed, our analysis shows that the classical interval can be interpreted as a set of candidate values for Y n+1 that are all sufficiently plausible based on the observed data Y n . This is important because the traditional notion of "confidence" is difficult for non-experts interpret and, arguably, not especially relevant to conclusions drawn in a particular application, so practitioners will interpret the prediction intervals in terms of plausibility anyway. Therefore, our analysis here gives a mathematical justification for the more intuitive interpretation of prediction intervals that practitioners actually use.
Second, readers familiar with the imprecise probability literature are sure to recognize some connections to the nonparametric predictive inference presented in, e.g., Augustin and Coolen (2004) and elsewhere. Based on what is referred to in this literature as Hill's assumption (Hill 1968(Hill , 1993, a pair of lower and upper probabilities for Y n+1 , given Y n = y n , are defined as, respectively, where B is some generic set in the support of Y n+1 and Y y n (v), for v ∈ I n+1 , is as in (18). Coolen (2006) interprets these in a "best of both worlds" sense: on one hand, these are genuine post-data lower and upper probabilities and, therefore, inherit certain coherence properties (e.g., de Finetti 1990;Walley 1991) which are relevant to one-off/personal decision making; on the other hand, through their connection to the underlying datagenerating process through Hill's assumption, they inherit certain calibration properties (e.g., Lawless and Fredette 2005) which are relevant to frequentist statistical inference. It turns out that these lower and upper probabilities can also be obtained from the above IM construction. Indeed, if we take S = {Ṽ }, forṼ ∼ Unif(I n+1 ), then the induced distribution of the data-dependent random set Y y n (S) generates the lower and upper probabilities defined above. Although the aforementioned lower and upper probability approach is a special case of that presented above, there are some key differences between the two kinds of output, resulting from our proposed random set S being nested and, therefore, structurally different from the non-nested singleton random sets that generate the (P, P) output. In particular, there is no useful sense in which once can assess the "plausibility" that the next observation, Y n+1 , will be equal to some specified y-the lower probability equals 0 and the upper probability equals (n + 1) −1 for every y. Compare this to the plausibility contour displayed in Figure 1(b), where it is clear that there is one y which is the "most plausible" value of Y n+1 based on the observed data y n . Having this plausibility contour makes reading off prediction intervals with desired coverage level straightforward. Of course, this same calibration property is embedded in the lower and upper probabilities above, but it is much more difficult, and unnatural, to extract when represented in that form. To our knowledge, given a lower probability like that above, the most direct way to extract a nominal 100(1 − α)% prediction interval is through the following calculation: In this case, with P as defined above, this calculation returns the classical prediction interval solution given in (20). If one is interested in determining which y are plausible values for Y n+1 based on the observed y n , the output produced by our proposed IM, with a nested random set, gives a more direct answer to the relevant questions.

One-dimensional prediction
Consider a random sample of size 50 of some univariate quantitative variable Y , whose histogram is shown in Figure 1(a). The goal is to construct a nonparametric IM for prediction of Y 51 such that the output plausibility function is centered on the median of the observed data, and equal-tailed prediction intervals can be obtained. Figure 1(b) shows, in red, the plausibility contour in (13) with and, in blue, the plausibility contour in (19). Recall that the blue plausibility contour is the one that leads to prediction plausibility intervals which coincide with the classic equal-tailed prediction intervals based on order statistics; see Section 4.2. For example, by thresholding it at 0.05 we obtain the 95% classic prediction interval [y (1) , y (50) ]. Note that the plausibility contour based on (13) is narrower than that based on (19), which indicates that the former is more efficient than the latter.

Two-dimensional prediction
An important feature of the nonparametric IM developed in Section 3 is its flexibility to deal with multivariate responses. Algorithm 1 is all that is needed to compute pl y n (y), regardless of the dimension. This is not the case for the classical prediction interval methods, discussed in Section 4.3, that depend on order statistics, since there is no natural ordering of multivariate data (Oja 2013).
Recall that the computation of pl y n (y) requires the specification of a non-conformity measure A. For multivariate responses, there are a number of options, for example, Lei et al. (2013) construct conformal prediction regions using a multivariate kernel density  (19) and (13), colored blue and red, respectively.
estimator. As an alternative, here we use data depth (e.g., Liu et al. 1999) to construct a non-conformity measure and, in turn, a prediction plausibility function as in (13). Roughly speaking, the concept of data depth amounts to assigning an appropriate ordering to the multivariate data, one with respect to a distance that measures how far away a point is from the center of a data cloud. This ordering is determined by a depth function and, among the variety of depth functions that appear in the literature, one of the most widely used is the so-called half-space depth proposed by Tukey (1975). The half-space depth, H(y | y n ), of a point y relative to a data set y n is determined as the smallest fraction of data points contained in a closed half-space with boundary through y. It ranges from 0-for the points that lie beyond the convex hull of the data-to its maximum value 1 2 attained at the Tukey median (Dyckerhoff and Mozharovskyi 2016). Computation of the half-space depth is discussed in Cuesta-Albertos and Nieto-Reyes (2008) and can be carried out using the depth.halfspace function in depth R package (Genest et al. 2019). Therefore, for multivariate prediction problems, this suggests the following non-conformity measure where H is the half-space depth function defined above.
As an illustration, consider the bivariate data y 60 = (y 1 , . . . , y 60 ), where each y i consists of the weight and engine displacement of one of 60 cars; these data are available in the R library rpart and are plotted in Figure 2(a). The goal is to predict Y 61 , the weight and engine displacement of the 61st car. Figures 2(b), 2(c) and 2(d) show, in different angles, the plausibility contour in (13)  prediction plausibility region for Y 61 . For further illustration, we consider a small simulation experiment. Assume we observe n = 100 bivariate iid observations, y 100 , either coming from a normal distribution or a bivariate Student-t distribution with three degrees of freedom, both with zero mean vector and covariance matrix composed by unit variances and a correlation coefficient equal to 0.5. The goal is to compare the performance of the nonparametric IM, again with non-conformity measure as in (22), with a parametric Bayes solution that assumes normality and assigns the Jeffreys prior to the parameters; see Yang and Berger (1998). For each distribution we evaluated the coverage probability and expected area of 95% prediction regions for Y 101 of each method, based on 1,000 Monte Carlo samples. The results are shown in Table 1. As expected, the nonparametric IM is valid regardless of the data distribution, and is almost as efficient as the parametric Bayesian solution when the assumed model is correct. When the assumed model is incorrect, the Bayesian 95% based on the two listed methods, for data coming from an bivariate normal or a bivariate Student-t distribution with 3 degrees of freedom, both with zero mean vector and covariance matrix composed by unit variances and a correlation coefficient equal to 0.5.
prediction ellipsoids end up being too small, leading to under-coverage.

Prediction with covariates
So far we have focused exclusively on the case where the data consists of a sequence of response variables, Y 1 , . . . , Y n , Y n+1 , but it is possible to extend all of the aforementioned developments to the case where the responses are coupled with corresponding covariates. That is, the data consists of an exchangeable sequence D 1 , D 2 , . . . , D n , D n+1 , where each D i = (X i , Y i ) is a covariate and response variable pair, where X i ∈ R p and Y i is a scalar. The goal is to predict the response Y n+1 at a given X n+1 , based on the observed data D n . Following the framework presented in Section 3, we consider a transformation of the data D n and the new D n+1 = (X n+1 , Y n+1 ) pair via a suitable non-conformity measure that is symmetric in its first vector argument. That is, As before, this transformation preserves the exchangeability of the original data sequence. For example, like in the A-step in Section 3.2, a reasonable choice of the A function is is the estimated mean function based on a model fit to d n+1 −i . Once the function A is specified by the data analyst, the same three-step IM construction in Section 3 can easily be followed: A-step. If r(·) denotes the (ascending order) ranking operator as before, then the dimensionreduced association in (11) becomes r(T n+1 ) = V , where V ∼ Unif{I n+1 }, which can be expressed with the v-indexed collection of sets P-step. Use the random set S in (12) to predict the unobserved auxiliary variable V .
C-step. Combine (23) and S above to get a new random set on the Y n+1 -space, The corresponding prediction plausibility contour for Y n+1 is given by Algorithm 2: Plausibility contour with covariates Initialize: data d n = {d 1 , . . . , d n } and x n+1 , non-conformity measure A, and a grid of y values; for each y value on the grid do Set d n+1 = (x n+1 , y) and write d n+1 = d n ∪ {d n+1 } ; Return pl d n (y | x n+1 ) for each y on the grid.
The computation of the prediction plausibility contour in (24) proceeds as described in Algorithm 2. Its validity can be easily demonstrated by an argument analogous to Theorem 1, and the plausibility region is a valid 100(1 − α)% prediction plausibility region, i.e., it satisfies where the infimum is over all distributions P for the data process. Despite being valid without any assumptions on P orμ, the efficiency of these intervals is depends on the quality of the estimatorμ of the underlying mean function µ (Lei et al. 2018).
To illustrate the practicality and flexibility of this nonparametric IM for prediction involving covariates, consider the following example. Let X i iid ∼ Unif(0, 1), i = 1, . . . , n, with n = 200, and let Y i = µ(X i ) + 0.1ε i , where µ(x) = sin 3 (2πx 3 ), and ε 1 , . . . , ε n iid ∼ t(5). Figure 3 displays the data, the true regression function µ(x) and the fitted regression curveμ(x) based on a B-spline with 12 degrees of freedom. A 95% prediction band is also displayed, derived by (25) and x n+1 taking values along the observed x 200 .

Conclusion
In this paper we have focused on a new "model-free" approach to the construction of an IM for prediction of a future response. Our key insight was that an association between observable data, the quantity of interest, and a set of unobsrvable auxiliary variables with a known distribution can be established based on suitable transformations of the original data as opposed to based on an explicit model/description of the data generating process. Our specific implementation of this insight had close connections to the conformal prediction methodology recently developed in the statistics and machine learning literature. This connection between conformal prediction and IMs is itself important because it helps to shed light on what the latter actually is, namely, a procedure for constructing a valid, non-additive plausibility function (Martin 2019b) to quantify uncertainty about Y n+1 , given Y n , one with frequentist prediction coverage guarantees, uniformly over all exchangeable distributions over the data sequence.
Exchangeability was crucial to our IM construction, that is, without exchangeability, we cannot establish the distribution of the auxiliary variables. While exchangeability is a relatively weak assumption compared to iid from a parametric family, there are, of course, situations where exchangeability is inappropriate, e.g., in time series or spatial data applications. Work to develop conformal prediction methods in non-exchangeable settings is an active area of current research-e.g., Mao et al. (2020) for spatial data-and the insights from the connection to IMs may prove to be useful in this endeavor. Another direction worth pursuing is the use of this form of generalized association in the context of statistical inference as opposed to prediction. A first step would be the development of IM-based solutions to the classical nonparametric problems, e.g., inference on quantiles without distributional assumptions, moving on to more modern/complex problems. Closely related to inference is that of assessing the quality or appropriateness of a posited statistical model, and here conformal-type methods have already proved to be useful (e.g., Lei et al. 2018;Rinaldo et al. 2020). That is, since a model is judged to be "good" only if it predicts the observed data reasonably well, methods that provide highquality prediction can be leveraged to provide equally high-quality model assessment. The work in the present paper provides a starting point for the construction of IM-based procedures for valid model assessment. More generally, while IMs had previously focused on classical problems with a fully specified statistical model, which limited its flexibility and applicability, the work presented in this paper reveals an opportunity for IMs play a role in solving modern problems where valid inference/prediction is required with as few model assumptions as possible.