Shrinkage estimation of Poisson means is considered when observations are given in the form of a two-way contingency table. Assuming a multiplicative Poisson model, estimators which shrink to the specified values or an order statistic in one dimension and in two dimensions are considered and are shown to dominate the maximum likelihood estimator (MLE) under normalized squared error loss

Chang, Yuan-Tsung

Shinozaki, Nobuo

English

Kyoto University Research Information Repository

RIGHT:URL:CITATION:AUTHOR(S):ISSUE DATE:TITLE:Simultaneous estimation of Poisson meansin two-way contingency tables undernormalized squared error loss inmultiplicative models (Bayesian approachesand statistical inference)Chang, Yuan-Tsung; Shinozaki, NobuoChang, Yuan-Tsung ...[et al]. Simultaneous estimation of Poisson means in two-way contingency tables undernormalized squared error loss in multiplicative models (Bayesian approaches and statistical inference). 数理解析研究所講究録 2022, 2221: 90-1002022-06http://hdl.handle.net/2433/27717190Simultaneous estimation of Poisson means in two-way contingency tables under normalized squared error loss in multiplicative models *Chang Yuan-Tsung, **Shinozaki Nobuo *Department of Social Information, Mejiro University **Faculty of Science and Technology, Keio University Abstract Shrinkage estimation of Poisson means is considered when observations are given in the form of a two-way contingency table. Assuming a multiplicative Poisson model, estimators which shrink to the specified values or an order statistic in one dimension and in two dimensions are considered and are shown to dominate the maximum likelihood estimator (MLE) under normalized squared error loss. 1 Introduction Consider a two-way multiplicative model where X;j, i = 1, ... , I, j = l, ... , J, are inde-pendent random Poisson random variables with means A;j = >.a;(Jj, i = 1, ... , I, j = l, ... , J, where a; 2: 0 and (Jj 2'. 0 satisfy I:{=1 a; = 1 and I:{=1 (Jj = 1, respectively. We denote the one-dimensional frequencies and the total frequency by J I I J Xi+= LXij, 'i = l, ... 'I, X+j = LXij, j = l, ... 'J, X++ = LLXij· j=l i=l i=l j=l As discussed in Hara and Takemura (2006) complete sufficient statistics are x 1 =(xi+, ... , xI+) and x 2 = (x+1, ... , X+J ). The MLE of A;j is if x++-=/ 0 if X++ = 0. They have given a class of improved estimators which shrink the MLE toward the origin under the normalized squared error loss. The simple one is 8:JT=Xi+X+j{l- d }, i=l, ... ,J, j=l, ... ,J,. X++ X++ + d 91where dis a positive constant, 0 < d::; 2(J + J - 2) and I, J ~ 2. The following lemma is a special case of Lemma 2.1 of Hara and Takemura (2006) and is useful to evaluate the risk of the shrinkage estimators when normalized squared error loss is concerned. Lemma 1.1. If g(x 1 , x 2 ) is a real-valued function satisfying Elg(x1 , x 2 ) I < oo and g(x1 , x 2 ) = 0 when xi+= 0 or X+j = 0, then { g(x1, x2)} { (x++ + 1) r J } E A;j = E (xi++ l)(x+i + 1)9(x1 + e;' X2 + ej) ' where e{ (ef) is Ix 1 (J x 1) unit vector with i-th (j-th) component 1. In Section 2 multiplicative Poisson model is assumed and shrinkage to a specified value or an order statistic is considered in one dimension and in two dimensions. In subsections 2.1 and 2.2 we give improved estimators which shrink xi+ in the MLE to an order statistic and a given non-negative constant, respectively. Further, in subsections 2.3 and 2.4 we are concerned with impoved estimators which shrink xi+ and X+j simultaneously. A discussion is given in subsection 2.5 and simulation results which show the performance of the proposed estimators are given in subsection 2.6. 2 Shrinkage estimators in the multiplicative Poisson model First, we consider one-dimensional shrinkage to an order statistic. 2.1. One-dimensional shrinkage to an order statistic Let X(f)+ be the £-th smallest observation among xl+, ... , XJ+- We assume that J ~ £ + 2 and consider the following estimator which shrinks xi+ toward X(t)+ when xi+ ~ X(t)+= ,(l) = x+i { . _ (W) (xi+ - X(f)+)+} u,J X++ x,+ cp W + d , i = 1, ... , I, j = 1, ... , J, where W = L{=1(x;+ - X(f)+)+, a+= max(0,a) and dis a positive constant. Then we have the following. Theorem 2.1. Suppose that cp(W) is a non-decreasing function satisfying 0 :::; cp(W) :::; 2(J - £ - 1) and that d ~ sup cp(W)/2. Then 5gl, i = 1, ... , I improves upon the MLE >.tf L, i = 1, ... ,I under the loss function L{=1 ( >.ii - A;j )2 / A;j for any j = 1, ... , J. Remark 2.1. Theorem 2.1 can be generalized directly to the case of Poisson multiplica-tive model for a multi-way contingency tables by using a lemma (Lemma 3.1 of Hara and Takemura (2006)) which is a generalization of Lemma 2.1. For example, consider the case of a 3-way contingency table Xijk, i = 1, ... , I, j = 1, ... , J, k = 1, ... , K where Xijk are independent Poisson random variables with means >.ijk· Let xi++, x+i+, x++k and X+++ 92denote the one-dimensional marginal frequencies and the total frequency. Let j and k be arbitrarily fixed and consider the simultaneous estimation of Aijk, ... , AJjk under the loss function I:{=1 (>.ijk - A;jk)2 / A;jk· Then, by adopting similar notations and conditions on cp(W) and d, we see that the estimator X+j+X++k { _ _ (W) (xi++ - X(f)++)+} 2 x,++ cp W + d , i = l, ... , I X+++ improves upon the MLE X;++X+i+X++k! x!++' i = 1, ... , I. Next we consider the estimators shrink >.ff L to a specified non-negative values, b;. 2.2. One-dimensional shrinkage to a specified point Let b; 2: 0, i = 1, ... , I be given numbers and we propose the following shrinkage estimator which shrinks xi+ to b; when xi+ 2: b;: t5(2) = x+i { _ _ (N W) (xi+ - b;)+} 'J x++ x,+ 1P ' W + d(N) ' i = 1, ... 'I, j = l, ... 'J, where W = I:{=1(xi+ - b;)+ and N = #{ilxi+ 2: b;}. Then we have the following. Theorem 2.2, Suppose that cp(N, W) is a non-decreasing function of W and satisfies 0 :S cp(N, W) :S 2(N - 1)+ for any O :SN :SI. Suppose that d(N) 2: supw cp(N, W)/2. Then 8&), i = 1, ... , I improves upon the MLE >.ff L, i = 1, ... , I under the loss function I:{=1(>.;j - A;j)2/A;j for any j = l, ... , J. It may be noticed that the shrinkage is made only when N 2: 2. Next we consider two-dimensional shrinkage to order statistics or to a specified point. 2.3. Two-dimensional shrinkage to order statistics Let X(f)+ and X+(m) be the £-th and m-th smallest observation among Xi+, ... , Xf+ and x+l, ... , X+J, respectively. We assume that I 2: £ + 2 and J 2: m + 2 and consider the estimator which shrinks xi+ toward X(t)+ when xi+ 2: X(f)+ in the first dimension and shrinks X+j toward X+(m) when X+j 2: X+(m) in the second dimension simultaneously. To improve upon the MLE >.ff L, we propose the following estimator : ,(3) = _1_{ . _ (W) (xi+ - X(f)+)+ }{ . _ (W) (x+j - X+(m))+} u,1 x,+ ip1 1 W + d X+1 1P2 2 W. + d , X++ 1 1 2 2 i = 1, ... , I, j = l, ... , J, where W1 = I:{=1(x;+-X(f)+)+ and W2 = I:f=1(x+i-X+(m))+ and d1 and d2 are positive constants. Then we have the following. Theorem 2.3. Suppose that cp1(W1) and cp2(W2) are non-decreasing functions sat-isfying O :S cp1 (W1) :S I - £ - l and O :S cp2(W2) :S J - m - 1, respectively. If d1 2: (J-£-1)/(J-£)supcp1(W1) andd2 2: (J-m-1)/(J-m)supcp2(W2). Then 93Jt), i = 1, ... , I, j = 1, .. . , J improves upon the MLE 5.ff L under the loss function I:{=1 I:f=l().ij - A;j)2/>...ij· Remark 2.3. Theorem 2.3 is directly generalized to the case of multi-way contingency tables. Since the notations and conditions are essentially the same, we only give a sketch of the result for the case of 3-way contingency table. We shrink xi++ toward X(e)++ when xi++ 2'.'. X(e)++ in the first dimension and shrink x+i+ toward X+(m)+ when x+i+ 2'.'. X+(m)+ in the second dimension. Under the loss function I J Ll)5.ijk -Aijk)2/>...ijk, i=l j=l where k = l, · · ·, K is arbitrarily fixed, the improved estimator is given by 2.4. Two-dimensional shrinkage to a specified point Let b; 2'.'. 0, i = 1, ... , I and Cj 2'.'. 0,j = 1, ... , J be given numbers. Assuming that I, J 2'.'. 2, we shrink xi+ to b; when X;+ 2'.'. b; and X+j to Cj when X+j 2'.'. Cj. To improve upon the MLE 5.f!L, we propose the following estimator where W1 = I:{=1(x;+ - b;)+, W2 = I:f=1(x+i - ci)+, N1 = #{ilxi+ 2:'. b;,i = 1, ... ,I} and N2 = #{Jlx+i 2:'. c;,j = 1, ... , J}. Although it may be natural to put the condition I:{=1 b; = I:f =l Cj, we do not need it in the following. Theorem 2.4. Suppose that cpe(Ne, We) is a non-decreasing function of Wt and satisfies 0 :S: <f!t(Nt, Wt) :S: (Nt-1)+ for any Nt 2'.'. 0, and that dt(Nt) 2'.'. (Ne-l)+ / Nt SUPwe <f!t(Nt, Wt), for any Ne 2'.'. 0, £ = 1, 2. Then J&) improves upon the MLE 5.ff L under the loss function I J ' 2 Li=l Lj=l(>...ij - A;j) />...ij· It may be noticed that the shrinkage in the £-th dimension is made only when Nt 2'.'. 2. 2.5. A discussion Here we mention the possibility of the two-dimensional shrinkage estimators other than Jt) and J&) given in subsections 2.3 and 2.4, respectively. We only give two alternative estimators for J&). The following estimator is the simple average of the one-dimensional 94shrinkage estimator 8g) and its counterpart which makes shrinkage in the second dimen-sion: xi+X+j 'P1(N1, W1) X+j (xi+ - bi)+ 'P2(N2, W2) xi+ (x+j - cj)+ ----x++ 2 x++ W1 +d1(N1) 2 x++ W2 + d2(N2)' where Wiand Ni, i = 1, 2, are defined in 2.4. It is easily shown that this estimator improves upon the MLE when <p(Ni, Wi) and di(Ni), i = 1, 2, satisfy the similar conditions as given in Theorem 2.2. We may pool W1 and W2 and consider the following estimator xi+x+j _ <p(N, W) X+j(xi+ - bi)++ xi+(x+j - cj)+ x++ 2 x++{W + d(N)} where W = (W1 + W2)/2 and N = N1 + N2 . Although this estimator will dominate the MLE under suitable conditions on <p(N, W) and d(N), we do not pursue it here further. Unfortunately, these two estimators do not give the estimates which belong to the parameter space of the multiplicative Poisson models, whereas the estimators 8i) and 8i~) do. 2.6. Simulation study Through simulation study we evaluate the risk performance of the following estimators: 8?j=Xi+X+j{1- d }, i=l, ... ,1,j=l, ... ,.l, X++ X++ + d s:min - X+j { - - d (xi+ - X(I)+)+ } uij - x,+ 1 , i = 1, ... , I, j = 1, ... , J, x++ ~k=l (xk+ - X(l)+) + d s:,bJ_ = X+j {x,·+ - d(N) I (xi+ - bi)+ } ' u i = 1, ... , I, j = 1, ... , J. x++ ~k=l (xk+ - bk)++ d(N) For an estimator 8ij of Aij, i = 1, ... , I, j = 1, ... J, let 8 j = ( 81j, ... , 8 I j), j = 1, ... , J and 8 = ( 81 , ... , 8 J). When we estimate Aij by 8ij, the expected loss is given by E{~(8 ->..--) 2 } =R(.X 8-) \ . . '] '] - ' '] ' /\iJ where A= {>..ij, i = 1, ... , I, j = 1, ... , .!}. The total risk is I J R( .X, 8) = L L R( .X, 8ij). i=l j=l We are also concerned with the columnwise risks I R(.X, 8j) = L R(.X, 8iJ, j = 1, ... , J. i=l 95Before describing the details of our simulation study, we state some fundamental properties of the total risk function. We consider a class of estimators 6ij of the form 7Pi(x 1)x+j/x++, i = 1, ... , I,j = 1, ... , J, where x 1 = (x1+, ... , XJ+)- We notice that all the estimators considered here are of this form. Applying Lemma 1.1, we have where e{ is the I x 1 unit vector with i-th component 1. Thus we have which is a function of (>.1+, ... , AJ+) and is independent of (/31 , ... , f3J). Therefore we see that the total risk R(>.., c5) is independent of (/31, ... , f3J ). Further, if 7Pi(x1) is of the form ¢(x++)xi+, we have Therefore the total risk depends only on>.=>-++· We may notice here that the total risk of 5,.ML (¢(x++ = 1) is expressed as which is a decreasing function of >.. Now we describe the details of the estimators used in our simulation study. We only consider the case I = 5 and J = 4. The constant din c5min is chosen as 2(5- 2) = 6, which is an upper bound of cp(W) given in Theorem 2.1. Although simulation experiments are performed also for the cased= 3, the performance for the cased= 6 is generally better. Similarly, d(N) in c5b is chosen as 2(N-1)+. The constant din c5° is set d = 2(5+4-2) = 14, since it is of interest and does not show largely different performance compared with the choice of d = 7. We have tried two choices of bin c5b: Although these choices use the true values of >.i+'s, it will make clear how much improve-ment is possible by choosing t5b pertinently. To evaluate the risk improvement we use the (total) relative improvement rate and also the columnwise relative improvement rate In Tables 2.1 and 2.2, based on 100,000 replications, the total relative improvement rates and the columnwise relative improvement rates are given respectively for the four estima-tors. 96a(l) = (0.2, 0.2, 0.2, 0.2, 0.2), /3(l) = (0.25, 0.25, 0.25, 0.25), a(2) = (0.3, 0.2, 0.2, 0.15, 0.15), 13(2) = (0.4, 0.3, 0.15, 0.15), a(3l = (0.4, 0.3, 0.1, 0.1, 0.1), 13(3) = (0.6, 0.2, 0.1, 0.1), 9 multiplicative models are considered in Table 2.1, where >. = 10, 50, 100 and a= a(1), a(2), a(3) (/3 = 13(1) is fixed for convenience). In Table 2.2, 9 multiplicative mod-els are considered, where >. = 10, 50, 100 and We summarize the simulation results as follows. 1. When we are concerned with the total risk, the relative improvement rate of Jmin is larger than that of JD for the case where a; = 0.2, i = 1, 2, ... , 5 which seems to be the most favorable case for Jmin. Although the relative improvement rate of Jmin gets smaller when the values of a;'s are more imbalanced, especially for larger values of >., it is still nearly identical to that of JD. 2. The relative improvement rate of JbC2J is the largest among those of all the estimators for all cases. This implies that if we choose b;, i = 1, ... , I pertinently based on the prior information on >.i+, i = 1, ... , I, then we are able to get large improvement. However, the choice b?) = >.i+, i = 1, ... , I is not pertinent and choosing smaller values such as max(0, >.i+ - A+), i = 1, ... , I is preferable. 3. The value of >. largely affects the performance of not only JD but also Jmin and Jb_ When the value of>. gets larger, the improvement rate of the estimators gets smaller. 4. When f31 = · · · = fJJ, the columnwise risks are all the same and the columnwise relative improvement rates are equal to the total relative improvement rate for any estimator. The imbalance among the values of fJj 's affects the columnwise risk performance of JD the most as we see from Tables 2.2. For the columns with relatively larger values of fJj, the columnwise relative improvement rates of JD get smaller and are sometimes negative. This implies that JD will not dominate the MLE in terms of the columnwise risks for large values of d. 97Table 2.1 Total risks and relative improvement rates Total risk 1 - R(>.., 6)/ R(>.., >,.ML) a Estimators d/d(N) >. = 10 >. = 50 >. = 100 >. = 10 >. = 50 >. = 100 >,.ML 9.181 8.237 8.123 60 14 5.219 7.721 7.966 0.432 0.063 0.019 a(ll 6min 6 5.055 6.370 6.752 0.449 0.227 0.169 6b(l) 2(N - 1)+ 5.816 6.630 6.934 0.367 0.195 0.146 6b(2) 2(N - 1)+ 4.898 6.168 6.620 0.466 0.251 0.185 >,.ML 9.174 8.238 8.121 60 14 5.217 7.720 7.963 0.431 0.063 0.019 a(2) 6min 6 5.161 6.747 7.223 0.437 0.181 0.111 6b(l) 2(N - 1)+ 5.988 6.688 6.922 0.347 0.188 0.148 6b(2) 2(N - 1)+ 4.970 6.180 6.626 0.458 0.250 0.184 >,.ML 9.180 8.258 8.107 60 14 5.222 7.735 7.951 0.431 0.063 0.019 a(3) 6min 6 5.429 7.249 7.577 0.409 0.122 0.065 6b(l) 2(N - 1)+ 5.675 6.571 6.853 0.382 0.204 0.155 6b(2) 2(N - 1)+ 4.598 6.165 6.540 0.499 0.253 0.193 98Table 2.2 Columnwise relative improvements rates a,(3 ,x_ Estimator 1st column 2nd column 3rd column 4th column total risk 10 50 0.316 0.388 0.541 0.544 0.431 t5min 0.415 0.436 0.482 0.483 0.449 Jb(l) 0.387 0.375 0.347 0.348 0.367 t:5b(2) 0.431 0.452 0.500 0.502 0.466 50 JU -0.020 0.031 0.147 0.150 0.063 a=a(ll t5min 0.237 0.230 0.216 0.216 0.226 (3=(3(2) Jb(l) 0.220 0.205 0.168 0.169 0.195 t:5b(2) 0.264 0.257 0.237 0.237 0.251 100 5u -0.033 -0.001 0.073 0.075 0.019 t5min 0.185 0.175 0.152 0.152 0.169 Ob(l) 0.169 0.155 0.122 0.122 0.146 t:5b(2) 0.204 0.193 0.165 0.165 0.185 10 50 0.218 0.480 0.615 0.616 0.431 t5min 0.386 0.463 0.503 0.504 0.449 t5b;1) 0.403 0.358 0.336 0.336 0.367 t:5b(2) 0.400 0.480 0.523 0.524 0.466 50 5u -0.083 0.100 0.205 0.208 0.062 a=a(ll t5min 0.246 0.221 0.207 0.208 0.226 (3=(3(3) Ob(l) 0.241 0.183 0.149 0.150 0.195 t:5b(2) 0.276 0.244 0.226 0.227 0.251 100 JU -0.069 0.044 0.111 0.111 0.020 t5min 0.199 0.161 0.139 0.140 0.169 Ob(l) 0.187 0.135 0.104 0.105 0.146 t:5b(2) 0.220 0.176 0.150 0.151 0.185 10 50 0.316 0.386 0.543 0.545 0.431 t5min 0.359 0.389 0.457 0.459 0.409 Ob(l) 0.404 0.390 0.360 0.360 0.382 t:5b(2) 0.466 0.486 0.531 0.533 0.499 50 JU -0.020 0.030 0.149 0.150 0.062 a=a.(3) t5min 0.103 0.114 0.142 0.142 0.122 (3=(3(2) Ob(l) 0.231 0.216 0.175 0.176 0.204 t:5b(2) 0.267 0.259 0.238 0.238 0.253 100 5u -0.030 -0.001 0.072 0.074 0.020 t5min 0.056 0.061 0.075 0.076 0.065 Ob(l) 0.179 0.164 0.128 0.128 0.154 t:5b(2) 0.214 0.200 0.170 0.171 0.193 99References [1] Chang, Y.-T. and Shinozaki, N. (2006). Estimation of ordered means of two Poisson distributions. Comm. Statist.-Theory and Methods-, 35, 1993-2003. [2] Chang, Y.-T. and Shinozaki, N. (2019). New types of shrinkage estimators of Poisson means under the normalized squared error loss. Comm. Statist.-Theory and Methods-, 48, 1108-1122. [3] Chou, J.P. (1991). Simultaneous estimation in discrete multivariate exponential fam-ilies. Ann. Statist., 19, 314-328. [4] Clevenson, M. and Zidek, J. (1975). Simultaneous estimation of the means of inde-pendent Poisson Laws. J. Amer. Statist. Assoc., 70, 698-705. [5] Ghosh, M. and Yang, M. C. (1988). Simultaneous estimation of Poisson means under entropy loss. Ann. Statist., 16, 278-291. [6] Ghosh, M., Hwang, J. T. and Tsui, K. W. (1983). Construction of improved estimators in multiparameter estimation for discrete exponential families. Ann. Statist., 11, 351-367. [7] Hamura, Y. and Kubokawa, T. (2019). Simultaneous estimation of parameters of Pois-son distributions with unbalanced sample sizes. Japanese J. Statist. and Data Analysis, https://doi.org/10.1007 /s42081-019-00039-x. [8] Hara, H. and Takemura, A. (2006). Simultaneous estimation of the means in some Poisson log linear models. J. Japan Statist. Soc., 36, 17-36. [9] Hwang, J. T. (1982). Improving upon standard estimators in discrete exponential families with applications to Poisson and negative binomial cases. Ann. Statist., 10, 857-867. [10] Komaki, F. (2004). Simultaneous prediction of independent Poisson observables. Ann. Statist., 32, 1744-1769. [11] Komaki, F. (2015). Simultaneous prediction for independent Poisson process with different durations. J. Multivariate Anal., 141, 35-48. [12] Peng, J. C. M. (1975). Simultaneous estimation of the parameters of independent Poisson distributions. Technical Reports 78. Department of Statistics, Stanford Uni-versity. [13] Tsui, K. W. (1984). Robustness of Clevenson-Zidek-type estimators. J. Amer. Statist. Assoc., 79, 152-157. [14] Tsui, K. W. (1986). Further developments on the robustness of Clevenson-Zidek-type means estimators. J. Amer. Statist. Assoc., 81, 176-180. [15] Tsui, K. W. and Press, S. J. (1982). Simultaneous estimation of several Poisson parameters under k-normalized squared error loss. Ann. Statist., 10, 93-100. 100Acknowledgement: This work was supported by the Research Institute for Mathemat-ical Sciences, an International Joint Usage/Research Center located in Kyoto University. This work is salso upported by Grant-in-Aid for Scientific Research (C) No. 18K11196. Two authors also thank editor in chief Professor Koike Ken-ichi, Nihon University for his many supports. *Chang Yuan-Tsung Mejiro University Tokyo 161-8539 JAPAN E-mail address: chogenso@gmail.com **Shinozaki Nobuo Keio University Kanagawa 223-8522 JAPAN 

Simultaneous estimation of Poisson means in two-way contingency tables under normalized squared error loss in multiplicative models (Bayesian approaches and statistical inference)

https://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/277171/1/2221-09.pdf

Simultaneous estimation of Poisson means in two-way contingency tables under normalized squared error loss in multiplicative models (Bayesian approaches and statistical inference)

Abstract

Similar works

Full text

Available Versions

Kyoto University Research Information Repository