Stein's method and duality of Markov processes

TL;DR

This paper introduces a novel method using duality and Markov jump processes to explicitly solve Stein equations for diffusion processes, simplifying the evaluation of solutions for distributional approximation.

Contribution

It presents a new approach leveraging the semi-group of simpler Markov jump processes to solve Stein equations for diffusions, avoiding complex semi-group calculations.

Findings

New explicit solutions for Stein equations in diffusion processes

Applicable to both univariate and multivariate distributions

Simplifies the process of distributional approximation

Abstract

One of the key ingredients to successfully apply Stein's method for distributional approximation are solutions to the Stein equations and their derivatives. Using Barbour's generator approach, one can solve for the solutions to the Stein equation in terms of the semi-group of a Markov process, which is typically a diffusion process if it is a continuous distribution. For an arbitrary diffusion it can a difficult task to evaluate the semi-group and its derivatives. In this paper, for polynomial test functions, instead of calculating the semi-group of a diffusion, via a duality argument, we instead utilise the semi-group of a much simpler Markov jump process. This approach yields a new method for explicitly solving for the solutions of Stein equations for diffusion processes. We present both the general idea of the approach and examples for both univariate and multivariate distributions.

Equations

$$\mathcal{A}f(x)=f^{\prime \prime }(x)-x{f}^{\prime }(x).$$

$$\mathcal{A}{f}_h(x)=h(x)-\text{mathbbm}Eh(Z),$$

$$f_h(x)=-{\int }_0^{\infty }[\text{mathbbm}Eh(Z_x(t))-\text{mathbbm}Eh(Z)]dt,$$

$$\text{mathbbm}{E}_{x}H(X(t),y)=\text{mathbbm}{E}^{y}H(x,Y(t)),$$

$$\text{mathbbm}{P}_x(X_0\le Y_t)=\text{mathbbm}{P}^y(X_t\le Y_0).$$

$$\mathcal{A}f(x)=f^{\prime \prime }(x)-x{f}^{\prime }(x),$$

$$\mathcal{B}f(y)=y(y-1)(f(y-2)-f(y)),$$

$$\text{mathbbm}E[X(t{)}^{Y(0)}]=\text{mathbbm}E\left[X(0{)}^{Y(t)}\exp \left({\int }_0^t(Y^2(u)-2Y(u))du\right)\right].$$

$$f_h(x)=-{\int }_0^{\infty }\left\{\text{mathbbm}E\left[X(0{)}^{Y(t)}\exp \left({\int }_0^t(Y^2(u)-2Y(u))du\right)\right]-\text{mathbbm}Eh(Z)\right\}dt.$$

$$f_h(x)$$

$$=-\text{mathbbm}E{\int }_0^{{\tau }_1}{x}^3\exp \left({\int }_0^t{Y}^2(u)-2Y(u)du\right)dt-\text{mathbbm}E{\int }_{{\tau }_1}^{\infty }x\exp \left({\int }_0^t{Y}^2(u)-2Y(u)du\right)dt$$

$$=-\text{mathbbm}E{\int }_0^{{\tau }_1}{x}^3\exp \left\{3t\right\}dt-\text{mathbbm}E{\int }_{{\tau }_1}^{\infty }x\exp \left(3{\tau }_1-(t-{\tau }_1)\right)dt$$

$$=-\frac{x^3}{3}\text{mathbbm}E(e^{3{\tau }_1}-1)-x\text{mathbbm}E{e}^{3{\tau }_1}$$

$$=-\frac{x^3}{3}-2x.$$

$$f_h(x)=-{\int }_0^{\infty }\left\{\text{mathbbm}Eh(e^{-s}x+\sqrt{1-e^{-2s}}Z)-\text{mathbbm}Eh(Z)\right\}ds.$$

$$f_h(x)$$

$$=-{\int }_0^{\infty }{x}^3{e}^{-3s}+3e^{-s}x(1-e^{-2s})ds$$

$$=-\frac{x^3}{3}-2x,$$

$$\mathcal{A}f(x)$$

$$=k(k-1)x^{k-2}-k{x}^k$$

$$=k((k-1)x^{k-2}-x^k).$$

$$\text{mathbbm}Eh(Z_x(t))=\text{mathbbm}E(Y_x(1,t)).$$

$$f_h(x)=-{\int }_0^{\infty }\text{mathbbm}Eh(Z_x(t))$$

