Efficiency of maximum likelihood estimation for a multinomial distribution with known probability sums

TL;DR

This paper analyzes how knowing the sum of certain categories' probabilities in a multinomial distribution affects the efficiency of maximum likelihood estimation, revealing limited benefits and potential risks with small samples.

Contribution

It provides an asymptotic expansion of the MLE risk in both full and submodels, highlighting the limited improvement and possible increased risk with small samples.

Findings

Risk reduction is small in some cases when using the submodel.

Using the submodel can increase risk with small sample sizes.

Asymptotic expansion of MLE risk for both models is derived.

Abstract

For a multinomial distribution, suppose that we have prior knowledge of the sum of the probabilities of some categories. This allows us to construct a submodel in a full (i.e., no-restriction) model. Maximum likelihood estimation (MLE) under this submodel is expected to have better estimation efficiency than MLE under the full model. This article presents the asymptotic expansion of the risk of MLE with respect to Kullback--Leibler divergence for both the full model and submodel. The results reveal that, using the submodel, the reduction of the risk is quite small in some cases. Furthermore, when the sample size is small, the use of the subomodel can increase the risk.

Tables

Table 1. Table 1 : Example 1-1—Risk and R.S.S.—

n	f.risk.app	s.risk.app	ratio.app	r.s.s.app	f.risk.sim	s.risk.sim	ratio.sim	r.s.s.sim
100	1.3283	1.3332	1.0037	100	1.004589(5.3e-6)	1.006799(5.3e-6)	1.0022	101
200	0.5808	0.5808	0.9999	200	0.571158(4.1e-6)	0.571305(4.1e-6)	1.0003	201
300	0.3687	0.3681	0.9985	300	0.384620(3.3e-6)	0.384268(3.3e-6)	0.9991	300
400	0.2696	0.2690	0.9977	399	0.283346(2.6e-6)	0.282871(2.6e-6)	0.9983	400
500	0.2123	0.2117	0.9972	499	0.221933(2.2e-6)	0.221424(2.2e-6)	0.9977	499
1000	0.1028	0.1024	0.9961	996	0.104093(1.1e-6)	0.103697(1.1e-6)	0.9962	996
2000	0.0506	0.0504	0.9955	1991	0.050657(5.1e-7)	0.050430(5.1e-7)	0.9955	1992

Table 2. Table 2 : Example 1-2—Risk and R.S.S.—

n	f.risk.app	s.risk.app	ratio.app	r.s.s.app	f.risk.sim	s.risk.sim	ratio.sim	r.s.s.sim
100	1.9209	1.9609	1.0208	101	0.800707(6.9e-6)	0.804596(5.7e-6)	1.0049	101
200	0.7290	0.7377	1.0120	202	0.458571(4.1e-6)	0.460303(3.9e-6)	1.0038	201
300	0.4345	0.4379	1.0077	302	0.322643(2.9e-6)	0.323684(2.9e-6)	1.0032	301
400	0.3066	0.3082	1.0051	402	0.250190(2.1e-6)	0.250864(2.3e-6)	1.0027	402
500	0.2360	0.2368	1.0034	501	0.205073(1.7e-6)	0.205388(1.9e-6)	1.0015	502
1000	0.1088	0.1087	0.9995	1000	0.107846(9,2e-7)	0.107870(1.0e-6)	1.0002	1000
2000	0.0521	0.0519	0.9974	1955	0.053458(5.1e-7)	0.053359(5.5e-7)	0.9981	1997

Table 3. Table 3 : Breast cancer classification

	30–39	40–49	50–59	60–69	70–79
1	0.025	0.063	0.088	0.060	0.014
2	0.060	0.168	0.137	0.084	0.004
3	0.042	0.084	0.112	0.056	0.004
Column sum	0.126	0.317	0.337	0.200	0.021

Table 4. Table 4 : Example 2—Risk and R.S.S.—

	$n=200$	$n=400$	$n=600$	$n=800$	$n=1000$
Full Model	0.0367(0.0361)	0.0179(0.0180)	0.0119(0.0119)	0.0089(0.0090)	0.0071(0.0071)
Submodel	0.0281(0.0277)	0.0133(0.0135)	0.0087(0.0088)	0.0064(0.0066)	0.0051(0.0052)
Risk Ratio	0.766(0.769)	0.741(0.752)	0.732(0.741)	0.728(0.734)	0.725(0.730)
R.S.S.	158(155)	302(305)	445(450)	588(591)	732(739)
R.S.S./$n$	0.790(0.775)	0.755(0.762)	0.742(0.750)	0.735(0.739)	0.732(0.739)

Table 5. Table 5 : Household classification

	H1	H2	H3	H4	H5	H6
Y1	0.00161	0.00232	0.00512	0.00395	0.00468	0.00066
Y2	0.00331	0.0081	0.00953	0.00783	0.01145	0.00278
Y3	0.00974	0.02109	0.02046	0.01499	0.02536	0.00494
Y4	0.00799	0.03519	0.03229	0.02017	0.0338	0.00708
Y5	0.00547	0.0376	0.04362	0.02442	0.02675	0.00398
Y6	0.00494	0.05082	0.09003	0.05772	0.03732	0.00452
Y7	0.00126	0.02106	0.0543	0.05531	0.01999	0.00234
Y8	0.00071	0.00961	0.0323	0.04043	0.0108	0.00122
Y9	0.00011	0.00201	0.01204	0.02184	0.00466	0.00052
Y10	0.00006	0.00139	0.00697	0.01582	0.00344	0.00022
Col. Sum	0.03520	0.18919	0.30666	0.26248	0.17825	0.02826

