Computing functions of random variables via reproducing kernel Hilbert space representations

We describe a method to perform functional operations on probability distributions of random variables. The method uses reproducing kernel Hilbert space representations of probability distributions, and it is applicable to all operations which can be applied to points drawn from the respective distributions. We refer to our approach as kernel probabilistic programming. We illustrate it on synthetic data and show how it can be used for nonparametric structural equation models, with an application to causal inference.


Introduction
Data types, derived structures, and associated operations play a crucial role for programming languages and the computations we carry out using them.Choosing a data type, such as Integer, Float, or String, determines the possible values, as well as the operations that an object permits.Operations typically return their results in the form of data types.Composite or derived data types may be constructed from simpler ones, along with specialised operations applicable to them.
The goal of the present paper is to propose a way to represent distributions over data types, and to generalize operations built originally for the data types to operations applicable to those distributions.Our approach is nonparametric and thus not concerned with what distribution models make sense on which statistical data type (e.g., binary, ordinal, categorical).It is also general, in the sense that in principle, it applies to all data types and functional operations.The price to pay for this generality is that • our approach will, in most cases, provide approximate results only; however, we include a statistical analysis of the convergence properties of our approximations, and • for each data type involved (as either input or output), we require a positive definite kernel capturing a notion of similarity between two objects of that type.Some of our results require, moreover, that the kernels be characteristic in the sense that they lead to injective mappings into associated Hilbert spaces.
In a nutshell, our approach represents distributions over objects as elements of a Hilbert space generated by a kernel, and describes how those elements are updated by operations available for sample points.If the kernel is trivial in the sense that each distinct object is only similar to itself, the method reduces to a Monte Carlo approach where the operations are applied to sample points which are propagated to the next step.Computationally, we represent the Hilbert space elements as finite weighted sets of objects, and all operations reduce to finite expansions in terms of kernel functions between objects.
The remainder of the present article is organized as follows.After describing the necessary preliminaries, we provide an exposition of our approach (Section 3).Section 4 analyses an application to the problem of cause-effect inference using structural equation models.We conclude with a limited set of experimental results.

Kernel Maps 2.1 Positive definite kernels
The concept of representing probability distributions in a reproducing kernel Hilbert space (RKHS) has recently attracted attention in statistical inference and machine learning [2,36].One of the advantages of this approach is that it allows us to apply RKHS methods to probability distributions, often with strong theoretical guarantees [44,45].It has been applied successfully in many domains such as graphical models [38,39], two-sample testing [11], domain adaptation [16,13,23], and supervised learning on probability distributions [24,48].We begin by briefly reviewing these methods, starting with some prerequisites.
We assume that our input data {x 1 , . . ., x m } live in a nonempty set X and are generated i.i.d. by a random experiment with Borel probability distribution p.By k, we denote a positive definite kernel on X × X , i.e., a symmetric function satisfying the following nonnegativity condition: for any m ∈ N, and a 1 , . . ., If equality in (3) implies that a 1 = • • • = a m = 0, the kernel is called strictly positive definite.

Kernel maps for points
Kernel methods in machine learning, such as Support Vector Machines or Kernel PCA, are based on mapping the data into a reproducing kernel Hilbert space (RKHS) H [3,34,35,15,47], where the feature map (or kernel map) Φ satisfies for all x, x ∈ X .One can show that every k taking the form ( 6) is positive definite, and every positive definite k allows the construction of H and Φ satisfying (6).The canonical feature map, which is what by default we think of whenever we write Φ, is with an inner product satisfying the reproducing kernel property Mapping observations x ∈ X into a Hilbert space is rather convenient in some cases.
If the original domain X has no linear structure to begin with (e.g., if the x are strings or graphs), then the Hilbert space representation provides us with the possibility to construct geometric algorithms by using the inner product of H .Moreover, even if X is a linear space in its own right, it can be helpful to use a nonlinear feature map in order to construct algorithms that are linear in H while corresponding to nonlinear methods in X .

