I found a nice problem recently, but could not come up with a solution:
Find all functions $f:\mathbb{Z}_{\geq 0}\to\mathbb{Z}_{\geq 0}$ such that for all $0< \lambda < 1$, if $X\sim G(\lambda)$, then there exists $0<\mu < 1$ with $f(X) \sim G(\mu)$.
Thanks to the answers below, I was able to understand how to solve it!
Spoiler: the solutions are all the functions $f_d:x\mapsto \lfloor x/d\rfloor$ for a fixed $d\geq 1$.
Let $f$ be a function such that for each $\lambda\in (0,1)$, if $X\sim\mathcal G(\lambda)$, then there exists $\mu\in (0,1)$ such that $f(X)\sim\mathcal G(\mu)$.
First of all, $f$ should satisfy $f\left(\mathbb Z_{\geqslant 0}\right)=\mathbb Z_{\geqslant 0}$ because $\mathbb P(f(X)=k)=\mathbb P(f^{-1}(\{k\}))$ is positive for each $k$ and $X$ satisfying a geometric distribution.
For $j,k\in\mathbb Z_{\geqslant 0}$, define $a_{k,j}=\begin{cases}1&\mbox{if }f(j)=k\\ 0&\mbox{otherwise} \end{cases}$, $a_{k,-1}=0$ and $I_k=f^{-1}(\{k\})$.
We observe that if $X\sim\mathcal G(\lambda)$, then $f(X)\sim\mathcal G(\mu(\lambda))$ where $$ \mu(\lambda)=\sum_{j\in I_0}(1-\lambda)^j\lambda. $$ Indeed, $\mu(\lambda)=\mathbb P(f(X)=0)=\mathbb P(X\in I_0)$.
It will be more convenient to write $\mu(1-\lambda)$ as a power series, namely, $$ \mu(1-\lambda)=\sum_{j\in I_0}\lambda^j-\sum_{j\in I_0}\lambda^{j+1} =\sum_{j=0}^{\infty}\left(a_{0,j}-a_{0,j-1}\right)\lambda^j. $$
By definition of geometric distribution with parameter $\mu(1-\lambda)$, we derive that for each $k\geqslant 0$ and each $\lambda\in (0,1)$,
$$\tag{*}
\left(1-\sum_{j=0}^{\infty}\left(a_{0,j}-a_{0,j-1}\right)\lambda^j\right)^k
\sum_{j=0}^{\infty}\left(a_{0,j}-a_{0,j-1}\right)\lambda^j
=\sum_{j=0}^{\infty}\left(a_{k,j}-a_{k,j-1}\right)\lambda^j.
$$
Note that all the involved series are convergent since $0\leqslant a_{k,j}\leqslant 1$ and $0<\lambda<1$.
Let us look first at the constant coefficient, which can be obtained by letting $\lambda\to 0$. We get that for each $k\geqslant 0$, $$ \left(1-a_{0,0}\right)^ka_{0,0}=a_{k,0} $$ If $a_{0,0}=0$, we would get that $a_{k,0}=0$ for each $k\geqslant 0$ but since the sets $I_k$ are disjoint and their union is $\mathbb Z_{\geqslant 0}$, we should have $\sum_{k\geqslant 0}a_{k,0}=1$ hence we get $a_{0,0}=1$ hence $f(0)=0$.
Let us now update $(*)$. Separating the index $j=0$ in the two sums of the left hand side, we derive that for each $k\geqslant 1$ and each $\lambda\in (0,1)$,
$$\tag{**}
\left(-\sum_{j=1}^{\infty}\left(a_{0,j}-a_{0,j-1}\right)\lambda^j\right)^k
\left(1+\sum_{j=1}^{\infty}\left(a_{0,j}-a_{0,j-1}\right)\lambda^j\right)
=\sum_{j=0}^{\infty}\left(a_{k,j}-a_{k,j-1}\right)\lambda^j.
$$
Notice that in the left hand side of $(**)$, the lowest degree in $\lambda$ is $k$ hence we should have for $0\leqslant j\leqslant k-1$ that $a_{k,j}-a_{k,j-1}=0$, and since $a_{j,-1}=0$, $a_{k,0}=a_{k,1}=\dots=a_{k,k-1}=0$ or in other words, $j<k$ implies that $f(j)\neq k$ hence $f(j)\leqslant j$.
Looking now at the coefficient in front of $\lambda^k$, we get that $(1-a_{0,1})^k=a_{k,k}$. If $f(1)=1$, then $f(k)=k$ for all $k$.
Suppose that $f(1)\neq 1$. By 6, we have $f(1)=0$. Let $d$ be the smallest natural number such that $f(d)\neq 0$, that is, $a_{0,d}=0$ and $a_{0,j }=1$ for $j<d$. Looking now at the coefficient in front of $\lambda^j$, $0\leqslant j\leqslant j_0k$, we find that $f(dk)=k$.
Define $b_{k,j}=a_{k,j}-a_{k,j-1}$ and $c_j$ by $c_0=0$ and $c_j=-b_{0,j}$ for $j \geqslant 1$. Since $-\sum_{j=1}^{\infty}\left(a_{0,j}-a_{0,j-1}\right)\lambda^j=\sum_{j=0}^\infty c_j\lambda^j$, for $k=1$ in $(**)$, recognizing a Cauchy product between power series, we get that \begin{align} b_{1,j}&=\sum_{i=0}^j b_{0,i}c_{j-i}\\ &=-\sum_{i=0}^{j-1} b_{0,i}b_{0,j-i}\\ &=-b_{0,j}-\sum_{i=1}^{j-1} b_{0,i}b_{0,j-i}\\ &=-b_{0,j}-\sum_{i=1}^{j-1} b_{0,i}\mathbf{1}_{i\geqslant d}b_{0,j-i}\mathbf{1}_{j-i\geqslant d} \end{align} where we used in the last equality that $b_{0,i}=0$ if $1\leqslant i<d$. We thus got $$\tag{1} b_{1,j}=-b_{0,j}-\mathbf{1}_{j\geqslant 2d}\sum_{i=d}^{j-d} b_{0,i}b_{0,j-i} $$
We will now prove by induction on $\ell\geqslant 1$ that $b_{0,d+\ell}=0$.
Since $$ b_{1,2d+1}+b_{0,2d+1}=\mathbf{1}_{f(2d+1)\in\{0,1\}}-\mathbf{1}_{f(2d)\in\{0,1\}}\in \{-1,0,1\}, $$ this forces $b_{0,d+1}$ to be equal to $0$.
Going back to (**), we can drastically simplify the left hand side by $$ \left(-\sum_{j=1}^{d}\left(a_{0,j}-a_{0,j-1}\right)\lambda^j\right)^k \left(1+\sum_{j=1}^{d}\left(a_{0,j}-a_{0,j-1}\right)\lambda^j\right) =\sum_{j=0}^{\infty}\left(a_{k,j}-a_{k,j-1}\right)\lambda^j. $$ and since $a_{0,j}-a_{j-1}=0$ for $1\leqslant j\leqslant d-1$, we get $$ \lambda^{dk} \left(1-\lambda^d\right) =\sum_{j=0}^{\infty}\left(a_{k,j}-a_{k,j-1}\right)\lambda^j $$ hence identifying the coefficient in $dk$ and $d(k+1)$ gives that $f(k)=\lfloor k/d\rfloor$, where $\lfloor x\rfloor$ is the unique integer for which $\lfloor x\rfloor\leqslant x<\lfloor x\rfloor+1$.
It remains to check that any function of the form $f(k)=\lfloor k/d\rfloor $ for $d\geqslant 1$ does the job. First, $I_0=\{0,\dots,d-1\}$ hence $\mu(\lambda)=1-(1-\lambda)^{d}$ and if $\ell$ is a non-negative integer, \begin{align} \mathbb P(f(X)=\ell)&=\mathbb P(d\ell\leqslant X\leqslant d\ell+d-1)\\ &=\sum_{j=d\ell}^{d\ell+d-1}(1-\lambda)^j\lambda&\\ &=(1-\lambda)^{d\ell}\sum_{j=0}^{ d-1}(1-\lambda)^j\lambda&\\ &=(1-\lambda)^{d\ell}\left(1-(1-\lambda)^{d}\right) \\ &=(1-\mu(\lambda))^{ \ell}\mu(\lambda)&. \end{align}
Below are only proofs of the known results from the question, which is a characterization of all non-decreasing $f$.
For $\lambda\in[0,1]$ let $G(\lambda)$ denote the geometric distribution, that is the distribution of the number of failures before the first success in a Bernoulli trial given by $$\mathbb P(X_\lambda=k)=(1-\lambda)^{k}\lambda$$ for $k\in\mathbb Z_{\ge 0}$, where $X_\lambda\sim G(\lambda)$. Let $\mathcal F=\{f:\mathbb Z_{\ge 0}\rightarrow\mathbb Z_{\ge 0}: \forall\lambda\in(0,1)\exists\mu\in(0,1)\,f(X_\lambda)\sim G(\mu)\}$. For $f\in\mathcal F$ and $\lambda\in(0,1)$ let $\mu_\lambda\in(0,1)$ be such that $f(X_\lambda)\sim G(\mu_\lambda)$. Notice that $\mu_\lambda=\mathbb P(f(X_\lambda)=0)=\sum_{k\in f^{-1}(0)}(1-\lambda)^k\lambda$ is determined by $f$ and $\lambda$.
Claim 1: Each $f\in\mathcal F$ is surjective.
Proof: For $k\in\mathbb Z_{\ge 0}$ we have $\mathbb P(X_\lambda\in f^{-1}(k))=\mathbb P(f(X_\lambda)=k)=\mathbb P(X_{\mu_\lambda}=k)>0$ and hence $f^{-1}(k)\neq\emptyset$.
Claim 2: For $f\in\mathcal F$ we have $f(0)=0$ and $\mu_\lambda\ge\lambda$, where $f(X_\lambda)\sim G(\mu_\lambda)$ for $\lambda\in(0,1)$.
Proof: We have $(1-\mu_\lambda)^{f(0)}\mu_\lambda=\mathbb P(f(X_\lambda)=f(0))\ge\mathbb P(X_\lambda=0)=\lambda$. For $a\in\mathbb R_{\ge 0}$ and $g:[0,1]\rightarrow[0,1]$, $x\mapsto (1-x)^ax$ we have $g'(x)=(1-x)^a-a(1-x)^{a-1}x$ and thus a unique maximum at $x=\frac{1}{a+1}$, yielding $(1-\frac{1}{f(0)+1})^{f(0)}\frac{1}{f(0)+1}\ge (1-\mu_\lambda)^{f(0)}\mu_\lambda\ge\lambda$. Since this holds for all $\lambda\in(0,1)$, we have $1\ge(1-\frac{1}{f(0)+1})^{f(0)}\ge f(0)+1$, thereby $f(0)\le 0$ and hence $f(0)=0$. Substituting this in the first inequality gives $\mu_\lambda\ge\lambda$.
Claim 3: Let $d=\inf\{k\in\mathbb Z_{\ge 0}:f(k)=1\}$ for $f\in\mathcal F$. Then we have $d\in\mathbb Z_{>0}$ and $\lim_{\lambda\rightarrow 1}\frac{1-\mu_\lambda}{(1-\lambda)^d}=1$.
Proof: We have $d<\infty$ by Claim 1 and $d>0$ by Claim 2, so $d\in\mathbb Z_{>0}$. Further, we have $(1-\mu_\lambda)\mu_\lambda=\mathbb P(f(X_\lambda)=1)=\mathbb P(X_\lambda\in f^{-1}(1))=(1-\lambda)^d\lambda\sum_{k\in f^{-1}(1)}(1-\lambda)^{k-d}$, which yields $$ \frac{1-\mu_\lambda}{(1-\lambda)^d}=\frac{\lambda}{\mu_\lambda}\sum_{k\in f^{-1}(1)}(1-\lambda)^{k-d}. $$ Using the geometric series we have $1\le\sum_{k\in f^{-1}(1)}(1-\lambda)^{k-d}\le\frac{1}{\lambda}$, and further $\lambda\le\mu_\lambda\le 1$ using Claim 2, so we indeed obtain $\lim_{\lambda\rightarrow 1}\frac{1-\mu_\lambda}{(1-\lambda)^d}=1$.
Claim 4: For $f\in\mathcal F$ and $n\in\mathbb Z_{\ge 0}$ we have $dn=\min\{k\in\mathbb Z_{\ge 0}:f(k)=n\}$.
Proof: Notice that $m(n)=\min\{k\in\mathbb Z_{\ge 0}:f(k)=n\}$ is well-defined by Claim 1, that we have $m(0)=0$ by Claim 2, and that $m(1)=d$ trivially holds, so let $n\in\mathbb Z_{>1}$. Notice that \begin{align*} (1-\mu_\lambda)^{n}\mu_\lambda &=\mathbb P(f(X_\lambda)=n)=\mathbb P(X_\lambda\in f^{-1}(n))=\sum_{k\in f^{-1}(n)}(1-\lambda)^k\lambda\\ &=(1-\lambda)^{m(n)}\lambda\sum_{k\in f^{-1}(n)}(1-\lambda)^{k-m(n)}. \end{align*} As before, we have $1\le\sum_{k\in f^{-1}(n)}(1-\lambda)^{k-m(n)}\le\frac{1}{\lambda}$, so with Claim 2 and Claim 3 we have $\lim_{\lambda\rightarrow 1}(1-\lambda)^{dn-m(n)}=1$ since \begin{align*} (1-\lambda)^{dn-m(n)} &=\left(\frac{(1-\lambda)^d}{1-\mu_\lambda}\right)^n \frac{\lambda\sum_{k\in f^{-1}(n)}(1-\lambda)^{k-m(n)}}{\mu_\lambda}\frac{(1-\mu_\lambda)^n\mu_\lambda}{(1-\lambda)^{m(n)}\lambda\sum_{k\in f^{-1}(n)}(1-\lambda)^{k-m(n)}}\\ &=\left(\frac{(1-\lambda)^d}{1-\mu_\lambda}\right)^n \frac{\lambda\sum_{k\in f^{-1}(n)}(1-\lambda)^{k-m(n)}}{\mu_\lambda}. \end{align*} Since $dn-m(n)$ does not depend on $\lambda$, this implies that $dn=m(n)$.
Now, we restrict to non-decreasing functions. Let $\mathcal F_+=\{f\in\mathcal F:\forall n\le n'\,f(n)\le f(n')\}$.
Claim 5: For $d\in\mathbb Z_{>0}$ let $f_d:\mathbb Z_{\ge 0}\rightarrow\mathbb Z_{\ge 0}$, $k\mapsto\lfloor\frac{k}{d}\rfloor$. We have $\mathcal F_+=\{f_d:d\in\mathbb Z_{>0}\}$.
Proof: Let $f\in\mathcal F_+$ and $d\in\mathbb Z_{>0}$ as defined in Claim 3. Since $f$ is non-decreasing, Claim 4 directly implies $f=f_d$. Conversely, for $d\in\mathbb Z_{>0}$ the map $f_d$ is clearly non-decreasing. We further have $\mu_\lambda=\sum_{k=0}^{d-1}(1-\lambda)^k\lambda=1-(1-\lambda)^d$. Further, we have $\mathbb P(f(X_\lambda)=n)=\sum_{k=dn}^{d(n+1)-1}(1-\lambda)^k\lambda=(1-\lambda)^{dn}\mu_\lambda=(1-\mu_\lambda)^n\mu_\lambda$ and thus $f(X_\lambda)\sim G(\mu_\lambda)$.
