# Limit behaviour of random walks on $\mathbb{Z}^{m}$ with two-sided membrane

Victor Bogdanskii111National Technical University of Ukraine “Igor Sikorsky
Kyiv Polytechnic Institute”, ave. Pobedy 37, Kiev 03056, Ukraine;
<EMAIL_ADDRESS>Ilya Pavlyukevich222Institute for Mathematics, Friedrich
Schiller University Jena, Ernst–Abbe–Platz 2, 07743 Jena, Germany;
<EMAIL_ADDRESS>and Andrey Pilipenko333Institute of Mathematics,
National Academy of Sciences of Ukraine, Tereshchenkivska Str. 3, 01601, Kiev,
Ukraine 444National Technical University of Ukraine “Igor Sikorsky Kyiv
Polytechnic Institute”, ave. Pobedy 37, Kiev 03056, Ukraine;
<EMAIL_ADDRESS>
(August 27, 2024)

###### Abstract

We study Markov chains on $\mathbb{Z}^{m}$, $m\geq 2$, that behave like a
standard symmetric random walk outside of the hyperplane (membrane)
$H=\\{0\\}\times\mathbb{Z}^{m-1}$. The transition probabilities on the
membrane $H$ are periodic and also depend on the incoming direction to $H$,
what makes the membrane $H$ two-sided. Moreover, sliding along the membrane is
allowed. We show that the natural scaling limit of such Markov chains is a
$m$-dimensional diffusion whose first coordinate is a skew Brownian motion and
the other $m-1$ coordinates is a Brownian motion with a singular drift
controlled by the local time of the first coordinate at $0$. In the proof we
utilize a martingale characterization of the Walsh Brownian motion and
determine the effective permeability and slide direction. Eventually, a
similar convergence theorem is established for the one-sided membrane without
slides and random iid transition probabilities.

Keywords: Skew Brownian motion; Walsh’s Brownian motion; perturbed random
walk; two-sided membrane; weak convergence; martingale characterization.

2020 Mathematics Subject Classification: 60F17∗, 60G42, 60G50, 60H10, 60J10,
60J55, 60K37.

##  1 Introduction

A multidimensional Brownian motion is a fundamental stochastic process that
describes an idealized mathematical model of a free physical diffusion in a
homogeneous medium. Having in mind the observations of pollen made my Brown or
a later kinetic theory of gases, we can interpret a diffusion as a collective
motion of independent random walkers whose distribution density in space obeys
the isotropic Gaussian distribution.

However in real physical or biological systems, the space is often separated
into compartments by _membranes_ that impede or facilitate the passage of the
walker and create an anisotropy in the walkers’ collective motion. From the
physical point of view, a membrane is a thin slice of a material whose
diffusivity is different from the diffusivity of the environment. Diffusions
through membranes are observed in biological tissues where they control the
transport of ions, water molecules and gases, or in porous and composite
materials.

A rigorous mathematical justification of the interpretation of a diffusion as
a limit of scaled random walks is given by the Functional Central Limit
Theorem (FCLT, the Donsker–Prokhorov invariance principle). For $m\geq 1$, let
$\\{\mathrm{e}_{1},\dots,\mathrm{e}_{m}\\}$ be a standard basis in
$\mathbb{R}^{m}$. Consider a symmetric random walk $Z=(Z(n))_{n\geq 0}$ on
$\mathbb{Z}^{m}$ defined by identical one-step transition probabilities

$\displaystyle\mathbf{P}\Big{(}Z(n+1)-Z(n)=\pm\mathrm{e}_{k}\Big{|}Z(0),\dots,Z(n)\Big{)}=\frac{1}{2m},\quad
k=1,\dots m,\quad n\geq 0.$ (1.1)

Then the FCLT yields the weak convergence in the uniform topology

$\displaystyle\Big{(}\frac{Z([nt])}{\sqrt{n}}\Big{)}_{t\geq
0}\Rightarrow\Big{(}\frac{1}{\sqrt{m}}W(t)\Big{)}_{t\geq 0},$ (1.2)

where $W$ is a standard $m$-dimensional Brownian motion.

In this paper we consider a novel class of scaling limits of random walks in
the presence of extended spatially non-homogeneous one- or two-sided barrier
(membrane, interface). More precisely, in the $m$-dimensional space we
consider an semi-permeable hyperplane that separates the space into two half-
space compartments. The random walker will perform the symmetric random walk
outside of the membrane, however the probability of the passage through the
membrane will depend on the hitting position of the membrane by the walker and
on its incoming direction. Thus the probability to penetrate the membrane will
depend of the fact whether the walker has reached it from the “right” of from
the “left”. Furthermore we will assume that the membrane has a spatially
regular structure according to one of the following setting. First, we assume
that the penetration probabilities are periodic in space, so that the membrane
reminds of a two-sided fabric. Alternatively, the membrane can imitate a
“random” perforated surface so that the penetration probabilities are iid
random variables. We will show the scaled limit of such Markov chains is a
skew Brownian motion in the direction perpendicular to the membrane. The other
coordinates converge to a standard Brownian motion, maybe, with a singular
drift given by the local time of the first coordinate at the origin.

The proofs of these results are purely probabilistic and employ methods of
homogenization and dynamics of singular differential equations. Our results
give a path-wise picture of the diffusion through a two-side semi-permeable
interface. In the physical language this corresponds to the
Langevin–Smoluchowsky approach to diffusions. It should be emphasized that
physical papers (see, e.g., Novikov et al. [24], Grebenkov et al. [7], Moutal
and Grebenkov [21] and references therein) devoted to similar problems use
analytical methods, mainly the analysis of the Fokker–Planck equation.

One-dimensional locally perturbed random walks were considered from different
points of view in Harrison and Shepp [10], Minlos and Zhizhina [20], Pilipenko
and Pryhod’ko [28], Pilipenko and Prikhod’ko [31], Pilipenko and Sakhanenko
[29], Ngo and Paigné [23]. In this paper, we use the multidimensional
martingale characterization approach previously considered in Iksanov and
Pilipenko[11] in dimension one.

It should be noted that if transition probabilities of a multidimensional
random walk are perturbed on a finite set or on a hyperplane of co-dimension
2, then under some natural assumptions its scaling limit is a Brownian motion,
see Szász and Telcs [35], Yarotskii [37], and Paulin and Szász [25].

We also refer the reader’s attention to the following related mathematical
works: the monograph by Portenko [33] and the research papers by Lejay [18],
Mandrekar and Pilipenko [19], Aryasova and Pilipenko [1], Iksanov et al. [12],
Pilipenko and Khomenko [30], Pilipenko [27], and Pilipenko and Prykhodko [32].

Notation. The weak convergence in the Skorokhod space
$D([0,\infty),\mathbb{R}^{m};J_{1})$ is denoted by $\Rightarrow$. It should be
noted, however, that all limit processes in this paper are continuous. The
convergence in distribution of random variables is denoted by
$\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}$.

Acknowledgements. I.P. and A.P. acknowledge support by the DFG project PA
2123/7-1. A.P. thanks the Institute of Mathematics of the FSU Jena for
hospitality. V.B. and A.P. acknowledge support by the National Research
Foundation of Ukraine (project 2020.02/0014 “Asymptotic regimes of perturbed
random walks: on the edge of modern and classical probability”).

##  2 Two-sided periodic membrane: the model and the main result

Let $m\geq 2$ and let $\\{\mathrm{e}_{1},\dots,\mathrm{e}_{m}\\}$ be a
standard basis in $\mathbb{R}^{m}$. Consider a Markov chain
$Z=\\{Z(n)\\}_{n\geq 0}$ on $\mathbb{Z}^{m}$ that behaves as a simple random
walk outside of the hyperplane $H:=\\{0\\}\times\mathbb{Z}^{m-1}$, i.e., for
all $k=1,\dots,m$

$\mathbf{P}\Big{(}Z({n+1})=z\pm\mathrm{e}_{k}\Big{|}Z(n)=z\Big{)}=\frac{1}{2m},\
z\notin H.$ (2.1)

For each $n\geq 1$, we denote the first coordinate of the process $Z$ by $X$,
and the other $(m-1)$ coordinates by $Y$ so that $Z=(X,Y)$.

We will interpret $H$ as a semipermeable two-sided membrane that may let a
particle pass from one half-space to another with certain probabilities that
can depend on the crossing direction. Moreover, the particle can “slide” along
the membrane.

The membrane has to be homogeneous, i.e., the transition probabilities are
periodic in space.

Notice that if the membrane is two-sided, then $Z$ is not a Markov chain,
generally. Indeed, its position upon leaving the membrane is determined by
both current the location on the membrane and on the particle’s incoming
direction. Hence, in order to introduce a Markov structure we have to enlarge
the state space by splitting the membrane $H$ into two parts $\mathcal{H}^{-}$
and $\mathcal{H}^{+}$ corresponding to its “right” and “left” sides.

To formalize the setting, we consider the Markov chain
$\mathcal{Z}=(\mathcal{Z}(n))_{n\geq 0}$ on the state space

$\displaystyle\\{\dots,-3,-2,-1,-0,+0,1,2,3,\dots\\}\times\mathbb{Z}^{m-1}.$
(2.2)

We also write $\mathcal{Z}=(\mathcal{X},Y)$, so that the process $Y$ coincides
with the process $Y$ of the original random walk $Z$ whereas the process
$\mathcal{X}$ is defined on the enlarged space
$\\{\dots,-3,-2,-1,-0,+0,1,2,3,\dots\\}$.

Denote the “right” and “left” hand side of the membrane by

$\displaystyle\mathcal{H}^{-}:=\\{-0\\}\times\mathbb{Z}^{m-1},\quad\mathcal{H}^{+}:=\\{+0\\}\times\mathbb{Z}^{m-1},$
(2.3)

and we set $\mathcal{H}:=\mathcal{H}^{-}\cup\mathcal{H}^{+}$.

We assume that transition probabilities of $\mathcal{Z}$ outside of
$\mathcal{H}$ satisfy (2.1), namely

$\mathbf{P}\Big{(}\mathcal{Z}({n+1})=z\pm\mathrm{e}_{k}\,\Big{|}\,\mathcal{Z}(n)=z\Big{)}=\frac{1}{2m},\
z\notin\mathcal{H},$ (2.4)

where we agree that for each $y\in\mathbb{Z}^{m-1}$

$\displaystyle(1,y)-\mathrm{e}_{1}=(+0,y)$
$\displaystyle\text{and}\quad(-1,y)+\mathrm{e}_{1}=(-0,y),$ (2.5)
$\displaystyle(+0,y)+\mathrm{e}_{1}=(1,y)$
$\displaystyle\text{and}\quad(-0,y)+\mathrm{e}_{1}=(1,y),$
$\displaystyle(+0,y)-\mathrm{e}_{1}=(-1,y)$
$\displaystyle\text{and}\quad(-0,y)-\mathrm{e}_{1}=(-1,y).$

Suppose also that the perturbed random walk does not stay on the membrane,
i.e.,

$\displaystyle\mathbf{P}\Big{(}\mathcal{Z}({n+1})\in\\{\pm
1\\}\times\mathbb{Z}^{m-1}\,\Big{|}\,\mathcal{Z}(n)=z\Big{)}=1,\quad
z\in\mathcal{H}.$ (2.6)

However it is allowed that upon hitting the membrane (from the left or from
the right) the particle can “slide” along the membrane, so that its
$y$-coordinate changes. We assume that transition probabilities from
$\mathcal{H}$ have periodic structure.

Notation. For $k_{2},\dots,k_{m}\geq 1$ fixed, let $U$ denote the box in
$\mathbb{Z}^{m-1}$ defined by

$\displaystyle U:=[0,k_{2}-1]\times\cdots\times[0,k_{m}-1].$ (2.7)

For $y\in\mathbb{Z}^{m-1}$, $y=(y_{2},\dots,y_{m})$, and $j\in U$,
$j=(j_{2},\dots,j_{m})$, we say that

$\displaystyle y\equiv j\quad\text{if}\quad y_{i}\ (\textrm{mod}\
k_{i})=j_{i},\ i=2,\dots,m.$ (2.8)

For each $j\in U$ we set

$\displaystyle H_{j}^{\pm}:=\\{(\pm 0,y)\in H^{\pm}\colon y\equiv j\\}.$ (2.9)

Without loss of generality let us assume that the random walk $\mathcal{Z}$
starts on the hyperplane $\mathcal{H}$, i.e., $\mathcal{X}(0)=\pm 0$.

The following are our key assumptions concerning the transition probabilities
of the random walk in the membrane.

A${}_{\text{periodic}}$. Periodicity of the transition probabilities. We
assume that there exist $k_{2},\dots,k_{m}\geq 1$ such that for all
$l_{2},\dots,l_{m}\in\mathbb{Z}$, for all $z_{0}\in\mathcal{H}$, for all
$z_{1}\in\\{-1,+1\\}\times\mathbb{Z}^{m-1}$

$\displaystyle\mathbf{P}\Big{(}\mathcal{Z}({n+1})=z_{1}$
$\displaystyle\Big{|}\,\mathcal{Z}(n)=z_{0}\Big{)}$ (2.10)
$\displaystyle=\mathbf{P}\Big{(}\mathcal{Z}(n+1)=z_{1}+k_{2}l_{2}\mathrm{e}_{2}+\cdots+k_{m}l_{m}\mathrm{e}_{m}\,\Big{|}\,\mathcal{Z}(n)=z_{0}+k_{2}l_{2}\mathrm{e}_{2}+\cdots+k_{m}l_{m}\mathrm{e}_{m}\Big{)}.$

Aγ. To describe the transitions through the membrane, we denote by
$0=\tau_{0}<\tau_{1}<\cdots$ the successive arrivals of $\mathcal{Z}$ to
$\mathcal{H}$ or, equivalently, of $\mathcal{X}$ to $\\{-0,+0\\}$.

On the finite state space $\\{-0,+0\\}\times U$ we consider an auxiliary
embedded process $\hat{\mathcal{Y}}=(\hat{\mathcal{Y}}(n))_{n\geq 0}$ as
follows. We set

$\displaystyle\hat{\mathcal{Y}}(n)=\begin{cases}(-0,j),\ \text{ if
}Y(\tau_{n})\equiv j\in U\text{ and }X(\tau_{n}+1)=-1,\\\ (+0,j),\ \text{ if
}Y(\tau_{n})\equiv j\in U\text{ and }X(\tau_{n}+1)=1.\end{cases}$ (2.11)

It follows from the strong Markov property of $\mathcal{Z}$ that
$\hat{\mathcal{Y}}$ is a Markov chain. It is easy to see that all the states
of the set $\\{-0,+0\\}\times U$ are connected for $\hat{\mathcal{Y}}$. Hence
there is a unique stationary distribution

$\displaystyle\pi=\\{\pi_{\\{-0\\}\times j},\pi_{\\{+0\\}\times j}\\}_{j\in
U}$ (2.12)

of $\hat{\mathcal{Y}}$ on $\\{-0,+0\\}\times U$. By the strong law of large
numbers for Markov chains, for each $j\in U$

$\lim_{n\to\infty}\frac{1}{n}\sum_{k=0}^{n}1\\!\\!\\!\;{\rm
I}(\hat{\mathcal{Y}}(k)=\\{\pm 0\\}\times j)=\pi_{\\{\pm 0\\}\times j}\ \
\text{a.s.}$ (2.13)

With the help of the stationary distribution $\pi$ we introduce the _effective
permeability_

$\displaystyle\gamma$ $\displaystyle:=\sum_{j\in U}(\pi_{\\{+0\\}\times
j}-\pi_{\\{-0\\}\times j})\in[-1,1].$ (2.14)

Ac. Finally, we describe the “slides” along the membrane. We denote

$\displaystyle\alpha_{\pm,j}$
$\displaystyle:=\mathbf{E}\Big{[}Y(1)-Y(0)\,\Big{|}\,\mathcal{X}(0)=\pm
0,Y(0)=j\Big{]}\in\mathbb{R}^{m-1},\quad j\in U,$ (2.15)

the mean slide sizes along the “right” of the “left” side of the membrane, and
assume that they are finite. Introduce the _effective slide_ as

$\displaystyle c$
$\displaystyle:=\mathbf{E}_{\pi}\Big{[}Y(1)-Y(0)\Big{]}=\sum_{j\in
U}\Big{(}\pi_{\\{+0\\}\times j}\cdot\alpha_{+,j}+\pi_{\\{-0\\}\times
j}\cdot\alpha_{-,j}\Big{)}\in\mathbb{R}^{m-1}.$ (2.16)

Now we are ready to formulate the main result of the paper. Let $W_{X}$ and
$W_{Y}$ be independent standard Brownian motions in $\mathbb{R}$ and
$\mathbb{R}^{m-1}$ respectively.

