Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
58,911 characters · 13 sections · 56 citation commands
Bayesian Analysis of Mixtures of Lognormal Distribution with an Unknown Number of Components from Grouped Data
Finding a distribution that fits the data well is one of the main challenges in the estimation of income distributions. However, we face the trade-off between the interpretation of parameters and the fit of the hypothetical distribution. To explore the fit of the distribution, several flexible hypothetical distributions are proposed, including: the generalized beta distribution of first and second kind M84; generalized beta distribution MX95; double Pareto-lognormal distribution RJ04; and $\kappa$-generalized distribution CGK07. Several of these support interpretations that are economically meaningful. \footnote{ To better fit income distribution models to empirical data, Bayesian Model Averaging (BMA), as explored by GCR05, has also been proposed as a viable approach. }<++>
Conversely, the mixture distribution models are also considered to fit the distribution to the data because the assumed underlying distributions are easy to interpret and the distribution fits better than single component models in many cases. The greater level of detail offered by mixture distribution models, such as a subgroups' information, is evident from the model's adoption in PvD98,GH12, among other studies.
Mixture distribution models have also considered the framework of household income distributions from a Bayesian point of view using Markov chain Monte Carlo (MCMC) methods. For example, in the case of lognormal distribution, LN16 considered a finite mixtures of lognormal (MLN) distribution model from individual data and determined the number of components by the marginal likelihood C95 and DIC SBCvdL02. The income inequality was then decomposed into between-subgroup and within-subgroup components. Moreover, it is also considered in gamma distribution cases. WIR01 considered the mixtures of gamma distribution model with a known and unknown number of components. CG08 examined the Canadian income data using two components' mixtures of gamma densities, which is the known number of components case in WIR01. However, with the exception of WIR01, the number of components were assumed in advance or determined after estimation in these studies, and they used individual or household data as mentioned above.
As with WIR01, there are two main approaches for dealing with an unknown number of components in a mixture model: one uses a Dirichlet process prior EW95, and the other uses a reversible jump MCMC algorithm RG97, which is used in WIR01. The reversible jump MCMC algorithm, which was first proposed by G95, is one of the most powerful tools in model determination. RG97 proposed the algorithm in the framework of the mixtures of normal distribution model. Subsequently, some scholars have proposed extensions to the multivariate normal distribution K09 and the mixtures of normal distribution with the same component means PI09. In addition, MH13 pointed out that the posterior from a Dirichlet process prior for the number of components was not consistent---unlike with the reversible jump MCMC, the Dirichlet process prior did not converge at the true number. Therefore, we consider the reversible jump MCMC algorithm in this study, because we are also interested in the number of components in the analysis of income distribution.
Although the availability of individual and household data has improved, it remains difficult to access, especially in developing countries. Alternatively, the grouped data, which partitions the sample space of observations into several non-overlapping groups, is widely available. Using this type of data, GdL14 considered the MCMC sampling scheme for finite mixtures of normal distribution.
This study extends their approach in two significant directions. First, we generalize the assumed distribution in GdL14 from the normal distribution to the lognormal distribution. This allows for a more realistic modeling of income data, which is typically skewed and strictly positive. Second, instead of fixing the number of components in advance, we adopt the reversible jump MCMC algorithm proposed by RG97, enabling us to estimate the number of components directly from the data. This extension enhances the model's flexibility and allows it to capture potential overdispersion in the structure of income subgroups, which is especially valuable when analyzing heterogeneous populations based on grouped data.
Exploring this model is worthwhile because the number of components provides information about population subgroups, as discussed in LN16. Therefore, if it is possible to determine the number of components from grouped data, this approach can be used for detailed comparisons of income inequalities in developing countries.
This study aims to develop a reversible jump MCMC method for the mixtures of lognormal (MLN) distribution model from grouped data to examine the income distributions and income inequalities in Japan. Our proposed algorithm is discussed using simulated data examples. From these, we can confirm that our proposed algorithm works well in terms of the accuracy of the parameters and in fitting the distribution. The data also suggests that the posterior distributions of the Gini coefficients are accurate. Hence, we applied it to real data in Japan in 2020 to examine the income distributions and inequalities. From the results, we identified two subgroups in both two-or-more person households and workers' households. We also observed that the Gini coefficient of two-or-more person households was larger than that of workers' households.
The rest of this paper is organized as follows. In Section (ref), we summarize the MLN distribution model using grouped data with its Gini coefficient and obtain a joint posterior distribution. Section (ref) discusses the computational strategy of the MCMC method. In Section (ref), our approach is illustrated using simulated datasets. Section (ref), examines the empirical examples using real datasets from Japan. Finally, brief conclusions are offered in Section (ref).
Let $x > 0$, which means the annual income of households or individuals, for example, follow any hypothetical distribution. Let $x_{i}$, $i = 1, 2, \ldots, n$ observations be sampled from the distribution. Then, the grouped data partitions the sample space of observations into $K > 1$ non-overlapping intervals of the forms $(t_{0}, t_{1}]$, $(t_{1}, t_{2}]$, $\ldots$, $(t_{K-1}, t_{K})$, where $t_{0} = 0$ and $t_{K} =\infty$. Moreover, only the number, $n_{k}$ of observations falling in each interval $(t_{k-1},t_{k}]$, $k = 1, 2,\ldots, K$, can be observed with $\displaystyle \sum_{k = 1}^{K} n_{k} = n$. It should be mentioned that the class income mean $\bar{x}_{k}$, which means the average of $x_{i}$ in the interval $(t_{k-1},t_{k}]$, is also available in many cases.
Let $\text{\boldmath{$\theta$}}$ be the vector of parameters of any underlying hypothetical distribution, which we assume in advance. Let $f(x|\text{\boldmath{$\theta$}})$ and $F(x|\text{\boldmath{$\theta$}})$ be the probability density function (PDF) and cumulative distribution function (CDF), respectively. Given the PDF and CDF, we define the likelihood function, which is based on the concept of selected order statistics, to estimate the parameters of the distribution. \footnote{ MR79 considered the likelihood based on the multinomial distribution, whereas NK11 considered the likelihood based on the selected order statistics. As is pointed out by EB21, the likelihood based on the multinomial distribution is applicable to the data with known fixed boundaries and random frequencies, while the likelihood based on the selected order statistics is applicable to the data with known random boundaries and fixed frequencies. In this study, we follow the likelihood based on NK11, because we are interested in the decile data, whose features are with known random boundaries and fixed frequencies. It should be mentioned that our approach merely treats the special case of DGP1 in EB21. } To explain the likelihood function, let $\mathbf{t} = (t_{1}, t_{2}, \ldots, t_{K-1})^{\prime}$ be the vector of the endpoints of the intervals and let $\mathbf{n} = (n_{1}, n_{2}, \ldots, n_{K})^{\prime}$ be the vector of frequencies, which fall in the intervals. Then, the likelihood function is defined as follows:
Once the parameter estimate for $\text{\boldmath{$\theta$}}$ is obtained from (ref) using maximum likelihood and so on, the Gini coefficient can be estimated by using
where $\mu$ is the mean of the distribution. \footnote{ In the numerical integration, we use the expression
because it is equivalent to (ref) D79 and easier than calculating (ref). }
In the empirical analysis we need to specify the hypothetical income distribution. First, we start with the lognormal (LN) distribution, following NK11, because the distribution fits to the Japanese data, which is also used in this empirical example. Although we could consider the other distributions, such as a gamma distribution and so on, we restrict our discussion on the LN distribution to focus on our empirical example. Let $x \sim \mathcal{LN}(\mu, \sigma^{2})$, which means $x$ follows LN distribution, where the PDF is expressed by
and the CDF is expressed by
where $\Phi(\cdot)$ is the CDF of the standard normal distribution. If we substitute (ref) and (ref) for (ref), it becomes the likelihood function for the LN distribution model and its Gini coefficient has a closed form, expressed by
To extend the above results, we consider the MLN distribution model with $R$ components. Let us begin with the fixed number of components model. Let $\text{\boldmath{$\pi$}} = (\pi_{1}, \pi_{2},\ldots, \pi_{R})^{\prime}$, $\text{\boldmath{$\theta$}}_{r} = (\mu_{r}, \sigma_{r}^{2})^{\prime}$, and $\mathbf{\Theta} = \left\{ \text{\boldmath{$\theta$}}_{r} \right\}_{r = 1}^{R}$, where $\displaystyle \sum_{r = 1}^{R} \pi_{r} = 1$. Then, the PDF of the MLN distribution with $R$ components is expressed by
and the CDF is expressed by
If we substitute (ref) and (ref) for (ref), it becomes the likelihood function for the MLN distribution model. However, its Gini coefficient does not have a closed form. Therefore, it is calculated from (ref). In the next section, we will consider the MLN distribution model with an unknown number of components, where $R$ is also treated as one of the parameters.
The likelihood function given in (ref) for the MLN distribution model is not particularly useful for Bayesian inference because its full conditional distributions are not the standard forms. In this study, we consider an alternative approach based on the framework by GdL14. They proposed augmenting the model with vectors of latent variables $\mathbf{x} = (x_{1}, x_{2}, \ldots, x_{n})^{\prime}$ and $x_{n_{k}^{*}}$ is set to $t_{k}$ for $k = 1,\ldots, K-1$, where $\displaystyle n_{k}^{*} = \sum_{j = 1}^{k}n_{j}$ and $\mathbf{z} = (z_{1}, z_{2}, \ldots, z_{n})^{\prime}$, where $z_{i} = r \in \{1,2, \ldots, R\}$. To complete this augmented likelihood, we also introduce a vector of observed variable $\displaystyle \mathbf{d} = (\underbrace{1,1,\ldots,1}_{n_{1}},\ldots, \underbrace{k,k,\ldots,k}_{n_{k}},\ldots,\underbrace{K,K,\ldots,K}_{n_{K}})^{\prime}$, instead of $\mathbf{n}$. Then, the joint likelihood of $(\mathbf{x}, \mathbf{d}, \mathbf{z}, \text{\boldmath{$\pi$}}, \mathbf{\Theta})$ can be specified as
where $n_{r} = \#\left\{ i: z_{i} = r \right\}$ and $I(A)$ denotes the indicator function of the event $A$.
As we adopt a Bayesian approach and extend the model to allow the number of components to change, we complete the model by specifying the following hierarchical prior distributions over the parameters ($R, \text{\boldmath{$\pi$}}, \mathbf{\Theta}$):
where $\mathcal{PO}(\lambda)$ means a Poisson distribution and $I(R \le R_{\max})$ imposes the preassigned upper limit $R_{\max}$ on the number of components. $\mathcal{D}(\alpha, \alpha, \ldots, \alpha)$ means a symmetric Dirichlet distribution and $\mathcal{G}(a,b)$ is a gamma distribution with scale and shape parameters $a$ and $b$, respectively.
It should be mentioned that the use of midpoint and range of data is recommended as the hyper-parameters in RG97. On the other hand, GdL14 used an improper prior. For the grouped data in income distribution, it is difficult to find a midpoint and range of data. Therefore, the hierarchical prior is assigned and they are also treated as parameters in the model to avoid it.
As the joint posterior distribution is much simplified, we can now use MCMC methods. The Markov chain sampling scheme can be constructed from birth-and-death process, split-or-combine process, and the full conditional distributions of $R, \text{\boldmath{$\pi$}}, \left\{ \mu_{r} \right\}, \left\{ \sigma_{r}^{2} \right\}, \left\{ z_{i} \right\}, \left\{ x_{i} \right\}_{i:i \ne n_{k}^{*}}, \mu, \tau^{2}, \beta$.
To implement a birth and death process, we first make a random choice between birth and death with the probability $b_{R}$ and $d_{R}$, where $d_{R} = 1 - b_{R} = 0.5$ except for $d_{1} = 0$ and $b_{R_{\max}} = 0$. For a birth process, a weight and parameters for the proposed new component are sampled from
where $\mathcal{B}(p,q)$ is a beta distribution. For a death process, a random choice is made between any empty components and the chosen component is deleted. Then, the acceptance probabilities $\min(1, A)$ and $\min(1, A^{-1})$ for birth and death are evaluated by
where $g_{p,q}$ denotes the $\mathcal{B}(p,q)$ density, $B(p,q)$ is a beta function and $R_{0}$ is the number of empty components RG98.
Using the same probabilities $b_{R}$ and $d_{R}$ as above, we make a random choice between attempting to split or combine, depending on $R$. Our combine proposal begins by choosing a pair of adjacent components ($r_{1}, r_{2}$) at random, which satisfies that there is no other $\mu_{r}$ in the interval $[\mu_{r_{1}}, \mu_{r_{2}}]$. Then, the combined component, here labeled $r^{*}$, is created as ($\pi_{r^{*}}, \mu_{r^{*}}, \sigma_{r^{*}}^{2}$), which satisfies the following equations:
To make a split proposal, we begin with choosing a component, here labeled $r^{*}$, and drawing a three-dimensional random variables as follows:
Then, the split proposal is made as follows:
where the adjacency condition that there is no other $\mu_{r}$ in the interval $[\mu_{r_{1}}, \mu_{r_{2}}]$ is satisfied.
Finally, the acceptance probability $\min(1, A)$ for split is evaluated by
where $l_{1}$ and $l_{2}$ are the numbers of observations proposed to be assigned to $r_{1}$ and $r_{2}$ and $P_{\text{alloc}}$ is the probability that this particular allocation is made. For the corresponding combine move, the acceptance probability is $\min(1, A^{-1})$, using the same expression for $A$, but some obvious differences in the substitutions.
Unsurprisingly, we can use the same acceptance probabilities with RG97 in spite of the fact that this model includes additional parameters ($\mu,\ \tau^{2}$) and a latent vector ($\mathbf{x}$). This is because the additional parameters and a latent vector are independent from the number of components. These two processes are explained in more detail in RG97.
With the exception of some hyper-parameters, the other parameters are easily sampled from the standard distributions following DR94 and GdL14. The full conditional distribution for $\text{\boldmath{$\pi$}}$ remains Dirichlet in form:
where we use `$|\cdots$' to denote conditioning on all other variables.
The full conditionals for $\left\{ \mu_{r} \right\}$ and $\left\{ \sigma_{r}^{2} \right\}$ are
where $\displaystyle \hat{\tau}_{r}^{2} = \left( \sigma_{r}^{-2} n_{r} + \tau^{-2} \right)^{-1}$, $\displaystyle \hat{\mu}_{r} = \hat{\tau}_{r}^{2} \left( \sigma_{r}^{-2} \sum_{i: z_{i} = r} \ln x_{i} + \tau^{-2} \mu \right)$, $\hat{\nu}_{r} = 0.5 n_{r} + \nu_{0}$ and $\displaystyle \hat{\beta}_{r} = 0.5 \sum_{i: z_{i} = r}\left( \ln x_{i} - \mu_{r} \right)^{2} + \beta$. Although LN16 considered several restrictions to avoid the label switching problem, we simply assume that $\mu_{1}<\mu_{2}\ldots<\mu_{R}$ to help remove the label switching problem DR94.
For the allocation variables, we have \footnote{ It should be mentioned that GdL14 derived the full conditional distribution of (B) $z_{i} | \text{\boldmath{$\pi$}}, \mathbf{\Theta}, d_{i}$, which is one without the condition on $x_{i}$. However, it is not required to decide the initial values of $z_{i}$ and we can choose any $z_{i}$. The simplest example is to start from $R=1$. Even if we start from any $R > 1$, (B) in Step 2 in GdL14 is not required. We can start from any random $z_{i}$. }
For the latent variables $x_{i}$, $i=1,2,\ldots,n$ except for $i = n_{k}^{*}$, $k = 1,2,\ldots, K-1$, the full conditional distributions are
For the hyper-parameters that we are not treating as fixed, $\mu$, $\tau^{2}$ and $\beta$, have
where $\hat{\tau}^{2} = (\tau^{-2}R + \tau_{0}^{-2})^{-1}$, $\displaystyle \hat{\mu} = \hat{\tau}^{2}(\tau^{-2}\sum_{r=1}^{R} \mu_{r} + \tau_{0}^{-2} \mu_{0})$, $\hat{g} = R \nu_{0} + g_{0}$, $\displaystyle \hat{h} = \sum_{r=1}^{R}\sigma_{r}^{-2} + h_{0}$.
To illustrate the Bayesian approach discussed in the previous section, we consider two simulated data examples. One is the case where the true data generating process (DGP) is the MLN distribution, and the other is the case where the true DGP is the generalized beta distribution of the second kind (GB2 distribution). In the first example, we examine the performance of our method and compare the GB2 distribution. In the second example, we explore the possibility of the MLN distribution model assuming that the true DGP is the GB2 distribution. The reason for choosing the GB2 distribution as the competing distribution against the MLN distribution is that the GB2 distribution is reported to fit the data well in the empirical analyses. Therefore, it is worthwhile to examine the fit of the MLN distribution when the true DGP is the GB2 distribution. All the results reported here were generated using Ox version 9.30 (macOS_64/Parallel) D13 and all the figures are drawn using R version 4.5.1 R25.
In the first simulated example, we examine the performance of our algorithm, and then compare the distribution with the GB2 distribution. The simulated data, wherein the DGP is the MLN distribution, is generated as follows. First, $x_{i}$, $i = 1,\ldots, 10,000$ were generated from the MLN distribution with three components ($R =3$) with parameters $\text{\boldmath{$\pi$}} = (0.2,\ 0.5,\ 0.3)^{\prime}$, $\text{\boldmath{$\mu$}} = (2.0,\ 3.0,\ 4.0)^{\prime}$ and $\text{\boldmath{$\sigma$}}^{2} = (0.3,\ 0.1,\ 0.2)^{\prime}$. The generated random numbers are sorted in ascending order, and $x_{n_{k}}$ corresponds to the $n_{k}$th observation, where $\displaystyle n_{k} = n \times \frac{k}{K}$. Then, $k=1,2,\ldots,K-1$ is picked up and $\mathbf{t} = (t_{1}, t_{2},\ldots,t_{K-1})^{\prime}$, where $t_{k} = x_{n_{k}}$, are collected. In this example, $K$ is set to $10$, which means that the dataset is decile data. Figure (ref) shows the true distribution and histogram, which is drawn from the simulated data. From the figure, we can observe the following features under this setting: (i) the first and second components' modes are found easily whereas the third one is not in the true distribution; (ii) the histogram looks like a unimodal distribution with heavy tail, that is, it is difficult to identify the first component as well from the grouped data.
To proceed with the Bayesian analysis, we need to set the hyper-parameters. For the prior distributions, we set the hyper-parameters as follows:
With the simulated data, we ran the MCMC algorithm using $500,000$ and discarding the first $100,000$ iterations.
Figure (ref) shows the posterior distribution of $R$. From the figure, we can confirm that the true number of components $R = 3$ shows the highest posterior mass. Therefore, we will focus on the result of the conditional posterior results on $R=3$ hereafter, if we report the result of a conditional one. We can also conclude that our algorithm can identify the exact number of components. Figure (ref) shows the unconditional predictive distribution and the predictive distribution conditioned on $R=3$. From the figure, we find that the unconditional predictive distribution and the predictive distribution conditioned on $R=3$ show the similar shape. In addition, they also show a similar shape with the true distribution, although there is a slight difference around the first component, where the first mode of the predictive distribution is lower than that of the true distribution. To see the difference between the predictive and true distributions, Table (ref) shows the conditional posterior estimates on $R=3$. From the table, we can see that all the estimates include the true values in the 95% credible intervals and most of the posterior means are close to the true values. However, it seems difficult to identify the first component because the standard deviations of the first component are larger than those of other components, and the 95% credible intervals of the first component are wider than those of others. Nevertheless, we can conclude that our method not only identifies the true number of components, but also estimates the parameters accurately.
Based on the favorable performance of the model and our algorithm in the simulated data, we can consider the goodness-of-fit of this model. BMM97,KN19 confirmed that the fit of the GB2 distribution proposed by M84 performed quite well in the empirical analyses. Therefore, we examine the goodness-of-fit of the GB2 distribution model when the true distribution is the MLN distribution model. The parameters of the GB2 distribution is estimated using the Tailored randomized block Metropolis-Hastings (TaRBMH) algorithm by KN19, which is first proposed by CR10 for estimating the DSGE model, and we examine the goodness-of-fit using the marginal likelihoods in line with KN19. Table (ref) shows the log of marginal likelihoods for these distributions, which are calculated by the harmonic mean estimator proposed by NR94 for its simplicity. \footnote{ The parameters and the harmonic mean estimate from GB2 distribution are estimated independently from the MLN distribution and we ran the MCMC algorithm using $40,000$ and discarding the first $10,000$ iterations. For the parameters $a$, $b$, $p$, $q$ in the GB2 distribution, we assume the following prior distributions:
and the hyper-parameters are set to $\alpha_{0} = \beta_{0} = \gamma_{0} = \delta_{0} = \epsilon_{0} = \zeta_{0} = \eta_{0} = \theta_{0} = 1.0$. } From the table, we can confirm that the unconditional MLN distribution model shows the highest and the GB2 distribution model shows the lowest log of marginal likelihood. The log of marginal likelihood of the conditional MLN distribution on $R=3$ lies between the unconditional MLN and GB2 distribution models and is much higher than that of GB2 distribution. Figure (ref) shows the unconditional predictive distributions for the MLN distribution, conditional predictive distribution for the MLN distribution on $R=3$, and predictive distribution for the GB2 distributions. From the figure, we can confirm that the unconditional predictive distribution for the MLN distribution and the conditional predictive distribution for the MLN distribution on $R=3$ seem to be similar to the true distribution, whereas the predictive distribution for the GB2 distribution seems to be different from the true distribution. Therefore, we can conclude that the unconditional predictive distribution for the MLN distribution model is a good fit for the income distribution, and, that the conditional predictive distribution, which includes the number of components, is useful for interpretation in an economically meaningful way.
Once the parameters are estimated, the Gini coefficient can be calculated from the parameters using (ref). To confirm the accuracy of the estimates, we calculated the Gini coefficient by means of the numerical integration. Figure (ref) shows the posterior distributions of the Gini coefficients for the unconditional MLN distribution, the conditional MLN distribution on $R=3$, and the GB2 distribution. The posterior summaries are also reported in Table (ref). To see the accuracy of the Gini coefficients, the nonparametric lower and upper bounds of the Gini coefficient ($0.4144, 0.4226$), which was proposed by G72, are shaded in the figure. From the results, we can confirm that not only all of the $95$% credible intervals are wider than the nonparametric bound, but also the posterior means are estimated outside the bound. Moreover, the posterior means for the unconditional MLN distribution model and the conditional MLN distribution on $R=3$ are similar, whereas the posterior mean for the GB2 distribution is farther than these means from the true value ($0.4196$). Simultaneously, we can see that the posterior modes for both the unconditional and conditional MLN distribution models on $R=3$ seem to be close to the true value. Therefore, we also calculated the posterior modes using the half sample mode estimator by RC74. From the results, we can confirm that the posterior modes for the unconditional MLN distribution are estimated accurately. Therefore, we can conclude that our algorithm can accurately estimate not only the number of components, but also the posterior estimates including the Gini coefficient, which is constructed as the function of the parameters.
In the previous subsection, we confirm that our algorithm can identify the true number of components and that the Gini coefficient can be calculated accurately. In this subsection, we will consider the contrary situation, wherein the true distribution is different from the MLN distribution. As we have found that when the MLN distribution is the true distribution, the GB2 distribution may not fit the distribution well or the calculation of the Gini coefficient may not be accurate, we explore the possibility of the MLN distribution model by considering the contrary case. The simulated data, wherein the DGP is the GB2 distribution, is generated as follows. First, $x_{i}$, $i = 1,\ldots, 10,000$ were generated from the GB2 distribution with parameters $a = 2.0$, $b=10.0$, $p=2.5$, and $q=1.5$. The generated random numbers are sorted in ascending order, and $x_{n_{k}}$ corresponds to the $n_{k}$th observation, where $\displaystyle n_{k} = n \times \frac{k}{K}$. Then, $k=1,2,\ldots,K-1$ is picked up and $\mathbf{t} = (t_{1}, t_{2},\ldots,t_{K-1})^{\prime}$, where $t_{k} = x_{n_{k}}$, are collected. In this example, $K$ is also set to $10$. Figure (ref) shows the true distribution and histogram, which is drawn from the simulated data. From the figure, we can confirm that the shape of the distribution exhibits the standard shape for income distribution, that is, it is unimodal and is right-skewed. Using the same hyper-parameters with the previous subsection, we ran the MCMC algorithm using $500,000$ and discarding the first $100,000$ iterations. \footnote{ For the case of the GB2 distribution, we ran the MCMC algorithm using $50,000$ and discarding the first $10,000$ iterations. }
To keep this paper focused on the results of the MLN distribution model, we do not report the estimation result of the GB2 distribution; however, it should be mentioned that it was estimated quite well. Figure (ref) shows the posterior distribution of $R$ and the number of components $R=2$ shows the highest posterior mass. Therefore, we will focus on the result of the conditional posterior result on $R=2$ hereafter. As with the previous subsection, we will examine the goodness-of-fit of the distributions to examine the possibility of the MLN distribution model. Table (ref) shows the log of marginal likelihoods for these distributions. From the table, we can confirm that the log of marginal likelihood of the GB2 distribution model shows the highest value. However, that of the unconditional MLN distribution is not so different if we take the standard error into account. Conversely, that of the conditional MLN distribution model on $R=2$ indicates a smaller value than these distribution models. Therefore, if we only focus on the fit of the distribution, the unconditional MLN distribution model becomes the alternative candidate against the GB2 distribution model. However, if we are interested in the economically meaningful interpretation, the marginal likelihood can distinguish which distribution is preferred, because the log of marginal likelihood of the conditional MLN distribution model is smaller than that of the GB2 distribution.
Figure (ref) shows the posterior predictive distributions for these distribution models. From the figure, we can see that all the predictive distributions are close to the true distribution. However, we can also confirm that the conditional posterior predictive distribution on $R=2$ exhibits a small difference from the true distribution and that difference may lead to the difference in the log of marginal likelihoods. Therefore, the choice of the distribution is very important and the marginal likelihoods play an important role to examine the goodness of fit.
Finally, we examine the accuracy of the Gini coefficients. Figure (ref) shows the posterior distributions of the Gini coefficients for the unconditional MLN distribution, conditional MLN distribution on $R=2$, and GB2 distribution. The posterior summaries are also reported in Table (ref). To see the accuracy of the Gini coefficients, the nonparametric lower and upper bounds of the Gini coefficients ($0.3343, 0.3439$) are shaded in the figure. From the results, we can confirm that all of the posterior means are included in the nonparametric bound, although all of the $95$% credible intervals are wider than the nonparametric bound. Moreover, the posterior mean for the GB2 distribution is closest to the true value ($0.3438$). However, it should be mentioned that posterior modes for the unconditional and conditional MLN distribution models are estimated outside the bound. Therefore, we can conclude that the Gini coefficients from both distributions can be calculated accurately. However, those form these distributions infer to the GB2 distribution in terms of accuracy if the true DGP is the GB2 distribution.
Using the Japanese household survey, Family Income and Expenditure Survey (FIES) in 2020, which is compiled by the Statistics Bureau of the Ministry of Internal Affairs and Communications, we will consider the income distributions and income inequalities in Japan. The survey offers data, presented in Table 3, for two types of households: Yearly Average of Monthly Receipts and Disbursements per Household by Yearly Income Quintile Group, and by Yearly Income Decile Group (Two-or-more-person Households). Yearly pre-tax income is surveyed for two types of households, depending on the occupation of the households' head: workers' households are employed as clerks or wage earners by public or private enterprises, such as government office, private companies, factories, schools, hospitals, shops, etc, whereas two-or-more person households include those other than workers' households, such as individual proprietors households and households whose heads are merchants, artisans or administrators of unincorporated enterprise. The sample size for each dataset is $n = 10,000$ and the dataset in decile form is utilized, therefore $n_{k} = 1,000$ for $k = 1,2,\ldots,K=10$. Thus, only $\mathbf{t}$ (unit: million yen) is different in each dataset.
Using the same hyper-parameters as in the previous section, we ran the MCMC algorithm using $500,000$ iteration and discarding the first $100,000$ iterations for each dataset. Figure (ref) shows the posterior distribution of $R$ and we can see that the posteriors for $R$ favor the model with 2 components both for two-or-more person households and for workers' households. Therefore, we proceed our discussion based on the results from $R=2$ both for two-or-more person households and for workers' households, when we examine the conditional ones. Table (ref) shows the log of marginal likelihoods for the unconditional MLN, conditional MLN on $R=2$, and GB2 distributions for both datasets. From the table, we can confirm that the unconditional and conditional MLN distributions are preferred to the GB2 distribution. Therefore, it allows us to interpret the parameters in an economically meaningful way using the results from the conditional MLN distribution on $R=2$. The results from two-or-more person households and for workers' households suggest that there are two groups: a lower and higher income group. In other words, Japanese households are divided into two groups both for two-or-more person households and for the workers' households, when the MLN distribution is assumed. However, the interpretation of the parameters is different in each dataset.
Table (ref) shows the posterior estimates both for two-or-more person households and for workers' households. From the table, we can make the following observations: that 35.5% households belong to the lower income group ($r=1$) in two-or-more person households from $\text{\boldmath{$\pi$}}$, whereas 64.8% households belong to the lower income group in workers' households; if we focus on the posterior estimates of $\text{\boldmath{$\mu$}}$, the posterior estimate $\mu_{2}$ of two-or-more person households is close to that of workers' households, whereas the posterior estimate $\mu_{1}$ of two-or-more person households is much smaller than that of workers' households; and the posterior estimate $\sigma_{1}^{2}$ of two-or-more person households is smaller than $\sigma_{2}^{2}$, whereas $\sigma_{1}^{2}$ of workers' households is larger than $\sigma_{2}^{2}$. As is described in the data explanation, workers' households are a subset of two-or-more person households. Therefore, we can guess that the households, which are not included in the workers' households, make the meaning of the components different. Even if the two components in workers' households are interpreted as lower and higher income class, the components wherein the households belong in workers' households might be different from the components wherein the households belong in two-or-more person households, because most of the estimates are different between two-or-more person and workers' households. In addition to these interpretations, these differences may lead to the difference in the shape of income distribution.
Figure (ref) shows the histogram from the data, predictive distributions for the unconditional and conditional MLN distributions, and GB2 distribution both for two-or-more person and workers' households. First, we can confirm that the predictive distributions for the unconditional and conditional MLN distributions both from two-or-more person and workers' households are very similar. They seem to fit to the histograms, whereas those from the GB2 distribution are different from those from the unconditional and conditional MLN distributions, especially in the case of two-or-more person households. Therefore, the goodness-of-fit of the MLN distribution model seems to be adequate for the Japanese income data.
Finally, as the predictive distributions for both datasets seem to fit to the histograms, we are also estimate the Gini coefficients as is discussed in the previous section. This is because the Gini coefficients are sometimes used for policy making and related matters. To examine the features of the Gini coefficients, the estimated Gini coefficients are shown in Table (ref) with the posterior distributions shown in Figure (ref). At first, we can confirm that the Gini coefficient from two-or-more person households is larger than that of workers' households. This might be caused by the additional households, which do not appear in the workers' households. Secondly, the 95% credible intervals do not overlap between the unconditional and conditional MLN distribution and GB2 distribution for two-or-more person households, whereas they overlap in the case of workers' households. This may suggest that the choice of the distribution in two-or-more person households is more pronounced than that of workers' households. Furthermore, if we assume the GB2 distribution as the hypothetical income distribution, the Gini coefficients are overestimated both in two-or-more person and workers' households. However, the log of marginal likelihood of the GB2 distribution is much smaller than that of the MLN distribution. Therefore, we can avoid such an overestimation, if the marginal likelihoods are appropriately utilized.
There is a strong argument for employing a reversible jump MCMC algorithm for the MLN distribution model with an unknown number of components from grouped data. Based on the simulated data examples, our proposed algorithm worked well in terms of fitting the distribution and enabled us to calculate the Gini coefficient accurately. The unconditional MLN distribution model is useful if we are interested in the fit of the income distribution, whereas the conditional MLN distribution model is useful if we are interested in an economically meaningful interpretation. A major strength of the reversible jump MCMC algorithm is that it can provide both results simultaneously in one estimation. This, along with the ability of the marginal likelihood to choose an appropriate distribution, makes it an algorithm of choice for estimating the MLN model.
The robust results support the case for using the MLN distribution model to compare other candidate distributions. Finally, using FIES datasets in 2020, the income distributions and inequalities in Japan were examined. The results indicated two subgroups, both in two-or-more person households and in workers' households. However, the meanings of the two subgroups might be different in each dataset. We also observed that the Gini coefficient of two-or-more person households are larger than that of workers' households. Moreover, if we calculate the Gini coefficients from the GB2 distribution, the Gini coefficients are overestimated.
Finally, we discuss the remaining issue. Although a reversible jump MCMC algorithm for grouped data is considered to determine the number of components, more sophisticated algorithms, which can determine the number of components, are proposed, for example, by MFG16. We need to examine more efficient algorithm, but our finding that a reversible jump MCMC algorithm can identify the number of components correctly even from grouped data, represents an interesting first step.