# How to simulate the following joint distribution?

I am working with a multi-armed bandit problem in which there are finitely many arms, say $$K$$ of them. Here, $$K\geq 2$$ (for the specific problem I am looking at, each arm is a homogeneous and ergodic Markov chain on a finite state space $$S$$. The state space is common across all the arms).

At each time instant $$t=0,1,2,\ldots$$, only one of the arms may be selected. A state evolution is then observed on the selected arm. The unobserved arms also undergo state evolution, and hence there is a certain 'delay' associated with observing each arm. For any given arm, its delay is a positive integer that denotes how long it has been since it was last observed. For example, if the delay of arm $$1$$ is $$5$$, it means that arm $$1$$ was selected previously $$5$$ time instants ago (with reference to the current time), and has not been selected since then. The smallest value of delay for any arm can be $$1$$, signifying that this particular arm was sampled in the immediately previous time instant.

Suppose that $$\lambda(a,d),\quad a=1,2,\ldots,K,\quad d\geq 1$$ denotes a joint probability distribution on arms and delays. At any given time $$t\in\{0,1,2,\ldots\}$$, only the arm which has been selected at $$t-1$$ will have delay of $$1$$, whereas the delays of the remaining arms increment from their previous value by $$1$$.

I am interested in simulating the above setup on MATLAB. To begin with, I select each of the $$K$$ arms once. After that, I need to carefully select the next arm based on the information of the arm delays. I would like my simulation to finally reflect that the empirical distribution of arms and delays matches with the above $$\lambda(a,d)$$ with increasing number of iterations.

I am stuck at what probability distribution I should use to choose an arm at any given time, based on the arm selection and delay information available currently.

It would be helpful if anyone can give me some ideas.

PS: I will share here some of the approaches I have already tried.

Approach 1: Suppose $$d_1(t),d_2(t),\ldots,d_K(t)$$ denote the arm delays at time $$t$$ (among these, exactly one of them is $$1$$). Then, I tried selecting the arm at time $$t+1$$ according to the distribution $$\frac{\lambda(a,d_a(t))}{\sum\limits_{b=1}^{K}\lambda(b,d_b(t))}.$$ However, this does not lead to the empirical distribution of arms and delays to match with $$\lambda(a,d)$$ in the long run.

Approach 2: I also tried the following approach: Select arm at time $$t+1$$ according to the formula $$\text{arm at time }(t+1)=\arg\min\limits_{a\in\{1,2,\ldots,K\}}||\lambda_t(a,\cdot)-\lambda(a,\cdot)||_{\infty},$$ where $$\lambda_t(a,d)$$ denotes the empirical probability (at time $$t$$) of arm $$a$$ being selected when its delay is $$d$$. This too does not lead to convergence of empirical distribution to $$\lambda(a,d)$$ in the long run.

• "Suppose that λ(a,d),a=1,2,…,K,d≥1 denotes a joint probability distribution on arms and delays." I have hard time deciphering the meaning of this notation. As written, it means that $\lambda$ (which, apparently, has the meaning of probability of something) depends on the arm $a$ and the delay $d$ only but you call it a "joint probability distribution on arms and delays", which IMO means that elementary events are described not by arm-delay pairs (a,d) but by vectors $(d_1,\dots,d_a)$ (arm delays). Perhaps you really mean just individual distributions for each arm separately? Can you clarify? Commented Aug 6, 2019 at 3:56
• The matrix entries $\lambda(a,d)$ are to be viewed as the probability of selecting arm $a$ when its delay is $d$, i.e., the probability that arm $a$ was selected at the current time and that it was not selected for $d$ time instances before. However, when any arm is selected, its delay immediately gets set to $1$ and the delays of all the other arms get incremented by one. Commented Aug 6, 2019 at 16:06
• You are right that the delays $d_1,\ldots,d_K$ are what matter when trying to select an arm at any given time based on the information of what arms have been selected in the past. Commented Aug 6, 2019 at 16:07
• I want to figure out a way of selecting arms sequentially, one at each time instant, using the information of delays available to me. All the approaches I've described above are not leading me to where I'd like to get to. Commented Aug 6, 2019 at 16:10
• "The matrix entries λ(a,d) are to be viewed as the probability of selecting arm a when its delay is d, i.e., the probability that arm a was selected at the current time and that it was not selected for d time instances before" Eh? As written, those are different probabilities: the first one is conditional and the second isn't. Which one do you mean? Commented Aug 6, 2019 at 17:09