For $\gamma\in[-1,1]$ defined in (2.14), we consider the skew Brownian motion
$X^{\gamma}$ that is a unique strong solution of the SDE

$\displaystyle X^{\gamma}(t)=\gamma L(t)+W_{X}(t),\quad t\geq 0.$ (2.17)

where $L$ is the symmetric two-sided local time of $X^{\gamma}$ at zero, see
Harrison and Shepp [10]. Other characterizations and properties of the skew
Brownian motion can be found in the review by Lejay [17].

Furthermore, let

$\displaystyle Y^{c}(t):=cL(t)+W_{Y}(t),\quad t\geq 0,$ (2.18)

i.e., the process $Y^{c}$ is a $(m-1)$-dimensional Brownian motion that slides
in the direction $c$ in the time instants when $X^{\gamma}$ touches zero. Note
that (2.18) is not a stochastic differential equation because the process $L$
is already determined in (2.17).

For the original process $Z=(X,Y)$, denote the rescaled continuous time
processes

$\displaystyle X_{n}(t):=\frac{X([nt])}{\sqrt{n}},\quad
Y_{n}(t):=\frac{Y([nt])}{\sqrt{n}},\quad t\geq 0.$ (2.19)

######  Theorem 2.1.

Let $Z=(X,Y)$ be a perturbed random walk satisfying the preceding assumptions.
Then for any initial value $Z(0)\in\mathbb{Z}^{m}$, the weak convergence holds
true:

$\displaystyle(X_{n},Y_{n})\Rightarrow\frac{1}{\sqrt{m}}(X^{\gamma},Y^{c}),\quad
n\to\infty,$ (2.20)

where the processes $X^{\gamma}$ and $Y^{c}$ are defined in (2.17) and (2.18),
and $\gamma$ and $c$ are defined in (2.14) and (2.15) respectively.

The crux of the Theorem is transparent. Away from the membrane $H$, the
limiting process $(X^{\gamma},Y^{c})$ coincides with the standard
$m$-dimensional Brownian motion $(W_{X},W_{Y})/\sqrt{m}$. The perturbation of
the transition probabilities on the two-sided membrane results in the
appearance of a singular drift in the direction $x$ perpendicular to the
membrane. Hence the $x$-coordinate of the limiting process becomes a skew
Brownian motion with the effective permeability parameter $\gamma\in[-1,1]$.
The $y$-coordinates also become a singular drift in the effective sliding
direction $c\in\mathbb{R}^{m-1}$. This drift equals zero as long as the
limiting process stays away from the membrane. However upon hitting the
membrane, the limiting process performs a singular “sliding” in the direction
$c$ controlled by the local time at zero of the $x$-coordinate. Note that in
the limit, the two-sided membrane structure disappears, and the limiting
process $(X^{\gamma},Y^{c})$ is Markov in $\mathbb{R}^{m}$.

We illustrate Theorem 2.1 by an example.

######  Example 2.2.

Consider a perturbed random walk $Z=(X,Y)$ on the plane $\mathbb{Z}^{2}$,
$m=2$, separated by a two-sided membrane $\mathcal{H}$. Assume that the
membrane is $2$-periodic, i.e., $k_{1}=2$ and the set $U=\\{0,1\\}$.

Assume that upon hitting the membrane from the “right” the random walk gets
reflected in the even points and surely penetrates the membrane in the odd
points. The penetration probabilities from the “left” are given at Fig. 1 a).

Let

$\displaystyle\tau=\inf\\{k\geq
0\colon\mathcal{Z}(n)=(\mathcal{X}(n),Y(n))\in\mathcal{H}\\}$ (2.21)

be the first hitting time of the Markov chain $\mathcal{Z}$ of the two-sided
membrane $\mathcal{H}$. Since away of the membrane $\mathcal{H}$, the
increments of $\mathcal{Z}$ coincide with the increments of a translation
invariant two-dimensional symmetric random walk $Z$, we can easily calculate
the probabilities of hitting the membrane in an even or an odd point:

$\displaystyle\alpha$ $\displaystyle:=\mathbf{P}_{(1,0)}(Y(\tau)\equiv
0)=\mathbf{P}_{(1,1)}(Y(\tau)\equiv 1)=\mathbf{P}_{(-1,0)}(Y(\tau)\equiv
0)=\mathbf{P}_{(-1,1)}(Y(\tau)\equiv 1)=2-\sqrt{2},$ (2.22) $\displaystyle
1-\alpha$ $\displaystyle:=\mathbf{P}_{(1,0)}(Y(\tau)\equiv
1)=\mathbf{P}_{(1,1)}(Y(\tau)\equiv 0)=\mathbf{P}_{(-1,0)}(Y(\tau)\equiv
1)=\mathbf{P}_{(-1,1)}(Y(\tau)\equiv 0)=\sqrt{2}-1.$

With the help of (2.22) we calculate the transition probabilities of the
embedded Markov chain $\hat{\mathcal{Y}}$ on $\\{-0,+0\\}\times\\{0,1\\}$:

$\displaystyle\mathbb{P}=\begin{array}[]{c|cccc}&(-0,0)&(-0,1)&(+0,0)&(+0,1)\\\
\hline\cr\vphantom{1^{|}}(-0,0)&(1-p)(1-\alpha)&(1-p)\alpha&p\alpha&p(1-\alpha)\\\
(-0,1)&p\alpha&p(1-\alpha)&(1-p)(1-\alpha)&(1-p)\alpha\\\
(+0,0)&0&0&\alpha&1-\alpha\\\ (+0,1)&1-\alpha&\alpha&0&0\end{array}$ (2.23)

Solving the forward Kolmogorov equation $(\mathbb{P}^{T}-\text{Id})\pi=0$ we
obtain the stationary law $\pi$ of $\hat{\mathcal{Y}}$:

$\displaystyle\pi_{(-0,0)}$
$\displaystyle=\frac{(1-\alpha)(1-\alpha+(2\alpha-1)p)}{(1-\alpha)(2\alpha+1)+(2\alpha-1)p^{2}},\quad\pi_{(-0,1)}=\frac{(1-\alpha)\alpha}{(1-\alpha)(2\alpha+1)+(2\alpha-1)p^{2}},\quad$
(2.24) $\displaystyle\pi_{(+0,0)}$
$\displaystyle=\frac{\alpha(1-\alpha+(2\alpha-1)p^{2})}{(1-\alpha)(2\alpha+1)+(2\alpha-1)p^{2}},\quad\pi_{(+0,1)}=\frac{(1-\alpha)(\alpha-(2\alpha-1)(1-p)p)}{(1-\alpha)(2\alpha+1)+(2\alpha-1)p^{2}},$

and calculate the effective permeability and the effective slide according to
(2.14) and (2.15):

$\displaystyle\gamma$
$\displaystyle=\frac{(2\alpha-1)(\alpha(2p-1)+(p-1)^{2})}{(1-\alpha)(2\alpha+1)+(2\alpha-1)p^{2}},$
(2.25) $\displaystyle c$
$\displaystyle=-\frac{(1-\alpha)((2\alpha-1)p+1)}{(1-\alpha)(2\alpha+1)+(2\alpha-1)p^{2}}.$

$1$$p$$1-p$$1$$1-p$$p$$-2$$-1$$0$$1$$2$$3$$-2$$-1$$-0$$+0$$1$$2$a)

$1$$1-p$$p$$1$$p$$1-p$$-2$$-1$$0$$1$$2$$3$$-2$$-1$$-0$$+0$$1$$2$b)

Figure 1: Two models with a two-sided 2-periodic semi-permeable membrane. The
left-hand sides of the membranes a) and b) are shifted by 1 along the
$y$-axis.

######  Example 2.3.

Consider the random walk from Example 2.2, where transition probabilities from
the “right” side of the membrane are the same, and transition probabilities
from the “left” side are shifted by 1 along the $y$-axis, see Fig. 1 b).

The results for this model are obtained immediately by interchanging $p$ and
$1-p$ in (2.25). Notice that parameters of effective permeability and
effective slide are different, i.e., they depends on mutual disposition of the
“left” and the “right” sides of the membrane.

The idea of the proof of the Theorem 2.1 consists in a decomposition of the
process $\mathcal{Z}$ into excursions starting and ending on the membrane
$\mathcal{H}$. The excursions have a probability law of excursions of the
symmetric random walk.

If the membrane is homogeneous ($k_{2}=\cdots=k_{m}=1$), then it is well know
that the limit of $x$ coordinate is a skew Brownian motion, see Harrison and
Shepp [10]. To control a slide in $y$ direction we only have to control the
number of visits to 0 of $x$-coordinate of the walk. However, if the membrane
is periodic with a non-trivial period, then sign of the $x$-coordinate of an
excursion is selected in accordance with the transition probabilities (2.10),
i.e., its sign depends on the position of the random walk at the last visit to
the membrane. The slide along $y$ depends on that position too.

Due to the periodicity assumption (2.10), there is $d=2|U|:=2k_{2}\cdots
k_{m}$ different types of excursions between consecutive visits of the
membrane. In order to treat this number of excursions at the same time we have
to consider a natural generalization of the skew Brownian motion, namely the
Walsh’s Brownian motion. Hence we will show that the family of $d$ one-sided
random walks converge to a Walsh’s Brownian motion. This will allow us to
derive the effective permeability and sliding parameters.

The paper is organized as follows. In the next Section we introduce Walsh’s
Brownian motion and give its convenient realization as a $d$-dimensional
stochastic process that takes values on the positive coordinate half-axes. We
will also formulate two martingale characterizations of Walsh’s Brownian
motion. Section 4 is devoted to the proof of Theorem 2.1. In the final Section
5 we show how our method can be applied to a model of a one-sided membrane
that possesses some ergodic properties.

##  3 Walsh’s Brownian motion and its martingale characterizations

Walsh’s Brownian motion (WBM) was introduced in the Epilogue of Walsh [36] as
a diffusion on $d$ rays on a two-dimensional plane with the common origin. On
each ray WBM, is a standard one-dimensional Brownian motion that however can
change the ray upon hitting the origin, i.e., each ray is characterized by a
weight $p_{i}>0$, $i=1,\dots,d$, $p_{1}+\cdots+p_{d}=1,$ that heuristically
can be understood as a probability to go on the ray number $i$. Hence, the
conventional WBM is a process $X$ on the plane expressed in polar coordinates
as $X=(R_{t},\theta_{t})$ where $R$ is the reflecting Brownian motion and
$\theta$ is a random process taking values on the set of $d$ angles on
$[0,2\pi)$ and being constant during each excursion of $R$ from $0$. This
representation has been used in various works on WBM including Barlow et al.
[2], Freidlin and Sheu [5], Hajri [8], Hajri and Touhami [9], and Karatzas and
Yan [15].

In this paper prefer to embed the WBM into a $d$-dimensional Euclidean space
as it was indicated in Walsh [36]. To this purpose, let $E$ be the union of
non-negative coordinate half-axes in $\mathbb{R}^{d}$, i.e.,

$\displaystyle E=\\{x\in\mathbb{R}^{d}\colon x_{i}\geq 0\text{ and
}x_{i}x_{j}=0,\ i\neq j,\ i,j=1,\dots,d\\}.$ (3.1)

We fix probabilities $p_{1},\dots,p_{d}>0$, $p_{1}+\cdots+p_{d}=1$, and also
denote $q_{i}=1-p_{i}$, $i=1,\dots,d$. Let also
$(\Omega,\mathscr{F},\mathbb{F},\mathbf{P})$ be a filtered probability space
satisfying the usual hypotheses. Then adopting the Markovian characterization
of the WBM from Barlow et al. [2], we say that WBM is a continuous Feller
Markov process $X=(X_{1},\dots,X_{d})$ on $E$ with the one-dimensional laws
given by

$\displaystyle\mathbf{E}_{0}\mathrm{e}^{\lambda_{1}X_{1}(t)+\cdots+\lambda_{d}X_{d}(t)}$
$\displaystyle=\sum_{k=1}^{d}p_{k}\mathbf{E}_{0}\mathrm{e}^{\lambda_{k}|W(t)|},$
(3.2)
$\displaystyle\mathbf{E}_{x}\mathrm{e}^{\lambda_{1}X_{1}(t)+\cdots+\lambda_{d}X_{d}(t)}$
$\displaystyle=\mathbf{E}_{x_{j}}\Big{[}1\\!\\!\\!\;{\rm
I}(t<\tau_{0})\mathrm{e}^{\lambda_{j}W(t)}\Big{]}+\sum_{k=1}^{d}p_{k}\mathbf{E}_{x_{j}}\Big{[}1\\!\\!\\!\;{\rm
I}(t\geq\tau_{0})\mathrm{e}^{\lambda_{k}|W(t)|}\Big{]}$
$\displaystyle=\mathbf{E}_{x_{j}}\Big{[}\sum_{k=1}^{d}p_{k}\mathrm{e}^{\lambda_{k}|W(t)|}\Big{]}+\mathbf{E}_{x_{j}}\Big{[}1\\!\\!\\!\;{\rm
I}(t<\tau_{0})\Big{(}\mathrm{e}^{\lambda_{j}W(t)}-\sum_{k=1}^{d}p_{k}\mathrm{e}^{\lambda_{k}W(t)}\Big{)}\Big{]},$
$\displaystyle x$ $\displaystyle=(0,\dots,x_{j},\dots,0),\
x_{j}>0,\lambda\in\mathbb{C}^{d},$

$W$ is a standard Brownian motion and $\tau_{0}=\inf\\{t\geq 0\colon
W_{t}=0\\}$. Notice that the last expectation can be considered as the
expectation of a killed Brownian motion.

In Barlow et al. [2], the authors also gave the martingale characterization of
the WBM realized as a process on the plane. In terms of the $d$-dimensional
realization $X$ of the WBM, their characterization takes the following form.

######  Theorem 3.1 (Propositon 3.1 and Theorem 3.2 in Barlow et al. [2]).

Let $X=(X_{1}(t),\dots,X_{d}(t))_{t\geq 0}$ be an adapted continuous process.
Then $X$ is a WBM with parameters $p_{1},\dots,p_{d}>0$ if and only if it
satisfies the following conditions:

  1. 1.

$X_{i}(t)\geq 0$ and $X_{i}(t)X_{j}(t)=0$ for all $i\neq j$ and $t\geq 0$;

  2. 2.

for each $i=1,\dots,d$ the process

$\displaystyle N_{i}(t):=q_{i}X_{i}(t)-p_{i}\sum_{j\neq i}X_{j}(t),\quad t\geq
0,$ (3.3)

is a continuous martingale with respect to $\mathbb{F}$;

  3. 3.

for each $i=1,\dots,d$ the process

$\displaystyle(N_{i}(t))^{2}-\int_{0}^{t}\Big{(}q_{i}1\\!\\!\\!\;{\rm
I}(X_{i}(s)>0)-p_{i}1\\!\\!\\!\;{\rm
I}(X_{i}(s)=0)\Big{)}^{2}\,\mathrm{d}s,\quad t\geq 0,$ (3.4)

is a continuous martingale with respect to $\mathbb{F}$.

Note that since the product of the indicator functions in (3.4) is identically
zero, we have

$\displaystyle\langle N_{i}\rangle_{t}$
$\displaystyle=\int_{0}^{t}\Big{(}q_{i}1\\!\\!\\!\;{\rm
I}(X_{i}(s)>0)-p_{i}1\\!\\!\\!\;{\rm I}(X_{i}(s)=0)\Big{)}^{2}\,\mathrm{d}s$
(3.5) $\displaystyle=\int_{0}^{t}\Big{(}q_{i}^{2}1\\!\\!\\!\;{\rm
I}(X_{i}(s)>0)+p_{i}^{2}1\\!\\!\\!\;{\rm I}(X_{i}(s)=0)\Big{)}\,\mathrm{d}s.$

Furthermore, it is clear, see Lemma 2.2 in Barlow et al. [2], that the radial
process

$\displaystyle R(t):=\sum_{i=1}^{d}X_{i}(t)=\max_{1\leq i\leq d}X_{i}(t)$
(3.6)

is a reflecting Brownian motion, and hence it has a local time at $0$ defined
by:

$\displaystyle L^{X}(t):=L^{R}(t):=\lim_{\varepsilon\to
0+}\frac{1}{2\varepsilon}\int_{0}^{t}1\\!\\!\\!\;{\rm I}\Big{(}\max_{1\leq
i\leq d}X_{i}(s)\leq\varepsilon\Big{)}\,\mathrm{d}s.$ (3.7)

In the following theorem we give an equivalent martingale characterization of
the WBM that better fits into the setting of this paper.

######  Theorem 3.2.

Let $X=(X_{1}(t),\dots,X_{d}(t))_{t\geq 0}$ and $\nu=(\nu(t))_{t\geq 0}$ be
adapted continuous processes. Then $X$ is a WBM with parameters
$p_{1},\dots,p_{d}>0$, and $\nu$ is the local time of $X$ at $0$ if and only
if they satisfy the following conditions:

  1. a)