Kernel maps for sets and distributions
One can generalize the map Φ to accept as inputs not only single points, but also sets of points or distributions.It was pointed out that the kernel map of a set of points corresponds to a kernel density estimator in the input domain [33,32], provided the kernel is nonnegative and integrates to 1.However, the map (10) can be applied for all positive definite kernels, including ones that take negative values or that are not normalized.Moreover, the fact that µ[X] lives in an RKHS and the use of the associated inner product and norm will have a number of subsequent advantages.For these reasons, it would be misleading to think of (10) simply as a kernel density estimator.
The kernel map of a distribution p can be defined as the expectation of the feature map [2,36,11], all of X (i.e., µ injective) Table 1: What information about a sample X does the kernel map µ[X] (see (10)) contain?
moments of p up to order n ∈ N k(x, x ) characteristic/universal all of p (i.e., µ injective) Table 2: What information about p does the kernel map µ[p] (see (11)) contain?For the notions of characteristic/universal kernels, see [46,8,7]; an example thereof is the Gaussian kernel (26).
where we overload the symbol µ and assume, here and below, that p is a Borel probability measure, and A sufficient condition for this to hold is the assumption that there exists an or equivalently k(x, x) ≤ M 2 , on the support of p. Kernel maps for sets of points or distributions are sometimes referred to as kernel mean maps to distinguish them from the original kernel map.Note, however, that they include the kernel map of a point as a special case, so there is some justification in using the same term.If p = p X is the law of a random variable X, we sometimes write µ In all cases it is important to understand what information we retain, and what we lose, when representing an object by its kernel map.We summarize the known results [47,8,36,11,45] in Tables 1 and 2.
We conclude this section with a discussion of how to use kernel mean maps.To this end, first assume that Φ is injective, which is the case if k is strictly positive definite (see Table 1) or characteristic/universal (see Table 2).Particular cases include the moment generating function of a RV with distribution p, which equals (11) for k(x, x ) = e x,x using (8).
We can use the map to test for equality of data sets, or distributions, Two applications of this idea lead to tests for homogeneity and independence.In the latter case, we estimate µ[p x p y ] − µ[p xy ] [1,12]; in the former case, we estimate µ[p] − µ[p ] [11].
Estimators for these applications can be constructed in terms of the empirical mean estimator (the kernel mean of the empirical distribution) where X = {x 1 , . . ., x m } is an i.i.d.sample from p (cf. (10)).As an aside, note that using ideas from James-Stein estimation [17], we can construct shrinkage estimators that improve upon the standard empirical estimator (see e.g., [25,26]).One can show that µ[ pm ] converges at rate m −1/2 (cf.[36] and [37,Theorem 27]): where K i j := k(x i , x j ).
Independent of the requirement of injectivity, µ can be used to compute expectations of arbitrary functions f living in the RKHS, using the identity which follows from the fact that k represents point evaluation in the RKHS, A small RKHS, such as the one spanned by the linear kernel may not contain the functions we are interested in.If, on the other hand, our RKHS is sufficiently rich (e.g., associated with a universal kernel [46]), we can use (19) to approximate, for instance, the probability of any interval (a, b) on a bounded domain, by approximating the indicator function I (a,b) as a kernel expansion ∑ n i=1 a i k(x i , .), and substituting the latter into (19).See [19] for further discussion.Alternatively, if p has a density, we can estimate it using methods such as reproducing kernel moment matching and combinations with kernel density estimation [41,19].
This shows that the map is not a one-way road: we can map our objects of interest into the RKHS, perform linear algebra and geometry on them [34], and at the end answer questions of interest.In the next section, we shall take this a step further, and discuss how to implement rather general operations in the RKHS.
Before doing so, we mention two additional applications of kernel maps.We can map conditional distributions and perform Bayesian updates [8,51,9], and we can connect kernel maps to Fourier optics, leading to a physical realization as Fraunhofer diffraction [14].
3 Computing Functions of Independent Random Variables

