Generating random data from a discrete multimodal distribution

I have a $$n \times n$$ distribution grid that represents the likelihood of finding a person within an area (see below). From this I would like to be able to create $$x$$ samples that fit this distribution.

Ideally $$50\lt x \lt200$$, and $$10 \lt n \lt 50$$. This is to simulate some search and rescue scenarios for a few path planning algorithms based on this distribution grid. • What's a nxn distribution? Are you generating matrices? What do you "have" a joint pdf? All stuff that needs to be included in your question. Nov 5 '20 at 13:33
• Also how big would n be for practical purpose? That might change the answer. Nov 5 '20 at 13:33
• @ThreeDiag sorry, you're right. I've edited my question. Nov 5 '20 at 13:58

If I translate in mathematical terms what I believe you said, we have that if we designate by by $$(a_{i,j})_{1\le i,j\le n}$$ the cases of your grid and $$X$$ a person, then you have the data of the probability that person $$X$$ pop up in case $$a_{i,j}$$. i.e., you have access to the $$n^2$$ numbers $$(p_{i,j})_{1\le i,j\le n}$$ such that $$\mathbb{P}(X \text{ is in the case } a_{i,j})=p_{i,j}.$$ and you want to be able to sample from this distribution ($$\mathbb{P}$$ stands for "the probability that").

This is a multinomial distribution (not multimodal, multimodal means another thing, see wikipedia). Just draw according to a multinomial distribution with parameters $$(p_{1,1},p_{1,2}, \dots,p_{1,n},p_{2,n},\dots,p_{n,n})$$, the fact that this is in 2D does not really matter, you can use the multinomial right away saying that if you draw a $$1$$ from the multinomial, your person is in case $$a_{1,1}$$, if you draw a $$n$$ he is in case $$a_{1,n}$$, if you draw a $$(n+1)$$ he is in case $$a_{1, 2}$$, if you draw $$(3\times n+6)$$, he is in case $$a_{6,3}$$... (you can enumerate the cases of your grid as you wish, I only give an example).

Drawing from a multinomial distribution is relatively easy using python or matlab.

• Hi yes that's perfect! that's something that I hadn't even considred. So if I understand what you're saying is that I can use $(x,y)$ as a sort of key. Nov 5 '20 at 14:56

Using what @TMat said, I got this code to work using numpy

img = Image.open("path/to/img")
img = np.array(img)[:,:,0]
prob = (np.array(img)/np.sum(img)).flatten()

x,y = np.meshgrid(np.arange(0,img.shape),np.arange(0,img.shape))
x,y = x.flatten(),y.flatten()
xy  = np.vstack((x,y)).T

xy_indices = np.arange(len(xy))

choices = np.random.choice(xy_indices, 200, p=prob)

points = xy[choices] • Yes it is exactly so. Your numpy implementation is great the only tricky part is to be sure that your first flatten (for the probability) gives the same enumeration as the flatten then vstack you used to construct xy. And I wonder if you do not have a problem there, using your code I think I would have used the transposed matrix prob = (np.array(img)/np.sum(img)).T.flatten() or something like that. As is, I beleive you exchanged the columns and the lines but in your case this has no consequence because the image seems to be symmetric.
– TMat
Nov 5 '20 at 15:32
• try to check that img[xy] corresponds to prob (up to normalization) when img is not symmetric (img.T != img). Maybe I am wrong, I didn't check myself.
– TMat
Nov 5 '20 at 15:35
• @TMat yup that's exactly what happened when I tested this with a more complex image. Worked a charm! Thanks! Nov 5 '20 at 15:35