$X_{i}(t)\geq 0$ and $X_{i}(t)X_{j}(t)=0$ for all $i\neq j$, $t\geq 0$;

  2. b)

$\nu(0)=0$, $\nu$ is nondecreasing a.s., $\int_{0}^{\infty}1\\!\\!\\!\;{\rm
I}(X(s)\neq 0)\,\mathrm{d}\nu(s)=0$ a.s.;

  3. c)

the processes $M_{1},\dots,M_{d}$ defined by

$\displaystyle M_{i}(t):=X_{i}(t)-p_{i}\nu(t),\quad t\geq 0,$ (3.8)

are continuous square integrable martingales with respect to $\mathbb{F}$ with
the predictable quadratic variations

$\displaystyle\langle M_{i}\rangle_{t}=\int_{0}^{t}1\\!\\!\\!\;{\rm
I}(X_{i}(s)>0)\,\mathrm{d}s.$ (3.9)

  4. d)

$\int_{0}^{\infty}1\\!\\!\\!\;{\rm I}(X(s)=0)\,\mathrm{d}s=0$ a.s.

######  Remark 3.3.

There is a typo in §3 in Barlow et al. [2] (p. 281). One should remove the
indicator $1\\!\\!\\!\;{\rm I}_{r>0}$ in the formula for $h_{i}(r,\theta)$,
otherwise the zero process is also a solution of equations (3.2) or (3.3) from
their paper. Note that condition (d) in our Theorem 3.2 prohibits $0$ to be a
sticky point of $X$.

###### Proof.

We will show that conditions 1–3 of Theorem 3.1 and a)–d) of Theorem 3.2 are
equivalent.

i) We show that a)–d) $\Rightarrow$ 1–3.

First, it is obvious that a) implies 1.

By c), each process $N_{i}$ is a martingale as a linear combination of
martingales:

$\displaystyle N_{i}(t)$ $\displaystyle=q_{i}X_{i}(t)-p_{i}\sum_{j\neq
i}X_{j}(t)$ (3.10) $\displaystyle=q_{i}(M_{i}(t)+p_{i}\nu(t))-p_{i}\sum_{j\neq
i}(M_{j}(t)+p_{j}\nu(t))$ $\displaystyle=q_{i}M_{i}(t)-p_{i}\sum_{j\neq
i}M_{j}(t).$

Observe that quadratic covariations $\langle M_{i},M_{j}\rangle$ vanish for
$i\neq j$. Indeed, applying the Cauchy-type inequality for quadratic
covariations (Proposition 15.10 in Kallenberg [14]) we get

$\displaystyle\langle M_{i},M_{j}\rangle_{t}$
$\displaystyle=\int_{0}^{t}\mathrm{d}\langle
M_{i},M_{j}\rangle_{s}=\int_{0}^{t}\Big{(}1\\!\\!\\!\;{\rm
I}(X_{i}(s)>0)+1\\!\\!\\!\;{\rm I}(X_{i}(s)=0)\Big{)}\,\mathrm{d}\langle
M_{i},M_{j}\rangle_{s}$ (3.11)
$\displaystyle\leq\Big{(}\int_{0}^{t}1\\!\\!\\!\;{\rm
I}(X_{i}(s)>0)\,\mathrm{d}\langle
M_{j}\rangle_{s}\cdot\int_{0}^{t}\mathrm{d}\langle
M_{i}\rangle_{s}\Big{)}^{1/2}+\Big{(}\int_{0}^{t}1\\!\\!\\!\;{\rm
I}(X_{i}(s)=0)\,\mathrm{d}\langle
M_{i}\rangle_{s}\cdot\int_{0}^{t}\mathrm{d}\langle
M_{j}\rangle_{s}\Big{)}^{1/2}$
$\displaystyle\leq\Big{(}\int_{0}^{t}1\\!\\!\\!\;{\rm
I}(X_{i}(s)>0)1\\!\\!\\!\;{\rm I}(X_{j}(s)>0)\,\mathrm{d}s\cdot
t\Big{)}^{1/2}+\Big{(}\int_{0}^{t}1\\!\\!\\!\;{\rm
I}(X_{i}(s)=0)1\\!\\!\\!\;{\rm I}(X_{i}(s)>0)\,\mathrm{d}s\cdot
t\Big{)}^{1/2}$ $\displaystyle=0$

Consequently with the help of c) we obtain the martingale property 3. Indeed,
by the Itô formula we get

$\displaystyle(N_{i}(t))^{2}$
$\displaystyle=(N_{i}(0))^{2}+2\int_{0}^{t}N_{i}(s)\,\mathrm{d}N_{i}(s)+\langle
N_{i}\rangle_{t}$ (3.12)
$\displaystyle=(N_{i}(0))^{2}+2q_{i}\int_{0}^{t}N_{i}(s)\,\mathrm{d}M_{i}(s)-2p_{i}\sum_{j\neq
i}\int_{0}^{t}N_{i}(s)\,\mathrm{d}M_{j}(s)$ $\displaystyle+q_{i}^{2}\langle
M_{i}\rangle_{t}+p_{i}^{2}\sum_{j\neq i}\langle M_{j}\rangle_{t}.$

We have

$\displaystyle\int_{0}^{t}\Big{(}q_{i}1\\!\\!\\!\;{\rm
I}(X_{i}(s)>0)-p_{i}1\\!\\!\\!\;{\rm I}(X_{i}(s)=0)\Big{)}^{2}\,\mathrm{d}s$
$\displaystyle=\int_{0}^{t}\Big{(}q_{i}^{2}1\\!\\!\\!\;{\rm
I}(X_{i}(s)>0)+p_{i}^{2}1\\!\\!\\!\;{\rm I}(X_{i}(s)=0)\Big{)}\,\mathrm{d}s$
(3.13) $\displaystyle=q_{i}^{2}\langle M_{i}\rangle_{t}+p_{i}^{2}\sum_{j\neq
i}\langle M_{j}\rangle_{t},$

where we used (d) in the last equality. Hence we get the property 3.

Show that the process $\nu$ is the local time of $X$ at zero. Let $L^{X}$ be
the local time of $X$ ar zero and let

$\displaystyle R(t)=\max_{1\leq i\leq d}X_{i}(t)=\sum_{i=1}^{d}X_{i}(t)$
(3.14)

be the radial part of the WBM $X$. By Lemma 2.2 in Barlow et al. [2], $R$ is a
reflecting Brownian motion with a local time $L^{R}$, and hence

$\displaystyle R(t)-L^{R}(t)$ (3.15)

is a Brownian motion. On the other hand

$\displaystyle R(t)-\nu(t)=\sum_{i=1}^{d}M_{i}(t)$ (3.16)

and $Z(t):=\sum_{i=1}^{d}M_{i}(t)$ is a continuous martingale with the bracket
(we use 3.)

$\displaystyle\langle Z\rangle_{t}=\sum_{i=1}^{d}\int_{0}^{t}1\\!\\!\\!\;{\rm
I}(X_{i}(s)>0)\,\mathrm{d}s=t,$ (3.17)

thus, $Z$ is a $\mathbb{F}$-Brownian motion. By the uniqueness of the
semimartingale decomposition, $\nu=L^{R}$. Notice that $L^{R}=L^{X}$.

ii) We show that 1–3 $\Rightarrow$ a)–d).

Let $X$ be a WBM. Denote $\nu:=L^{X}$ its local time at 0. The properties a),
b) and d) follow immediately.

For each $i=1,\dots,d$ consider the process

$\displaystyle Y_{i}(t)=X_{i}(t)-\sum_{j\neq i}X_{j}(t),\quad t\geq 0.$ (3.18)

and a rescaled martingale

$\displaystyle\tilde{N}_{i}(t):=\frac{1}{p_{i}q_{i}}N_{i}(t)=\frac{1}{p_{i}}X_{i}(t)-\frac{1}{q_{i}}\sum_{j\neq
i}X_{j}(t),\quad t\geq 0,$ (3.19)

with the bracket

$\displaystyle\langle\tilde{N}_{i}\rangle_{t}=\int_{0}^{t}\Big{(}\frac{1}{p_{i}^{2}}1\\!\\!\\!\;{\rm
I}(X_{i}(s)>0)+\frac{1}{q_{i}^{2}}1\\!\\!\\!\;{\rm
I}(X_{i}(s)=0)\Big{)}\,\mathrm{d}s.$ (3.20)

Then obviously

$\displaystyle 1\\!\\!\\!\;{\rm I}(N_{i}(s)>0)$
$\displaystyle=1\\!\\!\\!\;{\rm I}(\tilde{N}_{i}(s)>0)=1\\!\\!\\!\;{\rm
I}(X_{i}(s)>0)=1\\!\\!\\!\;{\rm I}(Y_{i}(s)>0),$ (3.21) $\displaystyle
1\\!\\!\\!\;{\rm I}(N_{i}(s)<0)$ $\displaystyle=1\\!\\!\\!\;{\rm
I}(\tilde{N}_{i}(s)<0)=1\\!\\!\\!\;{\rm I}(Y_{i}(s)<0).$

Consider the processes

$\displaystyle W_{i}(t)$
$\displaystyle:=\int_{0}^{t}\Big{(}p_{i}1\\!\\!\\!\;{\rm
I}(\tilde{N}_{i}(s)>0)+q_{i}1\\!\\!\\!\;{\rm
I}(\tilde{N}_{i}(s)<0)\Big{)}\,\mathrm{d}\tilde{N}_{i}(s).$ (3.22)

These process are continuous martingales with the bracket

$\displaystyle\langle W_{i}\rangle_{t}=t,$ (3.23)

hence they are Brownian motions. Denoting
$\sigma_{i}(x):=\frac{1}{p_{i}}1\\!\\!\\!\;{\rm
I}(x>0)+\frac{1}{q_{i}}1\\!\\!\\!\;{\rm I}(x<0)$ we get that $\tilde{N}_{i}$
satisfies the SDE

$\displaystyle\tilde{N}_{i}(t)=\tilde{N}_{0}+\int_{0}^{t}\sigma_{i}(\tilde{N}_{i}(s))\,\mathrm{d}W_{i}(s),$
(3.24)

and hence each $N_{i}$ is the so-called oscillating Brownian motion, see
Keilson and Wellner [16]. By Nakao’s theorem, see Nakao [22], this SDE has a
unique strong solution. Let $r_{i}(x)=p_{i}1\\!\\!\\!\;{\rm
I}(x>0)+q_{i}1\\!\\!\\!\;{\rm I}(x<0)$ so that

$\displaystyle Y_{i}=r_{i}(\tilde{N}_{i}).$ (3.25)

Repeating literally the calculations from Section 5.2 in Lejay [17] we get
that

$\displaystyle Y_{i}(t)=Y_{i}(0)+W_{i}(t)+(2p_{i}-1)L^{Y_{i}}(t),$ (3.26)

where $L^{Y_{i}}(t)$ is the symmetric local time of ${Y_{i}}$ at 0.

Hence $Y_{i}$ is a skew Brownian motion with parameter $(2p_{i}-1)$. Moreover
the radial process $R=|Y_{i}|=|X|$ is a reflected Brownian motion and

$\displaystyle
R(t)=|Y_{i}(t)|=R(0)+\int_{0}^{t}\operatorname{sgn}(Y_{i})\,\mathrm{d}W_{i}(s)+L^{Y_{i}}(t).$
(3.27)

Futhermore,

$\displaystyle R(t)=R(0)+W(t)+L^{R}(t)$ (3.28)

for some Brownian motion $W$. From the uniqueness of the decomposition of $R$
as a semimartingale we get that $L^{Y_{i}}=L^{R}$ for all $i=1,\dots,d$.

It follows from [34, Theorem 1.7, Chapter VI] that the right and left local
times of $Y_{i}$ at 0 are equal to

$\displaystyle L^{Y_{i}}_{\text{left}}(t)=2(1-p_{i})L^{Y_{i}}(t),\ \
L^{Y_{i}}_{\text{right}}(t)=2p_{i}L^{Y_{i}}(t).$ (3.29)

Finally Tanaka’s formula yields

$\displaystyle X_{i}(t)$ $\displaystyle=\max\\{Y_{i}(t),0\\}=$ (3.30)
$\displaystyle=\max\\{Y_{i}(0),0\\}+\int_{0}^{t}1\\!\\!\\!\;{\rm
I}(Y_{i}(s)>0)\,\mathrm{d}W_{i}(s)+\frac{1}{2}L^{Y_{i}}_{\text{right}}(t)$
$\displaystyle=\max\\{Y_{i}(0),0\\}+\int_{0}^{t}1\\!\\!\\!\;{\rm
I}(Y_{i}(s)>0)\,\mathrm{d}W_{i}(s)+p_{i}L^{R}(t),$

and thus the process $M_{i}=X_{i}-p_{i}\nu$ is a continuous martingale with
the bracket (3.9). ∎

######  Corollary 3.4.

Let $X=(X_{1}(t),\dots,X_{d}(t))_{t\geq 0}$ be a WBM with parameters
$p_{1},\dots,p_{d}>0$ and let $I\subseteq\\{1,\dots,d\\}$. Let

$\displaystyle\gamma=\sum_{i\in I}p_{i}-\sum_{j\in I^{c}}p_{j}=2\sum_{i\in
I}p_{i}-1\in[-1,1].$ (3.31)

Then the process

$\displaystyle X^{\gamma}(t):=\sum_{i\in I}X_{i}(t)-\sum_{j\in
I^{c}}X_{j}(t),\quad t\geq 0,$ (3.32)

is a skew Brownian motion with the parameter $\gamma$.

###### Proof.

Without loss of generality assume that $I=\\{1,\dots,k\\}$ for some $0\leq
k\leq d$. By Theorem 3.2, the process

$\displaystyle W(t)$
$\displaystyle=\sum_{i=1}^{k}M_{i}(t)-\sum_{j=k+1}^{d}M_{j}(t)$ (3.33)

is a continuous martingale. Taking into account (3.8) we get

$\displaystyle W(t)$ $\displaystyle=X^{\gamma}(t)-\gamma\nu(t).$ (3.34)

Since the local times at $0$ of $X$ and $X^{\gamma}$ coincide, we get

$\displaystyle W(t)=X^{\gamma}(t)-\gamma L(t)$ (3.35)

where $L$ is the local time of $X^{\gamma}$ at $0$. The bracket of the
martingale $W$ equals

$\displaystyle\langle W\rangle_{t}=\sum_{i=1}^{d}\int_{0}^{t}1\\!\\!\\!\;{\rm
I}(X_{i}(s)>0)\,\mathrm{d}s=t,$ (3.36)

and thus $W$ is a Brownian motion. In other words, $X^{\gamma}$ satisfies the
SDE $X^{\gamma}(t)=W(t)+\gamma L(t)$ and thus is a skew Brownian motion, see
Harrison and Shepp [10]. ∎

##  4 Proof of Theorem 2.1

Consider the Markov chain $\mathcal{Z}=(\mathcal{X},Y)$ on the enlarged state
space $\\{\pm 0,\pm 1,\dots\\}\times\mathbb{Z}^{m-1}$. Let us decompose the
process $\mathcal{X}$ into $2|U|$ non-negative excursions parameterized by the
elements of the set $\\{-0,+0\\}\times U$. By $\sigma_{n}$, $n\geq 0$, denote
the time instant of the last visit of the random walk $\mathcal{Z}$ to $H_{0}$
before time $n$, i.e.,

$\displaystyle\sigma_{n}:=\max\Big{\\{}k\leq
n\colon\mathcal{X}(k)\in\\{-0,+0\\}\Big{\\}},\quad\sigma_{0}=0.$ (4.1)

Consider $2|U|$ processes

$\displaystyle X^{j,+}(n)$ $\displaystyle:=X(n)1\\!\\!\\!\;{\rm
I}(Y(\sigma_{n})\equiv j,X(n)>0),$ (4.2) $\displaystyle X^{j,-}(n)$
$\displaystyle:=|X(n)|1\\!\\!\\!\;{\rm I}(Y(\sigma_{n})\equiv j,X(n)<0),\quad
j\in U,$

and introduce a $2|U|$-dimensional Markov chain
$\\{(X^{j,+}(n),X^{j,-}(n))_{j\in U}\\}_{n\geq 0}$ on $\mathbb{N}_{0}^{2|U|}$.
Observe that all its coordinates are non-negative, only one coordinate of the
vector $(X^{j,+}(n),X^{j,-}(n))_{j\in U}$ may be non-zero, and
$X^{j,\pm}(\sigma_{n})=0$.