Introduction and Earlier Work
A random variable (RV) is a measurable function mapping possible outcomes of an underlying random experiment to a set E (often, E ⊂ R d , but our approach will be more general).The probability measure of the random experiment induces the distribution of the random variable.We will below not deal with the underlying probability space explicitly, and instead directly start from random variables X,Y with distributions p X , p Y and values in X , Y .Suppose we have access to (data from) p X and p Y , and we want to compute the distribution of the random variable f (X,Y ), where f is a measurable function defined on X × Y .For instance, if our operation is addition f (X,Y ) = X + Y , and the distributions p X and p Y have densities, we can compute the density of the distribution of f (X,Y ) by convolving those densities.If the distributions of X and Y belong to some parametric class, such as a class of distributions with Gaussian density functions, and if the arithmetic expression is elementary, then closed form solutions for certain favorable combinations exist.At the other end of the spectrum, we can resort to numerical integration or sampling to approximate f (X,Y ).
Arithmetic operations on random variables are abundant in science and engineering.Measurements in real-world systems are subject to uncertainty, and thus subsequent arithmetic operations on these measurements are operations on random variables.An example due to [42] is signal amplification.Consider a set of n amplifiers connected together in a serial fashion.If the amplification of the i-th amplifier is denoted by X i , then the total amplification, denoted by Y

, a product of n random variables.
A well-established framework for arithmetic operation on independent random variables (iRVs) relies on integral transform methods [42].The above example of addition already suggests that Fourier transforms may help, and indeed, people have used transforms such as the ones due to Fourier and Mellin to derive the distribution function of either the sum, difference, product, or quotient of iRVs [5,43,30,42].[49] proposes an approximation using Laguerre polynomials, and a notion of envelopes bounding the cumulative distribution function.This framework also allows for the treatment of dependent random variables, but the bounds can become very loose after repeated operations.[21] approximate the probability distributions of the input random variables as mixture models (using uniform and Gaussian distributions), and apply the computations to all mixture components.
[18] considers a numerical approach to implement arithmetic operations on iRVs, representing the distributions using piecewise Chebyshev approximations.This lends itself well to the use of approximation methods that perform well as long as the functions are well-behaved.Finally, Monte Carlo approaches can be used as well, and are popular in scientific applications (see e.g., [6]).
The goal of the present paper is to develop a derived data type representing a distribution over another data type, and to generalize the available computational operations to this data type, at least approximately.This would allow us to conveniently handle error propagation as in the example discussed earlier.It would also help us perform inference involving conditional distributions of such variables given observed data.The latter is the main topic of a subject area that has recently begun to attract attention, probabilistic programming [10].A variety of probabilistic programming languages has been proposed [50,27,4].To emphasize the central role that kernel maps play in our approach, we refer to it as kernel probabilistic programming (KPP).

Computing Functions of Independent Random Variables using Kernel Maps
The key idea of KPP is to provide a consistent estimator of the kernel map of an expression involving operations on random variables.This is done by applying the expression to the sample points, and showing that the resulting kernel expansion has the desired property.Operations involving more than one RV will increase the size of the expansion, but we can resort to existing RKHS approximation methods to keep the complexity of the resulting expansion limited, which is advisable in particular if we want to use it as a basis for further operations.The benefits of KPP are three-fold.First, we do not make parametric assumptions on the distributions associated with the random variables.Second, our approach applies not only to real-valued random variables, but also to multivariate random variables, structured data, functional data, and other domains, as long as positive definite kernels can be defined on the data.Finally, it does not require explicit density estimation as an intermediate step, which is difficult in high dimensions.
We begin by describing the basic idea.Let f be a function of two independent RVs X,Y taking values in the sets X , Y , and suppose we are given i.i.d.m-samples x 1 , . . ., x m and y 1 , . . ., y m .We are interested in the distribution of f (X,Y ), and seek to estimate its representation µ Although x 1 , . . ., x m ∼ p X and y 1 , . . ., y m ∼ p Y are i.i.d.observations, this does not imply that the { f (x i , y j )|i, j = 1, . . ., m} form an i.i.d.m 2 -sample from f (X,Y ), sinceloosely speaking -each x i (and each y j ) leaves a footprint in m of the observations, leading to a (possibly weak) dependency.Therefore, Theorem 1 does not imply that ( 22) is consistent.We need to do some additional work: Theorem 2 Given two independent random variables X,Y with values in X , Y , mutually independent i.i.d.samples x 1 , . . ., x m and y 1 , . . ., y n , a measurable function f : X × Y → Z , and a positive definite kernel on Z × Z with RKHS map Φ, then is an unbiased and consistent estimator of µ[ f (X,Y )].
Moreover, we have convergence in probability As an aside, note that ( 23) is an RKHS valued two-sample U-statistic.
Proof For any i, j, we have The convergence (24) can be obtained as a corollary to Theorem 3, and the proof is omitted here.