Table 6. Table 6 : Example 3—Risk and R.S.S.—

	$n=1000$	$n=1500$	$n=2000$	$n=2500$	$n=3000$
Full Model	0.0372(0.0295)	0.0231(0.0197)	0.0167(0.0148)	0.0130(0.0119)	0.0107(0.0099)
Submodel	0.0351(0.0275)	0.0216(0.0183)	0.0155(0.0137)	0.0121(0.0109)	0.0099(0.0091)
Risk Ratio	0.945(0.931)	0.937(0.926)	0.932(0.924)	0.929(0.923)	0.927(0.922)
R.S.S.	955(934)	1419(1392)	1880(1848)	2339(2305)	2799(2773)
R.S.S./$n$	0.955(0.934)	0.946(0.928)	0.940(0.924)	0.936(0.922)	0.933(0.924)

Equations

$$D_{\Psi }[\hat{\theta }|{\theta }^{\ast }]=D_{\Psi }[\hat{\theta }|{\hat{\theta }}^{\ast }]+D_{\Psi }[{\hat{\theta }}^{\ast }|{\theta }^{\ast }],$$

$${\hat{\theta }}^{\ast }=\underset{\theta \in {\mathcal{M}}^{\ast }}{\mathrm{argmin}}{D}_{\Psi }[\hat{\theta }|\theta ].$$

$$E[D_{\Psi }[\hat{\theta }|{\theta }^{\ast }]]\ge E[D_{\Psi }[{\hat{\theta }}^{\ast }|{\theta }^{\ast }]],$$

$$m\triangleq (m_0,\dots ,m_p),m_i\triangleq P(X\in C_i)>0,i=0,1,\dots ,p.$$

$$\sum _{i=0}^p{m}_i=1,$$

$$m_{i_1}+\cdots +m_{i_s}=c,0\le i_1<i_2<\dots <i_s\le p$$

$$\hat{m}\triangleq ({\hat{m}}_0,\dots ,{\hat{m}}_p),$$

$$D[\hat{m}:m]=\sum _{i=0}^p{\hat{m}}_i\log \frac{{\hat{m}}_i}{m_i}$$

$$ED=ED[\widehat{m}:m]\triangleq E\bigl{[}D[\widehat{m}:m]\bigr{]}.$$

$$\bar{\mathcal{M}}\triangleq\bigl{\{}m=(m_{0},\ldots,m_{p})\,|\,m_{i}>0,\ i=0,\ldots,p\bigr{\}},$$

$${\theta }^i\triangleq \log m_i+1,i=0,\dots ,p$$

$$\Psi (\theta )\triangleq \sum _{i=0}^p{e}^{{\theta }^i-1}=\sum _{i=0}^p{m}_i.$$

$$\frac{\partial }{\partial {\theta }^i}\Psi (\theta )=e^{{\theta }^i-1}=m_i,i=0,\dots ,p,$$

$$\begin{array}{rl}{D}_{\Psi }[{\theta }_1|{\theta }_2]& =\Psi ({\theta }_1)-\Psi ({\theta }_2)-\sum _{j=0}^p{m}_j^2({\theta }_1^j-{\theta }_2^j)\\ & =\sum _{i=0}^p{m}_i^1-\sum _{i=0}^p{m}_i^2+\sum _{i=0}^p{m}_i^2\log (m_i^2/m_i^1)\end{array}$$

$$D_{\Psi }[{\theta }_1|{\theta }_2]=\sum _{j=0}^p{m}_j^2\log (m_j^2/m_j^1).$$

$$D[{\theta }_2|{\theta }_1]=D_{\Psi }[{\theta }_1|{\theta }_2].$$

$$ED=\frac{p}{2n}+\frac{1}{12n^2}(M-1)+o(n^{-2}),$$

$$M\triangleq \sum _{i=0}^p{m}_i^{-1}.$$

$$m_{ij}\triangleq P(X\in C_{ij})(>0),m_{i\cdot }\triangleq \sum _{j=1}^{J_i}{m}_{ij},p_{ij}\triangleq \frac{m_{ij}}{m_{i\cdot }}$$

$$D[X_{1}:X_{2}]=D[U_{1}:U_{2}]+E_{U_{1}}\bigl{[}D[(X|U)_{1}:(X|U)_{2}]\bigr{]},$$

$$x_{i\cdot }=\sum _{j=1}^{J_i}{x}_{ij},\sum _{i=1}^I{x}_{i\cdot }=n.$$

$$m_f\triangleq (m_{1\cdot },\dots ,m_{I\cdot })$$

$$i(1\le i\le I)$$

$${\hat{m}}_{ij}={\hat{m}}_{i\cdot }{\hat{p}}_{ij},i=1,\dots ,I,j=1,\dots ,J_i,$$

