I'm trying to write a simulation that models population dynamics a set of organisms occupying a specific niche. Let's say I have 3 different species each with populations $n_1, n_2, n_3$. Together these populations sum to $N=\Sigma n_i$ which will be held constant. This system models these dynamics using discrete time steps. At each time there are 3 things that can happen. For each organism:
- With probability $\mu$ it will be replaced by an individual from any species with equal probability (including potentially it's own species, effectively resulting in nothing happening). This effectively models the organism randomly switching places with another one.
- With probability $d$ it will be replaced by an organism from species $i$ with probability $\phi_i$. Where $$\phi_i=\frac{(1+s_i)\frac{n_i}{N}}{1+\Sigma s_j \frac{n_j}{N}}$$ Here, the $s_i$ are called the selection coefficients (usually $\approx\pm .01$) and they are a measure of the organism's relative fitness. The $\phi_i$ thus models the probability that the organism will be replaced by a member of another species taking that species fitness and relative population size into account. Note that $\Sigma \phi_i=1$.
- If neither of the above happens then nothing else happens.
This occurs for every organism at each time step simultaneously, thus $\phi_i$ will not change in the middle of a time step. Also, because whenever one organism disappears it is immediately replaced by another, the total population $N$ remains constant. Typical values of $\mu$ are about $.01$, typical values of $d$ are about $.05$.
Now for my question. I know that I could simulate this with a brute force method by calculating these probabilities for each organism and then updating the variables at each time step. The problem is that I want to run this simulation for thousands of organisms across thousands of times steps and repeat this simulation thousands of times. The brute force method is slow and I want to speed it up.
The way that I think I can do this is by using multinomial distributions. At each time step, you have some set of $\{n_i\}$ summing to $N$ that have some probability of becoming a new set of $\{n'_i\}$ summing to $N$. If I could compute the probability distribution for each $\{n'_i\}$ I could speed this up considerably by removing the need to iterate through each organism and roll for possible outcomes.
I think that multinomial distributions hold the answer I seek, but I don't understand them all that well yet entirely so I'm not sure they could help here or how I would implement it in this case.
I wound up figuring out how to do this. Each step of the simulation can just be computed with a simple multinomial distributions. The first step is just
$$ \text{Split}_i:=multi(n_i,[\mu,d,1-\mu-d)]=(\text{mutating}_i,\text{dying}_i,\text{unchanged}_i)$$ And then this returns a 3-tuple containing the number of mutating, dying and unchanged organism. The next step is to draw a multinomial distribution for the mutating and dying population. for the mutating population I draw from $$\text{mutations}_i=multi(\text{mutating}_k,[\alpha_{1k},...,\alpha_{ik},...,\alpha_{Nk}])=(\text{mutations}_i^1,...,\text{mutations}_i^k,...,\text{mutations}_i^N)$$ where each $\alpha_{ik}=\frac{1-\delta_{ik}}{N-1}$. Where the delta function enforces that the probability of a species mutating to itself is $0$, while the probability of mutating to any other species is uniform. Then for replaced population we have $$\text{replaced}_i=multi(\text{dying}_i,[\phi_1,...,\phi_k,...,\phi_N])=(\text{replaced}_i^1,...,\text{replaced}_i^k,...,\text{replaced}_i^N)$$