Approximating Expansions
To keep computational cost limited, we need to use approximations when performing multi-step operations.If for instance, the outcome of the first step takes the form (23), then we already have m • n terms, and subsequent steps would further increase the number of terms, thus quickly becoming computationally prohibitive.
We can do so by using the methods described in Chapter 18 of [34].They fall in two categories.In reduced set selection methods, we provide a set of expansion points (e.g., all points f (x i , y j ) in ( 23)), and the approximation method sparsifies the vector of expansion coefficients.This can be for instance done by solving eigenvalue problems or linear programs.Reduced set construction methods, on the other hand, construct new expansion points.In the simplest case, they proceed by sequentially finding approximate pre-images of RKHS elements.They tend to be computationally more demanding and suffer from local minima; however, they can lead to sparser expansions.
Either way, we will end up with an approximation of (23), where usually p m • n.Here, the z k are either a subset of the f (x i , y j ), or other points from Z .
It is instructive to consider some special cases.For simplicity, assume that Z = R d .If we use a Gaussian kernel whose bandwidth σ is much smaller than the closest pair of sample points, then the points mapped into the RKHS will be almost orthogonal and there is no way to sparsify a kernel expansion such as (23) without incurring a large RKHS error.In this case, we can identify the estimator with the sample itself, and KPP reduces to a Monte Carlo method.If, on the other hand, we use a linear kernel k(z, z ) = z, z on Z = R d , then Φ is the identity map and the expansion (23) collapses to one real number, i.e., we would effectively represent f (X,Y ) by its mean for any further processing.By choosing kernels that lie 'in between' these two extremes, we retain a varying amount of information which we can thus tune to our wishes, see Table 1.

Computing Functions of RKHS Approximations
More generally, consider approximations of kernel means µ[X] and µ In our case, we think of ( 27) as RKHS-norm approximations of the outcome of previous operations performed on random variables.Such approximations typically have coefficients α ∈ R n and β ∈ R m that are not uniform, that may not sum to one, and that may take negative values [34], e.g., for conditional mean maps [40,9].We propose to approximate the kernel mean µ[ f (X,Y )] by the estimator where the feature map Φ z defined on Z , the range of f , may be different from both Φ x and Φ y .The expansion has m • n terms, which we can subsequently approximate more compactly in the form (25), ready for the next operation.Note that ( 28) contains ( 23) as a special case.One of the advantages of our approach is that ( 23) and ( 28) apply for general data types.In other words, X , Y , Z need not be vector spaces -they may be arbitrary nonempty sets, as long as positive definite kernels can be defined on them.

