https://doi.org/10.53656/math2025-5-1-nml

Резюме. In the present work we consider a special case of estimation in branching stochastic processes with random migration. The aim is to estimate the distributions of the offspring, of the emigration and of the immigration components, as well as the migration events probabilities, using information about the entire family tree. For this purpose, the nonparametric maximum likelihood method is chosen. Additionally, estimators of the probability mass functions are derived and evaluated through simulations and computational analysis.

Ключови думи: branching processes; migration; nonparametric estimators; simulation; maximum likelihood

1. Introduction

At some point in life, everyone becomes curious about how things evolve or disintegrate. For instance, we might wonder why certain species go extinct, or whether we can predict how a particular virus will spread. An even more interesting question is whether we can somehow control the dynamics of a population. These questions, and many more, can, in one way or another, be explored using the tools of branching processes (BP) (for example (CaronLormier et al. 2006) and (Haccou et al. 2005)).

In our work, motivated by (Nitcheva & Yanev 2000), we examine a BP that allows random migration. A significant early contribution to the theory of such processes, was made independently by (Nagaev & Khan 1981) and (Yanev & Mitov 1980). These models were developed as a special case of controlled branching processes, where the number of offspring in each generation can be influenced by random external factors, including the migration of individuals into or out of the system. Their work extended the classical Galton–Watson model by incorporating random migration terms, enriching the model’s applications to biological and ecological systems where such external movement is realistic and significant.

1.1. The general model

The model we simulate and derive estimators for is closely based on the one defined (Nitcheva & Yanev 2000) as follows.

On a probability space (\(\Omega, \mathcal{A}, \mathbb{P}\) ) , we define the following sets: \(X=\left\{X_{n, i}\right\}\), \(\eta=\left\{\left(\eta_{n, 1}, \eta_{n, 2}\right)\right\}\) and \(I=\left\{\left(I_{n}, I_{n}^{0}\right)\right\}\). Each set consists of independent nonnegative integer valued and identically distributed (i.i.d.) random variables (r.v.), and the sets themselves are mutually independent. The branching process with random migration (BPRM) \(\left\{Y_{n}\right\}\) is then defined by the recursive equations

(1) \[ Y_{0} \geq 0, \quad Y_{n+1}=\left(\sum_{i=1}^{Y_{n}} X_{n, i}+M_{n}\right)^{+}, \quad n=0,1, \ldots, \]

where

(2) \[ M_{n}= \begin{cases}-\left(\sum_{i=1}^{\eta_{n, 1}} X_{n, i}+\eta_{n, 2}\right), & \text { with probability } p \\ 0, & \text { with probability } q \\ I_{n} \mathbf{1}_{\left\{Y_{n} \gt 0\right\}}+I_{n}^{0} \mathbf{1}_{\left\{Y_{n}=0\right\}}, & \text { with probability } r\end{cases} \]

Here \(p+q+r=1\) and \(Y_{0}\) is independent of \(X, \eta\) and \(I\). As usual \(a^{+}=\) \(\max (0, a)\).

The interpretation of the components involved in the above definition relates to the branching and migration mechanisms within a population. The r.v. \(X_{n, i}\) denotes the number of offspring in the \((n+1)\)-st generation of the \(i\)-th i individual (called in some biologically driven situations the \(i\)-th particle) which exists in the \(n\)-th generation. The r.v. \(\eta_{n, 1}\) represents family emigration, that is, the number of families that emigrate before reproduction. In contrast, \(\eta_{n, 2}\) corresponds to individual emigration, meaning that \(\eta_{n, 2}\) individuals, randomly chosen from different families, are removed after reproduction. Finally, \(I_{n}\) and \(I_{n}^{0}\) represent the state-dependent immigration in the non-zero state and the zero state, respectively. The emigration, immigration or absence of migration in a given generation occur with the probabilities \(p, r\), r, and \(q=1-p-r\) r respectively. Note that there is a generalization of this model to the model with time non-homogeneous migration \(\left(p=p_{n}, q=q_{n}, r=r_{n}\right)\), investigated in (Yanev & Mitov 1985).

This definition incorporates many important cases. The process with \(q=1\) is the classical Bienayme-Galton-Watson (BGW) branching process. The process with \(r=1\) and \(I_{n}=I_{n}^{0}\) a.s. is the classical BGW process with immigration. The process with \(p=1\) (with pure emigration) was studied in the papers (Vatutin 1977), (Kaverin 1990), and (Grey 1988).