To study the limit behavior of the scaled process $X_{n}(\cdot)$ defined in
(2.19) we will prove that the properly scaled $2|U|$-dimensional process
$(X^{j,+}(n),X^{j,-}(n))_{j\in U}$ converges to a WBM with parameters
$\\{\pi_{\\{+0\\}\times j},\pi_{\\{-0\\}\times j}\\}_{j\in U}$. To this end,
for each $j\in U$ we decompose the processes $X^{j,+}$ and $X^{j,-}$ into a
sum of a martingale and a “local time in 0”, namely we set

$\displaystyle X^{j,\pm}(n)$
$\displaystyle=\sum_{k=1}^{n}(X^{j,\pm}(k)-X^{j,\pm}(k-1))1\\!\\!\\!\;{\rm
I}(X^{j,\pm}(k-1)>0)$ (4.3) $\displaystyle+\sum_{k=1}^{n}1\\!\\!\\!\;{\rm
I}(X^{j,\pm}(k-1)=0,X^{j,\pm}(k)=1)$
$\displaystyle=:M^{j,\pm}(n)+L^{j,\pm}(n),\qquad M^{j,\pm}(0)=L^{j,\pm}(0)=0.$

Since the original process $Z$ is a symmetric random walk outside of $H$, it
follows from the construction that

$\displaystyle\mathbf{E}\Big{[}(X^{j,\pm}(k)$
$\displaystyle-X^{j,\pm}(k-1))1\\!\\!\\!\;{\rm
I}(X^{j,\pm}(k-1)>0)\Big{|}\mathscr{F}_{k-1}\Big{]}$ (4.4)
$\displaystyle=1\\!\\!\\!\;{\rm
I}(X^{j,\pm}(k-1)>0)\cdot\mathbf{E}\Big{[}X^{j,\pm}(k)-X^{j,\pm}(k-1)\Big{|}\mathscr{F}_{k-1}\Big{]}=0$

and the sequences $(M^{j,+}(n))_{n\geq 0}$ and $(M^{j,-}(n))_{n\geq 0}$ are
martingales for any $j\in U$ with respect to filtration
$(\mathscr{F}_{n})_{n\geq 0}$ generated by $Z$. Moreover, since

$\displaystyle\langle M^{j,\pm}\rangle_{n}$
$\displaystyle=\sum_{k=1}^{n}\mathbf{E}\Big{[}(X^{j,\pm}(k)-X^{j,\pm}(k-1))^{2}\cdot
1\\!\\!\\!\;{\rm I}(X^{j,\pm}(k-1)>0)\Big{|}\mathscr{F}_{k-1}\Big{]}$ (4.5)
$\displaystyle=\frac{1}{m}\sum_{k=1}^{n}1\\!\\!\\!\;{\rm
I}(X^{j,\pm}(k-1)>0),$

the sequences

$\displaystyle(M^{j,\pm}(n))^{2}-\frac{1}{m}\sum_{k=1}^{n}1\\!\\!\\!\;{\rm
I}(X^{j,\pm}(k-1)>0),\quad n\geq 0,$ (4.6)

are martingales too.

Set

$\displaystyle X_{n}^{j,\pm}(t):=\frac{X^{j,\pm}([nt])}{\sqrt{n}},\quad
M_{n}^{j,\pm}(t):=\frac{M^{j,\pm}([nt])}{\sqrt{n}},\quad\nu_{n}^{j,\pm}(t):=\frac{L^{j,\pm}([nt])}{\sqrt{n}}.$
(4.7)

######  Proposition 4.1.

There are continuous processes $\\{X_{j}^{\pm}\\}_{j\in U}$, continuous
martingales $\\{M_{j}^{\pm}\\}_{j\in U}$ and a nondecreasing process $\nu$
such that

$\displaystyle(X_{n}^{j,\pm},M_{n}^{j,\pm},\nu_{n}^{j,\pm})_{j\in
U}\Rightarrow(X_{j}^{\pm},M_{j}^{\pm},\pi^{\pm}_{j}\nu)_{j\in U},\quad
n\to\infty.$ (4.8)

Moreover, the process
$\sqrt{m}(X_{j}^{\pm},M_{j}^{\pm},\pi^{\pm}_{j}\nu)_{j\in U}$ satisfies
conditions of Theorem 3.2.

To prepare the proof of Proposition 4.1 we notice that in problems of this
type it is often helpful to start with the study of the radial process

$\displaystyle S(n):=\sum_{j\in U}(X^{j,-}(n)+X^{j,+}(n)),\ n\geq 0.$ (4.9)

The process $S=\\{S(n)\\}_{n\geq 0}$ is a Markov chain on
$\mathbb{N}\cup\\{0\\}$ with reflection at $0$ whose steps in $\mathbb{N}$
have the distribution

$\displaystyle\mathbf{P}\Big{(}S(n+1)=j-1\Big{|}S(n)=j\Big{)}=\mathbf{P}\Big{(}S(n+1)=j+1\Big{|}S(n)=j\Big{)}=\frac{1}{2m},$
(4.10)
$\displaystyle\mathbf{P}\Big{(}S(n+1)=j\Big{|}S(n)=j\Big{)}=1-\frac{1}{m},\quad
j\geq 1,$ $\displaystyle\mathbf{P}\Big{(}S(n+1)=1\Big{|}S(n)=0\Big{)}=1.$

Set

$L(n):=\sum_{k=0}^{n-1}1\\!\\!\\!\;{\rm
I}(X(k)=0)=\sum_{k=0}^{n-1}1\\!\\!\\!\;{\rm I}(S(k)=0)=\sum_{j\in
U}(L^{+}_{j}(n)+L^{-}_{j}(n)),\quad n\geq 0,$ (4.11)

$L(n)$ being the number of visits of $S$ to the origin up to time $n$.

######  Lemma 4.2.

We have the weak convergence

$\displaystyle\Big{(}\frac{S([n\cdot])}{\sqrt{n}},\frac{L([n\cdot])}{\sqrt{n}}\Big{)}\Rightarrow\frac{1}{\sqrt{m}}\Big{(}|B(\cdot)|,L^{0}(\cdot)\Big{)},\
n\to\infty,$ (4.12)

where $B$ is a standard Brownian motion, $L^{0}(\cdot)$ is a local time of $B$
at 0 defined by (3.7).

###### Proof.

Whereas the convergence of the reflected random walks to a reflected Brownian
motion is straightforward, certain work should be done to ensure the
convergence of the local times. A similar result in a more general setting can
be found in Section 2 of Pilipenko and Prykhodko [32] so that here we just
briefly outline the argument.

First we construct an auxiliary slowed-down one-dimensional symmetric random
walk $Q$ on $\mathbb{Z}$ with steps $\pm 1$ and $0$, namely we set

$\displaystyle\mathbf{P}\Big{(}Q(n+1)$
$\displaystyle=j-1\Big{|}Q(n)=j\Big{)}=\mathbf{P}\Big{(}Q(n+1)=j+1\Big{|}Q(n)=j\Big{)}=\frac{1}{2m},$
(4.13) $\displaystyle\mathbf{P}\Big{(}Q(n+1)$
$\displaystyle=j\Big{|}Q(n)=j\Big{)}=1-\frac{1}{m},\quad j\in\mathbb{Z}.$

Clearly, by the functional central limit theorem,

$\displaystyle Q([n\cdot])/\sqrt{n}\Rightarrow B(\cdot)/\sqrt{m},\quad
n\to\infty.$ (4.14)

Then we construct a copy of the reflected random walk $S$ on the basis of the
random walk $Q$. To this end, we write the solution of the Skorokhod
reflection problem for $Q$,

$\displaystyle\tilde{Q}(n)=Q(n)+L^{Q}(n),$ (4.15)

where

$\displaystyle L^{Q}(n)=-\min_{0\leq k\leq n}(Q(k)\wedge 0).$ (4.16)

Then $\tilde{Q}(n)\geq 0$. Since $\tilde{Q}$ may spend several consecutive
steps at the origin, we perform a random time transformation with the help of
the local time $L^{Q}$ to obtain a copy of a Markov chain $S$. Essentially we
push the the process $\tilde{Q}$ from zero with probability one. Since the
number of visits of zero on the time interval $[0,n]$ is of the order
$\sqrt{n}$ this random time transformation (properly rescaled) converges to
identity. Since the reflection mapping is continuous, from the convergence
(4.14) we obtain convergence of the rescaled constructed copy of $S$ and the
local time $L^{Q}$ to the solution of the Skorokhod reflection problem for
$B/\sqrt{m}$, namely to $(\max_{s\leq t}B_{s}-B_{t},\max_{s\leq
t}B_{s})/\sqrt{m}$, which has the same law as
$(|B(\cdot)|,L^{0}(\cdot))/\sqrt{m}$ by Lévy’s theorem, see, e.g. [14,
Corollary 19.3]. ∎

######  Lemma 4.3.

The sequence
$\\{(X^{j,\pm}_{n}(\cdot),M^{j,\pm}_{n}(\cdot),\nu^{j,\pm}_{n}(\cdot))_{j\in
U}\\}_{n\geq 1}$ is weakly relatively compact in $D([0,T],\mathbb{R}^{6|U|})$
and any limit point is a continuous processes.

###### Proof.

Notice that the modulus of continuity (in the uniform topology) of any
$X^{j,\pm}_{n}$ is dominated by doubled modulus of continuity of
$\\{S(n\cdot)/\sqrt{n}\\}$, and the modulus of continuity of any
$M^{j,\pm}_{n}$ is dominated by doubled modulus of continuity of
$X^{j,\pm}_{n}$. The third coordinate is the difference of the first two.
Hence the proof of the Lemma follows from Lemma 4.2. ∎

_Proof of Proposition 4.1_. To show convergence $\\{(X_{n}^{j,\pm})_{j\in
U}\\}_{n\geq 0}$ to the WBM it suffices to verify that for any subsequence
$\\{(X_{n_{k}}^{j,\pm})_{j\in U}\\}_{k\geq 0}$ there is a subsubsequence that
converges to the WBM. Due to Lemma 4.3 without loss of generality we will
assume that the sequence
$\Big{\\{}\Big{(}X^{j,\pm}_{n}(\cdot),M^{j,\pm}_{n}(\cdot),\nu^{j,\pm}_{n}(\cdot)\Big{)}_{j\in
U}\Big{\\}}_{n\geq 1}$ converges in distribution to a continuous process
$\Big{(}X^{j,\pm}(\cdot),M^{j,\pm}(\cdot),\nu^{j,\pm}(\cdot)\Big{)}_{j\in U}$.

Let us check the conditions a)–d) of Theorem 3.2 for the process
$\sqrt{m}\Big{(}X^{j,\pm}(\cdot),M^{j,\pm}(\cdot),\nu^{j,\pm}(\cdot)\Big{)}_{j\in
U}$.

a) It follows from the construction that $X^{j,\pm}(t)\geq 0$, $t\in[0,T]$.
Moreover, only one of these processes may be positive at any fixed time.

b) and d) The processes $\nu^{j,\pm}(\cdot)$ are non-decreasing a.s. and
$\nu^{j,\pm}(0)=0$. Lemma 4.2 yields that

$\displaystyle\Big{(}|X(\cdot)|,\nu(\cdot)\Big{)}\stackrel{{\scriptstyle\mathrm{d}}}{{=}}\frac{1}{\sqrt{m}}\Big{(}|B(\cdot)|,L^{0}(\cdot)\Big{)},$
(4.17)

where

$\displaystyle\nu(t)=\sum_{j\in U}(\nu^{j,-}(t)+\nu^{j,+}(t)).$ (4.18)

Since $\int_{0}^{T}1\\!\\!\\!\;{\rm I}(B(t)=0)\,\mathrm{d}t=0$ and
$\int_{0}^{T}1\\!\\!\\!\;{\rm I}(|B(t)|>0)\,\mathrm{d}L^{0}(t)=0$ a.s. for any
$T>0$, we have

$\displaystyle\int_{0}^{T}1\\!\\!\\!\;{\rm I}(|X(t)|=0)\,\mathrm{d}t=0,$
(4.19) $\displaystyle\int_{0}^{T}1\\!\\!\\!\;{\rm
I}(|X(t)|>0)\,\mathrm{d}\nu(t)=0$

almost surely.

c) It follows from the construction that

$\displaystyle X^{j,\pm}(t)=M^{j,\pm}(t)+\nu^{j,\pm}(t)\ \text{a.s.}$ (4.20)

for all $j\in U,t\geq 0$. To show that

$\displaystyle\nu^{j,\pm}(t)=\pi^{\pm}_{j}\nu(t)\ \text{a.s.}$ (4.21)

we recall the strong law of large numbers (2.13) for Markov chains. For any
$t>0$ the process $L([nt])$ defined in (4.11) increases to $+\infty$ a.s. as
$n\to\infty$. Hence

$\displaystyle\frac{\nu_{n}^{j,\pm}(t)}{\sum_{k\in
U}\nu_{n}^{k,\pm}(t)}=\frac{1}{L([nt])}\sum_{k=0}^{L([nt])}1\\!\\!\\!\;{\rm
I}(\hat{\mathcal{Y}}(k)=\\{\pm 0\\}\times j)\to\pi^{\pm}_{j}\ \text{a.s.}$
(4.22)

where $\hat{\mathcal{Y}}$ is defined in (2.11).

The processes $M_{n}^{j,\pm}(\cdot)$, $j\in U$, are local martingales with
respect to the filtration generated by
$\\{X_{n}^{j,\pm}(\cdot),M_{n}^{j,\pm}(\cdot)\\}_{j\in U}$. Since the jumps of
each $M_{n}^{j,\pm}(\cdot)$ are uniformly bounded, the limits $M^{j,\pm}$ are
local martingales with respect to filtration generated by
$\\{X^{j,\pm},M^{j,\pm}\\}_{k\in U}$ due to Lemma 1.17 in Chapter IX of Jacod
and Shiryaev [13]. Moreover, the limit processes are continuous due to Lemma
4.3.

It is left to show that

$\displaystyle\langle
M^{j,\pm}\rangle_{t}=\frac{1}{m}\int_{0}^{t}1\\!\\!\\!\;{\rm
I}(X^{j,\pm}(s)>0)\,\mathrm{d}s\ \text{a.s.\ for }j\in U,\ t\geq 0.$ (4.23)

By Skorokhod’s representation theorem there is a probability space and the
copies

$\displaystyle\Big{\\{}\Big{(}\tilde{X}^{j,\pm}_{n}(\cdot),\tilde{M}^{j,\pm}_{n}(\cdot),\tilde{\nu}^{j,\pm}_{n}(\cdot)\Big{)}_{j\in
U}\Big{\\}}_{n\geq 1}\text{ and
}\Big{\\{}\Big{(}\tilde{X}^{j,\pm}(\cdot),\tilde{M}^{j,\pm}(\cdot),\tilde{\nu}^{j,\pm}(\cdot)\Big{)}_{j\in
U}\Big{\\}}$ (4.24)

of

$\displaystyle\Big{\\{}\Big{(}X^{j,\pm}_{n}(\cdot),M^{j,\pm}_{n}(\cdot),\nu^{j,\pm}_{n}(\cdot)\Big{)}_{j\in
U}\Big{\\}}_{n\geq 1}\text{ and
}\Big{\\{}\Big{(}X^{j,\pm}(\cdot),M^{j,\pm}(\cdot),\nu^{j,\pm}(\cdot)\Big{)}_{j\in
U}\Big{\\}}$ (4.25)

such that on any interval $[0,T]$ we have a.s. uniform convergence

$\Big{(}\tilde{X}^{j,\pm}_{n}(\cdot),\tilde{M}^{j,\pm}_{n}(\cdot),\tilde{\nu}^{j,\pm}_{n}(\cdot)\Big{)}_{j\in
U}\to\Big{(}\tilde{X}^{j,\pm}(\cdot),\tilde{M}^{j,\pm}(\cdot),\tilde{\nu}^{j,\pm}(\cdot)\Big{)}_{j\in
U},\quad n\to\infty.$ (4.26)

To prove (4.23) it suffices to verify that with probability 1 the sequence

$\displaystyle\Big{(}\tilde{M}^{j,\pm}_{n}(t)\Big{)}^{2}-\frac{1}{m}\int_{0}^{[nt]/n}1\\!\\!\\!\;{\rm
I}(\tilde{X}^{j,\pm}_{n}(s)>0)\,\mathrm{d}s$ (4.27)

converges uniformly over $t\in[0,T]$ to

$\displaystyle(\tilde{M}^{j,\pm}(t))^{2}-\frac{1}{m}\int_{0}^{t}1\\!\\!\\!\;{\rm
I}(\tilde{X}^{j,\pm}(s)>0)\,\mathrm{d}s.$ (4.28)

Here we again use Lemma 1.17 in Chapter IX of Jacod and Shiryaev [13] and a
localization procedure.

It follows from (4.26) that we have to prove the convergence of the integrals
only.