Convergence analysis in an idealized setting
We analyze the convergence of (28) under the assumption that the expansion points are actually samples x 1 , . . ., x m from X and y 1 , . . ., y n from Y , which is for instance the case if the expansions (27) are the result of reduced set selection methods (cf.Section 3.3).Moreover, we assume that the expansion coefficients α 1 , . . ., α m and β 1 , . . ., β n are constants, i.e., independent of the samples.
The following proposition gives a sufficient condition for the approximations in (27) to converge.Note that below, the coefficients α 1 , . . ., α m depend on the sample size m, but for simplicity we refrain from writing them as α 1,m , . . ., α m,m ; likewise, for β 1 , . . ., β n .We make this remark to ensure that readers are not puzzled by the below statement that ∑ m i=1 α 2 i → 0 as m → ∞.Proposition 1 Let x 1 , . . ., x m be an i.i.d.sample and (α i ) m i=1 be constants with where X and X are independent copies of x i .Then, the convergence Proof From the expansion the assertion is straightforward.
The next result shows that if our approximations (27) converge in the sense of Proposition 1, then the estimator (28) (with expansion coefficients summing to 1) is consistent.
Theorem 3 Let x 1 , . . ., x m and y 1 , . . ., y n be mutually independent i.i.d.samples, and the constants Proof By expanding and taking expectations, one can see that which implies the norm in the assertion of the theorem converges to zero at under the assumptions on α i and β j .Here ( X, Ỹ ) is an independent copy of (X,Y ).This concludes the proof.Note that in the simplest case, where α i = 1/m and β j = 1/n, we have ∑ i α 2 i = 1/m and ∑ j β 2 j = 1/n, which proves Theorem 2. It is also easy to see from the proof that we do not strictly need ∑ i α i = ∑ j β j = 1 -for the estimator to be consistent, it suffices if the sums converge to 1.For a sufficient condition for this convergence, see [19].
More general expansion sets To conclude our discussion of the estimator (28), we turn to the case where the expansions (27) are computed by reduced set construction, i.e., they are not necessarily expressed in terms of samples from X and Y .This is more difficult, and we do not provide a formal result, but just a qualitative discussion.
To this end, suppose the approximations (27) n ∑ j=1 β j = 1 and for all j, β j > 0, (30) and we approximate µ[ f (X,Y )] by the quantity (28).We assume that ( 27) are good approximations of the kernel means of two unknown random variables X and Y ; we also assume that f and the kernel mean map along with its inverse are continuous.We have no samples from X and Y , but we can turn (27) into sample estimates based on artificial samples X and Y, for which we can then appeal to our estimator from Theorem 2.
To this end, denote by X = (x 1 , . . ., x m ) and Y = (y 1 , . . ., y n ) the expansion points in (27).We construct a sample X = (x 1 , x 2 , . . . ) whose kernel mean is close to ∑ m i=1 α i Φ x (x i ) as follows: for each i, the point x i appears in X with multiplicity m•α i , i.e., the largest integer not exceeding m • α i .This leads to a sample of size at most m.Note, moreover, that the multiplicity of x i , divided by m, differs from α i by at most 1/m, so effectively we have quantized the α i coefficients to this accuracy.Since m is constant, this implies that for any ε > 0, we can choose m large enough to ensure that 1 We may thus work with 1 m ∑ m 1=1 Φ x (x i ), which for strictly positive definite kernels corresponds uniquely to the sample X = (x 1 , . . ., x m ).By the same argument, we obtain a sample Y = (y 1 , . . ., y n ) approximating the second expansion.Substituting both samples into the estimator from Theorem 2 leads to where αi = m • α i /m, and βi = n • β i /n.By choosing sufficiently large m, n, this becomes an arbitrarily good approximation (in the RKHS norm) of the proposed estimator (28).Note, however, that we cannot claim based on this argument that this estimator is consistent, not the least since Theorem 2 in the stated form requires i.i.d.samples.
Larger sets of random variables Without analysis, we include the estimator for the case of more than two variables: Let g be a measurable function of jointly independent RVs U j ( j = 1, . . ., p).Given i.i.d.observations u j 1 , . . ., u j m ∼ U j , we have in probability.Here, in order to keep notation simple, we have assumed that the samples sizes for each RV are identical.
As above, we note that (i) g need not be real-valued, it can take values in some set Z for which we have a (possibly characteristic) positive definite kernel; (ii) we can extend this to general kernel expansions like (28); and (iii) if we use Gaussian kernels with width tending to 0, we can think of the above as a sampling method.