$$=-\text{mathbbm}E{\int }_0^{{\tau }_1}{x}^{3}dt-\text{mathbbm}E{\int }_{{\tau }_1}^{{\tau }_2}2xdt$$

$$=-x^3/3-2x,$$

$$\mathcal{A}f(x)=\sum _{i=1}^k{a}_i[b_i{x}^{k-i}-x^k],$$

$$f_h(x)=-{\int }_0^{\infty }\text{mathbbm}E{Y}_x(1,t)-Eh(Z)dt,$$

$$\mathcal{B}x(g(a,k))=\sum _{i=1}^k{a}_i[g_x(a{b}_i,k-i)-g_x(a,k)].$$

$$\text{mathbbm}Eh(X_x(t))-h(x)$$

$$={\int }_0^t\sum _{i=1}^k\text{mathbbm}E{a}_i[b_i{X}_x(s{)}^{k-i}-X_x(s{)}^k]ds,$$

$$\frac{d}{dt}\text{mathbbm}Eh(X_x(t))=\text{mathbbm}E\sum _{i=1}^k{a}_i[b_i{X}_x(t{)}^{k-i}-X_x(t{)}^k].$$

$$\mathcal{A}f(x)$$

$$=k(k-1)x^{k-2}-k{x}^k$$

$$=k((k-1)x^{k-2}-x^k),$$

$$f_h(x)$$

$$=-\text{mathbbm}E\sum _{i=0}^k{\int }_{{\tau }_i}^{{\tau }_{i+1}}\frac{(2k)!!}{(2k-2i)!!}{x}^{2k+1-2i}dt$$

$$=-\sum _{i=0}^k\frac{(2k)!!}{(2k-2i)!!(2k-2i+1)}{x}^{2k+1-2i},$$

$$f_h(x)$$

Full text

1 Introduction

Stein’s method using generator approach, pioneered in Barbour, [1988, 1990], involves characterising a target distribution as the stationary distribution of an (infinitesimal) generator of a stochastic process. Such a generator is then used as the Stein operator or identity. As an example, consider the following generator,

[TABLE]

This is the generator of an Ornstein-Uhlenbeck process, and notably the stationary distribution of such a process is the standard normal distribution. The general idea of the generator method can be summarised as follows.

Let $\mathcal{A}f$ be a generator, and assume that its associated stationary distribution $\pi$ exists and is unique. Then for any function $h$ from a suitably rich family of test functions $\mathcal{F}$, let $f_{h}$ be the solution to

[TABLE]

where $Z\sim\pi$. For any random variable $W$, our aim is to bound $|\mathbbm{E}h(W)-\mathbbm{E}h(Z)|$ for all $h\in\mathcal{F}$. To achieve this, we substitute $x=W$ into Equation (1.1), and take expectations of both sides. The problem is now transferred to taking expectations of the generator $\mathcal{A}f_{h}(W)$. It turns out that properties of the function $f_{h}$ are essential, in particular supremum bounds on the function and its derivatives. Under the right conditions, (for example if $\mathcal{A}$ generates a Feller semi-group), it can be shown that for any $h\in\mathcal{F}$,

[TABLE]

where $Z_{x}(t)$ is a stochastic process following generator $\mathcal{A}$ such that $Z_{x}(0)=x$. The main difficulty with using such an approach is that to evaluate (1.2), one typically needs to have a very strong understanding of the process $Z_{x}(t)$. Our approach in this paper is to bypass this problem by evaluating $\mathbbm{E}h(Z_{x}(t))$ for a diffusion by instead calculating the semi-group for a much simpler Markov jump process via a duality argument.

We begin by defining the canonical definition of duality of two Markov processes with respect to a duality function.

Definition 1.1**.**

Let $X(t)$ and $Y(t)$ be two Markov processes with state spaces $E_{1}$ and $E_{2}$, and $H$ be a (bounded measurable) function on $(E_{1},E_{2})$. Then $X(t)$ is dual to $Y(t)$ with respect to a function $H$ if for all $x\in E_{1},y\in E_{2}$ and $t\geq 0$,

[TABLE]

where $\mathbbm{E}_{x}$ indicates $X(0)=x$ and $\mathbbm{E}^{y}$ indicates $Y(0)=y$.

For a summary of Markov dual processes, see the survey article Jansen and Kurt, [2014]. Note that the restriction that $H$ be bounded and measurable can be weakened, and the equality above takes a slightly different form, see [Ethier and Kurtz,, 1986, Section 4.4] for full details.