$${\hat{m}}_{i\cdot }=m_{i\cdot }(\hat{\tau }),{\hat{p}}_{ij}=p_{ij}({\hat{\theta }}^{(i)}),i=1,\dots ,I,j=1,\dots ,J_i,$$

$$\log f(x;m)$$

$$\displaystyle=\sum_{i=1}^{I}\Bigl{(}\sum_{j=1}^{J_{i}}x_{ij}\bigl{(}\log m_{ij}/m_{i\cdot}+\log m_{i\cdot}\bigr{)}\Bigr{)}+c(x)$$

$$\displaystyle=\sum_{i=1}^{I}\Bigl{(}\sum_{j=1}^{J_{i}}x_{ij}\bigl{(}\log p_{ij}(\theta^{(i)})+\log m_{i\cdot}(\tau)\bigr{)}\Bigr{)}+c(x)$$

$$=\sum _{i=1}^I{x}_{i\cdot }\log m_{i\cdot }(\tau )+\sum _{i=1}^I\sum _{j=1}^{J_i}{x}_{ij}\log p_{ij}({\theta }^{(i)})+c(x)$$

$$\displaystyle\Bigl{(}\triangleq\log f(x;\tau,\theta),\quad\theta=(\theta^{(1)},\ldots,\theta^{(I)})\Bigr{)}$$

$$\underset{\tau ,\theta }{\max }f(x;\tau ,\theta )=\underset{\tau }{\max }\sum _{i=1}^I{x}_{i\cdot }\log m_{i\cdot }(\tau )+\sum _{i=1}^I\underset{{\theta }^{(i)}}{\max }\sum _{j=1}^{J_i}{x}_{ij}\log p_{ij}({\theta }^{(i)}).$$

$$\sum _{i=1}^I{x}_{i\cdot }\log m_{i\cdot }(\tau )$$

$$\sum _{j=1}^{J_i}{x}_{ij}\log p_{ij}({\theta }^{(i)})$$

Full text

MSC(2010) Subject Classification: Primary 60F99; Secondary 62F12

Key words and phrases: Kullback-Leibler divergence, asymptotic risk, asymptotic expansion, multinomial distribution

1 Introduction

As Eguchi [4] , [5] showed, the geometrical properties of the manifold of a parametric distribution model can be induced by a divergence on it. Hence, if we make an analysis for a statistical model that is characterized by some certain geometrical properties (especially flatness), it seems natural that we use the corresponding divergence.

Let’s consider this issue in view of the submodel. Consider an estimation problem in a parametric model $\mathcal{M}$ with the parameter $\theta$. If we have prior knowledge on $\theta$, it gives a submodel, $\mathcal{M^{*}}$. We usually project an estimator (say $\hat{\theta}$) in $\mathcal{M}$ onto $\mathcal{M^{*}}$, and use the projected point as an alternative estimator, $\hat{\theta}^{*}$. Naturally we conjecture that $\hat{\theta}^{*}$ performs better than $\hat{\theta}$ because of the use of the prior knowledge. This is theoretically guaranteed under some conditions; For example, if we use Bregman divergence induced from a convex function $\Psi(\theta)$ (say $D_{\Psi}[\cdot\,|\,\cdot]$) and $\mathcal{M^{*}}$ is the affine subspace with respect to the dual parameter $\eta$, then the generalized Pythagorean theorem (see Amari [1]) says

[TABLE]

where $\theta^{*}$ is the parameter of the true distribution, and $\hat{\theta}^{*}$ is the projected point with respect to $D_{\Psi}[\hat{\theta}\,|\,\cdot]$, namely

[TABLE]

Consequently the risk of $\hat{\theta}^{*}$ is smaller than that $\hat{\theta}$, i.e.

[TABLE]

where the expectation is taken under the true distribution. The point is pairing of $D_{\Psi}[\cdot\,|\,\cdot]$ and the submodel that is affine with respect to its dual coordinate.

This paper investigates the risk of the submodel for the case the pair is mismatched. The framework of the paper is given as follows:

The object of the estimation is the finite dimensional discrete distribution, i.e. the multinomial distribution. 2. 2.

The prior knowledge is completely correct and is given as the sum of the probabilities of some bins (categories). 3. 3.

The estimation method is MLE and we use Kullback-Leibler divergence (K-L divergence) for the estimation evaluation.

The statistical problems under this framework are often observed in the surveys where a prior survey is available. It is an important issue whether we should adjust the result of the present survey according to the prior survey. We show two concrete examples:

– Household Survey –

   From a household sample survey, we have formulated a two-way contingency table by categorizing each household in view of the yearly income level and the age group of the main earner. The relative frequency of each category to the sample size gives us MLE estimation of the corresponding two-dimensional discrete distribution. Suppose that a national census was carried out in the previous year, then we have the highly reliable knowledge on the distribution of the age group. This means we have the prior knowledge on the sum of the probabilities for each row (or column) in the two-way table. The natural question arises whether we should modify the survey result so that the row (column) sums are equal to those of the census.

