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<x<1$$
to
$$0<y<a$$
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.