Let $\omega\in\Omega$ be such that (4.26) holds. If $s\in[0,T]$ is such that
$\tilde{X}^{k,\mathfrak{s}}(s)>0$ for some $1\leq k\leq|U|$ and
$\mathfrak{s}\in\\{-,+\\}$ then $\tilde{X}^{k,\mathfrak{s}}_{n}(s)>0$ for
large $n$. Since only one of the processes
$\\{\tilde{X}^{j,\mathfrak{s}}(s)\\}_{j\in U,\mathfrak{s}\in\\{-,+\\}}$ and
only one of $\\{\tilde{X}^{j,\mathfrak{s}}_{n}(s)\\}_{j\in
U,\mathfrak{s}\in\\{-,+\\}}$ may be non-zero we have convergence of the
indicators for all $j\in U$:

$\lim_{n\to\infty}1\\!\\!\\!\;{\rm
I}(\tilde{X}^{j,\pm}_{n}(s)>0)=1\\!\\!\\!\;{\rm I}(\tilde{X}^{j,\pm}(s)>0).$
(4.29)

The process $\tilde{X}=(\tilde{X}^{j,+},\tilde{X}^{j,-})_{j\in U}$ spends zero
time in 0 with probability 1 because
$\sqrt{m}\sum_{j=1}^{|U|}(\tilde{X}^{j,+}+\tilde{X}^{j,-})$ is a reflected
Brownian motion, see Lemma 4.2. Therefore for a.a. $\omega$ and a.a.
$s\in[0,T]$ there is $k$ and $\mathfrak{s}$ such that
$\tilde{X}^{k,\mathfrak{s}}(s)>0$ and we have (4.29) for any index $(j,\pm)$,
$j\in U$. So by the Fubini theorem and by the Lebesgue dominated convergence
theorem for a.a. $\omega$ and all $j\in U$ we have convergence of the
integrals

$\displaystyle\lim_{n\to\infty}\int_{0}^{t}1\\!\\!\\!\;{\rm
I}(\tilde{X}^{j,\pm}_{n}(s)>0)\,\mathrm{d}s=\int_{0}^{t}1\\!\\!\\!\;{\rm
I}(\tilde{X}^{j,\pm}(s)>0)\,\mathrm{d}s.$ (4.30)

This completes the proof of Proposition 4.1. $\Box$

To treat convergence of $Y$, similarly to the representation (4.3) for
$\\{X^{j,\pm}(n)\\}$ we decompose the sequence $Y=\\{Y(n)\\}_{n\geq 0}$ into
the sum

$\displaystyle Y(n)$
$\displaystyle=\sum_{k=0}^{n-1}\Big{(}Y(k+1)-Y(k)\Big{)}1\\!\\!\\!\;{\rm
I}(X(k)\neq 0)$ (4.31) $\displaystyle+\sum_{j\in
U}\sum_{k=0}^{n-1}\Big{(}Y(k+1)-Y(k)\Big{)}1\\!\\!\\!\;{\rm
I}(\mathcal{X}(k)=\pm 0,Y(k)\equiv j)$ $\displaystyle=M^{Y}(n)+D^{Y}(n),\quad
n\geq 0,$

and we define

$\displaystyle M^{Y}_{n}(t):=\frac{M^{Y}([nt])}{\sqrt{n}},\quad
D^{Y}_{n}(t):=\frac{D^{Y}([nt])}{\sqrt{n}}.$ (4.32)

The following Proposition is proven analogously to the previous reasoning.

######  Proposition 4.4.

Let $X$ and $\nu$ be as in Proposition 4.1. Then

$\displaystyle\Big{(}X_{n}(\cdot),\nu_{n}(\cdot),M^{Y}_{n}(\cdot),D^{Y}_{n}(\cdot)\Big{)}\Rightarrow\Big{(}X(\cdot),\nu(\cdot),\frac{1}{\sqrt{m}}W_{Y}(\cdot),c\nu(\cdot)\Big{)},$
(4.33)

where $W_{Y}$ is a $(m-1)$-dimensional Brownian motion independent of $X$, and
$c$ is defined in (2.15).

_Proof of Theorem 2.1._ We combine Propositions 4.1 and 4.4 together with
Corollary 3.4. $\Box$

##  5 One-sided membrane with ergodic properties

The same method of decomposition of the perturbed Markov chain into a sum of
excursions combined with the strong law of large numbers (2.13) can be applied
for the analysis of a one-sided membrane that has ergodic properties.

As in Section 2, let $m\geq 2$ and let
$\\{\mathrm{e}_{1},\dots,\mathrm{e}_{m}\\}$ be a standard basis in
$\mathbb{R}^{m}$. Consider a Markov chain $Z=(X,Y)$ on $\mathbb{Z}^{m}$ that
behaves as a simple random walk outside of the hyperplane
$H:=\\{0\\}\times\mathbb{Z}^{m-1}$, i.e., (2.1) holds true.

Let $\\{p_{y}\\}_{y\in H}\subset[0,1]$. Now we interpret $H$ as a
semipermeable non-homogeneous membrane that may let a particle into one half-
space with probabilities $\\{p_{y}\\}_{y\in H}$. More precisely, we assume
that for each $z=(0,y)\in H$:

$\displaystyle
p_{y}=\mathbf{P}\Big{(}Z(n+1)=z+\mathrm{e}_{1}\Big{|}Z(n)=z\Big{)}=1-\mathbf{P}\Big{(}Z(n+1)=z-\mathrm{e}_{1}\,\Big{|}\,Z(n)=z\Big{)}.$
(5.1)

Note that the particle leaves the membrane in the direction orthogonal to $H$,
i.e., $Y(n+1)=Y(n)$ for $Z(n)\in H$ with probability 1, and hence there is no
slide along the membrane.

We assume that the membrane has the following ergodic property:

A${}_{\text{SLLN}}(\beta)$: there is $\beta>0$ such that

$\displaystyle\lim_{A\to\infty}\frac{1}{|V(A,o)|}\sum_{y\in
V(A,o)}p_{y}=:\bar{p}\in[0,1],$ (5.2)

where the limit is taken over all cubes $V(A,o)\subseteq H$ of volume
$|V(A,o)|$ with side size larger than $A$ and whose centre $o$ is within
distance $A^{\beta}$ from the origin.

We give two clarifying examples for the assumption A${}_{\text{SLLN}}(\beta)$.

######  Example 5.1.

Assume that the family $\\{p_{y}\\}_{y\in H}$ has a periodic structure: there
are $k_{2},\dots,k_{m}\geq 1$ such that for all
$l_{2},\dots,l_{m}\in\mathbb{Z}$ and for all $y\in H$

$\displaystyle
p_{y}=p_{y+k_{2}l_{2}\mathrm{e}_{2}+\cdots+k_{m}l_{m}\mathrm{e}_{m}}.$ (5.3)

Then $\\{p_{y}\\}$ clearly satisfy assumption A${}_{\text{SLLN}}(\beta)$ for
any $\beta>0$ with

$\displaystyle\bar{p}=\frac{1}{k_{2}\cdots
k_{m}}\sum_{i_{2},\dots,i_{m}=1}^{k_{2},\dots,k_{m}}p_{(i_{2},\dots,i_{m})}.$
(5.4)

######  Example 5.2.

Let $\\{p_{y}\\}_{y\in H}$ be i.i.d. random variables with values in $[0,1]$
defined on a probability space
$(\Omega^{\prime},\mathscr{F}^{\prime},\mathbf{P}^{\prime})$. Then for each
fixed $\omega^{\prime}\in\Omega^{\prime}$ the family
$\\{p_{y}(\omega^{\prime})\\}$ defines a random “environment”. Then the
assumption A${}_{\text{SLLN}}(\beta)$ with $\beta>0$ is satified with

$\displaystyle\bar{p}=\mathbf{E}^{\prime}p_{y}.$ (5.5)

To see this, let $\beta>0$, and let $V(A,o)$ denote a cube with the size
$A\in\mathbb{N}$ and the centre at $o\in\mathbb{Z}^{m-1}$. Then

$\displaystyle\mathbf{P}\Big{(}\frac{1}{|V(A,o)|}\sum_{y\in
V(A,o)}p_{y}\not\to\bar{p}\Big{)}=\mathbf{P}\Big{(}\bigcup_{m=1}^{\infty}\bigcap_{k=1}^{\infty}\bigcup_{A=k}^{\infty}\bigcup_{|o|\leq
L^{\beta}}\Big{\\{}\frac{1}{|V(A,o)|}\Big{|}\sum_{y\in
V(A,o)}(p_{y}-\bar{p})\Big{|}>\frac{1}{m}\Big{\\}}\Big{)}.$ (5.6)

For each $m\geq 1$ and $A\geq 1$ we apply Hoeffding’s inequality, see Chapter
III, §5.8 in Petrov [26]:

$\displaystyle\mathbf{P}\Big{(}\bigcup_{|o|\leq
A^{\beta}}\Big{\\{}\frac{1}{|V(A,o)|}$ $\displaystyle\Big{|}\sum_{y\in
V(A,o)}(p_{y}-\bar{p})\Big{|}>\frac{1}{m}\Big{\\}}\Big{)}\leq(2A^{\beta})^{m-1}\mathbf{P}\Big{(}\frac{1}{|V(A,o)|}\Big{|}\sum_{y\in
V(A,o)}(p_{y}-\bar{p})\Big{|}>\frac{1}{m}\Big{)}$ (5.7) $\displaystyle\leq
2(2A^{\beta})^{m-1}\mathrm{e}^{-2|V(A,o)|/m^{2}}=2(2A^{\beta})^{m-1}\mathrm{e}^{-2A^{m-1}/m^{2}}.$

Hence, the probability in (5.6) equals to 0.

######  Theorem 5.3.

Let assumption _A ${}_{\text{SLLN}}(\beta)$_ holds true for $\beta>1$. Then
for any initial value $Z(0)\in\mathbb{Z}^{m}$ the weak convergence holds true:

$\displaystyle(X_{n},Y_{n})\Rightarrow\frac{1}{\sqrt{m}}(X^{\gamma},Y^{0}),\quad
n\to\infty,$ (5.8)

where $X^{\gamma}$ and $Y^{0}$ are defined in (2.17) and (2.18), with
$\gamma=2\bar{p}-1$, and $c=0$.

###### Proof.

1\. Without loss of generality assume that $Z(0)=0$. As in Section 4 we
decompose the Markov chain $Z=(X,Y)$ into the “left” and the “right”
excursions. Since the membrane is one-sided, and there is no need of
introducing the set $U$, the notation of Section 4 simplifies significantly.
Similarly to (4.2) we define

$\displaystyle X^{+}(n)$ $\displaystyle:=X(n)\cdot 1\\!\\!\\!\;{\rm
I}(X(n)>0),$ (5.9) $\displaystyle X^{-}(n)$ $\displaystyle:=|X(n)|\cdot
1\\!\\!\\!\;{\rm I}(X(n)<0),$

so that

$\displaystyle X(n)=X^{+}(n)-X^{-}(n),\quad n\geq 0.$ (5.10)

Then we decompose these processes similarly to (4.3) as

$\displaystyle X^{\pm}(n)$
$\displaystyle=\sum_{k=1}^{n}\Big{(}X^{\pm}(k)-X^{\pm}(k-1)\Big{)}\cdot
1\\!\\!\\!\;{\rm I}(X^{\pm}(k-1)>0)+\sum_{k=1}^{n}1\\!\\!\\!\;{\rm
I}(X^{\pm}(k-1)=0,X^{\pm}(k)=1)$ (5.11)
$\displaystyle=:M^{\pm}(n)+L^{\pm}(n),\qquad M^{\pm}(0)=L^{\pm}(0)=0.$

The processes $M^{\pm}$ are martingales, and $L^{\pm}$ are non-decreasing
processes. Recall the processes $X_{n}(\cdot)$ and $Y_{n}(\cdot)$ defined in
(2.19), and define additionally the scaled processes $X^{\pm}_{n}$,
$M^{\pm}_{n}$, $\nu^{\pm}$ similarly to (4.7) (omitting the index $j$), so
that

$\displaystyle X_{n}(t)=X^{+}(t)-X^{-}(t),\quad
M_{n}(t)=M^{+}(t)-M^{-}(t),\quad\nu_{n}(t):=\nu^{+}_{n}(t)+\nu^{-}_{n}(t).$
(5.12)

Analogously to reasoning of the previous section (Propositions 4.1 and 4.4),
we have that the sequence

$\displaystyle\\{(X^{+}_{n},X^{-}_{n},X_{n},M^{+}_{n},M^{-}_{n},M_{n},\nu^{+}_{n},\nu^{-}_{n},\nu_{n},Y_{n})\\}_{n\geq
0}$ (5.13)

is weakly relatively compact in $D(\mathbb{R}_{+},\mathbb{R}^{10})$ each its
limit point $(X^{+},X^{-},X,M^{+},M^{-},M,\nu^{+},\nu^{-},\nu,Y)$ is
continuous, and

$\displaystyle X(t)=X^{+}(t)-X^{-}(t),\quad
M(t)=M^{+}(t)-M^{-}(t),\quad\nu(t)=\nu^{+}(t)+\nu^{-}(t),$ (5.14)

where $M^{+}$ and $M^{-}$ are local martingales with the brackets

$\displaystyle\langle
M^{\pm}\rangle_{t}=\frac{1}{m}\int_{0}^{t}1\\!\\!\\!\;{\rm
I}(X^{\pm}(s)>0)\,\mathrm{d}s.$ (5.15)

Moreover, the process $\sqrt{m}|X|=\sqrt{m}(X^{+}+X^{-})$ is a standard
reflected Brownian motion, $\sqrt{m}\nu$ is its local time at 0, and the
process $\sqrt{m}Y$ is a standard $(m-1)$-dimensional Brownian motion
independent of $X$.

Hence to prove Theorem 5.3 it remains to show that

$\displaystyle\frac{\nu^{+}(t)}{\nu(t)}=\bar{p}\quad\text{a.s.\ for all
$t>0$}.$ (5.16)

To verify this, it suffices to prove that

$\displaystyle\lim_{n\to\infty}\frac{\nu^{+}_{n}(t)}{\nu_{n}(t)}=\lim_{n\to\infty}\frac{L^{+}([nt])}{L([nt])}=\bar{p}\quad\text{a.s.\
for all $t>0$}.$ (5.17)

2\. Let $\tau_{k}$ be the moment of $k$th visit of $X$ to 0, $k\geq 0$,
$\tau_{0}=0$, so that $Z(\tau_{k})=(0,Y(\tau_{k}))$. By the strong Markov
property we have

$\displaystyle\operatorname{Law}(\tau_{k}-(\tau_{k-1}+1))=\operatorname{Law}(\tau|X(0)=1),\quad
k\geq 1,$ (5.18)

and

$\displaystyle\operatorname{Law}(Y(\tau_{k})-Y(\tau_{k-1}))=\operatorname{Law}(Y(\tau)|X(0)=1,Y(0)=0),\quad
k\geq 1,$ (5.19)

where $\tau$ is the first return time to zero of $X$,

$\displaystyle\tau=\inf\\{k\geq 1\colon X(k)=0\\}.$ (5.20)

It is well known that $\tau_{k}<\infty$ a.s. and $\tau_{k}\to+\infty$ a.s. as
$k\to\infty$. Hence

$\displaystyle
Y(\tau_{n})=\sum_{k=1}^{n}\Big{(}Y(\tau_{k})-Y(\tau_{k-1})\Big{)},\ \ {n\geq
0},$ (5.21)

is a random walk on $\mathbb{Z}^{m-1}$ whose jumps have the probability
distribution $\operatorname{Law}(Y(\tau)|X(0)=1,Y(0)=0)$.

We claim that

$\displaystyle\frac{Y(\tau_{n})}{n}\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}S,\quad
n\to\infty,$ (5.22)

where $S$ is a $1$-stable random variable on $\mathbb{R}^{m-1}$ with the
characteristic function $\mathbf{E}\mathrm{e}^{\mathrm{i}\langle
u,S\rangle}=\mathrm{e}^{-|u|/\sqrt{m}}$, $u\in\mathbb{R}^{m-1}$. Indeed, for
each $n\geq 1$ consider the symmetric random walk
$\tilde{Z}=(\tilde{X},\tilde{Y})$ on $\mathbb{Z}^{m}$ starting at
$\tilde{Z}(0)=(n,0,\dots,0)$. Denote

$\displaystyle\tilde{\tau}_{0}$ $\displaystyle=0,$ (5.23)
$\displaystyle\tilde{\tau}_{1}$ $\displaystyle=\inf\\{k\geq
1\colon\tilde{X}(k)=n-1\\},$ $\displaystyle\cdots$
$\displaystyle\tilde{\tau}_{n}$ $\displaystyle=\inf\\{k\geq
1\colon\tilde{X}(n)=0\\}.$