– Land Survey –

   We are interested in the state of land in a certain district. Land is categorized into three kinds of areas: non-forest land area, commercial forest land area, noncommercial forest land area. We will make a field investigation from randomly chosen points and make an estimation of the ratio among the above three areas. Prior to investigation, we already have obtained an aerial photo of the district, which informs us about the ratio between the forest and non-forest land areas. Hence, we have the prior information on the probability of one category (non-forest area) of the multinomial distribution. The same question as above arises.

Now we formulate the estimation problem. First we consider the multinomial distribution as follows; In each independent trial, the random variable $X$ takes a value which belongs to one of the categories $C_{i}\ (i=0,1,\ldots,p)$ with the probability

[TABLE]

Due to the restriction

[TABLE]

the dimension of this probability model is equal to $p$. The model with no other restrictions than (3) will be called ”full model.”

Suppose we have prior information on the sum of the probabilities over some categories, namely

[TABLE]

with some known constant $c$. If we have a collection of such restriction equations as (5), it formulates a submodel (hereafter we call it as ”$m$-aggregation submodel”). If each parameter $m_{i}$ appears no more than once in all restriction equations, we call the model “non-overlapping.” If some $m_{i}$’s appear multiple times in the restriction, the model is referred to as “overlapping.” For example, if we only know the column sums of a two-way contingency table, the submodel is non-overlapping. If both row- and column-wise sums are known, an overlapping submodel is formulated. Especially if we know the probability of one category, that is, $s=1$ in (5), we call the restriction (or information) is “solid.”

Hereafter we denote the MLE of the parameter $m$ (in the full model or submodel according to the context) by

[TABLE]

and measure the discrepancy between the true distribution and the predictive distribution using K-L divergence,

[TABLE]

Note that the MLE for the submodel is given by the projection onto the submodel of the MLE of the full model. (See [6] for the use of Kullback–Leibler divergence in the inference for the multinomial distribution.)

We evaluate the performance of MLE $\widehat{m}$ through the risk, that is

[TABLE]

Here we show that the $m$-aggregation submodel and K-L divergence are “mismatched”. Treat $m$ as an element of the set of $p+1$-dimensional positive measures, i.e.

[TABLE]

where the sum is not necessarily equal to one. Let

[TABLE]

and

[TABLE]

Then $\Psi(\theta)$ is the convex function of $\theta=(\theta^{0},\ldots,\theta^{p})$, and as

[TABLE]

$m$ is the dual coordinate of $\theta$. The divergence induced from $\Psi(\theta)$ on $\bar{\mathcal{M}}$ equals

[TABLE]

for $\theta_{i}=(\theta_{i}^{0},\ldots,\theta_{i}^{p}),$ and their dual coordinates $m^{i}=(m^{i}_{0},\ldots,m^{i}_{p}),\ i=1,2$. If the two points $\theta_{1}$ and $\theta_{2}$ are probabilities, i.e. $\sum_{j=0}^{p}m^{i}_{j}=1$ for $i=1,2$, then

[TABLE]

This is the conjugate of K-L divergence, and not K-L divergence itself. Actually

[TABLE]

Since the $m$-aggregation submodel is an affine subspace in $\bar{\mathcal{M}}$ with respect to $m$, if we use $D_{\Psi}[\cdot\,|\,\cdot]$, Pythagorean theorem guarantees the smaller risk of $\hat{\theta}_{*}$ than $\hat{\theta}$. However if we use K-L divergence, this is not theoretically guaranteed.

In this paper, the asymptotic expansion of risk is derived for non-overlapping $m$-aggregation submodel with respect to the sample size $n$ up to the second-order term (Section 2.1). The first-order term ($n^{-1}$-order term) and the second-order term ($n^{-2}$-order term) provide information on the asymptotic efficiency of the MLE for the submodel.

Sheena [7] derived the asymptotic expansion of the risk of MLE with regard to $\alpha$-divergence for a general parametric model that is well-specified (namely the model that contains the true distribution). (See Theorem 1 of [7].) As an application to the full model multinomial distribution, it demonstrated (see (42) of [7]) that

[TABLE]

where

[TABLE]

The first-order term is proportional to the ratio of the model’s dimension to the sample size (“$p-n$ ratio”), which indicates (relative) complexity of the model. Since this holds true for any parametric model, the first-order term for any well-specified submodel is smaller than that of the full model, which is the benefit of the dimension reduction for the submodel. Consequently when the sample size is large enough, the risk of the submodel is smaller than that of the full model, and as the sample size goes to infinity, the risk ratio between the sub and full models converges to the dimensions ratio.

However as the result of this paper (Corollary 1) shows, the second-order term of the non-overlapping $m$-aggregation submodel is always larger than that of the full model unless the prior information is solid. This indicates if the sample size $n$ is not large enough to neglect the second-order term, the risk reduction of the submodel does not become as small as we expect from the dimension reduction, or even anomaly, reversed risk order, could happen. In Section 3, simulation studies demonstrate that these phenomena actually happen.

Improvement of the inference efficiency using prior knowledge is an important theme also in hypothesis testing. For this issue, we refer to Trosset et. al. [9]. They treated Marden’s conjecture that the additional information on the alternative hypothesis increase the power of the likelihood test. They showed that this conjecture fails in some concrete distributions even when the information reduces the dimension in the alternative hypothesis.