In this paper we will be focusing on using a form of moment duality for the convergence of moments of distributions. We note that it can be trivially extended to polynomial test functions. Convergence of $k$-th moments, is generally not a sufficient condition to ensure convergence in distribution. If in addition to convergence of all moments, the limiting distribution has finite support, or finite moment generating function in a neighborhood of 0, then this yields convergence in distribution Billingsley, [2013]. The examples we present will usually satisfy at least one of these two conditions.

Our goal is to use a dual process to calculate the semi-group for test functions that correspond to moments of random variables. The corresponding duality function for moments is therefore functions $H$ of the form $H(x,y)=x^{y}$. It remains an interesting open problem whether other forms of duality can also yield fruitful results. For example a common form of duality known as Siegmund duality, takes the form

[TABLE]

Such a duality would give an alternative approach to calculating bounds in Kolmogorov distance.

We begin with a simple motivating example for the Ornstein-Uhlenbeck process and Normal approximation, to illustrate how we can utilise duality to explicitly calculate solutions to the Stein equation for moment test functions. Recall

[TABLE]

as the generator of the Ornstein-Uhlenbeck process. Furthermore define the following

[TABLE]

as the generator of a Markov jump process on non-negative integers such that when it is in state $y$, it waits an exponential $y(y-1)$ amount of time before it jumps down in step sizes of 2, and this process continues until it is absorbed at $1$ or [math]. It can be shown that these two processes are dual to each other with respect to the moment duality function $H(x,y)=x^{y}$, see [Ethier and Kurtz,, 1986, Example 4.4.9] for full details. Set $X(t)$ and $Y(t)$ to be processes with generators $\mathcal{A}$ and $\mathcal{B}$ respectively. Then it can be shown that

[TABLE]

This slightly more complicated form of duality compared to (1.3) due to the fact that our dual function is unbounded. Recall the solution to the Stein equation (1.2) takes the form $f_{h}(x)=-\int_{0}^{\infty}[\mathbbm{E}h(Z_{x}(t))-\mathbbm{E}h(Z)]dt$, where $Z\sim\mathrm{N}(0,1)$. Set $h(x)=x^{k}$, $X(0)=x$, and $Y(0)=k$, then (1.4) gives us an alternative approach to evaluating $\mathbbm{E}h(Z_{x}(t))$,

[TABLE]

The primary advantage of such an approach is that our solution rather than being in terms of the Ornstein-Uhlenbeck diffusion, is now written purely in terms of a relatively simple Markov jump process.

We compute a simple example case where $k=3$ to illustrate our general approach and also to confirm that this yields the correct solution. Let $\tau_{1}$ be an exponential random variable with parameter 6. Then,

[TABLE]

For the Ornstein-Uhlenbeck process and Normal approximation, a closed form for the solution to the Stein equation is well known Barbour, [1990]. We can thus use this to confirm that our approach above is successful. For a general function $h$, the solution to the Stein equation is

[TABLE]

Hence for $h(x)=x^{3}$,

[TABLE]

which is the desired answer. As the example indicates, by finding an appropriate dual process we can certainly explicitly solve for solutions to the Stein equation for certain test functions. The primary question then becomes, how does one find an appropriate dual process? The existence of dual processes is often unclear, let alone a systematic approach for finding such a process, as discussed in the introduction of Jansen and Kurt, [2014]. To address this issue, in this paper we use a modified version of duality, that will be relatively easy to find.

If $f(x)=x^{k}$, then the generator of the Ornstein-Uhlenbeck process can be rewritten in the following form,

[TABLE]

This suggests an alternative duality with a Markov jump process on the space of functions. The process begins with the function $x^{k}$, and at a rate of $k$ it transitions to the function $(k-1)x^{k-2}$. Then at a rate of $k-2$ it transitions to the function $(k-1)(k-3)x^{k-4}$ and so forth until it reaches its absorbing state. Depending on whether $k$ is even or odd, this process will either be absorbed at $(k-1)!!$ or 0 respectively, which importantly are the moments of the standard Normal distribution. Note if $f(x)=\alpha x$, $\mathcal{A}f(x)=-\alpha x=1(0-\alpha x)$, so can be interpreted as transitioning from $\alpha x$ to 0 at a rate of 1. Denote by $Y_{x}(1,t)$ to be the process on functions as described, such that $Y_{x}(1,0)=x^{k}$. (The argument 1 corresponds to the coefficient of $x^{k}$ in our test function.) We now continue our simple example where $k=3$.