Then clearly

$\displaystyle(\tau_{1},\dots,\tau_{n},Y(\tau_{1}),\dots,Y(\tau_{n}))\stackrel{{\scriptstyle\mathrm{d}}}{{=}}(\tilde{\tau}_{1},\dots,\tilde{\tau}_{n},\tilde{Y}(\tilde{\tau}_{1}),\dots,\tilde{Y}(\tilde{\tau}_{n}))$
(5.24)

By the functional central limit theorem,

$\displaystyle\operatorname{Law}\Big{(}\sqrt{m}\frac{\tilde{Z}([n^{2}\,\cdot])}{n}\Big{|}\tilde{Z}(0)=(n,0,\dots,0)\Big{)}\Rightarrow\operatorname{Law}\Big{(}W(\cdot)\Big{|}W(0)=(1,0\dots,0)\Big{)},$
(5.25)

where $W$ is a standard $m$-dimensional Brownian motion, and thus

$\displaystyle\Big{(}n^{2}\tilde{\tau}_{n},n^{-1}\sqrt{m}\tilde{Y}(\tilde{\tau}_{n})\Big{)}\stackrel{{\scriptstyle\mathrm{d}}}{{\to}}\Big{(}\tau^{W},\big{(}W_{2}(\tau^{W}),\dots,W_{m}(\tau^{W})\big{)}\Big{)},$
(5.26)

where $\tau^{W}=\inf\\{t\geq 0\colon W_{1}(t)=0\\}$. It is well known that
$(W_{2}(\tau^{W}),\dots,W_{m}(\tau^{W}))$ is a 1-stable random vector, see,
e.g., Theorem II.1.16 in Bass [3].

3\. With (5.22) in hand, we apply Theorem 1.3 from Dolgopyat et al. [4], that
states that under the assumption A${}_{\text{SLLN}}(\beta)$ with $\beta>1$

$\displaystyle\frac{1}{n}\sum_{k=0}^{n-1}p_{Y(\tau_{k})}\to\bar{p}\ \
\text{a.s.},\quad n\to\infty.$ (5.27)

4\. Now we are able to finish the proof of (5.16) Then

$\displaystyle L^{+}(n)$ $\displaystyle=\sum_{k=1}^{n}1\\!\\!\\!\;{\rm
I}(X(k-1)=0,X(k)=1)=\sum_{k\colon\tau_{k}\leq n-1}1\\!\\!\\!\;{\rm
I}(X(\tau_{k}+1)=1),$ (5.28) $\displaystyle L(n)$
$\displaystyle=\sum_{k=1}^{n}1\\!\\!\\!\;{\rm I}(X(k-1)=0)=\max\\{k\geq
1\colon\tau_{k}\leq n-1\\}.$

Since $\tau_{k}\to+\infty$ a.s., $L(n)\to\infty$ a.s. as well as at least one
of the processes $L^{+}(\cdot)$ and $L^{-}(\cdot)$.

Consider the process

$\displaystyle Q^{+}(n):=\sum_{k=0}^{n-1}1\\!\\!\\!\;{\rm
I}(X(\tau_{k}+1)=1),$ (5.29)

and note that

$\displaystyle\lim_{n\to\infty}\frac{L^{+}(n)}{L(n)}=\lim_{n\to\infty}\frac{Q^{+}(n)}{n}$
$\displaystyle=\lim_{n\to\infty}\frac{1}{n}\sum_{k=0}^{n-1}\Big{(}1\\!\\!\\!\;{\rm
I}(X(\tau_{k}+1)=1)-p_{Y(\tau_{k})}\Big{)}+\lim_{n\to\infty}\frac{1}{n}\sum_{k=0}^{n-1}p_{Y(\tau_{k})}$
(5.30)
$\displaystyle=\lim_{n\to\infty}\frac{1}{n}\sum_{k=0}^{n-1}\Big{(}1\\!\\!\\!\;{\rm
I}(X(\tau_{k}+1)=1)-p_{Y(\tau_{k})}\Big{)}+\bar{p}.$

The process

$\displaystyle V(n):=\sum_{k=0}^{n-1}\Big{(}1\\!\\!\\!\;{\rm
I}(X(\tau_{k}+1)=1)-p_{Y(\tau_{k})}\Big{)}$ (5.31)

is a martingale difference with

$\displaystyle\mathbf{E}\Big{[}1\\!\\!\\!\;{\rm
I}(X(\tau_{k}+1)=1)-p_{Y(\tau_{k})}\Big{]}=0\quad\text{and}\quad\mathbf{E}\Big{[}1\\!\\!\\!\;{\rm
I}(X(\tau_{k}+1)=1)-p_{Y(\tau_{k})}\Big{]}^{2}\leq 2.$ (5.32)

Hence by the strong law of large numbers for martingales (Theorem 8b in
Chapter II, §3 of Gikhman and Skorokhod [6]) we have convergence

$\displaystyle\frac{V(n)}{n}\to 0\ \ \text{a.s.},\quad n\to\infty,$ (5.33)

what finishes the proof of Theorem 5.3. ∎

## References

  * [1] O. V. Aryasova and A. Yu. Pilipenko.  On Brownian motion on the plane with membranes on rays with a common endpoint.  Random Operators and Stochastic Equations, 17(2):139–157, 2009\.
  * [2] M. T. Barlow, J. Pitman, and M. Yor.  On Walsh’s Brownian motions.  Séminaire de probabilités de Strasbourg, 23:275–293, 1989\.
  * [3] R. F Bass.  Probabilistic Techniques in Analysis.  Probability and its Applications. Springer-Verlag, New York, 1995.
  * [4] D. Dolgopyat, M. Lenci, and P. Nándori.  Global observables for RW: Law of large numbers.  57(1):94–115, 2021.
  * [5] M. Freidlin and S.-J. Sheu.  Diffusion processes on graphs: stochastic differential equations, large deviation principle.  Probability Theory and Related Fields, 116(2):181–220, 2000.
  * [6] I. I. Gikhman and A. V. Skorokhod.  The Theory of Stochastic Processes I, volume 210 of Grundlehren der mathematischen Wissenschaften.  Springer-Verlag, Berlin, 1974.
  * [7] D. S. Grebenkov, D. Van Nguyen, and J.-R. Li.  Exploring diffusion across permeable barriers at high gradients. I. Narrow pulse approximation.  Journal of Magnetic Resonance, 248:153–163, 2014.
  * [8] H. Hajri.  Stochastic flows related to Walsh Brownian motion.  Electronic Journal of Probability, 16:1563–1599, 2011.
  * [9] H. Hajri and W. Touhami.  Itô’s formula for Walsh’s Brownian motion and applications.  Statistics & Probability Letters, 87:48–53, 2014.
  * [10] J. M. Harrison and L. A. Shepp.  On skew Brownian motion.  The Annals of Probability, 9(2):309–313, 1981.
  * [11] A. Iksanov and A. Pilipenko.  A functional limit theorem for locally perturbed random walks.  Probability and Mathematical Statistics, 36(2):353–368, 2016.
  * [12] A. Iksanov, A. Pilipenko, and O. Prykhodko.  Functional limit theorems for random walks perturbed by positive alpha-stable jumps.  arXiv preprint arXiv:2107.00760, 2021.
  * [13] J. Jacod and A. N. Shiryaev.  Limit Theorems for Stochastic Processes, volume 288 of Grundlehren der Mathematischen Wissenschaften.  Springer, Berlin, second edition, 2003.
  * [14] O. Kallenberg.  Foundations of Modern Probability.  Probability and its Applications. Springer, New York, second edition, 2002.
  * [15] I. Karatzas and M. Yan.  Semimartingales on rays, Walsh diffusions, and related problems of control and stopping.  Stochastic Processes and their Applications, 129(6):1921–1963, 2019\.
  * [16] J. Keilson and J. A. Wellner.  Oscillating Brownian motion.  Journal of Applied Probability, 15(2):300–310, 1978.
  * [17] A. Lejay.  On the constructions of the skew Brownian motion.  Probability Surveys, 3:413–466, 2006.
  * [18] A. Lejay.  The snapping out Brownian motion.  The Annals of Applied Probability, 26(3):1727–1742, 2016.
  * [19] V. Mandrekar and A. Pilipenko.  On a Brownian motion with a hard membrane.  Statistics & Probability Letters, 113:62–70, 2016.
  * [20] R. A. Minlos and E. A. Zhizhina.  Limit diffusion process for a non-homogeneous random walk on a one-dimensional lattice.  Russian Mathematical Surveys, 52(2):327–340, 1997.
  * [21] N. Moutal and D. Grebenkov.  Diffusion across semi-permeable barriers: spectral properties, efficient computation, and applications.  Journal of Scientific Computing, 81(3):1630–1654, 2019.
  * [22] S. Nakao.  On the pathwise uniqueness of solutions of one-dimensional stochastic differential equations.  Osaka Journal of Mathematics, 9(3):513–518, 1972.
  * [23] H.-L. Ngo and M. Peigné.  Limit theorem for perturbed random walks.  Theory of Stochastic Processes, 24(40):61–78, 2019.
  * [24] D. S. Novikov, E. Fieremans, J. H. Jensen, and J. A. Helpern.  Random walks with barriers.  Nature Physics, 7(6):508–514, 2011.
  * [25] D. Paulin and D. Szász.  Locally perturbed random walks with unbounded jumps.  Journal of Statistical Physics, 141(6):1116–1130, 2010.
  * [26] V. V. Petrov.  Sums of Independent Random Variables, volume 82 of Ergebnisse der Mathematik und ihrer Grenzgebiete.  Springer-Verlag, New York, 1975.
  * [27] A. Pilipenko.  A functional limit theorem for excited random walks.  Electronic Communications in Probability, 22, 2017.
  * [28] A. Pilipenko and Yu. Pryhod’ko.  Limit behavior of symmetric random walks with a membrane.  Theory of Probability and Mathematical Statistics, 85:93–105, 2012\.
  * [29] A. Y. Pilipenko and L. A. Sakhanenko.  On a limit behavior of a one-dimensional random walk with non-integrable impurity.  Theory of Stochastic Processes, 20(2):97–104, 2015.
  * [30] A. Yu. Pilipenko and V. Khomenko.  On a limit behavior of a random walk with modifications upon each visit to zero.  Theory of Stochastic Processes, 22(1):71–80, 2017.
  * [31] A. Yu. Pilipenko and Yu. E. Prikhod’ko.  On the limit behavior of a sequence of Markov processes perturbed in a neighborhood of the singular point.  Ukrainian Mathematical Journal, 67(4):564–583, 2015.
  * [32] A. Yu. Pilipenko and Yu. E. Prykhodko.  Limit behavior of a simple random walk with non-integrable jump from a barrier.  Theory of Stochastic Processes, 19(1):52–61, 2014.
  * [33] N. I. Portenko.  Generalized Diffusion Processes, volume 83 of Translations of Mathematical Monographs.  American Mathematical Society, Providence, RI, 1990.
  * [34] D. Revuz and M. Yor.  Continuous Martingales and Brownian Motion, volume 293 of Grundlehren der Mathematischen Wissenschaften.  Springer, Berlin, corrected third printing of the third edition, 2005\.
  * [35] D. Szász and A. Telcs.  Random walk in an inhomogeneous medium with local impurities.  Journal of Statistical Physics, 26(3):527–537, 1981.
  * [36] J. B. Walsh.  A diffusion with a discontinuous local time.  Astérisque, 52(53):37–45, 1978.
  * [37] D. A. Yarotskii.  Central limit theorem for a class of nonhomogeneous random walks.  Mathematical Notes, 69(5):690–695, 2001.


longtable

# Tight relative estimation in the mean of Bernoulli random variables

Mark Huber

#### Abstract

Given a stream of Bernoulli random variables, consider the problem of
estimating the mean of the random variable within a specified relative error
with a specified probability of failure. Until now, the Gamma Bernoulli
Approximation Scheme (GBAS) was the method that accomplished this goal using
the smallest number of average samples. In this work, a new method is
introduced that is faster when the mean is bounded away from zero. The process
uses a two-stage process together with some simple inequalities to get
rigorous bounds on the error probability.

##  1 Introduction

Suppose that $B_{1},B_{2},\ldots$ are a stream of independent, identically
distributed (iid) random variables that have a Bernoulli distribution. That
is, for each $i$, ${\mathbb{P}}(B_{i}=1)=p$ and ${\mathbb{P}}(B_{i}=0)=1-p$.
Consider the problem of estimating $p$.

This is one of the oldest problems in probability, and was the inspiration for
the first Law of Large Numbers of Jacob Bernoulli [4]. There are numerous
modern applications, including network analysis [2], testing sequence
association [3], testing if a population is in Hardy-Weinburg equilibrium
[10], ray tracing [15], and many others [7, 1, 13, 8, 14].

The basic algorithm simply generates a fixed number of $B_{i}$ values, say
$n$, and then estimates $p$ using the sample average
$\hat{p}_{n}=(B_{1}+\cdots+B_{n})/n$. While both simple and unbiased, this
method fails to have any guarantees on the relative error of the estimate.
Taking a random number of samples can help create an estimate with such a
guarantee.

The goal here will be, given parameters $\epsilon>0$ and $\delta\in(0,1)$,
find an integrable stopping time $T$ with respect to the natural filtration
and an estimate

$\hat{p}=\hat{p}(B_{1},B_{2},\ldots,B_{T})$ (1)

such that the relative error is greater than $\epsilon$ with probability at
most $\delta$. That is,

${\mathbb{P}}\left(\left|\frac{\hat{p}}{p}-1\right|>\epsilon\right)\leq\delta.$
(2)

This is also known as a _randomized approximation scheme_ , or _ras_ for
short.

The time to find the estimate is dominated by the time needed to generate
$B_{1},\ldots,B_{T}$. Hence the running time of the procedure is taken to be
${\mathbb{E}}(T)$, the average number of samples needed by the ras.

Returning to the sample average for the moment, the variance of $\hat{p}_{n}$
is $v_{n}=p(1-p)/n$. The standard deviation is then $\sqrt{v_{n}}$. To get the
standard deviation to be $\epsilon p$, $n=(1-p)\epsilon^{-2}p^{-1}$ samples
are necessary. Using this together with the mean of medians technique [11,
12], it is possible to get an ras for this sample average using

$\Theta((1-p)p^{-1}\epsilon^{-2}\ln(\delta^{-1}))$ (3)

samples. Of course, here the Big-Theta notation ($\Theta$) is hiding the
constant involved.

While the Central Limit Theorem (because it is only an asymptotic result) does
not yield an ras, it does give a rough estimate of the number of samples
required. A CLT approximation indicates that about
$2(1-p)p^{-1}\epsilon^{-2}\ln(2\delta^{-1})$ samples should be necessary and
sufficient to obtain an $(\epsilon,\delta)$-ras.

This is a problem, because the number of samples needed has a factor of
$(1-p)/p$, and $p$ is the very thing being estimated! Dagum, Karp, Luby, and
Ross [5] found an elegant solution to this difficulty through an algorithm
that will be called DKLR throughout this work.

Their idea was to use the Bernoulli samples to create iid geometric random
variables. Say that $G\sim\textsf{Geo}(p)$ if for all $i\in\\{1,2,\ldots\\}$,
${\mathbb{P}}(G=i)=p(1-p)^{i-1}$. Then it is straightforward to see that for

$R=\inf\\{r:B_{1}+\cdots+B_{r}=1\\},$ (4)

then $R\sim\textsf{Geo}(p)$.

The mean of $R$ is $1/p$, and the variance is $(1-p)/p^{2}$. Hence for the
sample average of $k$ iid draws $R_{1},\ldots,R_{k}$ with the same
distribution as $R$,

$\frac{R_{1}+\cdots+R_{k}}{k}$ (5)

has mean $1/p$ and standard deviation $(1-p)^{1/2}/[k^{1/2}p]$.

That is, the standard deviation is already of the same order as the mean being
estimated! The DKLR algorithm then runs as follows.

#### DKLR

  1. 1.

Given $\epsilon$ and $\delta$, set $k$ equal to

$\lceil 1+(1+\epsilon)4(e-2)\epsilon^{-2}\ln(2\delta^{-1})\rceil.$

  2. 2.

Use the Bernoulli stream to generate $R_{1},\ldots,R_{k}$ iid
$\textsf{Geo}(p)$.

  3. 3.

Return

$\hat{p}_{\textsf{DKLR}}=\frac{k-1}{R_{1}+\cdots+R_{k}}.$

The choice of $k$ in the first line comes from [5], and is based on their
analysis of the value needed to obtain an $(\epsilon,\delta)$-ras. The use of
$k-1$ in the numerator of the final line is a change from their algorithm,
which used $k$. The $k-1$ used here makes this algorithm more easily
comparable to the improvement on the algorithm, GBAS, that is described later
on in this section.