The properties of the model (1) are considered in the papers (Yanev & Mitov 1980), (Yanev & Mitov 1981), (Yanev & Mitov 1985), (Yanev & Yanev 1991), (Yanev & Yanev 1995), (Yanev & Yanev 1996) and (Yanev & Yanev 1997). In these papers additional parameters of criticality are introduced:

\[ M=\mathbb{E}\left[M_{n} \mid Y_{n-1} \gt 0\right]=r \mathbb{E}\left[I_{1}\right]-p \mathbb{E}\left[\eta_{1,1}+\eta_{1,2}\right] \] and

\[ \theta=\mathbb{E}\left[M_{n} \mid Y_{n-1} \gt 0\right] /\left(\tfrac{1}{2} \operatorname{Var} X_{1,1}\right)=M /\left(\tfrac{1}{2} \operatorname{Var} X_{1,1}\right) . \] They play an important role for the asymptotic behaviour of the process, i.e. in the critical situation when \(\mathbb{E}\left[X_{1,1}\right]\) the limiting behaviour depends on whether \(\theta\) is smaller, equal to or greater than zero. Hence, for this type of processes the standard classification of subcritical (\(\mathbb{E}\left[X_{1,1}\right] \lt 1\) ), critical \(\left(\mathbb{E}\left[X_{1,1}\right]=1\right)\) and supercritical \(\left(\mathbb{E}\left[X_{1,1}\right] \gt 1\right)\) process must be additionally specified by the parameter \(M\) or \(\theta\). From statistical point of view, if \(M \lt 0\), on the average, the emigration predominates the immigration, and if \(M \gt 0\), the immigration overtakes the effect of emigration (Yanev & Yanev 2000).

Note that the state zero is a reflecting barrier for the regenerative process \(\left\{Y_{n}\right\}\), which is an aperiodic and irreducible Markov chain.

1.2. The model under consideration

In the next sections we use two main simplifications of the above model. First, we exclude family emigration, migration before reproduction, but permit only individual emigration occurring after reproduction. Second, unlike the original model that employs two different immigration distributions (depending on whether the previous generation size is zero or positive), we use a single immigration distribution regardless of the population size in the previous generation:

(3) \[ M_{n}= \begin{cases}-\eta_{n}, & \text { with probability } p \\ 0, & \text { with probability } q \\ I_{n}, & \text { with probability } r\end{cases} \]

and \(p+q+r=1\). Hence, the three exploited sets of r.v. become \(X=\left\{X_{n, i}\right\}\), \(\eta=\left\{\eta_{n}=\eta_{n, 2}\right\}\) and \(I=\left\{I_{n}\right\}\). Additionally, our estimations would be based on the information about the entire family tree, rather than relying only on the size of each generation. Specifically, we assume that for each generation, we observe the migration events, including the number of individuals added or removed.

The paper is organized as follows. In Section 2 we consider the formularandom migration component and full information about the branching tree. Section 3 discusses some statistical properties of the estimators. Simulations and computational results are presented in Section 4.

2. The nonparametric maximum likelihood estimation

To generate the \((n+1)\)-st generation from the \(n\)-th, we proceed as follows. First, we observe the reproduction of the individuals in generation \(n\) through the r.v. \(\left\{X_{n, i}\right\}\), which are the number of offspring produced by the \(i\)-th individual, existing in the \(n\)-th generation. After reproduction, a migration event occurs with one of the following outcomes: emigration with probability \(p\), no migration with probability \(q\), or immigration with probability \(r\) and only one of the three can occur at each step. Also, \(a^{+}=\max \{0, a\}\), which guarantees that a generation cannot have a negative size.

For the estimations that follow, we use the method of maximum likelihood, assuming that complete information about the entire family tree is available. We will also introduce some notation that will be used throughout the paper:

1) \(m=\mathbb{E}\left[X_{n, i}\right], e=\mathbb{E}\left[\eta_{n}\right]\) and \(\lambda=\mathbb{E}\left[I_{n}\right]\) : the means of the offspring, emigration, and immigration distributions, respectively. It is a well-known fact (see (Harris 1963), among others) that, according to the value of \(m\), branching processes are categorized into three types: subcritical \((m \lt 1)\), critical (\(m=1\) ) and supercritical (\(m \gt 1\) ), as mentioned in the introduction.

