# Infinite discounted sum of betas

Let $$0 \leq \gamma < 1$$, $$X_i \sim \text{Beta}(\alpha, \beta)$$, and

$$Y \sim \sum_{i = 0}^\infty \gamma^i X_i$$

What is the distribution of $$Y$$? Does it have a closed form? Can it be sampled efficiently? We know that \begin{align} \operatorname{supp}(Y) &= \left[0, \frac{1}{1-\gamma}\right) \\ \operatorname{E}(Y) &= \frac{\operatorname{E}(X)}{1-\gamma} \\ \operatorname{Var}(Y) &= \frac{\operatorname{Var}(X)}{1-\gamma^2} \\ \end{align}

Also, note that \begin{align} Y = X + \gamma Y' \end{align} where $$Y$$ and $$Y'$$ have the same distribution. Therefore, denoting the characteristic function of a random variable $$A$$ as $$\phi_A$$, \begin{align} \phi_Y(t) &= \phi_{X + \gamma Y'}(t) \\ &= (\phi_X \cdot \phi_{\gamma Y'})(t) \\ &= \phi_X(t) \cdot \phi_{\gamma Y'}(t) \\ &= \phi_X(t) \cdot \phi_{Y'}(\gamma t) \\ &= \phi_X(t) \cdot \phi_Y(\gamma t) \\ \frac{\phi_Y(t)}{\phi_Y(\gamma t)} &= \phi_X(t) \\ &= {}_1F_1(\alpha; \alpha + \beta; i t) \end{align}

where the last RHS is the characteristic function of the beta distribution. Is there a way to solve for $$\phi_Y(t)$$ or at least sample efficiently from $$Y$$'s distribution?

• Is the $\gamma$ random too, or is it a fixed vector of values? – jkm Dec 3 '19 at 13:00
• @jkm It’s a fixed number. – user76284 Dec 3 '19 at 18:14
• This is another case of the situation considered by Vincent Granville at stats.stackexchange.com/questions/432136. The connection is this: let $\phi:\mathbb{R}_{+}\to\mathbb{R}_{+}$ be $\phi(x)=\sqrt{x+1}-1$ for him and $\phi(x)=\gamma x$ for you, so that in both cases $0$ is the sole attractor of $\phi.$ Take a sequence of iid positive random variables $X_1,X_2,X_3,\ldots$ and consider the distribution of $\phi(X_1+\phi(X_2+\phi(X_3+\cdots))).$ These two situations are thereby quite alike, indicating an analysis of one will suggest an effective analysis of the other. – whuber Dec 4 '19 at 4:55
• I noticed, thanks. I had also worked out the relationship between the cfs, but stopped there because it doesn't appear to get us any closer to a solution. You could pursue it a little further to deduce the cumulants of $Y$--but then, when you explore the possibility of writing a closed formula even for the cumulants of a Beta distribution, you realize that's not going to be useful either. – whuber Dec 4 '19 at 15:24

Note I've used slightly different notation to you.

Firstly, let's consider the underlying random variables we are dealing with. We have that all

$$X_{i}\sim\text{Beta}(\alpha,\beta)$$

where all random variables are independent of each other. If we now consider some scaled form of these random variables:

$$Y=a X$$

where we have that $$a<1$$ because $$\gamma<1$$. Now, for a continuous random variable

\begin{align} f_{Y}(y)&=\frac{f_{X}(y/a)}{a}\\ &=\frac{(y/a)^{\alpha-1}(1-y/a)^{\beta-1}}{a\cdot\text{B}(\alpha,\beta)} \end{align}

Notice that the support has changed from $$0 to $$0

If we define

$$\text{B}'(\alpha,\beta,a)=\int_{0}^{a}t^{\alpha-1}(a-t)^{\beta-1}\,dt$$

as a variant on the usual Beta function (this function equals the Beta function for $$a=1$$), then we can show that

$$\text{B}'(\alpha,\beta,a)=a^{\alpha+\beta-1}\cdot\text{B}(\alpha,\beta)$$

This means that