While this was the first ras for this problem, there are some drawbacks.
First, the constant of $4(e-2)\approx 2.873\ldots$ is not close to the
constant of 2 predicted by the CLT. Second, it does not take advantage of
second order effects, that is, the tails of the geometric random variable
decline in a similar fashion to normals, allowing for second order reductions
in the run time that are significant. Third, it is biased.

To solve these issues, the author introduced the _Gamma Bernoulli
Approximation Scheme_ (GBAS) [9]. This algorithm took the reduction of
Bernoullis to Geometrics one step further. An easy to show fact about
Geometric random variables is the following. Suppose that $A$ is an
exponentially distributed random variable with rate $\lambda>0$ so that for
all $a>0$, it holds that ${\mathbb{P}}(A>a)=\exp(-\lambda a)$. Write
$A\sim\textsf{Exp}(\lambda)$.

######  Theorem 1.1.

If $G\sim\textsf{Geo}(p)$ and $A_{1},\ldots,A_{G}$ are iid exponential random
variables of rate 1, then $A_{1}+\cdots+A_{G}\sim\textsf{Exp}(p)$.

(See, for instance, [9] for a proof.)

For $S\sim\textsf{Exp}(p)$, the mean of $S$ is $1/p$ and standard deviation is
also $1/p$. Moreover, for $S_{1},\ldots,S_{k}$ iid $\textsf{Exp}(p)$,
$M=S_{1}+\cdots+S_{k}$ has a gamma distribution with shape parameter $k$ and
rate parameter $p$. Write $M\sim\textsf{Gamma}(k,p)$.

The density of both $M$ and $1/M$ (which has an inverse gamma distribution)
are known precisely, and the mean of $1/M$ is $p/(k-1)$. So

$\hat{p}_{\textsf{GBAS}}=(k-1)/M$ (6)

is an unbiased estimator for $p$.

In addition, gamma random variables are scalable. Multiplying a gamma
distributed random variable by a constant results in a new gamma distributed
random variable where the new rate constant is the old constant divided by the
multiplication factor. Hence

$\frac{p}{\hat{p}_{\textsf{GBAS}}}=\frac{pM}{k-1}\sim\textsf{Gamma}(k,k-1).$
(7)

That is to say, the relative error of the GBAS estimate does not depend on the
value $p$ that is being estimated! This makes it a straightforward numerical
calculation to find the minimum value of $k$ such that for
$X\sim\textsf{Gamma}(k,k-1)$,

${\mathbb{P}}\left(\left|\frac{1}{X}-1\right|>\epsilon\right)\leq\delta.$ (8)

This allows the GBAS algorithm to immediately give an $(\epsilon,\delta)$-ras.

#### GBAS

  1. 1.

Given $\epsilon$ and $\delta$, set $k$ equal to the smallest value such that a
$\textsf{Gamma}(k,k-1)$ random variable lies outside the interval
$[1/(1+\epsilon),1/(1-\epsilon)]$ with probability at most $\delta$.

  2. 2.

Use the Bernoulli stream to generate $S_{1},\ldots,S_{k}$ iid
$\textsf{Exp}(p)$.

  3. 3.

Return

$\hat{p}_{\textsf{GBAS}}=\frac{k-1}{S_{1}+\cdots+S_{k}}.$

This corrects most of the three deficiencies noted earlier for the DKLR
algorithm. First, the value of $k$ can be shown to be at most
$2\epsilon^{-2}\ln(2\delta^{-1})$ for $\epsilon$ and $\delta$ both less than
1. This follows directly from Theorem 1 of [6]. Second, the tails of a gamma
decline much like those of a normal, and so the second order effects are
retained. Third, the result is unbiased.

Unfortunately, one new deficiency has been added. For a stream of $R_{i}$ that
are iid $\textsf{Geo}(p)$, the mean of $R_{1}+\cdots+R_{k}$ is $1/p$ while the
standard deviation is the slightly smaller $(1-p)^{1/2}/p$. For $S_{i}$ iid
$\textsf{Exp}(p)$ the mean is $1/p$ while the standard deviation is exactly
$1/p$.

From a Monte Carlo point of view, this means that the DKLR algorithm (if $p$
was known ahead of time) should take about $1-p$ times the number of samples
that the GBAS algorithm does. While this does not affect much if $p$ is close
to 0, for $p$ in the middle of the interval $[0,1]$ this can have a
significant effect on the running time. Of course, $p$ is not known ahead of
time, and so there has been no way of getting this $1-p$ factor in the run
time exactly.

This paper for the first time presents a way of running DKLR close to
optimally in order to regain a factor that is asymptotically close for small
$\epsilon$ and $\delta$ to the $1-p$ factor in the run time. The rest of the
paper is organized as follows. The next section introduces the two-stage
method of running and proves the veracity of the estimate. The following
section covers how to gain most of the advantage of the $1-p$ factor while
still obtaining an unbiased estimate. The last section then concludes.

##  2 The two-stage algorithm

Suppose that $U_{1},U_{2},\ldots$ forms an iid sequence of random variables
that are uniform over $[0,1]$. Suppose $B_{i}={\mathbbm{1}}(U_{i}\leq p)$,
where ${\mathbbm{1}}$ is the indicator function that evaluates to 0 if the
argument is false and 1 if it is true. Then the $B_{1},B_{2},\ldots$ form an
iid Bernoulli sequence with mean $p$. Let

$T_{k}(p)=\inf\\{t:{\mathbbm{1}}(U_{1}<p)+{\mathbbm{1}}(U_{2}<p)+\cdots{\mathbbm{1}}(U_{t}<p)=k\\}.$

That is, $T_{k}(p)$ is the smallest time $t$ such that in the first $t$
members of the Bernoulli iid sequence with mean $p$, there were exactly $k$ of
the ${\mathbbm{1}}(U_{i}<p)$ values equal to 1 and the rest are equal to 0. In
particular, $T_{1}(p)\sim\textsf{Geo}(p)$. Then $T_{k}$ is the sum of $k$
independent geometric random variables, and so has a _negative binomial_
distribution with parameters $k$ and $p$. Write
$T_{k}\sim\textsf{NegBin}(k,p)$.

######  Theorem 2.1.

Let $p\leq p^{\prime}$, $k\in\\{1,2,\ldots\\}$,
$T_{k}(p)\sim\textsf{NegBin}(k,p)$ and
$T_{k}(p^{\prime})\sim\textsf{NegBin}(k,p^{\prime})$. Then for all $a$,
${\mathbb{P}}(T_{k}(p)\leq a)\leq{\mathbb{P}}(T_{k}(p^{\prime})\leq a)$.

###### Proof.

Let $U_{1},U_{2},\ldots$ be iid uniforms over $[0,1]$ used to form $T_{k}(p)$
and $T_{k}(p^{\prime})$.

Then the Bernoulli random variables ${\mathbbm{1}}(U_{i}<p)$ are monotonically
increasing in $p$, so $T_{k}(p)$ is monotonically decreasing in $p$. That is,
for $p\leq p^{\prime}$, $T_{k}(p)\geq T_{k}(p^{\prime})$. This means that for
all $a$, ${\mathbb{P}}(T_{k}(p)\leq a)\leq{\mathbb{P}}(T_{k}(p^{\prime})\leq
a)$.

This holds for negative binomials coupled using the same uniforms, but the
theorem statement is a statement about the marginal probabilities, which holds
or fails regardless of the coupling. Since it holds for one coupling, it holds
for all. ∎

What does this do for us? In the DKLR algorithm, let
$T_{k}(p)=R_{1}+\cdots+R_{k}$. Then failure occurs when
$\hat{p}_{\textsf{DKLR}}=k/T_{k}(p)$ is either too large or too small. Suppose
that

$p_{1}\leq p\leq p_{2}$

where $p$ is unknown to the user but $p_{1}$ and $p_{2}$ are both known. Then
the probability that $T_{k}(p_{1})$ is too large is an upper bound on the
probability that $T_{k}(p)$ is too large for all $p\in[p_{1},p_{2}]$.
Similarly, the probability that $T_{k}(p_{2})$ is too small is an upper bound
on the probability that $T_{k}(p)$ is too small for all $p\in[p_{1},p_{2}]$.

To bound the upper tail for relative error, note that

$\displaystyle{\mathbb{P}}\left(\frac{k-1}{T_{k}(p)}>(1+\epsilon)p\right)$
$\displaystyle=$
$\displaystyle{\mathbb{P}}\left(T_{k}(p)<\frac{k-1}{(1+\epsilon)p}\right)$ (9)
$\displaystyle\leq$
$\displaystyle{\mathbb{P}}\left(T_{k}(p)<\frac{k-1}{(1+\epsilon)p_{1}}\right)$
(10) $\displaystyle\leq$
$\displaystyle{\mathbb{P}}\left(T_{k}(p_{2})<\frac{k-1}{(1+\epsilon)p_{1}}\right)$
(11)

Similarly,

${\mathbb{P}}\left(\frac{k-1}{T_{k}(p)}<(1-\epsilon)p\right)\leq{\mathbb{P}}\left(T_{k}(p_{1})>\frac{k-1}{(1-\epsilon)p_{2}}\right)$
(12)

At this point the upper and lower tails for the estimate with $p$ (the size of
which are unknown) have been bounded in terms of an upper tail for $p_{1}$ and
a lower tail for $p_{2}$. These bounds can be used to find a value for $k$ to
use in DKLR that gives an $(\epsilon,\delta)$-ras for all $p\in[p_{1},p_{2}]$.

#### find_k

  1. 1.

Given $\epsilon$, $\delta$, and $a<b$, return the smallest value of $k$ such
that

${\mathbb{P}}\left(T_{k}(a)>\frac{k-1}{(1-\epsilon)b}\right)+{\mathbb{P}}\left(T_{k}(b)<\frac{k-1}{(1+\epsilon)a}\right)\leq\delta.$

The output of this algorithm is, for all $p\in[a,b]$, a value of $k$ that
yields an $(\epsilon,\delta)$-ras using DKLR. The problem is, if the length of
the interval $b-a$ is very large, the value of $k$ this gives could be very
much larger than optimal. To solve this issue, break $[a,b]$ into smaller
intervals and find $k$ for each one of the smaller intervals. Then the largest
of these $k$ values found will work for all intervals taken individually, and
hence works also for the union of the smaller intervals, which is just
$[a,b]$.

In particular, if all that is known is that $p\in[a,1]$, then it is possible
to partition the large interval $[a,1]$ into $m$ small intervals using $a_{i}$
satisfying

$a=a_{0}<a_{1}<a_{2}<\cdots<a_{m}=1.$

The probability that the $T$ is too large or too small for $p\in[a,1]$ is
bounded above by the maximum value that $T$ is too large or too small for each
of the subintervals $[a_{i},a_{i+1}]$. This leads to the following algorithm.

#### two_stage_dklr

  1. 1.

Given $\epsilon$ and $\delta$, run GBAS with $\sqrt{\epsilon}$ and $\delta/2$.
This gives an estimate $\hat{p}_{1}$. Then for $a=\hat{p}_{1}/(1+\epsilon)$,
${\mathbb{P}}(p<a)\leq\delta/2$.

  2. 2.

Break the interval $[a,1]$ into subintervals of width at most
$\epsilon(1-a)/10$.

  3. 3.

For each subinterval, calculate the smallest value of $k$ (using find_k) that
yields for $\epsilon$ error at most $\delta/2$.

  4. 4.

Draw $R_{1},\ldots,R_{k}$ iid $\textsf{Geo}(p)$ using the iid stream of
Bernoulli random variables.

  5. 5.

Return $(k-1)/(R_{1}+\cdots+R_{k})$.

From the previous results, the following holds.

######  Theorem 2.2.

Let $\hat{p}$ be the output of the algorithm two_stage_dklr. Then $\hat{p}$
satisfies

${\mathbb{P}}(|(\hat{p}/p)-1|>\epsilon)\leq\delta).$ (13)

The choice of $\epsilon(1-a)/10$ for the subinterval width ensures that the
$\epsilon$ relative error is much larger than the widths, and so errors
arising from the width choice are negligible. The choice to use 10 here is not
important, larger values like 100 are still valid, but will take longer to
compute the error probabilities.

###  2.1 Tilting

In [6], the idea of _tilting_ the resulting estimate was introduced. Tilting
takes advantage of the fact that one tail of the sum of geometrics is
unbounded and goes down slightly more slowly, while the other is bounded by
$k$ and so goes down slightly more quickly. By adjusting the estimate by a
multiplicative factor that is close to 1, the asymptotic behavior of the tails
can be balanced to go down as rapidly as possible.

In the case of both DKLR and GBAS, balancing the tails requires dividing the
estimate from earlier by

$c_{\text{tilt}}=\frac{2\epsilon}{(1-\epsilon^{2})\ln(1+2\epsilon/(1-\epsilon))}.$

This is very close to 1. More precisely,
$c_{\text{tilt}}=1+(2/3)\epsilon^{2}+O(\epsilon^{4})$.

###  2.2 Running time

How does the two_stage_dklr with tilting perform in practice? Generally
speaking, it is better than both the original DKLR and GBAS, especially when
$p$ is close to 1.

The running time of all these algorithms is determined by $k/p$. The following
table shows what value of $k$ is needed for the first and second stages in
two_stage_dklr. The first stage sets the error at $\epsilon^{-1/2}$, and since
the GBAS method uses $\Theta(\epsilon^{-2})$ samples, the first stage uses
approximately

$(\epsilon^{1/2})^{-2}=\epsilon^{-1}/(1-p)$ (14)

times the overall optimal number of samples.

Now consider the second stage. The second stage runs slower the smaller the
estimate $\hat{p}_{1}$ is. The smallest it can be (assuming stage one was
succesful) is $(1-\sqrt{\epsilon})p$. It is divided by $1+\sqrt{\epsilon}$ to
get a lower bound on $p$. Taken together, the gap between 1 and the lower
bound on $p$ is

$1-p(1-\sqrt{\epsilon})/(1+\sqrt{\epsilon}).$ (15)

Therefore, this factor of samples is gained, but $1/(1-p)$ is still lost. So
stage two uses about

$\frac{1-p(1-\sqrt{\epsilon})/(1+\sqrt{\epsilon})+\epsilon}{1-p}$ (16)

times the optimal number of DKLR samples, which itself has a speedup of
$1/(1-p)$ over GBAS.

A final effect to consider is that GBAS has $\delta$ as the error bound, while
each stage of two stage DKLR has $\delta/2$ as the error bound so that the
total probability of error is at most $\delta$ by the union bound. Combining
this effect with the earlier calculations, the speedup from GBAS to two stage
DKLR is about

$\rho=\frac{\ln(\delta^{-1})}{\ln(2)+\ln(\delta^{-1})}\cdot\frac{1}{1-p(1-\sqrt{\epsilon})/(1+\sqrt{\epsilon})+\epsilon}.$
(17)

This is just an estimate of the speedup, but as Table 1 shows, this estimate
is fairly accurate. As $\epsilon$ and $\delta$ approach zero, this estimate
$\rho$ of the speedup approaches the optimal speedup of $1/(1-p)$.

Table 1: Determination of $k$ for various values of $(p,\epsilon,\delta)$. The speedup is the k for GBAS divided by the sum of the $k$ values for Stage 1 and Stage 2. The absolute best this could be for any algorithm is $1/(1-p)$. The estimate of the speedup for the two stage algorithm is $\rho=\ln(\delta^{-1})/[(2+\ln(\delta^{-1}))/(1-p(1-\sqrt{\epsilon})/(1+\sqrt{\epsilon})+\epsilon)]$. $p$ |  $\epsilon$ |  $\delta$ |  $k$ GBAS |  $k$ Stage 1 |  $k$ Stage 2 |  Speedup |  $\rho$ |  $1/(1-p)$
---|---|---|---|---|---|---|---|---
0.9 |  0.10 |  1e-02 |  661 |  76 |  413 |  1.35 |  1.37 |  10.00
0.9 |  0.10 |  1e-06 |  2380 |  239 |  1317 |  1.53 |  1.51 |  10.00
0.9 |  0.01 |  1e-06 |  239268 |  2513 |  66203 |  3.48 |  3.48 |  10.00
0.5 |  0.10 |  1e-02 |  661 |  76 |  551 |  1.05 |  1.03 |  2.00
0.5 |  0.10 |  1e-06 |  2380 |  239 |  1760 |  1.19 |  1.13 |  2.00
0.5 |  0.01 |  1e-06 |  239268 |  2513 |  145055 |  1.62 |  1.58 |  2.00
0.1 |  0.10 |  1e-02 |  661 |  76 |  595 |  0.99 |  0.83 |  1.11
0.1 |  0.10 |  1e-06 |  2380 |  239 |  1901 |  1.11 |  0.91 |  1.11
0.1 |  0.01 |  1e-06 |  239268 |  2513 |  191853 |  1.23 |  1.03 |  1.11