2) \(M=\mathbb{E}\left[M_{n}\right]\) : the migration mean.

3) \(N_{l}^{o f f}(u)\) : the number of individuals in the \(l\)-th generation that have \(u\) offspring in the (\(l+1\) )-st generation.

4) \(N_{l}^{e m}(u)\) : the indicator of the event that emigration occurs in the \(l\)-th generation (after reproduction, i.e., during individual emigration), with exactly \(u\) emigrants removed from the population. These individuals do not participate in the (\(l+1\) )-st generation.

5) \(N_{l}^{i m m}(u)\) : the indicator of the event that immigration occurs in the \(l\)th generation (after reproduction), and exactly \(u\) immigrants are added to the population. This means that \(N_{l}^{i m m}(i)=1\), when \(i=u\) and \(N_{l}^{i m m}(i)=0\) for \(i \neq u\). These imm(i) = 0 for i = u. These immigrants then participate in the (\(l+1\) )st generation and reproduce according to the offspring distribution of the process.

6) \(N_{l}^{\text {em }}, N_{l}^{\text {null }}, N_{l}^{\text {imm }}\) : indicator random variables showing whether there

7) \(p_{u}^{o f f}, p_{u}^{e m}, p_{u}^{i m m}\) : the probability mass functions (p.m.f.s) of the number of offspring, emigrants, and immigrants, respectively.

Theorem 1. Let \(\left\{Y_{n}\right\}_{n \geq 0}\) be a BPRM defined with (1) and (3). The maximum likelihood estimators (MLEs) for \(p, q, r\) and \(p_{u}^{\tau} \gt 0\) for \(u=0,1,2, \ldots\) an \(\tau \in\{o f f, e m, i m m\}\), based on the entire family tree up to the \(n\)-th generation, are given by:

(4) \[ \hat{p}=\tfrac{\sum_{l=0}^{n-1} N_{l}^{e m}}{n}, \quad \hat{q}=\tfrac{\sum_{l=0}^{n-1} N_{l}^{\text {null }}}{n}, \quad \hat{r}=\tfrac{\sum_{l=0}^{n-1} N_{l}^{\text {imm }}}{n} \]

(5) \[ \hat{p}_{u}^{\tau}=\tfrac{\sum_{l=0}^{n-1} N_{l}^{\tau}(u)}{\sum_{l=0}^{n-1} \sum_{v=0}^{\infty} N_{l}^{\tau}(v)} \]

Proof. Throughout the proof, we use \(\tau \in\{o f f, e m, i m m\}\) , representing either the offspring, emigration, or immigration distribution. Unless otherwise specified, all expressions involving \(\tau\) are assumed to range over this set.

Using the fact that the individual reproduction is independent and the Markov property of the process, the likelihood function can be given by the following expression:

(6) \[ L\left(\mathbf{N} \mid p, q, r,\left\{p_{u}^{\tau}\right\}\right)=\prod_{l=0}^{n-1} p^{N_{l}^{e m}} q^{N_{l}^{n u l l}} r^{N_{l}^{i m m}} \prod_{u=0}^{\infty} \prod_{\tau}\left(p_{u}^{\tau}\right)^{N_{l}^{\tau}(u)}, \]

where

\[ \mathbf{N}=\left(N_{l}^{e m}, N_{l}^{n u l l}, N_{l}^{i m m}, N_{l}^{\tau}(u)\right), \] which gives us the information about the entire family history. For clarity, we will omit the parameters for the function in what follows and simply write \(L\).

In the above expression, the index \(u\) may range over different sets depending on the distribution of \(\tau\), since the support of each corresponding probability mass function can vary. For notational simplicity, we write the likelihood as a single expression, even though some values of \(u\) may not appear in the data for certain distributions. In such cases, the corresponding terms are implicitly excluded from the product, but the statements are still true.

It is well known that it is easier to inspect the logarithm of the likelihood function, since it reaches its maximum at the same points as the original. Moreover, because \(p+q+r=1\) and \(\sum p_{u}^{\tau}=1\), we introduce Lagrange multipliers \(\lambda^{M}\) and \(\lambda^{\tau}\). Therefore, the function to be maximized is:

(7) \[ \mathcal{L}=\log L+\lambda^{M}(1-p-q-r)+\sum_{\tau} \lambda^{\tau}\left(1-\sum_{u=0}^{\infty} p_{u}^{\tau}\right) \]

Taking into account that the \(\log L\) l looks like this:

(8) \[ \log L=\sum_{l=0}^{n-1}\left[N_{l}^{e m} \log p+N_{l}^{n u l l} \log q+N_{l}^{i m m} \log r\right]+\sum_{l=0}^{n-1} \sum_{u=0}^{\infty} \sum_{\tau} N_{l}^{\tau}(u) \log p_{u}^{\tau}, \] we now differentiate the function (7) with respect to each parameter, set the resulting expressions to zero, and solve for the estimates.

First,

(9) \[ \tfrac{\partial \mathcal{L}}{\partial p}=\tfrac{\sum_{l} N_{l}^{e m}}{p}-\lambda^{M}=0 \Rightarrow \hat{p}=\tfrac{\sum_{l} N_{l}^{e m}}{\lambda^{M}} \]

Analogously, for \(q\) and \(r\) we have

(10) \[ \hat{q}=\tfrac{\sum_{l} N_{l}^{\text {null }}}{\lambda^{M}} \quad \text { and } \quad \hat{r}=\tfrac{\sum_{l} N_{l}^{i m m}}{\lambda^{M}} \]

Differentiating with respect to \(\lambda^{M}\), we obtain the normalization constrain \(p+q+r=1\), from which it follows that

(11) \[ \hat{\lambda}^{M}=\sum_{l=0}^{n-1}\left(N_{l}^{e m}+N_{l}^{n u l l}+N_{l}^{i m m}\right)=n \]

Substituting (11) into (9) and (10), we derive the desired estimators for \(p, q\) and \(r\).

Next, for the estimators of the p.m.f.s from (7) we get

(12) \[ p_{u}^{\tau}=\tfrac{\sum_{l=0}^{n-1} N_{l}^{\tau}(u)}{\lambda^{\tau}} \]

Differentiating again but with respect to \(\lambda^{\tau}\) and solving the resulting equation, we find

(13) \[ \hat{\lambda}^{\tau}=\sum_{l=0}^{n-1} \sum_{v=0}^{\infty} N_{l}^{\tau}(v) \]

Substituting (13) into (12), we complete the proof of the theorem.

Corollary 1. Let \(\left\{Y_{n}\right\}_{n \geq 0}\) be a BPRM defined with (1)-(3). Using the MLEs obtained in (4) -(5) the MLEs for the offspring mean, \(m\), and the migration mean, \(M\), are given by:

(14) \[ \hat{m}=\sum_{u=0}^{\infty} u \hat{p}_{u}^{o f f} \quad \text { and } \quad \hat{M}=\hat{\lambda} \hat{r}-\hat{e} \hat{p} \]

where \(\hat{\lambda}=\sum_{u=0}^{\infty} u \hat{p}_{u}^{i m m}\) and \(\hat{e}=\sum_{u=0}^{\infty} u \hat{p}_{u}^{e m}\).

Proof. This follows from the invariance property of MLEs (see (Casella & Berger 2002), Theorem 7.2.10) and by using the induced likelihood function.

Comment 1. Note that the estimators of the offspring probabilities \(p_{u}^{o f f}\) have the same form as the estimators of the offspring probabilities for the general controlled branching process, proposed and studied in (Gonzalez et al. 2016). This is expected, because the BPRM process defined by (1) and (3) can be considered as a special case of a controlled branching process with random control function, defined by the migration component. In our setting, however, using the likelihood over the BPRM process, one can directly estimate the probabilities for migration and the distributions of the immigration and the emigration components.

3. On the statistical properties of the nonparametric maximum likelihood estimators

Note that the obtained estimators are correspondingly the relative frequencies of the following events: emigration occurs, immigration occurs or neither of these; \(u\) emigrants or immigrants in a generation are present, a particle has exactly \(u\) offspring. From the classical statistical theory of independent and identically distributed random observations these relative frequencies are “good” estimators of the theoretical probabilities (i.e. they are unbiased, consistent and asymptotically normal). We expect that the estimators from the previous section will also behave well.

Proposition 1. The estimators (4) − (5) of the migration probabilities \(p\), \(q\) and \(r\) and of the distribution of the immigration and of the emigration components in the BPRM process defined by (1) and (3) are unbiased and, conditionally upon nonextinction, consistent.