\begin{align} f_{Y}(y)&=\frac{(y/a)^{\alpha-1}(1-y/a)^{\beta-1}}{a\cdot\text{B}(\alpha,\beta)}\\ &=\frac{y^{\alpha-1}(a-y)^{\beta-1}}{\text{B}'(\alpha,\beta,a)}\cdot\frac{\text{B}'(\alpha,\beta,a)}{a^{\alpha+\beta-1}\cdot\text{B}(\alpha,\beta)}\\ &=\frac{y^{\alpha-1}(a-y)^{\beta-1}}{\text{B}'(\alpha,\beta,a)} \end{align}

For the purposes of this answer, let's denote this random variable as

$$Y\sim\text{Beta}'(\alpha,\beta,a)$$ This is essentially a scaled Beta distribution, where the support has been limited to $$[0,a]$$. The important point here is that random values generated from this random variable will be bound between $$[0,a]$$.

Now, in your summation, each Beta distribution $$X_{i}$$ is scaled by $$\gamma^{i}$$ for increasing $$i$$. This means that each summand will have an increasingly reduced support $$[0,\gamma^{i}]$$ and each subsequent summand's probability mass will be increasingly narrow.

We know this series will converge for $$0\leq\gamma<1$$ because, for each $$i$$, the maximum value that each subsequent random value can take is $$\gamma^{i}$$. So, at most, the series is

$$1+\gamma+\gamma^{2}+\gamma^{3}+\gamma^{4}+\ldots$$

which is a simple geometric series converging to

$$\frac{1}{1-\gamma}$$

Interestingly, this tells us the upper bound of the support of $$Y$$ for a given $$\gamma$$. If we define

$$Z=\sum_{i=0}^{\infty}Y_{i}$$

then the support for $$Z$$ would then be $$[0,1/(1-\gamma)]$$.

Now, you can derive the moment generating function of $$Y$$ as

$$M_{Y_{i}}(t)=1+\sum_{k=1}^{\infty}\Bigg(\prod_{r=0}^{k-1}\bigg(\frac{\alpha+r}{\alpha+\beta+r}\bigg)\frac{(\gamma^{i}\cdot t)^{k}}{k!}\Bigg)$$

Then the moment generating function of $$Z$$ is

\begin{align} M_{Z}(t)&=\prod_{i=0}^{\infty}M_{Y_{i}}(t)\\ &=\Bigg[1+\sum_{k=1}^{\infty}\Bigg(\prod_{r=0}^{k-1}\bigg(\frac{\alpha+r}{\alpha+\beta+r}\bigg)\frac{t^{k}}{k!}\Bigg)\Bigg]\cdot \\ &\Bigg[1+\sum_{k=1}^{\infty}\Bigg(\prod_{r=0}^{k-1}\bigg(\frac{\alpha+r}{\alpha+\beta+r}\bigg)\frac{(\gamma\cdot t)^{k}}{k!}\Bigg)\Bigg]\cdots\\ \end{align}

I'm not sure how useful this is. But it's useful to note that for increasing $$i$$ the subsequent terms will decay in size so that $$M_{Y_{i}}(t)\rightarrow 1$$ for large $$i$$.

If we look at some simulations of this distribution, there are some interesting results. First we'll define some functions for $$Y$$.

dsbeta=function(x, alpha, beta, gamma) {
x[x<0]=0
x[x>gamma]=gamma
return(x^(alpha-1)*(gamma-x)^(beta-1)/(gamma^(alpha+beta-1)*beta(alpha, beta)))
}

psbeta=function(q, alpha, beta, gamma) {
q[q<0]=0
q[q>gamma]=gamma
return(pbeta(q/gamma, alpha, beta))
}

qsbeta=function(p, alpha, beta, gamma) {
p[p<0]=0
p[p>1]=1
return(gamma*qbeta(p, alpha, beta))
}

rsbeta=function(n, alpha, beta, gamma) {
return(gamma*rbeta(n, alpha, beta))
}

#Perform some simulations:
a=3
b=2
gamma=0.2 #Scaling of each X
N=10000   #Number of samples of Z
M=100     #Number of Y in each Z

#Generate random values of Z:
r=as.vector(matrix(rbeta(N*M, a, b), N, M) %*% gamma^seq(0,M-1))

#Final simulations:
hist(r, freq = FALSE, breaks = 50)

#Perform fit (here we provide the upper bound gamma as the theoretical upper limit):
f=fitdistrplus::fitdist(r, "sbeta", start = list(alpha = 2, beta = 3), fix.arg = list(gamma = 1/(1-gamma)))

#Visualize the fit:
x=seq(min(r), max(r), length.out = 1000)
lines(x, dsbeta(x, f$$estimate[1], f$$estimate[2], 1/(1-gamma)), lwd = 2, col = "blue")

#GOF tests:
ks.test(r, "psbeta", f$$estimate[1], f$$estimate[2], f$$fix.arg$$gamma)
goftest::cvm.test(r, "psbeta", f$$estimate[1], f$$estimate[2], f$$fix.arg$$gamma)$$p.value goftest::ad.test(r, "psbeta", f$$estimate[1], f$$estimate[2], f$$fix.arg$$gamma)$$p.value

#qq plot:
plot(sort(r), qsbeta(ppoints(length(r)), f$$estimate[1], f$$estimate[2], f$$fix.arg$$gamma), ylab = "Theoretical quantiles", xlab = "Empirical quantiles", col = "blue", type = "l", lwd = 2)
abline(0, 1, lty = 2, lwd = 2, col = "red")


It seems from the simulations that $$Z$$ can be well approximated by the scaled beta

$$Z\sim\text{Beta}'(\cdot,\cdot,1/(1-\gamma))$$

but as per @whuber's comments, for small $$\alpha,\beta$$ and $$\gamma\rightarrow 0$$ this breaks down. I think it has something to do with the bi-modality for $$\alpha,\beta<1$$ and the fewer 'non-degenerate' random variables due to $$\gamma$$ small.

• I think you will get more insight into the underlying issues by considering cases where at least one of $\alpha,\beta$ is close to zero (preferably both). This will more clearly show the limitations of approximating $Z$ as a scaled Beta distribution. Focus on relatively small $\gamma,$ too: for large $\gamma \approx 1$ there's a CLT effect. – whuber Dec 5 '19 at 14:15
• @whuber I'm only seeing a poor scaled Beta fit for $\alpha\approx\beta\approx 1$ and $\gamma\rightarrow 0$ where the simulated data looks quite uniform. – Ed P Dec 5 '19 at 22:28
• @whuber I wonder whether it can be proven that it approaches a scaled beta as $\gamma \rightarrow 1$ (or as $\alpha, \beta \rightarrow \infty$). – user76284 Dec 8 '19 at 20:15
• Under those conditions, properly standardized, the sum ought to approach a Normal distribution. – whuber Dec 8 '19 at 21:48