2 Risk of non-overlapping $m$-aggregation submodel

2.1 Two-stage multinomial distribution models

For the risk analysis of a non-overlapping submodel, it is convenient to decompose each trial (experiment) into two stages. In this subsection, we consider the “two-stage multinomial distribution model’.

We begin by defining the two-stage multinomial distribution. Suppose that the random variable $X$ takes values that belong to one of the categories $C_{ij}$ and that

[TABLE]

for $i=1,\ldots,I,\ j=1,\ldots,J_{i}$. (To eliminate a trivial case, we suppose that $I\geq 2$.) Let $C_{i\cdot}\triangleq\bigcup\limits_{j}C_{ij}$. Then, the first-stage model is given by focusing on which $C_{i\cdot}$ the value of $X$ belongs to. Its parameters are given by $m_{i\cdot}\ (i=1,\ldots,I)$. The $i$th model ($i=1,\ldots,I$) in the second stage is given by the categories $C_{ij}$ and the corresponding probabilities $p_{ij}$ ($j=1,\ldots,J_{i}$) under the condition $X\in C_{i\cdot}$.

If the restriction for the submodel is non-overlapping, the parameters $m_{ij}$ can be grouped by the restriction equation in which the parameter appears (the parameters that do not appear in any restrictions form one group). The full model is decomposed into two stages according to this grouping, and the full model can be considered as a two-stage model.

Kullback–Leibler divergence satisfies the so-called “chain rule.” Let $(X|U)$ denote the conditional distribution of $X$ when $U$ is given. Suppose a pair of the distributions $(X|U)_{i}$ and $U_{i}$ defines the distribution of $X_{i}$ for $i=1,2$. Then, the following relationship (chain rule) holds.

[TABLE]

where $E_{U_{1}}[\cdot]$ is the expectation under the condition $U\stackrel{{\scriptstyle d}}{{=}}U_{1}$. (For more details on divergences, see Vajda [10] and Amari and Nagaoka [2].)

Thanks to this property, we can decompose the MLE risk into several parts, each of which corresponds to the first- and second-stage distributions.

When a sample of size $n$ is taken from the two-stage multinomial distribution, let $x_{ij}$ and $x_{i\cdot}$ $(i=1,\ldots,I,j=1,\ldots,J_{i})$ denote the number of individuals that belong to $C_{ij}$ and $C_{i\cdot}$, respectively, hence

[TABLE]

The notations of $X_{ij}$ and $X_{i\cdot}$ are similarly defined as the corresponding random variables.

In this section, we always assume that $m_{i\cdot}$ and $p_{ij}$ are parameterized independently, that is,

[TABLE]

We obtain the following result.

Theorem 1**.**

The MLE for the two-stage multinomial distribution model, $\widehat{m}=(\widehat{m}_{ij})$, is given by

[TABLE]

with

[TABLE]

where $\widehat{\tau}$ is the MLE for the first-stage model based on $x_{1\cdot},\ldots,x_{I\cdot}$, and for each $i(i=1,\ldots,I)$, $\widehat{\theta}^{(i)}$ is the MLE for the second-stage model based on $x_{i1},\ldots,x_{iJ_{i}}$.

Proof.

For the sample $x=(x_{ij})$, the log-likelihood is expressed as

[TABLE]

We obtain

[TABLE]

We notice that

[TABLE]

is the log-likelihood function of $\tau$ up to constant for the first-stage model based on $x_{1\cdot},\ldots,x_{I\cdot}$, whereas, for each $i(1\leq i\leq I)$, the log-likelihood function of $\theta^{(i)}$ for the second-stage model is given by

[TABLE]

up to constant with the sample $x_{i1},\ldots,x_{iJ_{i}}$. Since $m_{ij}=m_{i\cdot}(\tau)p_{ij}(\theta^{(i)})$, all the results are obtained. ∎

The following decomposition of $D[\widehat{m}:m]$ holds;

[TABLE]

where

[TABLE]

When the first-stage model is full, the MLE of $m_{i\cdot}$ is given by $\widehat{m}_{i\cdot}=X_{i\cdot}/n,\ i=1,\ldots,I$, and the asymptotic expansion of the risk of the MLE for the first-stage model (denoted by $ED_{f}$) with the sample size $n$ is given by (see (12))

[TABLE]

where

[TABLE]

For the $i$-th model ($i=1,\ldots,I$) at the second stage (not necessarily a full model), let $\overset{-1}{A}_{(i)}/(24n^{2})$ denote the second-order term of $ED_{(i)}\triangleq E[D[\widehat{p}_{i}:p_{i}]]$ with the sample size $n$. Namely, we have

[TABLE]

Because the dimension of the $i$-th model at the second stage equals $s_{i}$ due to (17), the first-order term equals $s_{i}/(2n)$ (see Theorem 1 of [7]).