Proof. In the following proof, we assume that \(l=0, \ldots, n-1\) and \(u=0,1, \ldots\).

First, we will prove the stated properties for the estimators of the migration probabilities.

Since \(N_{l}^{e m}\) is an indicator random variable that takes the value 1 if we observe emigration in the \(l\)-th generation and 0 otherwise, we have that \(\mathbb{E}\left[N_{l}^{e m}\right]=p\) and \(\mathbb{D}\left[N_{l}^{e m}\right]=p(1-p)\). Using this, we obtain:

(15) \[ \mathbb{E}[\hat{p}(n)]=\tfrac{1}{n} \sum_{l=0}^{n-1} \mathbb{E}\left[N_{l}^{e m}\right]=\tfrac{n p}{n}=p \]

and

(16) \[ \mathbb{D}[\hat{p}(n)]=\tfrac{1}{n^{2}} \sum_{l=0}^{n-1} \mathbb{D}\left[N_{l}^{e m}\right]=\tfrac{p(1-p)}{n} \]

From (15) directly follows that \(\hat{p}\) is an unbiased estimator for \(p\).

Using the Chebyshev’s inequality, we have:

(17) \[ \mathbb{P}(|\hat{p}-\mathbb{E}[\hat{p}(n)]| \gt \varepsilon) \leq \tfrac{\mathbb{D}[\hat{p}(n)]}{\varepsilon^{2}} \]

Taking into account (15)–(17), we can conclude that \(\hat{p}(n)\) is also a weakly consistent estimator, because:

(18) \[ \lim _{n \rightarrow \infty} \mathbb{P}(|\hat{p}-p| \gt \varepsilon) \leq \lim _{n \rightarrow \infty} \tfrac{p(1-p)}{n \varepsilon^{2}}=0 \]

In fact, it is also strongly consistent, based on the law of large numbers (LLN). Similarly, the same results hold for \(\hat{q}\) and \(\hat{r}\) as well.

Now, we will consider the estimators for \(p_{k}^{i m m}\).

Note that \(\sum_{l} \sum_{v} N_{l}^{i m m}(v)\) is the total number of generations up to the \(n\)-th one, when immigration is observed. Hence

\[ P\left(\sum_{l} \sum_{v} N_{l}^{i m m}(v)=u_{i m m}\right)=\binom{n}{u_{i m m}} r^{u_{i m m}}(1-r)^{n-u_{i m m}} . \]

The sum \(\sum_{l} N_{l}^{i m m}(u)\) is the total number of generations in which immigration occurs and the number of the observed immigrants is exactly \(u . N_{l}^{i m m}(u)\) are independent for different generations \(l\) and in a fixed generation, say \(l\)-th, \(\sum_{u} N_{l}^{i m m}(u)=1\), if immigration is observed, and \(\sum_{u} N_{l}^{i m m}(u)=0\), if no immigration happened.

The random variable

\[ \left[\sum_{l} N_{l}^{i m m}(u) \mid \sum_{l} \sum_{v} N_{l}^{i m m}(v)=u_{i m m}\right] \] is interpreted as the number of generations, in which \(u\) immigrants are included, given that \(u_{i m m}\) times immigration is observed, hence

\[ \left[\sum_{l} N_{l}^{i m m}(u) \mid \sum_{l} \sum_{v} N_{l}^{i m m}(v)=u_{i m m}\right] \in B i\left(u_{i m m}, p_{u}^{i m m}\right) \] Using the law of total expectation and substituting the obtained estimator for \(\hat{p}_{u}^{i m m}\) we have that