Dependent RVs and Structural Equation Models
For dependent RVs, the proposed estimators are not applicable.One way to handle dependent RVs is to appeal to the fact that any joint distribution of random variables can be written as a structural equation model with independent noises.This leads to an interesting application of our method to the field of causal inference.
We consider a model X i = f i (PA i ,U i ), for i = 1, . . ., p, with jointly independent noise terms U 1 , . . .,U p .Such models arise for instance in causal inference [28].Each random variable X i is computed as a function f i of its noise term U i and its parents PA i in an underlying directed acyclic graph (DAG).Every graphical model w.r.t. a DAG can be expressed as such a structural equation model with suitable functions and noise terms (e.g., [29]).
If we recursively substitute the parent equations, we can express each X i as a function of only the independent noise terms U 1 , . . .,U p , Since we know how to compute functions of independent RVs, we can try to test such a model (assuming knowledge of all involved quantities) by estimating the distance between RKHS images, using the estimator described in (33) (we discuss the bivariate case in Theorem 4).
It may be unrealistic to assume we have access to all quantities.However, there is a special case where this is conveivable, which we will presently discuss.This is the case of additive noise models [29] Such models are of interest for cause-effect inference since it is known [29] that in the generic case, a model (36) can only be fit in one direction, i.e., if (36) holds true, then we cannot simultaneously have To measure how well (36) fits the data, we define an estimator Analogously, we define the estimator in the backward direction Here, we assume that we are given the conditional mean functions f : In practice, we would apply the following procedure: we are given a sample (x 1 , y 1 ), . . ., (x m , y m ).We estimate the function f as well as the residual noise terms u 1 , . . ., u m by regression, and likewise for the backward function g [29].Strictly speaking, we need to use separate subsamples to estimate function and noise terms, respectively, see [20].
Below, we show that ∆ emp converges to 0 for additive noise models (36).For the incorrect model (37), however, ∆ bw emp will in the generic case not converge to zero.We can thus use the comparison of both values for deciding causal direction.
Theorem 4 Suppose x 1 , . . ., x m and u 1 , . . ., u m are mutually independent i.i.d.samples, and y i = f (x i ) + u i .Assume further that the direction of the additive noise model is identifiable [29] and the kernel for x is characteristic.We then have ∆ bw emp → 0 (41) in probability as m → ∞.
Proof Equation (40 → 0 (all convergences in this proof are in probability).To prove (41), we assume that ∆ bw emp → 0 which implies The key idea is to introduce a random variable Ṽ that has the same distribution as V but is independent of Y and to consider the following decomposition of the sum appearing in (42): where the index for v is interpreted modulo m, for instance, v m+3 := v 3 .Since v i+k = x i+k − g(y i+k ) is independent of y i , it further follows from Theorem 2 that A m − Together with (42) this implies Since the kernel is characteristic, this implies with Y ⊥ ⊥ Ṽ , which contradicts the identifiability of the additive noise model.As an aside, note that Theorem 4 would not hold if in (39) we were to estimate 5 Experiments

Synthetic data
We consider basic arithmetic expressions that involve multiplication X • Y , division X/Y , and exponentiation X Y on two independent scalar RVs X and Y .Letting p X = N (3, 0.5) and p Y = N (4, 0.5), we draw i.i.d.samples X = {x 1 , . . ., x m } and Y = {y 1 , . . ., y m } from p X and p Y .
In the experiment, we are interested in the convergence (in RKHS norm) of our estimators to µ[ f (X,Y )].Since we do not have access to the latter, we use an independent sample to construct a proxy μ[ f (X,Y )] = (1/ 2 ) ∑ i, j=1 Φ z ( f (x i , y j )).We found that = 100 led to a sufficiently good approximation.
Next, we compare the three estimators, referred to as µ 1 , µ 2 and µ 3 below, for sample sizes m ranging from 10 to 50: 1.The sample-based estimator (23) 2. The estimator (28) based on approximations of the kernel means, taking the form μ and µ[Y ], respectively.We used the simplest possible reduced set selection method: we randomly subsampled subsets of size m ≈ 0.4 • m from X and Y, and optimized the coefficients {α 1 , . . ., α m } and {β 1 , . . ., β m } to best approximate the original kernel means (based on X and Y) in the RKHS norm [34,Section 18.3].
3. Analogously to the case of one variable (17), we may also look at the estimator μ3 [X, Y] := (1/m) ∑ m i=1 Φ z ( f (x i , y i )), which sums only over m mutually independent terms, i.e., a small fraction of all terms of (23).
For i = 1, 2, 3, we evaluate the estimates μi [ f (X,Y )] using the error measure We use (6) to evaluate L in terms of kernels.In all cases, we employ a Gaussian RBF kernel (26) whose bandwidth parameter is chosen using the median heuristic, setting σ to the median of the pairwise distances of distinct data points [12].Figure 1 depicts the error (43) as a function of sample size m.For all operations, the error decreases as sample size increases.Note that Φ z is different across the three operations, resulting in different scales of the average error in Figure 1.

Causal discovery via functions of kernel means
We also apply our KPP approach to bivariate causal inference problem (cf.Section 4).That is, given a pair of real-valued random variables X and Y with joint distribution p XY , we are interested in identifying whether X causes Y (denote as X → Y ) or Y causes X (denote as Y → X) based on observational data.We assume an additive noise model E = f (C) +U with C ⊥ ⊥ U where C, E,U denote cause, effect, and residual (or "unexplained") variable, respectively.Below we present a preliminary result on the CauseEffectPairs benchmark data set [22].The error reported is an average of 30 repetitions of the simulations.The expensive estimator μ1 (see ( 23)) performs best.The approximation μ2 (see (28)) works well as sample sizes increase.
For each causal pair (X, Y) = {(x 1 , y 1 ), . . ., (x m , y m )}, we estimate functions y ≈ f (x) and x ≈ g(y) as least-squares fits using degree 4 polynomials.We illustrate one example in Figure 2a.Next, we compute the residuals in both directions as u i = y i − f (x i ) and v j = x j − g(y j ). 1  Following Theorem 4, we can use the comparison between ∆ X→Y and ∆ Y →X to infer the causal direction.Specifically, we decide that X → Y if ∆ X→Y < ∆ Y →X , and that Y → X otherwise.In this experiment, we also use a Gaussian RBF kernel whose bandwidth parameter is chosen using the median heuristic.To speed up the computation of ∆ X→Y and ∆ Y →X , we adopted a finite approximation of the feature map using 100 random Fourier features (see [31] for details).We allow the method to abstain whenever the two values are closer than δ > 0. By increasing δ , we can compute the method's accuracy as a function of a decision rate (i.e., the fraction of decisions that our method is forced to make) ranging from 100% to 0%. Figure 2b depicts the accuracy versus the decision rate for the 81 pairs in the CauseEffectPairs benchmark collection.The method achieves an accuracy of 80%, which is significantly better than random guessing, when forced to infer the causal direction of all 81 pairs.

Conclusions
We have developed a kernel-based approach to compute functional operations on random variables taking values in arbitrary domains.We have proposed estimators for RKHS representations of those quantities, evaluated the approach on synthetic data, and showed how it can be used for cause-effect inference.While the results are encouraging, the material presented in this article only describes the main ideas, and much remains to be done.We believe there is significant potential for a unified perspective on probabilistic programming based on the described methods, and hope that some of the open problems will be addressed in future work.

Figure 1 :
Figure 1: Error of the proposed estimators for three arithmetic operations -multiplication X • Y , division X/Y , and exponentiation X Y -as a function of sample size m.The error reported is an average of 30 repetitions of the simulations.The expensive estimator μ1 (see (23)) performs best.The approximation μ2 (see(28)) works well as sample sizes increase.

Figure 2 :
Figure 2: (a) Scatter plot of the data of causal pair 78 in the CauseEffectPairs benchmarks, along with the forward and backward function fits, y = f (x) and x = g(y).(b) Accuracy of cause-effect decisions on all the 81 pairs in the CauseEffectPairs benchmarks..