It should be noted that it is possible that $X_{i\cdot}=0$. Hence, we can encounter a situation in which we are unable to estimate $p_{ij},\ j=1,\ldots,J_{i}$, the parameter of the $i$-th second-stage model, because there is no available sample. We overcome this problem by making it a rule to discard such samples with no estimation. In the following theorems, all the expectations are conditional on the state $X_{i\cdot}\neq 0,\ 1\leq\forall i\leq I$. However, as Lemma 1 in the Appendix shows, the conditional expectations of Kullback–Leibler divergence and $(X_{i\cdot}-nm_{i\cdot})^{k}(k=1,2,\ldots)$ differ from those that are unconditional by $o(n^{-s})$ for any $s>0$. Therefore, we use the same notation as that for the unconditional distribution.(Avoiding zero probability estimates is a practically important issue. See Darsheid et.al. [3].)

The next theorem is Theorem 1 of Sheena [8]. It provides the decomposition of $ED$ when the first-stage model is a full model.

Theorem 2**.**

If the first-stage model is a full model, then the risk of MLE is equal to

[TABLE]

where

[TABLE]

2.2 Asymptotic risk expansion for the submodel

Suppose that we have prior knowledge on $\tau$ in (16), equivalently $m_{f}$. The next theorem gives the risk of MLE for this situation.

Theorem 3**.**

If the first-stage model is a full model and $\tau_{i}\ (i=1,\ldots,I)$ are all known, then the risk of MLE is

[TABLE]

Especially, if all the second-stage models are full models, then the risk of MLE is

[TABLE]

Proof.

Since $\widehat{m}_{i\cdot}=m_{i\cdot},\ i=1,\ldots,I$, that is, $\widehat{m}_{f}=m_{f}$, we have

[TABLE]

(see (19)).

[TABLE]

where

[TABLE]

From (22), substituting $n$ with $X_{i}$, we have

[TABLE]

where

[TABLE]

Therefore we have

[TABLE]

Taylor expansions of $1/X_{i\cdot}$ and $1/X^{2}_{i\cdot}$ are given by

[TABLE]

with suitable choice of $m^{*}_{i\cdot}\leq 1$ for each expansion. Using $E[X_{i\cdot}-nm_{i\cdot}]=0$ and $E[(X_{i\cdot}-nm_{i\cdot})^{2}]=nm_{i\cdot}(1-m_{i\cdot})$ and the asymptotic distribution

[TABLE]

we have

[TABLE]

and

[TABLE]

From these evaluations, we have

[TABLE]

Consequently, we have

[TABLE]

When all the second-stage models are full,

[TABLE]

If we insert these results into (24), we have

[TABLE]

and

[TABLE]

∎

The prior knowledge on the values of $m_{f}$ is equivalent to the $m$-aggregation submodel

[TABLE]

where $c_{i}$’s are known constants. It is clear that any non-overlapping $m$-aggregation submodel can be treated as a two-stage model in which every second-stage model is a full model. Consequently, the risk for the non-overlapping $m$-aggregation submodel is given by (25).

We now compare the risk of (23) and that of the (24). Notice that the difference between Theorem 2 and Theorem 3 is the existence of the prior information on $M_{f}$. We call the models in Theorem 2 and 3 “M1” and “M’2,” respectively. If the second-stage models are all full, then M1 becomes the full model and M2 becomes the $m$-aggregation submodel.

If we neglect the $o(n^{-2})$ terms, the difference between (23) and (24) is equal to

[TABLE]

especially, when the second-stage models are all full models,

[TABLE]

Regarding the first-order terms, that of M2 is always smaller than that of M1. The risk ratio between the full and submodels is close to the dimension ratio

[TABLE]

under a large enough sample size.

However, the second-order term of M2 can be larger than that of M1. Suppose that $s_{i}\geq 1(i=1,\ldots,I)$, then

[TABLE]

the second inequality holds due to the rule that the arithmetic mean is not less than the harmonic mean

[TABLE]

When the second-stage models are all full, $s_{i}=0$ is equivalent to $J_{i}=1$, namely

[TABLE]

in (30). We call the restriction equation of this type “solid” as it provides solid information on the probability of a particular cell. The aforementioned result is applied to the $m$-aggregation model as a corollary.

Corollary 1**.**

If none of the restriction equations of the $m$-aggregation submodel are solid, then the second-order term of the submodel is larger than that of the full model.

This means that (32) is negative for some small values of $n$ when the restriction equations are all non-solid. To simplify the explanation, we roughly classify the sample size as follows: “small sample size” when (32) is negative; “medium sample size” when (32) is positive, but the second-order term of (32) is still non-negligible; “large sample size” when the second-order term of (32) is negligible.

Our conclusion is as follows:

• •

When $n$ is “small,” the rather pathological situation may occur that the submodel (i.e., the prior information) increases the risk, and hence, it is better not to use the submodel.

• •

When $n$ is “medium,” the submodel has an advantage over the full model. However, the submodel loses the estimation efficiency gained by the dimension reduction due to the larger second-order term. That is, the ratio of the risks between the sub and full models is not as small as the dimension ratio.

• •

When $n$ is “large,” the risk ratio between the sub and full models is close to the dimension ratio. It should be noted that the dimension ratio is close to one when the number of restriction equations in the submodel ($=I$) is quite small compared with the dimension of the full model ($=\sum_{i=1}^{I}J_{i}-1$).

3 Numeric analysis of some examples