(19) \[ \begin{aligned} & \mathbb{E}\left(\tfrac{\sum_{l} N_{l}^{i m m}(u)}{\sum_{l} \sum_{v} N_{l}^{i m m}(v)}\right) \\ & =\sum_{u_{i m m}=0}^{n-1} \mathbb{E}\left(\left.\tfrac{\sum_{l} N_{l}^{i m m}(u)}{\sum_{l} \sum_{v} N_{l}^{i m m}(v)} \right\rvert\, \sum_{l} \sum_{v} N_{l}^{i m m}(v)=u_{i m m}\right) \\ & \quad \times \mathbb{P}\left(\sum_{l} \sum_{v} N_{l}^{i m m}(v)=u_{i m m}\right) \\ & =\sum_{u_{i m m}=0}^{n-1} \tfrac{1}{u_{i m m}} \underbrace{\mathbb{E}\left(\sum_{l} N_{l}^{i m m}(u) \mid \sum_{l} \sum_{v} N_{l}^{i m m}(v)=u_{i m m}\right)}_{u_{i m m} p_{u}^{i m m}} \\ & \quad \times \mathbb{P}\left(\sum_{l} \sum_{v} N_{l}^{i m m}(v)=u_{i m m}\right) \\ & =p_{u}^{i m m} \sum_{u_{i m m}=0}^{n-1} \mathbb{P}\left(\sum_{l} \sum_{v} N_{l}^{i m m}(v)=u_{i m m}\right)=p_{u}^{i m m} \end{aligned} \]

which shows that these estimators are unbiased.

Let us now prove that \(\hat{p}_{u}^{i m m}\) is consistent. To do this, we divide the numerator and denominator by \(n\), and analyze each part separately. The LLN is again used to justify the limits, i.e., the convergences in probability.

For the numerator, we have

(20) \[ \tfrac{\sum_{l=0}^{n-1} N_{l}^{i m m}(u)}{n} \xrightarrow[n \rightarrow \infty]{P} \mathbb{E}\left[N_{l}^{i m m}(u)\right]=r \cdot p_{u}^{i m m} \]

and for the denominator,

(21) \[ \tfrac{\sum_{l=0}^{n-1} \sum_{v=0}^{\infty} N_{l}^{i m m}(v)}{n} \xrightarrow[n \rightarrow \infty]{P} \sum_{v=0}^{\infty} \mathbb{E}\left[N_{l}^{i m m}(v)\right]=r \]

From (20) –(21) , and as a direct consequence of the continuous mapping theorem, we conclude that

(22) \[ \hat{p}_{u}^{i m m} \xrightarrow[n \rightarrow \infty]{P} p_{u}^{i m m}, \]

which means that the estimators are consistent.

In case we observe the exact number of needed emigrants in the population, as we suppose, the above reasoning can be applied to the estimators of the emigration distribution, hence their unbiasedness and consistency can be proved. Such sampling scheme can be achieved, for instance, in a situation when one knows that a given number of individuals is needed to be removed from the population. In case, when the known number of individuals is greater than the population size, the available number of individuals in the generation, smaller than the required, is removed.

Comment 2. Let us comment on the estimation of the offspring probabilities. Although the number of offspring of the different individuals is independent and the likelihood function is expressed in a form of products of p.m.f.’s of independent r.v., we may expect that the estimators have properties similar to the classical MLE’s. However, one must take into account that the number of observations in each generation is a random variable (it is the total number of particles in the previous generation). This may lead to misleading inferences in some cases.

Note that \(N_{l}^{o f f}(u)\), as def (u) , as defined, is the number of particles in generation \(l\), that have exactly \(u\) offspring in the next \((l+1)\)-st generation. This offspring is of all particles available in the \(l\)-th generation, which may be “native” and immigrants. One has that \(\sum_{l} N_{l}^{\text {off }}(u)\) is the number of particles up to the \((n-1)\)-st generation with \(u\) offspring, and the r.v. Tot \(=\sum_{l} \sum_{v} N_{l}^{\text {off }}(v)\) represents the total number of all particles up to the \((n-1)\)-st generation, which is in fact the sum of the sizes of the first \((n-1)\) generations: T ot \(=\) \(\sum_{l=0}^{n-1} Y_{l}\). Hence the conditional joint distribution of the number of particles with a given number of offspring for the \(l\)-th generation given this generation is multinomial and therefore

\[ \left[N_{l}^{o f f}(u) \mid Y_{l}\right] \in B i\left(Y_{l}, p_{u}^{o f f}\right) \] Note that the following relations hold:

\[ \sum_{u} N_{l}^{o f f}(u)=Y_{l} \quad \text { and } \quad \sum_{u} u N_{l}^{o f f}(u)=Y_{l+1}+\eta_{n} . \] Because of these restrictions we cannot directly apply the steps in the proof above to the nonparametric maximum likelihood estimators for the individual probabilities.