Suppose the following duality equation were true: Let $h(x)=x^{3}$, and $Z_{x}(t)$ be the Ornstein-Uhlenbeck process started at $x$, and

[TABLE]

Then,

[TABLE]

where we used the fact that $\tau_{1}\sim\mathop{\mathrm{Exp}}(3)$ and $\tau_{2}-\tau_{1}\sim\mathop{\mathrm{Exp}}(1)$.

Our approach in this paper will be to use dual processes of the above form. To successfully apply Stein’s method, one typically requires uniform bounds on the supremum of not only the solutions to the Stein equation, but also their derivatives. If the support of our target distribution is unbounded, our approach will explicitly solve for the Stein solution, but will usually not yield uniform supremum bounds. However, these bounds can still be used if more care is taken, see for example Gaunt, [2015]. In the case where the support of the distribution is finite, then one can use our approach to find uniform bounds on derivatives of the Stein equation. We shall see this in a few examples.

The remainder of the paper will be structured as follows. In Section 2, we will focus on univariate distributions, formally define our notion of duality, justify its construction and give the general approach for finding a suitable dual process. It will include various examples, and their results will mostly coincide exactly with known results derived from alternative methods. Section 3 will be analogous to Section 2, but for multivariate distributions. The final section will conclude with a discussion of the results and further open questions.

2 Univariate solutions to Stein equations from duality

Stein’s method for univariate distributions has a rich and well developed theory, see for example the survey articles Ley et al., [2017]; Ross, [2011]. In this section we will focus on solutions to Stein equations for univariate distributions using a dual process. We begin by briefly explaining intuition behind the main idea of the duality. The generator of a Markov process has two arguments, a test function $f$ and a point/position in the sample space $x$. Loosely speaking, typically one interprets a generator by first fixing the test function $f$, and the generator in a sense describes how $x$ evolves in time. For our dual process, we reverse this interpretation, we instead fix $x$ and study how the test function $f$ evolves in time. Hence if the generator can be interpreted as or written in the form of a Markov jump process on functions, we can consider the dual process to be this jump process on functions. Now due to the fact that the semi-group of a Markov process can be characterised by its generator, given both the original and dual processes have the same generator, they must also have the same semi-group. We will confirm this fact by comparing our dual process solutions to known solutions using the original process when possible. The following lemma formalises the above argument.

Lemma 2.1**.**

Let $\mathcal{A}$ be a generator of a one dimensional diffusion process and let $\mathcal{D}(\mathcal{A})$ be its domain. Suppose $\mathcal{A}$ generates a Feller semi-group on $\mathcal{D}(\mathcal{A})$ and has unique stationary distribution $\mu$. If for $h(x)=x^{k}\in\mathcal{D}(\mathcal{A})$ where $n\in\{1,2,\ldots\}$, and there exist constants $\{a_{i}\},\{b_{i}\}$, such that $\mathcal{A}f(x)$ can be written as

[TABLE]

*then there exists a Markov jump process $Y(x,t)$ that is dual to $X(x,t)$ where $X(x,t)$ if a diffusion process with generator $\mathcal{A}$ and $X(x,0)=x$, such that the solution $f_{h}(x)$ to the Stein equation (1.1) satisfies *

[TABLE]

*where $Z\sim\mu$, and $Y_{x}(1,t)$ is a Markov jump process on functions of the form $g_{x}(a,k)=ax^{k}$ where $Y_{x}(1,0)=x^{k}$ and generator *

[TABLE]

Proof.

The condition that we require a Feller semi-group ensures the existence of the solution to the Stein equation [Ethier and Kurtz,, 1986, Proposition 1.5(a)]. Essentially it remains to show that the two semi-groups of $\mathcal{A}$ and $\mathcal{B}$ are the same, that is for all $t>0$, $\mathbbm{E}h(X_{x}(t))=\mathbbm{E}Y_{x}(1,t)$. [Ethier and Kurtz,, 1986, Proposition 1.5 (c)] yields

[TABLE]

hence

[TABLE]

However it follows from the definition of $\mathcal{B}$ and the equality of the two generators, that the process $Y_{x}(1,t)$ will satisfy exactly the same differential form. Given they both have the same initial value, the two semi-groups must be identical. ∎

2.1 Univariate examples

We first define some common notation which we will be using for all the examples. $\mathcal{A}$ will denote the characterising diffusion operator, $Y_{x}(1,t)$ the corresponding dual process and $Z$ a random variable following the target distribution. Our test functions will always be of the form $h(x)=x^{k}$ unless otherwise specified. Note that this can easily be extended to polynomial test functions.