(Skip ahead to Solution starts here to skip the buildup.)

One thing to keep in mind is that the joint probability distribution $$\lambda$$ needs to satisfy some constraints apart from the usual conditions on probability distributions.

The distribution of delays for a given arm $$P(d|a)$$ will directly affect $$P(a)$$. To see this intuitively, note that for a given arm $$a_i$$, if $$\lambda(a_i, 1)$$ is non-zero and $$\lambda(a_i, j) = 0$$ for all $$j>1$$, then $$P(a_i)$$ is can only be either 0 or 1. Either $$a_i$$ can never be picked, or must be picked at every step after the first time it is picked.

The actual condition on $$\lambda$$... (I'll assume that every arm has a non-zero probability of being picked, i.e., that each arm $$a$$ has at least one delay $$d$$ such that $$\lambda(a, d) > 0$$)

Expected delay between two instances of picking an arm $$a_i$$ is given by: $$\langle d\rangle_{a_i} = \frac{\sum\limits_{d}~d~\lambda(a_i, d)}{\sum\limits_{d}~\lambda(a_i, d)}$$ The expected number of times arm $$a_i$$ appears in the long run = Number of steps / $$\langle d\rangle_{a_i}$$. Or $$P(a_i) = 1/\langle d\rangle_{a_i}$$. This can be written as $$\sum\limits_{d}~\lambda(a_i, d) = \frac{\sum\limits_{d}~\lambda(a_i, d)}{\sum\limits_{d}~d~\lambda(a_i, d)} \Rightarrow \sum\limits_{d}~d~\lambda(a_i, d) = 1$$

So, $$\lambda(a, d)$$ needs to satisfy the following two conditions: $$\sum\limits_{a, d} \lambda(a, d) = 1$$ $$\sum\limits_{d}~d~\lambda(a, d) = 1~~~~ \forall a$$

I can't see just yet whether or not either of your approaches will work if these constraints are satisfied by $$\lambda$$. I'll get back to you soon...

EDIT: The following should work to simulate arm selection:

The intuition for what follows is to look the current state as $$(D_{a_1}, \dots, D_{a_K})$$ which are the times it has been since $$(a_1, \dots,a_K)$$ were last picked. At each step, the system can be thought of as moving from one state to another.

The probability distribution of delay $$d$$ after which a given arm $$a$$ will be chosen is given by $$p(d|a) = \frac{\lambda(a, d)}{\sum\limits_{d'} \lambda (a, d')}$$ Let $$\eta(a, D)$$ be the probability that asymptotically, at a randomly chosen time (before making the selection for that step), it has been $$D$$ steps since $$a$$ was last chosen. It is given by (using the form of $$p(d|a)$$) $$\eta(a, D) = \left[\sum\limits_{d=1}^{\infty}\lambda(a, d)\right]\times \left[1 - \frac{\sum\limits_{d=1}^{D-1}\lambda(a, d)}{\sum\limits_{d=1}^{\infty}\lambda(a, d)}\right] = \sum\limits_{d=1}^{\infty}\lambda(a, d) - \sum\limits_{d=1}^{D-1}\lambda(a, d)$$ You can verify that the condition on $$\lambda$$ mentioned above ensures that $$\sum_D \eta(a, D) = 1$$ for all $$a$$. To clarify how $$D$$ is counted...

$$\eta(a, 1)$$ is the probability that $$a$$ was picked in the previous step

$$\eta(a, 2)$$ is the probability that $$a$$ was picked two steps ago, and wasn't picked in the last step, and so on.

If $$\eta(a, D)$$ is satisfied by the selection process for each $$a$$, then $$\lambda(a, d)$$ will be satisfied.

Solution starts here

Let $$\mathbb{P}(D_{a_1},\dots, D_{a_K})$$ be the probability of a given state of the system (asymptotically at a randomly chosen time) represented by the number of steps $$D_{a_i}$$ it has been since each arm $$a_i$$ was chosen. There are many distributions, $$\mathbb{P}$$, that'll work to reproduce $$\eta$$, but we can pick

$$\mathbb{P}(D_1, \dots, D_K) = \prod\limits_{a_i}\eta(a_i, D_{a_i}) = \prod\limits_{a_i}\left[\sum\limits_{d=1}^{\infty}\lambda(a_i, d) - \sum\limits_{d=1}^{D_{a_i}-1}\lambda(a_i, d)\right]$$

To simulate the next arm selection, see what state each of the $$K$$ selections will lead to ($$D$$ will become 1 for the chosen arm, and will increment by 1 for all others arms). At each step, choose arm $$a$$ with probability $$\mathrm{Prob}(a) = \frac{\mathbb{P}(\text{state resulting from choosing arm } a)}{\sum\limits_b{\mathbb{P}(\text{state resulting from choosing arm }b)}}$$ with $$\mathbb{P}$$ given above.

How to make the computation of $$\mathrm{Prob}(a)$$ computationally tractable:

If a closed form doesn't exist for the infinite sum, for each $$a_i$$ we can pre-compute $$\Lambda^*(a_i) = \sum\limits_{d=1}^{N}\lambda(a_i, d)$$ for some large $$N$$. During each generative step, update this by adding some more terms to the sum ($$O(K)$$ operations in each step, since there are $$K$$ such sums to be updated). In the asymptotic limit, $$\Lambda^*(a_i)$$ will become a good approximation to the infinite sum.

Similarly, at each step keep a table of the sum upto $$D_{a_i}-1$$ for each $$a_i$$. The possible changes to this sum caused by choosing any arm in the next step can then be computed in $$O(K)$$ operations (either add an extra term or reset to 0. Again, there are $$K$$ sums to be updated).

Having computed the possible changes to the summation terms in O(K) time, it is possible to calculate $$\mathbb{P}$$ for the possible states resulting from choosing any of the $$K$$ arms in $$O(K)$$ steps. The key here is to note that most of the terms in the product will be shared by all possible new states. In fact, it is possible to simplify $$\mathrm{Prob}(a)$$ by pulling these common terms out in the numerator and denominator and cancelling them, thereby eliminating the product computation entirely. Doing this will lead to the following:

Better solution

If the current state is $$(D_{a_1}, \dots, D_{a_K})$$, choose $$a$$ with probability proportional to $$\mathrm{Prob}(a) \sim \frac{\sum\limits_{d=1}^{\infty} \lambda(a, d)}{\sum\limits_{d=1}^{\infty} \lambda(a, d) - \sum\limits_{d=1}^{D_a} \lambda(a, d)}$$ normalized so that the sum of $$\mathrm{Prob}$$ over all arms is 1. The sums in this can be kept track of, as mentioned above, in $$O(K)$$ operations per step.

PS: It's been nearly 7 years since I was on this site, and this question made my morning. Thanks! And hopefully this solution works (my first attempt (which I didn't post) was wrong, so there's that)

• Thank you so much for this suggestion. I was struggling with this for quite some time now. I will try your suggestion. Once again, thanks a lot. Commented Aug 9, 2019 at 3:15
• I was aware of the conditions on $\lambda(a,d)$ that you have mentioned. In fact, these are the only conditions i believe. But what I was struggling with was a way to pick an arm given the current delay states. You have come up with a very clever joint distribution on the delays and seem to have exploited this to pick an arm. +1 for this idea. Commented Aug 9, 2019 at 3:20
• @S Prashanth Do you have any strong intuition as to why your proposed method would work but not the ones that I have mentioned in my question? Also, I would like to know if there are any references you can point me to on the kind of ideas you have in your answer. Commented Aug 9, 2019 at 7:14
• @Karthik The method in the answer and your approach 1 are based on the result on limiting distribution of Markov chains given here: en.wikipedia.org/wiki/… You can make the distribution of states work by choosing appropriate transition probabilities between states. The issue with your approach 1 is it treats just the next arm (and delay) as the state of the system. But that can't be right since transition probabilities are based on the current delays of all arms. So the state description should also include that... Commented Aug 12, 2019 at 17:15
• Sure. I'll let you know. Commented Aug 12, 2019 at 17:19