For the controlled branching processes in (Gonzalez et al. 2016) under specified conditions on the process and the control function the consistency, the asymptotic unbiasedness and the asymptotic normality of the estimators \(\widehat{p}_{u}{ }^{o f f}\) are proved. Our simulation results show that these properties apply to the BPRM process as well.

4. Simulations

Although suitable empirical data have not yet been identified at the time of writing, we maintain that the proposed model has clear potential for application to practical problems involving migration mechanisms. A natural direction for future work is to seek such applications, particularly in fields such as biology and ecology, where migration plays a fundamental role in population dynamics. In the present study, however, we restrict our focus to simulations of the process, implemented in Python, which allow us to investigate the behavior of the proposed estimators and to assess their potential effectiveness in practice.

For the simulation we have used the following Binomial distributions: \(X_{n, i} \in \operatorname{Bi}(2,0.55), \eta_{n} \in \operatorname{Bi}(4,0.2)\) and \(I_{n} \in \operatorname{Bi}(5,0.4)\).

From the above distributions, we can see that each individual in the population can have up to 2 offspring; in the case of emigration, up to 4 individuals can be removed, and in the case of immigration, up to 5 individuals can be added. Also, we can calculate the actual means of the given distributions, which are: \(m=1.1\) (supercritical process), \(e=0.8\) and \(\lambda=2\).

For the migration probabilities, we used a randomly generated vector from the Dirichlet distribution, resulting in the values: \(p=0.5797, q=0.3530\) and \(r=0.0673\). Based on these values and the previously calculated means, the resulting migration mean is \(M=-0.3292\), which is to be expected, because the most probable event is the emigration.

Using the derived estimators in Theorem 1, the above predefined distributions and probabilities, we get the following results for \(50,100,150\) and 200 generations \({ }^{1}\).

The estimates for \(p, q\) and \(r\) (Table 1) are moving closer to their true values as the number of generations increases, which is expected due to consistency of MLEs. An exception is \(r\), possibly because immigration occurs with low probability, so more generations are needed to accurately estimate it.

The MLEs for the offspring distribution (Table 2) demonstrate excellent performance. Even with just 50 generations the estimated probabilities are really close to the true values. This suggests that the offspring process is relatively stable, likely due to its central role in each generation’s evolution.

The MLEs for the emigration distribution (Table 3) provide a good approximation of the true probabilities, particularly for the more frequent emigration outcomes. The rare event of observing 4 emigrants has not been captured in any of the simulations, including those with the largest sample size, due to its extremely low probability. This highlights the inherent difficulty in estimating very low-probability outcomes using MLEs, as such events may not occur even in moderately large samples.

Based on the results shown in the above tables and Corollary 1, we also present the MLEs for the means (Table 5) of the process defined in Definition (1), using (3) for \(M_{n}\). In addition, we include results for the migration mean, denoted by \(M_{N Y}\), estimated by Nitcheva and Yanev using the weighted conditional least squares method in (Nitcheva & Yanev 2000). As mentioned previously, their work served as an inspiration for our study.

From the estimators for the migration mean, we observe that both yield very similar results for large numbers of generations. However, this is not the case for smaller ones. Therefore, we present, in Table 6, only the final estimated values of \(M\) for 20 and 30 generations (the real value is the same as in Table 5).

These results suggest that our estimator performs better with smaller samples, likely due to the additional information available about the family tree.

The following graphics, Figure 1 and Figure 2, display the dynamics of the population, giving us more visual representation of the process. We observe that at the beginning the process has more fluctuations, but with the increase of the generations, the population size begins to grow exponentially.

In conclusion, while the estimation of offspring probabilities in the BPRM process shares structural similarities with classical MLE settings, the randomness of generation sizes introduces subtleties that may affect inference. The constraints on the offspring counts and their dependence on random generation sizes mean that standard estimation techniques require careful adaptation. Nevertheless, our simulations suggest that desirable properties such as unbiasedness and consistency hold under the BPRM framework. Future work will delve deeper into the asymptotic behaviour of the estimators and explore more complex scenarios, such as cases with insufficient individuals for removal and corresponding adjustments to the offspring distribution.

Acknowledgments

The study is financed by the European Union-NextGenerationEU, through the National Recovery and Resilience Plan of the Republic of Bulgaria, project SUMMIT BG-RRP-2.004-0008-C01.

NOTES

1. In the tables, a “−” denotes that the given outcome did not occur in the observed data.