Example 2.2** (Normal distribution).**

Recall from the introduction

[TABLE]

and the description of the dual process therein. Note that $\mathbbm{E}(Z^{2k})=(2k-1)!!$ and $\mathbbm{E}(Z^{2k+1})=0$. Then for $h(x)=x^{2k+1}$,

[TABLE]

where $\tau_{0}=0$ and $\tau_{i+1}-\tau_{i}\sim\mathop{\mathrm{Exp}}(2k-2i+1)$. Similarly, for $h(x)=x^{2k}$,

[TABLE]

where $\tau_{0}=0$ and $\tau_{i+1}-\tau_{i}\sim\mathop{\mathrm{Exp}}(2k-2i)$ and noting $Y_{x}(1,t)=(2k-1)!!$ for $t\geq\tau_{k}$. Routine calculations can confirm that these solutions coincide with (1.6).

Example 2.3** (Gamma distribution).**

Luk, [1994] developed Stein’s method for the Gamma distribution using the generator method and we follow the same parameterisation. The Gamma distribution with shape parameter $r>0$ and scale parameter $\lambda>0$ has density function,

[TABLE]

Furthermore if $X\sim\mathrm{Gamma}(r,\lambda)$, then $\mathbbm{E}X^{k}=\frac{\Gamma(r+k)}{\lambda^{k}\Gamma(r)}$. The generator with the Gamma distribution as its stationary distribution is

[TABLE]

Hence for $f(x)=x^{k}$,

[TABLE]

The dual process can be therefore described as follows. Given the process is in state $ax^{k}$, it waits an $\mathop{\mathrm{Exp}}(\lambda n)$ distributed amount of time, before transitioning to $\frac{a(k-1+r)}{\lambda}x^{k-1}$, and the process will be eventually absorbed at the constant $\frac{\Gamma(r+k)}{\lambda^{k}\Gamma(r)}$.

Let $\tau_{0}=0$ and $\tau_{i+1}-\tau_{i}\sim\mathop{\mathrm{Exp}}(\lambda(k-i))$. Then for $h(x)=x^{k}$,

[TABLE]

Remark 2.4**.**

[Luk,, 1994, Theorem 2.6] shows that for a smooth test function with bounded derivatives $h$, the $j$-th derivative of the solution to the Stein equation satisfies

[TABLE]

where $\|\cdot\|$ denotes the sup norm and $f^{(j)}$ denotes the $j$-th derivative. Our test functions do not have bounded derivatives, with exception of the case where $h(x)=x$ and for the first derivative. In which case our bounds yield $f_{h}^{(1)}(x)=\frac{1}{\lambda}$, so the two bounds are identical.

Example 2.5** (Beta distribution).**

Döbler, [2015]; Goldstein and Reinert, [2013]; Gan et al., [2017] developed Stein’s method for the Beta$(\alpha,\beta)$ distribution. The generator presented is as follows:

[TABLE]

This is the generator of (one of the many forms of) a neutral Wright-Fisher diffusion. There are a variety of dual processes associated with various forms of Wright-Fisher diffusion, see for example Kingman, [1982]; Griffiths, [1980]; Griffiths et al., [2016] for just a small sample. Our form of duality will be similar in flavour and analogous to the ancestral/looking backwards in time Kingman’s coalescent. Set $f(x)=x^{k}$,

[TABLE]

Our dual process can therefore be described in the following manner. For the function $h(x)=ax^{k}$,

• •

At a rate of $k(k-1)$, $ax^{k}$ transitions to $ax^{k-1}$,

• •

At a rate of $k(\alpha+\beta)$, $ax^{k}$ transitions to $\frac{\alpha}{\alpha+\beta}ax^{k-1}$.

Notably the net transition rates of this process ($k(k-1+\alpha+\beta)$) are equivalent to Kingman’s coalescent.

Let $\tau_{0}=0$ and $\tau_{i+1}-\tau_{i}\sim\mathop{\mathrm{Exp}}((k-i)(k-i-1+\alpha+\beta))$. Then given $Y_{x}(1,0)=x^{k}$,

[TABLE]

Similarly, for $1\leq j\leq k$,

[TABLE]

Notably, $\mathbbm{E}Y_{x}(1,\tau_{k})=\prod_{r=0}^{k-1}\frac{\alpha+r}{\alpha+\beta+r}=\mathbbm{E}Z^{k}$, as expected and required.