##  3 Unbased approach

While two_stage_dklr gives the fewest number of samples, it does suffer from
the fact that the estimate is biased. To deal with this issue while
maintaining most of the benefit from the $1-p$ factor, a _shifted grid
approach_ can be used.

For positive integer $n$, consider a grid of points
$x=(0,1/n,2/n,\ldots,(n-1)/n)$. If $U$ is uniform over $[0,1]$, then
$x+U\text{ mod }1$ consists of a sets of points that are not independent, but
which are identically distributed. That is, each point will be uniform over
$[0,1]$. As usual, mod 1 means that for any value in $[1,2)$, 1 is subtracted
from the value. For $U\sim\textsf{Unif}([0,1])$, set

$(U_{1},\ldots,U_{n})=\left(\frac{0+U}{n},\frac{1+U}{n},\frac{2+U}{n},\ldots\frac{n-1+U}{n}\right)\text{
mod }1.$ (18)

The GBAS method first uses the Bernoulli values to draw $M$ as the sum of $k$
iid geometric random variables with mean $1/p$. Then given $M$, draw $G$ as
gamma with shape parameter $M$ and rate parameter 1. The result is a gamma
distribution with parameters $k$ and $p$.

Given a uniform $U$ over $[0,1]$, the _inverse transform method_ (ITM) can be
used to draw from a given one dimensional distribution. For a random variable
$X$ with cumulative distribution function $F_{X}$, the pseudoinverse is
defined as

$F_{X}^{-1}(a)=\inf\\{b:F_{X}(b)\geq a\\}.$ (19)

The ITM method relies on the fact that for $U$ uniform over $[0,1]$, $X$ and
$F_{X}^{-1}(U)$ have the same distribution.

This implies that for all $i\in\\{1,\ldots,n\\}$,

$\hat{p}_{i}=\frac{k-1}{F_{[G\mid M]}^{-1}(U_{i})}$ (20)

is an unbiased estimator of $p$, and the sample average of the $\hat{p}_{i}$
will be as well.

This sample average $\hat{p}_{i}$ tends to be close to $1/M$ (which is the
DKLR estimate), as seen next.

###  3.1 Non-interval centered unbiased estimate

By using the new DKLR method, a confidence interval of width at most
$2\epsilon p$ is created with the biased estimate at its center. In the
process, a negative binomial $M$ with parameters $k$ and $p$ is created. At
the same time, using the shifted grid method an unbiased estimator can also be
found at the same time from $M$.

The unbiased estimator will most likely be very close to the estimate
$(k-1)/M$. This is because the relative error between
$[G|M]\sim\textsf{Gamma}(M,1)$ and $M$ is typically small. Let

$\hat{p}_{\text{b}}=(k-1)/M$ (21)

be the biased estimate and for $U_{i}$ drawn using (18),

$\hat{p}_{\text{unb}}=(k-1)\frac{1}{n}\sum_{i=1}^{n}1/F^{-1}_{[G\mid
M]}(U_{i})$ (22)

be the unbiased estimate.

In calculating $\hat{p}_{\text{unb}}/\hat{p}_{\text{b}}$, the $k-1$ factor
cancels out, leaving us with

$r=\frac{\hat{p}_{\text{unb}}}{\hat{p}_{\text{b}}}=M\frac{1}{n}\sum_{i=1}^{n}1/F^{-1}_{[G\mid
M]}(U_{i}).$ (23)

When $r=1$, the unbiased and biased entries are the same, so $|r-1|$ measures
the relative error between the two. With probability $\delta_{1}$,
$U\in[\delta_{1}/2,1-\delta_{1}/2]$. Plugging these interval endpoints for $U$
into (23) gives an upper bound on $|r-1|$.

Table 2 shows the bounds on the distribution of $r$ for some other values of
the variables. As seen in this table, the chance of even a small relative
difference between the biased and unbiased estimates becomes very small very
fast as $M$ grows. It should also be noted that for $n=10^{3}$, the unbiased
estimate can be computed in a fraction of a second.

Table 2: Bounds on the relative error between the unbiased and biased estimate. Here $r$ is the unbiased estimate divided by the biased estimate. $M$ |  $n$ |  $\delta_{1}$ |  $x$ such that ${\mathbb{P}}(|r-1|\leq x)=1-\delta_{1}$
---|---|---|---
$10^{4}$ |  $10^{3}$ |  $10^{-6}$ |  $0.0001497$
$10^{4}$ |  $10^{4}$ |  $10^{-6}$ |  $0.0001049$
$10^{4}$ |  $10^{3}$ |  $10^{-8}$ |  $0.0001597$
$10^{4}$ |  $10^{2}$ |  $10^{-8}$ |  $0.0006899$
$10^{5}$ |  $10^{3}$ |  $10^{-8}$ |  $0.0000283$

##  4 Conclusion

Although the DKLR algorithm provides an $(\epsilon,\delta)$-ras for estimating
the mean of a Bernoulli random variable from a stream of such variables, it
has not been widely used. This is partially due to loose bounds on the number
of samples needed. The best speedup over existing methods that could be
achieved is $1/(1-p)$. This work gives a new version of the algorithm that
uses a number of samples that approaches the optimal speedup. Roughly, the
speedup for this two stage algorithm is

$\frac{\ln(\delta^{-1})}{\ln(2)+\ln(\delta^{-1})}\cdot\frac{1}{1-p(1-\sqrt{\epsilon})/(1+\sqrt{\epsilon})+\epsilon}.$

By using a randomly shifted grid and the inverse transform method, this
approach can also yield an unbiased estimate that is very close to the biased
estimate. Therefore the interval surrounding the biased estimate can also be
used to give an interval around the unbiased estimate of only slightly smaller
size.

## References

  * [1] Arvind Arasu, Michaela Götz, and Raghav Kaushik.  On active learning of record matching packages.  In Proceedings of the 2010 ACM SIGMOD International Conference on Management of data, pages 783–794, 2010.
  * [2] Susan Athey, Dean Eckles, and Guido W Imbens.  Exact $p$-values for network interference.  Journal of the American Statistical Association, 113(521):230–240, 2018.
  * [3] Roger Bakeman, Byron F Robinson, and Vicenç Quera.  Testing sequential association: Estimating exact $p$ values using sampled permutations.  Psychological methods, 1(1):4, 1996.
  * [4] Jakob Bernoulli.  Ars Conjectandi.  Impensis Thurnisiorum, fratrum, 1713.
  * [5] P. Dagum, R. Karp, M. Luby, and S. Ross.  An optimal algorithm for Monte Carlo estimation.  Siam. J. Comput., 29(5):1484–1496, 2000.
  * [6] J. Feng, M. Huber, and Y. Ruan.  Monte Carlo with user-specified relative error.  In P. W. Glynn and A. Owen, editors, Proceedings in Mathematics & Statistics: Monte Carlo and Quasi-Monte Carlo methods, volume 241, chapter 12. Springer, 2018.
  * [7] Haipeng Guo and William Hsu.  A survey of algorithms for real-time Bayesian network inference.  In Join Workshop on Real Time Decision Support and Diagnosis Systems, pages 1–12, 2002.
  * [8] Jiewen Huang, Lyublena Antova, Christoph Koch, and Dan Olteanu.  Maybms: a probabilistic database management system.  In Proceedings of the 2009 ACM SIGMOD International Conference on Management of data, pages 1071–1074, 2009.
  * [9] M. Huber.  A Bernoulli mean estimate with known relative error distribution.  Random Structures Algorithms, 50:173–182, 2017.  arXiv:1309.5413.
  * [10] Mark Huber, Yuguo Chen, Ian Dinwoodie, Adrian Dobra, and Mike Nicholas.  Monte Carlo algorithms for Hardy–Weinberg proportions.  Biometrics, 62(1):49–53, 2006.
  * [11] M. Jerrum, L. Valiant, and V. Vazirani.  Random generation of combinatorial structures from a uniform distribution.  Theoret. Comput. Sci., 43:169–188, 1986.
  * [12] Mark Jerrum and Alistair Sinclair.  Approximating the permanent.  SIAM journal on computing, 18(6):1149–1178, 1989.
  * [13] Christoph Koch and Dan Olteanu.  Conditioning probabilistic databases.  arXiv preprint arXiv:0803.2212, 2008.
  * [14] Dan Olteanu, Jiewen Huang, and Christoph Koch.  Approximate confidence computation in probabilistic databases.  In 2010 IEEE 26th International Conference on Data Engineering (ICDE 2010), pages 145–156. IEEE, 2010.
  * [15] Dongeek Shin, Ahmed Kirmani, Vivek K Goyal, and Jeffrey H Shapiro.  Computational 3d and reflectivity imaging with high photon efficiency.  In 2014 IEEE International Conference on Image Processing (ICIP), pages 46–50. IEEE, 2014.

##  Appendix A Code for the paper

The following code produces the tables in the main text. First load in the
needed libraries.

library(kableExtra)

library(tidyverse)

###  A.1 Preliminary functions

The first function, tilt_value, calculates the amount to divide the unbiased
estimate in GBAS to make the tails go down at the same rate.

tilt_value <- function(epsilon)

2 * epsilon / (1 - epsilon^2) / log(1 + 2 * epsilon /

(1 - epsilon))

This next function calculates the minimum possible $k$ needed for a given
$\epsilon$ and $\delta$ for GBAS. Note that once the function has been run
once, the parameter min_k can be set appropriately to make it run more quickly
on a second run.

find_k_GBAS <- function(epsilon, delta, tilt = FALSE, min_k = 1) {

tilt_constant <- 1 - tilt + tilt_value(epsilon) * tilt

k <- min_k

accept <- FALSE

while (!accept) {

k <- k + 1

high_error <- 1 - pgamma(1 / (1 - epsilon),

shape = k,

rate = (k - 1) / tilt_constant)

low_error <- pgamma(1 / (1 + epsilon),

shape = k,

rate = (k - 1) / tilt_constant)

accept <- (high_error + low_error) < delta

}

return(k)

}

It is helpful to have a vector version of this function.

find_k_GBAS_vector <- function(epsilon, delta, tilt = FALSE, min_k = 1) {

A <- cbind(epsilon, delta, tilt, min_k)

return(apply(A, 1, function(x) find_k_GBAS(x[1], x[2], x[3], x[4])))

}

For a minimum value of $p$, the find_error_dklr_small_num function finds the
largest probability of the DKLR estimate being too small and too large, and
returns the results as a length two vector.

find_error_dklr_small_num <- function(minp, epsilon, k, tilt = FALSE) {

change_p <- (1 - minp) * epsilon / 100

tilt <- tilt_value(epsilon)

seqp <- seq(minp, 1, by = change_p)

lower <- seqp[-length(seqp)]

upper <- seqp[-1]

lower_error <- 1 - pnbinom((k - 1) / (lower * (1 - epsilon)) / tilt - k, k,
upper)

upper_error <- pnbinom((k - 1) / (upper * (1 + epsilon)) / tilt - k, k, lower)

A <- cbind(seqp[-length(seqp)],

seqp[-1],

lower_error,

upper_error)

return(c(max(A[,3], na.rm = TRUE), max(A[,4], na.rm = TRUE)))

}

find_k_DKLR <- function(min_p, epsilon, delta, tilt = FALSE, min_k = 1) {

tilt_constant <- 1 - tilt + tilt_value(epsilon) * tilt

k <- min_k

accept <- FALSE

while (!accept) {

k <- k + 1

error <- sum(find_error_dklr_small_num(min_p, epsilon, k, tilt))

accept <- error < delta

if ((k - min_k) %% 1000 == 0) print(k)

}

return(k)

}

It is helpful to have a version of this function that can take vectors as
input and return a vector of outputs.

find_k_DKLR_vector <-

function(min_p, epsilon, delta, tilt = FALSE, min_k = 1)

apply(cbind(min_p, epsilon, delta, tilt, min_k), 1,

function(x) find_k_DKLR(x[1], x[2], x[3], x[4], x[5]))

###  A.2 Making Table 1

Set up the values needed for Table 1.

p <- c(rep(0.9, 3), rep(0.5, 3), rep(0.1, 3))

epsilon <- rep(c(0.1, 0.1, 0.01), 3)

delta <- rep(c(10^(-2), 10^(-6), 10^(-6)), 3)

Reset the variable k_min1 and k_min2 after the first run to speed up
subsequent runs of the code.

# Set up for GBAS

k_gbas <- rep(1, length(p))

# Set up for stage one

k_min1 <- rep(1, length(p))

k_min1 <- c(74, 237, 2511, 74, 237, 2511, 74, 237, 2511)

# Set up for stage two

k_min2 <- rep(1, length(p))

k_min2 <- c(411, 1315, 66200, 549, 1758, 145053, 593, 1899, 191851)

table1 <- tibble(p, epsilon, delta, tilt = TRUE, k_min1) |>

mutate(sq_ep = sqrt(epsilon),

ha_delta = delta / 2,

low_p = p * (1 - sq_ep) / (1 + sq_ep)) |>

mutate(`k GBAS` =

find_k_GBAS_vector(epsilon, delta, tilt, k_gbas),

.before = `tilt`) |>

mutate(`k Stage 1` =

find_k_GBAS_vector(sq_ep, ha_delta, tilt, k_min1),

.before = `tilt`) |>

mutate(`k Stage 2` =

find_k_DKLR_vector(low_p, epsilon, ha_delta, TRUE, k_min2),

.before = tilt) |>

mutate(speedup = round(`k GBAS` / (`k Stage 1` + `k Stage 2`), 2),

rho = round(

log(delta^(-1)) / (log(2 * delta^(-1))) /

(1 - p * (1 - sq_ep) / (1 + sq_ep) + epsilon), 2),

geo_speed = round(1 / (1 - p), 2),

.before = tilt) |>

select(1:9)

table1 |>

kable(booktabs = TRUE) |>

kable_styling(latex_options = c("scale_down"))

p | epsilon | delta | k GBAS | k Stage 1 | k Stage 2 | speedup | rho | geo_speed
---|---|---|---|---|---|---|---|---
0.9 | 0.10 | 1e-02 | 661 | 76 | 413 | 1.35 | 1.37 | 10.00
0.9 | 0.10 | 1e-06 | 2380 | 239 | 1317 | 1.53 | 1.51 | 10.00
0.9 | 0.01 | 1e-06 | 239268 | 2513 | 66203 | 3.48 | 3.48 | 10.00
0.5 | 0.10 | 1e-02 | 661 | 76 | 551 | 1.05 | 1.03 | 2.00
0.5 | 0.10 | 1e-06 | 2380 | 239 | 1760 | 1.19 | 1.13 | 2.00
0.5 | 0.01 | 1e-06 | 239268 | 2513 | 145055 | 1.62 | 1.58 | 2.00
0.1 | 0.10 | 1e-02 | 661 | 76 | 595 | 0.99 | 0.83 | 1.11
0.1 | 0.10 | 1e-06 | 2380 | 239 | 1901 | 1.11 | 0.91 | 1.11
0.1 | 0.01 | 1e-06 | 239268 | 2513 | 191853 | 1.23 | 1.03 | 1.11

###  A.3 Making Table 2

Given a value of $M$, $n$, and $\delta_{1}$, find the maximum relative error
between the unbiased and the biased estimate.

find_error_shifted_grid <- function(M, n , delta1) {

# create a shifted grid on [0, 1] with n elements

grid <- seq(0, 1 - 1 / n , by = 1 / n)

# Do the 0 part first

lower <- abs(1 -

M * mean(1 / qgamma(grid + delta1 / 2, M, 1)))

# Do the high part second

high <- abs(1 -

M * mean(1 / qgamma(grid + 1 / n - delta1 / 2, M, 1)))

return(max(lower, high))

}

find_error_shifted_grid_vector <- function(M, n , delta1)

apply(cbind(M, n, delta1), 1,

function (x) find_error_shifted_grid(

x[1], x[2], x[3]))

table2 <-

tibble(

M = c(10^4, 10^4, 10^4, 10^4, 10^5),

n = c(10^3, 10^4, 10^3, 10^2, 10^3),

delta1 = c(10^(-6), 10^(-6), 10^(-8), 10^(-8), 10^(-8)),

final = find_error_shifted_grid_vector(M, n, delta1)

)

table2 |>

kable(booktabs = TRUE, digits = 8) |> kable_styling()

M | n | delta1 | final
---|---|---|---
1e+04 | 1000 | 1e-06 | 0.00014967
1e+04 | 10000 | 1e-06 | 0.00010491
1e+04 | 1000 | 1e-08 | 0.00015871
1e+04 | 100 | 1e-08 | 0.00068990
1e+05 | 1000 | 1e-08 | 0.00002826