Transition probabilities in multi dimensional birth-death process

Question

I've got an urn problem that I believe can be nicely modeled as a birth-death process (I'm very new to markov chains, so maybe this is simply the wrong approach). Suppose we have $k$ urns each of which can hold a finite number of $n$ distinct balls labelled $1, \cdots, n$. The number of urns can increases without bound and there are no empty urns allowed, however, the state where there are no urns ($k = 0$) is allowed.

1. With some probability $q_k$ choose an urn uniformly at random

   1.1 With probability $1/2$ add a ball whose label is not already in the urn, otherwise remove a ball.

   • if the urn is full, do nothing.
   • if we remove the last ball from the urn, remove the urn all together and $k \rightarrow k-1$.

2. With probability $1-q_k$, add a new urn with a single ball chosen uniformly at random. ($k \rightarrow k+1$)

If $K$ is a random variable representing the number of existing urns and $U_i$ a random variable representing the number of balls currently in urn $i$ (number of balls in each urn is assumed to be independent), the goal is for

$$P\left(K = k\right) = \frac{1}{k!}e^{-1}$$ and

$$ P(U) = \prod\limits_{i=1}^k P(U_i = n_i) = \frac{1}{\prod\limits_{i=1}^k{n \choose n_i}} = \prod\limits_{i=1}^k\frac{n_i!\left(n-n_i\right)!}{n!} $$

where $1\leq n_i\leq n$ is the number of balls in urn $i$.

i.e. the number of balls in the urns should be uniformly distributed.

My goal is to get an expression for $q_k$.

I have a solution for the case where urns can only hold one ball:

This what I believe the markov chain looks like for this case, where each number in the circle represents the number of urns.

Let $\alpha_i = \pi_{i,i+1}$ and $\beta_i = \pi_{i,i-1}$, where p_i = P(K=i)

Detailed balance requires that: $p_i\alpha_i = p_{i+1}\beta_{i+1}$ for $\sum_i p_i = 1$ and $p_j = \sum_i p_i\pi_{ij}$

When $i = 0$ $$ p_0 = \beta_1p_1 +(1-\alpha_0)p_0 $$

and in general

$$ p_i = \alpha_{i-1}p_{i-1} + \beta_{i+1}p_{i+1} + (1-\alpha_i-\beta_i)p_i \tag{1} $$

the base case is such that $$ \alpha_0p_0 = \beta_1p_1 $$ assuming $\alpha_{i-1}p_{i-1} = \beta_ip_i$ then we can substitute this into Equation (1) and get

$$ \alpha_i p_i = \beta_{i+1}p_{i+1} $$

meaning the detailed balance can be satisfied pairwise

now suppose $\alpha_i = \alpha$ and $\beta_{i+1} = (i+1)\beta$ for all $i$, then

$$ p_1 = \frac{\alpha}{\beta}p_0 $$ and we can do induction on $i$ and get

$$ p_i = \frac{\left(\frac{\alpha}{\beta}\right)^i}{i!}p_0 $$

normalizing this, implies $p_0 = \exp{\left(-\frac{\alpha}{\beta}\right)}$ and letting $\frac{\alpha}{\beta} = 1$

$$ p_i = \frac{1}{i!}e^{-1} $$

which is the distribution I want, assuming each urn can only hold $1$ ball and I can set $q_i = \alpha/(i+1)$ for any $\alpha < 1/2$.

I'm having difficulty, however, following the line of thinking above to generalize this to urns that can hold more than $1$ ball.

I have a solution where I made a rather ad-hoc order $1/n$ approximation and get that $q_i = \frac{1}{2n(i+1)}$ where $i$ is the current number of urns and $n$ is the maximum number of balls an urn can hold (I verfied with some simulations that this does indeed sample from the distributions that I am looking for but I will spare the derivation to avoid making this post longer than it already is). But I feel like I there should be a better way than this approximation I made to solve this problem, something similar to what I did for the $n =1$ case. Is there a way to solve this problem following such a line of thinking?

Answer 1

Firstly, I think the derivation in the 1D case doesn't match with the description of the problem given as it says that $q_i$ decreases with $i$ which means you are likely to keep on adding urns which doesn't make sense with steady state probabilities decreasing with $i$. Assuming you have a solution for the 1D case for the given description, the following might work

One possibility is to change your state space to also include n. Specifically, if you take the state space to be $\mathcal{S} = \mathbb{N}^n$ where the markov chain being in a state $\vec{x} = (x_1 , \ldots , x_n )$ represents that currently you have $x_i$ urns with $i$ balls (I am guessing this is the multi in the multi-dimensional asked in the title). Now you can define the transitions in a more refined way.

Let $S(\vec{x})$ denote the sum of all the components of $\vec{x}$. Then from what I understand in the question, we require that the chain be irreducible and the steady state probabilities satisfy $$\Sigma_{\{\vec{x} | S(\vec{x})= k\}} \pi_{\vec{x}} = \frac{1}{k!} e^{-1} = \pi_k $$

From state $\vec{x}$ and denoting $S(\vec{x})$ by $k$, with probability $1-q_k$ we add a new urn $\textit{i.e}$ $$(x_1, \ldots, x_n) \xrightarrow{1-q_k} (x_1+1, \ldots, x_n)$$

Next in the remaining probability $q_k$ we choose an urn uniformly at random and either add or remove a ball at random. This implies that $$ \vec{x} \longrightarrow \begin{cases} (x_1-1, x_2, \ldots, x_n) & \textrm{w.p } (x_1/2k).q_k \\ (x_1-1, x_2+1, \ldots, x_n) & \textrm{w.p } (x_1/2k).q_k \\ (x_1, x_2, \ldots, x_{i-1}+1,x_{i}-1, \ldots, x_n) & \textrm{w.p } (x_i/2k).q_k \; \forall \, 2\leqslant i \leqslant n \\ (x_1, x_2, \ldots, x_{i}-1,x_{i+1}+1, \ldots, x_n) & \textrm{w.p } (x_i/2k).q_k \; \forall \, 2\leqslant i < n \\ \vec{x} & \textrm{w.p } (x_n/2k).q_k \end{cases} $$

It should be clear that the chain defined this way is irreducible and aperiodic over $\mathcal{S}$. Now defining $$\pi_{\vec{x}} = \frac{\binom{k}{x_1, \ldots, x_n}}{n^k} \pi_k$$ I think will satisfy the detailed balance equations