[TABLE]

Given the Beta distribution has finite support on $[0,1]$ we can use the above solution to find uniform bounds for $f_{h}(x)$ and its derivatives. For example, noting that $x\leq 1$, the first derivative satisfies

[TABLE]

Similar computations show for the $n$-th derivative,

[TABLE]

Notably, these bounds correspond exactly to the bounds in [Gan et al.,, 2017, Theorem 5]. The bounds in Gan et al., [2017] are more general, they essentially cover all smooth test functions, where as the approach in this paper will only cover polynomial test functions. However, given polynomials are dense we have not lost too much, and this proof is far shorter and less complicated than the proof presented in the aforementioned work which involves couplings and coalescent theory.

3 Multivariate solutions to Stein equations from Duality

Results for Stein’s method for multivariate distributions remain relatively limited. Barbour et al., 2018a ; Barbour et al., 2018b  provide a framework for multivariate discrete distributions that is applicable in many scenarios. Outside a few examples such as multivariate Normal and Dirichlet, Stein’s method for multivariate continuous distributions is underdeveloped. The generator method is the primary approach multivariate setting as nothing fundamentally changes when considering multivariate Markov processes compared to the univariate case, whereas often more purely analytic approaches become intractable. In a similar vein, our dual process approach will be equally effective in the multivariate setting without too much additional complexity. We omit a multivariate version of Lemma 2.1 as it is completely analogous and not particularly enlightening, and we instead focus on examples.

3.1 Multivariate examples

In the multivariate setting we need to study not only the convergence of the moments of the marginals, but also convergence of all mixed moments. Hence we will be considering test functions of the form $f(\textnormal{{x}})=\prod_{i=1}^{p}w_{i}^{k_{i}}$ for some $p$ and $k_{i}$.

Example 3.1** (Multivariate Normal).**

For $p$-dimensional Normal approximation with mean [math] and covariances $\{\sigma_{ij}\}^{p}_{i,j=1}$, the $p$-dimensional Ornstein-Uhlenbeck generator can be written in the form

[TABLE]

For $f(\textnormal{{x}})=\prod_{i=1}^{p}w_{i}^{k_{i}}$, and since the (multivariate) marginal distributions of the Normal distribution are still (multivariate) Normally distributed, without loss of generality assume that $k_{i}\geq 1$ for all $i$, and denote $f^{\{ij\}}(\textnormal{{x}})=\frac{f(\textnormal{{x}})}{w_{i}w_{j}}$. Then,

[TABLE]

where $\delta_{ij}$ denotes the Kronecker delta function. Given the above form, the process can be described as follows. For all $1\leq i\leq p$,

• •

At rate $k_{i}(k_{i}-1)$, $f(\textnormal{{x}})$ transitions to $\sigma_{ii}f^{\{ii\}}(\textnormal{{x}})$,

• •

At rate $k_{i}k_{j}$, $f(\textnormal{{x}})$ transitions to $\sigma_{ij}f^{\{ij\}}(\textnormal{{x}})$.

Analogously to Example 2.2 we can then use this dual process to solve for $f_{h}(x)$. We omit the calculation for the solution to the Stein equation as it is similar to Example 2.2 but rather long and tedious, and furthermore there already exists a well known nice closed form solution,

[TABLE]

where Z is a random vector having standard multivariate normal distribution of dimension $p$.

Example 3.2** (Dirichlet distribution).**

The Dirichlet distribution can be viewed as a multivariate extension of the Beta distribution. We follow the parameterisation and generator used in Gan et al., [2017]. The Dirichlet distribution with parameters $\textnormal{{a}}=(a_{1},\dots,a_{K})$, where $a_{1}>0,\ldots,a_{K}>0$, is supported on the $(K-1)$-dimensional open simplex, which we parameterise as

[TABLE]

On $\Delta_{K}$, the Dirichlet distribution with parameter a has density

[TABLE]

where $s=\sum_{i=1}^{K}a_{i}$, and where we set $x_{K}=1-\sum_{i=1}^{K-1}x_{i}$, as we shall often do whenever considering vectors taking values in $\Delta_{K}$.

The Stein operator used is the generator of the neutral Wright-Fisher diffusion with parent independent mutation,

[TABLE]