In the previous subsection, we observed the risk difference between the full model and the non-overlapping $m$-aggregation submodel through the approximated risks, that is, the asymptotic expansion of the risks up to the second order. In this subsection we confirm several results using simulation with three examples. Every example is a two-way contingency table and the non-overlapping $m$-aggregation submodel is given as the one in which the column sums are all known. As a two-stage model, the first-stage model consists of the multinomial distribution over the columns and each second-stage model is the distribution over the rows within a given column.

To compare the submodel and full model, we also use the indicator, “the required sample size (r.s.s) of the submodel to the full model under the condition $n_{0}$,” which is defined by the solution $n^{*}$ of the following equation

[TABLE]

where $ED_{s}(n)$ and $ED_{f}(n)$ are the risks of the sub and full model considered as the function of the sample size, respectively. This reveals the required sample size for the submodel is equal in terms of risk with the full model of the sample size $n_{0}$. The risks in the equation above are calculated by approximation or simulation.

The following abbreviations are commonly used in the three examples;

“f.risk.sim(app)” is the risk of the full model obtained by simulation(approximation).

“s.risk.sim(app)” is the risk of the submodel obtained by simulation(approximation).

“ratio.sim(app)” is the risk ratio between the submodel and the full model based on the simulated(approximated) risks.

“r.s.s.sim(app)” is the r.s.s obtained by simulation(approximation).

–Example 1–

The first example is an artificial setting to confirm some theoretical results in the previous subsection. We consider a 100 by 2 contingency table for the two cases of contrasting distributions. For both cases, $I=2,\ J_{1}=J_{2}=100,$ where in Example 1-1,

[TABLE]

and in Example 1-2,

[TABLE]

These distributions are chosen from the following consideration; Since the second-order term in the equation (32) equals

[TABLE]

which is maximized when $m_{1\cdot}=m_{2\cdot}=1/2$ (Example 1-1) and gets smaller when $m_{1\cdot}=1/10,\ m_{2\cdot}=9/10$ (Example 1-2).

The risks and r.s.s. under several values of $n$ are given in Table 1 for Example 1-1 and Table 2 for Example 1-2. The number in the parenthesis after the simulated risk is the standard deviation. The r.s.s is under the condition $n_{0}=n$. For the calculation of the simulated risk, we used $10^{4}$ sets of samples and took the average over these.

We observe the following results:

• •

For $n=100\text{ or }200$ in Ex.1-1 and all $n$’s except $n=2000$ in Ex.1-2, the use of the submodel increases the risk.

• •

Since the dimension ratio ($=198/199\doteqdot 0.9950$) is close to one, the use of the submodel causes little risk reduction. In view of r.s.s., the contribution of the submodel is almost negligible.

• •

The approximated values for the risk are somewhat deviated from the simulated values when $n$ is small; however, those for the risk ratio or r.s.s. are quite close to the simulated values even when $n$ is small.

–Example 2–