REFERENCES

Caron-Lormier, G., Masson, J. P., Menard, N., Pierre, J. S., (2006). A branching process, its application in biology: Influence of demographic parameters on the social structure in mammal groups. Journal of Theoretical Biology, 238 (3), pp. 564 – 574. issn: 0022-5193. DOI: 10.1016/j.jtbi.2005.06.010.

Casella, G., Berger, R., (2002). Statistical Inference. 2nd. Pacific Grove, CA: Duxbury.

Gonzalez, M., Minuesa, C., Puerto, I., (2016). Maximum likelihood estimation and expectation-maximization algorithm for controlled branching processes. Computational Statistics and Data Analysis, 93, pp. 209 — 227. DOI: 10.1016/j.csda.2015.01.015.

Grey, D. R., (1988). Supercritical branching processes with density independent catastrophes. Mathematical Proceedings of the Cambridge Philosophical Society, 42, pp. 1773 — 1776. DOI:

Haccou, P., Jagers, P., Vatutin, V., (2005). Branching Processes: Variation, Growth, and Extinction of Populations. Cambridge Studies in Adaptive Dynamics. Cambridge University Press. DOI: 10.1017/CBO9780511629136.

Harris, T. E., (1963). Theory of Branching Processes. Springer, New York.

Kaverin, S.V., (1990). A refinement of limit theorems for critical branching processes with an emigration. Theory of Probability & Its Applications, 35, pp. 574 — 580. DOI: 10.1137/1135080.

Nagaev, S. V., Khan, L. V., (1981). Limit Theorems for a Critical Galton–Watson Branching Process with Migration. Theory of Probability & Its Applications, 25 (3), pp. 514 — 525. DOI: 10.1137/1125063.

Nitcheva, D., Yanev, N., (2000). A system for simulation and estimation of branching processes. PLISKA. Studia Mathematica Bulgarica, 13. Available at: http://www.math.bas.bg/pliska/Pliska-13/Pliska-132000-173-178.pdf.

Vatutin, V. A., (1977). A critical Galton-Watson branching process with emigration. Theory of Probability & Its Applications, 22, pp. 465 — 481. DOI: 10.1137/1122058.

Yanev, G., Yanev, N., (1991). On a new model of branching migration processes. Comptes rendus de l’Académie bulgare des sciences: sciences mathématiques et naturelles, 44, pp. 19 -– 22. Available at: researchgate.net/publication/266519315_On_a_new_ model_of_branching_migration_processes.

Yanev, G., Yanev, N., (1995). Critical branching processes with random migration. Lecture Notes in Statistics, 99, pp. 36 — 46. DOI: 10.1007/978-1-4612-2558-4_5.

Yanev, G., Yanev, N., (1996). Branching Processes with Two Types Emigration and State-Dependent Immigration. In: Heyde, C. C., Prohorov, Y.V., Pyke, R., Rachev, S.T. (eds) Athens Conference on Applied Probability and Time Series Analysis. Lecture Notes in Statistics, 114, pp. 323 — 336. DOI: 10.1007/978-1-4612-0749-8_15.

Yanev, G., Yanev, N., (1997). Limit theorems for branching processes with random migration stopped at zero. In: Classical and Modern Branching Processes, eds. K.B.Athreya and P.Jagers, IMA Volumes in Mathematics and its Applications, 84, pp. 323 – 336. DOI: 10.1007/978-1-46121862-3_25.

Yanev, G., Yanev, N., (2000). Limit theorems for branching processes with random migration components. PLISKA. Studia Mathematica Bulgarica, 13. Available at: https://www.math.bas.bg/pliska/Pliska13/Pliska-13-2000-199-205.pdf.

Yanev, N., Mitov, K., (1980). Controlled branching processes: The case of random migration. Comptes rendus de l’Académie bulgare des sciences: sciences mathématiques et naturelles, 33, pp. 473 — 475.

Yanev, N., Mitov, K., (1981). Critical branching migration processes. (in Russian). Mathematics and Mathematical education, Proceedings of the 10th Spring Conference of the Union of Bulgarian Mathematicians, pp. 321 — 328.

Yanev, N., Mitov, K., (1985). Critical Branching Processes with Nonhomogeneous Migration. The Annals of Probability, 13 , pp. 923 — 933. Available at: http://www.jstor.org/stable/2243719.