where $\delta_{ij}$ denotes the Kronecker delta function. Similar the multivariate Normal distribution, the multivariate marginal distributions of the Dirichlet distribution are still Dirichlet distributed, hence we can without loss of generality assume $f(x)=\sum_{i=1}^{K-1}x_{i}^{k_{i}}$ where $k_{i}\geq 1$ for all $i$. As in the previous example set $f^{\{i\}}(\textnormal{{x}})=\frac{f(\textnormal{{x}})}{x_{i}}$. Then

[TABLE]

The dual process is unsurprisingly similar to the dual process in Example 2.5, with exception of an additional ‘interaction’ term in its generator. We describe the process as follows: For all $i,j$,

• •

At rate $k_{i}(k_{i}-1)$, $f(\textnormal{{x}})$ transitions to $f^{\{i\}}(\textnormal{{x}})$,

• •

At rate $sk_{i}$, $f(\textnormal{{x}})$ transitions to $\frac{a_{i}}{s}f^{\{i\}}(\textnormal{{x}})$,

• •

At rate $k_{i}k_{j}$, $f(\textnormal{{x}})$ transitions to 0.

For this example we will take a different approach, rather than explicitly solving for the solution to the Stein equation and then taking derivatives of the solution, we will instead use a coupling argument. Furthermore, we will only focus on bounding the derivative of the solution to the Stein equation rather than the solution itself, as derivative bounds are typically what are required in applications. Let $\varepsilon_{i}$ be the unit vector in the $i$-th direction, then without loss of generality set $i=1$.

[TABLE]

To bound this quantity, we first observe that for $\textnormal{{x}}\in\Delta_{K}$,

[TABLE]

We couple the transition times of the two processes $Y_{\textnormal{{x}}+\varepsilon e_{1}}(1,t)$ and $Y_{\textnormal{{x}}}(1,t)$ exactly, and we can hence focus predominantly upon only the $x_{1}$ term. Furthermore it suffices to ignore the third class of transition where $f(\textnormal{{x}})$ transitions to 0 as an upper bound as this will only increase the coupling time of the processes $(Y_{\textnormal{{x}}+\varepsilon e_{1}}(1,t),Y_{\textnormal{{x}}}(1,t)$), and we apply this to avoid interactions between the different dimensions in our process. Given we only focus on the one dimension, we can ignore all the transitions that involve the other dimensions, and use an identical calculation to (2.2),

[TABLE]

Hence

[TABLE]

As in Example 2.5 this bound corresponds to the bounds in Gan et al., [2017]. We leave the higher order derivatives as an exercise for interested readers. Interestingly this suggests that there may exist an opportunity to improve upon these bounds given we were slightly crude in ignoring the transitions of the type $f(\textnormal{{x}})\to 0$.

Example 3.3** (Poisson-Dirichlet distribution and Dirichlet process).**

Our final example is an application for an infinite dimensional distribution. The original motivation for this work was for the formulation of Stein’s method for Poisson-Dirichlet and Dirichlet process approximation. For the following, full details can be found in a forthcoming preprint Gan and Ross, [2019].

We begin by describing the Poisson-Dirichlet distribution. Define a random interval partition of the interval $[0,1]$. We encode the space of interval partitions as

[TABLE]

and define $\textrm{PD}(\theta)$, the Poisson-Dirichlet distribution on $\nabla_{\infty}$ with parameter $\theta>0$. The Poisson-Dirichlet distribution can be described/constructed in a number of ways, such as the limiting distribution of proportions in integer partitions, the limit as the dimension goes to infinity, of appropriate Dirichlet distributions, the distribution of the normalised points of a Poisson process with intensity measure $\theta e^{-x}/x,x>0$ Kingman, [1975] and via stick breaking schemes and the GEM distribution; see Pitman, [2006] for an overview.

We define a Dirichlet process with parameters $\theta>0$ and $\pi$ a probability measure on a compact metric space $E$, denoted $\textrm{DP}(\theta,\pi)$, to be the random element of $M_{1}(E)$, the space of the probabilities measures on $E$ equipped with the weak topology. We say $\mu\sim\textrm{DP}(\theta,\pi)$ if

[TABLE]

where $(P_{1},P_{2},\ldots)\sim\textrm{PD}(\theta)$ is independent of $\xi_{1},\xi_{2},\ldots$, which are i.i.d. with distribution $\pi$. Thus we think of $\textrm{DP}(\theta,\pi)$ as the Poisson-Dirichlet distribution with independent $\pi$-distributed labels.

The Stein operator used for Dirichlet process approximation is a particular class of Fleming-Viot processes. Define an operator $A$ by

[TABLE]