We use real data on breast cancer taken from the ”UCI machine learning repository” (https://archive.ics.uci.edu/ml/datasets/Breast+Cancer). We created a cross-tabulation table from the variable “the degree of malignancy” (3 levels:1,2,3) and “age group” (5 groups: 30–39,…,70–79) excluding one person in his or her twenties from the original data set. We gained the relative frequency by dividing each cell by the total number of individuals ,285(see Table 3). We assume that this is the true probability $m_{ij},\ i=1,\ldots,3,j=1,\ldots,5$ for each category. We consider the submodel according to the situation in which we have prior knowledge of the sums of each column, that is, the distribution over the age groups.

Table 4 reveals the approximated and simulated risk (in parentheses) for the full and submodels for several cases of sample size $n$. It also shows the approximated and simulated r.s.s. of the submodel to the full model under the condition $n_{0}=n$ (the simulated r.s.s is in parentheses).

We summarize the result as follows:

• •

For this example, when $n\leq 28$, the value of (31) (equivalently (32)) becomes negative. However, the simulation reveals that the risk of the submodel is always smaller than that of the full model. In fact, even if $n=10$, the former equals 0.6236 and the latter equals 0.6543.

• •

The dimension ratio is $10/14\doteqdot 0.714$. Since the table is small, the knowledge of the column sums significantly reduces the risk. However, we notice that the effect of the dimension reduction is lessened by the larger second-order term of the submodel as the “Risk Ratio” or “R.S.S./$n$” is larger than 0.714.

• •

The approximation method for risk calculation is effective under the sample sizes in the table.

–Example 3–

We use data from the “2014 National Survey of Family Income and Expenditure” by the Statistics Bureau in Japan (https://www.stat.go.jp/english/data/zensho/index.html). Table 5 is the relative frequency obtained from the classification of 100,006 households according to “Yearly income group (Y1,…,Y10)” and “Household age group (H1,…,H6).” We use this relative frequency as the population parameter $m_{ij},\ 1\leq i\leq 10,\ 1\leq j\leq 6$. We consider the submodel that is based on the prior knowledge of all the column sums. Table 6 presents the results.

We summarize the results as follows:

• •

For this example, when $n\leq 133$, the value of (31) (equivalently (32)) becomes negative. The simulation result (not presented in Table 6 ) shows that the pathological phenomenon (the reversed order in the risks) actually occurs, but only when $n$ is as small as 20.

• •

The dimension ratio is $54/59\doteqdot 0.915$. The table is larger than that of Example 2, and the knowledge of the column sums is not as useful in terms of reducing the dimension. The larger second-order term of the submodel is also a burden on risk reduction, which increases the “Risk Ratio” or “R.S.S./$n$” to over 0.915.

• •

The approximation method for risk calculation works effectively under the sample sizes in the table.

4 Submodel and Divergence

In this closing section, some complementary issues related to the pair of submodel and divergence are discussed. As explained in Introduction, it is natural to use the conjugate K-L divergence (11) when we treat the affine subspace of $\bar{\mathcal{M}}$ with respect to $m$. Let’s call the submodel defined by some linear restrictions in $m$, “$m$-submodel”. $m$-submodel is always a $m$-affince subspace of $\bar{\mathcal{M}}$, since the restriction $\sum_{i=0}^{p}m_{i}=1$ is also linear in $m$.

Now let $\mathcal{M}^{*}$ denote a $m$-submodel, and let $\hat{m}$ be MLE for the full model, then we can consider two projections of $\hat{m}$ onto the $\mathcal{M}^{*}$: One is the projection by K-L divergence

[TABLE]

The other is the projection by the conjugate K-L divergence

[TABLE]

In general, $\hat{m}^{+}$ is not equal to $\hat{m}^{*}$, however in the case of the non-overlapping $m$-aggregation submodel, $\hat{m}^{+}=\hat{m}^{*}$. This can be proved as follows; Suppose $\mathcal{M}^{*}$ is defined by the following $I$ restriction equations (including the equation $m_{0}+\cdots+m_{p}=1$),

[TABLE]

Then by the method of Lagrange multiplier, $\hat{m}^{+}$ is the solution of

[TABLE]

or equivalently

[TABLE]

For $i\neq j$, if

[TABLE]

then

[TABLE]

In the case of the non-overlapping $m$-aggregation model, for any $m_{i}$ and $m_{j}$ in the same restriction equation (5), (34) holds, hence (35) holds. Therefore $\hat{m}^{+}_{i}$ that appears in the restriction (5) is given by

[TABLE]

Similarly $\hat{m}^{*}$ is the solution of

[TABLE]

or equivalently

[TABLE]

For the non-overlapping m-aggregation submodel, the same argument as $\hat{m}^{+}$ holds, and $\hat{m}^{*}_{i}$ that appears in the restriction (5) is also given by

[TABLE]

For the overlapping submodel, $\hat{m}^{+}$ and $\hat{m}^{*}$ do not necessarily coincide, and are not simply expressed like (37) or (36). Besides, its risk (the asymptotic expansion up to the second-order) needs to be calculated by another method than the two-stage model.

Another important class of the submodels is the model given by the linear restriction in $\theta$ of (9) (say “$e$-submodel”). For example, in the two-way contingency table, if we assume each row (or column) is independently distributed, then this model is an $e$-submodel. In the two by two contingency table, the fixed odds-ratio produces an $e$-submodel. Notice that $e$-submodel is not an affine subspace of $\bar{\mathcal{M}}$, since the restriction $\sum_{i=0}^{p}m_{i}=1$ is not the linear restriction in $\theta$. Hence the use of the projection theorem and the related Pythagorean theorem is not straightforward.

Acknowledgment

This work was supported by JSPS KAKENHI Grant Number JP20K11706.

Declarations

On behalf of all authors, the corresponding author states that there is no conflict of interest.

Appendix

Lemma 1**.**

Let $X=(X_{0},\ldots,X_{p})$ be the random vector whose distribution is defined as the multinomial distribution with (3) and the sample size $n$. The distribution under the condition $X_{i}\neq 0,0\leq\forall i\leq p$ is considered. Let the unconditional and conditional expectations of a random variable $Y(X)$ be denoted by $E[Y(X)]$ and $E^{*}[Y(X)]$, respectively. If

[TABLE]

holds with some nonnegative numbers $a,b,c$, the difference between the two expectations decreases to zero with exponential speed as $n$ goes to infinity, namely for any $s>0$,

[TABLE]

In the special case, the following equations hold for any $s>0$;

[TABLE]

where $\widehat{m}(X)$ is the MLE of $m$.

Proof.

For $i(=1,\ldots,p+1)$ and $s(>0)$, because the following equivalence relationship holds

[TABLE]

we have

[TABLE]

Let $Z_{+}^{p+1}$ be the set of all $p+1$-dimensional vectors whose elements are nonnegative integers.

[TABLE]

Notice that

[TABLE]

which means $m^{*}=o(n^{-s})$. Because

[TABLE]

we have

[TABLE]

Choose arbitrary $s(>0)$. Because $\mathcal{X}$ and $\mathcal{X}^{*}$ are both finite sets, for some nonnegative constants $a^{\prime},a^{\prime\prime},b^{\prime},b^{\prime\prime}$,

[TABLE]

Because $m^{*}=o(n^{-t})$ for any $t(>0)$,

[TABLE]

as $n\to\infty$, which shows (38). The rest is obvious from the fact $D[\widehat{m}(X):m]$ is bounded (notice that $x\log x\to 0$ as $x\to 0$) and

[TABLE]

∎