where the domain of $A$ is the set of continuous functions $C(E)$. We define the generator $\mathcal{A}$ of our Fleming-Viot Markov process; for $F\in\mathcal{D}$ and $\mu\in M_{1}(E)$,

[TABLE]

where $\delta_{x}(\cdot)$ is Dirac measure, and for any measure $\nu\in M_{1}$ we define the directional (Gâteaux) derivatives by,

[TABLE]

with higher order derivatives analogously defined.

For $f:E\mapsto\mathbbm{R}$, let ${\langle f,\mu\rangle}$ denote the expectation of $f$ with respect to $\mu$, that is ${\langle f,\mu\rangle}=\int_{E}f(x)\mu(dx)$. Furthermore, denote let $\mu^{k}$ denote the product measure of $k$ independent copies of $\mu$, and for $f:E^{k}\mapsto\mathbbm{R}$ analogously define ${\langle f,\mu^{k}\rangle}$. The domain of the generator (3.4) we wish to consider is

[TABLE]

where $\mathrm{B}(E^{k})$ is the space of all bounded functions on $E^{k}$. To see how this space of test functions is related to moments, for example set $\varphi(x,y)=\mathbf{1}_{x=y}$, then $F(\mu)=\sum_{i=1}^{\infty}P_{i}^{2}$, where the $P_{i}$ are as in (3.2).

For $F={\langle\varphi,\mu^{k}\rangle}$, simple computations show that our generator (3.4) can be reduced to the form

[TABLE]

where $\Phi_{ij}^{(k)}\varphi(x_{1},\ldots,x_{k-1})=\varphi(x_{1},\ldots,x_{j-1},x_{i},x_{j},\ldots,x_{k-1})$. From the above formulation, we follow [Ethier and Kurtz,, 1993, Section 3] (in turn following Dawson and Hochberg, [1982]), and define our dual Markov jump process on $\cup_{j\geq 1}\mathrm{B}(E^{j})$ with $Y(0)=\varphi$. Informally, there are two types of transition for our jump process on functions

• •

At a rate of 1, for each $1\leq i,j\leq k$, the argument $x_{i}$ in $\varphi(\textnormal{{x}})$ is replaced with $x_{j}$,

• •

At a rate of $\frac{\theta}{2}$, for each $1\leq i\leq k$, the argument $x_{i}$ in $\varphi(\textnormal{{x}})$ is replaced with a $\pi$ distributed random variable.

As stated earlier, for full details about this formulation and how to use this to calculate bounds for derivatives of the Stein solution, see Gan and Ross, [2019].

4 Discussion

As seen by the numerous examples, this duality approach appears to give results that are on par with existing techniques for calculating Stein solutions. To the author’s knowledge, these are the first results where such a duality is used to calculate solutions to the Stein equation. We believe that there is potential for the ideas presented to be pushed and extended in a variety of interesting directions. In particular we present a few open questions that are of interest to the author.

Our focus in this paper was on test functions for moments. Are there other classes of test functions that can be approached using a dual process? For example is it possible to find a duality argument for test functions that correspond with Kolmogorov or Wasserstein distances? 2. 2.

In this paper, we convert difficult problems on diffusions to simpler problems on Markov jump processes. Is it possible to reverse this? Brown and Xia, [2001] give a flexible framework for Stein’s method for Markov birth-death jump processes that is widely applicable, but there are examples such as the compound Poisson distribution that have multiple births and/or death events that remain intractable. Is it possible that problems involving complicated Markov jump processes could be instead approached using a simpler diffusion or perhaps even a simpler jump process? 3. 3.

The key term in our Stein operators that enables the duality argument is the $xf^{\prime}(x)$ since this yields $xf^{\prime}(x)=kf(x)$ for $f(x)=x^{k}$. This is therefore similar to size bias distributions. It would be interesting to explore if there is a direct connection between size biasing and moment duality. 4. 4.

Not all Stein operators have the key $xf^{\prime}(x)$ term, for example the generalised Gamma distribution Peköz et al., [2016] has Stein operator $\mathcal{A}f(x)=f^{\prime\prime}(x)-\left(\frac{\alpha-1}{x}-\beta x^{\beta-1}\right)f^{\prime}(x)$, for some $\alpha,\beta>0$. Does there still exist some form of moment duality outside of the case $\beta=1$? Does there still exist a form of moment duality when we do not have an $xf^{\prime}(x)$ term?

5 Acknowledgments

The author would like to thank Nathan Ross for his contributions to Example 3.3 and for helpful discussions.