A Posterior p-Value for Homogeneity Testing of the Three Sample Problem

1. Introduction

In this paper we focus on the model proposed by Hosmer [1], which is used to study the halibut data. There are two different sources of halibut data. One is from the research cruises, where the sex, age and length of the halibut are available while another comes from the commercial catch where only age and length can be obtained since the fish have been cleaned before the boats returned to the port. The length distribution of an age class of halibut is closely approximated by a mixture of two normal distributions, which is

$$\begin{array}{cc}\hfill & X_i\stackrel{iid}{\sim }{f}_1\left(y\right),i=1,\dots ,n_1\hfill \\ \hfill & Y_j\stackrel{iid}{\sim }{f}_2\left(y\right),j=1,\dots ,n_2\hfill \\ \hfill & Z_k\stackrel{iid}{\sim }\lambda f_1\left(y\right)+(1-\lambda )f_2\left(y\right),k=1,\dots ,n_3,\hfill \end{array}$$

(1)

where $f_1$ and $f_2$ are the probability density functions of the normal distributions and $\lambda$ is the proportion of the male halibut in the commercial catches. Hosmer [1] estimated the parameters of the two distributions using the iterative maximum likelihood estimate method. Murray and Titterington [2] summarized the techniques for maximum likelihood estimation and Bayesian analysis. Anderson [3] proposed a semiparametric modeling assumption known as the exponential tilt mixture model. The estimating of the proportion is done by a general method based on direct estimation of the likelihood ratio. The estimation of the model is further studied by Qin [4], who extended Owen ’s [5] empirical likelihood to the semiparametric model. Qin gave the asymptotic variance formula for the maximum semiparametric likelihood estimation. Later, Zou et al. [6] showed that the partial profile empirical likelihood also worked well with realistic sample size. An EM algorithm for this semiparametric model is further given by Zhang [7]. Furthermore, Inagaki and Komaki [8] and Tan [9] respectively modified the profile likelihood function and provided better estimators for the parameters.

Except for the estimation of parameters, another important issue is to test the homogeneity of the model. Thus the null hypothesis is

$$H_0:f_1=f_2.$$

To test the null hypothesis, Liang and Rathouz [10] proposed a score test and applied it to genetic linkage. Lemdani and Pons [11] suggested using likelihood ratio tests(LRT) under modified likelihood function since classical results on LRT may be invalid. Chen et al. [12,13] proposed a modified LRT with a general parametric kernel distribution family and they proved that the modified LRT has a $\chi$-type of null limiting distribution and is asymptotically most powerful under local alternatives. Chen and Li [14] proposed an EM approach for the normal mixture models. Li et al. [15] used a high-order expansion to establish a nonstandard convergence rate $N^{-1/4}$ for the odds ratio parameter estimator and solved the problem of degeneration of the Fisher information. These works are applied to a variety of genetic mixture model in real applications. For example, genetic imprinting and quantitative trait locus mapping, see Li et al. [16] and Liu et al. [17].

Most mixture models discribed above mainly consider the case when $f_1$ and $f_2$ are normal distributions. In this paper, we want to extend the conclusion to more general cases. Similar question has been researched by Ren et al. [18]. In their paper a two-block Gibbs sampling method are proposed to obtain the samples of the generalized pivot quantities of the parameters. They studied both cases when $f_1$ and $f_2$ are normal and logistic distributions. In our paper, we assume that $f_1$ and $f_2$ are in a specified location-scale family with location parameter $\mu$ and scale parameter $\sigma$. We propose a posterior p-value based on the posterior distribution to test the homogeneity. To sample from the posterior distirbutions, we propose to use approximate Bayesian computation (ABC) method for the case when $f_1$ and $f_2$ are normal distributions, which is different from the cases when $f_1$ and $f_2$ are general distributions. This is because the posterior distirbution of the normal case can be regard as using the information contained in the first two samples as prior distribution and updating it via the third one without loss of information. We find in our simulation that this method is promising and efficient even we use the simplest reject-sampling. For the general case, since ABC method is no longer available, we use MCMC method such as Metropolis-Hastings sampling method proposed by Hannig et al. [19] and the two-block Gibbs sampling proposed by Ren et al. [18] to sample from the posterior distribution.

The paper is organized as follow. In Section 2 we first define the regular location-scale family and give some properties of the family. Then we propose our posterior p-value for testing the homogeneity. We further introduce the sampling method for different cases. A real data of the halibut is studied in Section 3 to illustrate the validity of our method. The simulation study is given in Section 4 while the conclusion is given in Section 5.

2. Test Procedure

In this section we consider model (1) where the distributions are in a certain regular location-scale family. Thus, we first give the definition in the following subsection.

2.1. Regular Location-Scale Family

In this section we first give the definition of the regular location-scale family.

Definition 1 (regular location-scale family).

Let $f\left(x\right)$ be a probability density function. If $f\left(x\right)$ satisfies

(1) 

$f\left(x\right)>0$, $-\infty <x<\infty$;

(2) 

$f^{″}\left(x\right)$ is continous;

(3) 

$\underset{x\to -\infty }{lim}{x}^2{f}^{\prime }\left(x\right)=\underset{x\to +\infty }{lim}{x}^2{f}^{\prime }\left(x\right)=0$;

(4) 

${\int }_{-\infty }^{+\infty }{x}^2\frac{{\left[f^{\prime }\left(x\right)\right]}^2}{f\left(x\right)}dx<\infty$.

Then $f\left(x\right)$ is defined as a regular density function, and

$${\mathcal{R}}_f=\left\{\frac{1}{\sigma }f\left(\frac{x-\mu }{\sigma }\right);\mu \in (-\infty ,+\infty ),\sigma \in (0,+\infty )\right\}$$

is defined as the regular location-scale family.

It is easy to verify that many families of distributions are regular location-scale families. For example, let

$$f_1\left(x\right)=\frac{1}{\sqrt{2\pi }}{e}^{-\frac{x^2}{2}},f_2\left(x\right)=\frac{e^{-x}}{{(1+e^{-x})}^2}.$$

Then $f_1\left(x\right)$ and $f_2\left(x\right)$ are regular density functions. The families of distributions that constructed by $f_1\left(x\right)$ and $f_2\left(x\right)$ are regular, which are the families of normal distributions and logistic distibutions, respectively. The two families of distributions are concerned in the paper later.

The following lemma highlights some properties of this family.

Lemma 1.

If $f\left(x\right)$ is a regular density function, then we have

(1) 

$\underset{x\to -\infty }{lim}xf\left(x\right)=\underset{x\to +\infty }{lim}xf\left(x\right)=0$;

(2) 

${\int }_{-\infty }^{+\infty }{f}^{\prime }\left(x\right)dx=0$;

(3) 

${\int }_{-\infty }^{+\infty }x{f}^{\prime }\left(x\right)dx=-1$;

(4) 

${\int }_{-\infty }^{+\infty }{f}^{″}\left(x\right)dx=0$;

(5) 

${\int }_{-\infty }^{+\infty }x{f}^{″}\left(x\right)dx=0$;

(6) 

${\int }_{-\infty }^{+\infty }{x}^2{f}^{″}\left(x\right)dx=0$.

The proof of this lemma is given in Appendix.

We further calculate the Fisher information matrix of the regular location-scale family with the following proposition.

Proposition 1.

Assume that $f(x;\xi )=\frac{1}{\sigma }f\left(\frac{x-\mu }{\sigma }\right)$ is in the regular location-scale family, where $\xi ={(\mu ,\sigma )}^{\top }$. The parameter space $\Omega =\left\{\right(\mu ,\sigma ):-\infty <\mu <\infty ,\sigma >0\}$. Let $l(\xi ,X)=logf(X;\xi )$. Then

(1) 

The score function satisfies

$$E_{\xi }\left[\frac{\partial l(X;\xi )}{\xi }\right]=0;$$

(2) 

The Fisher information matrix satisfies

$$0<I_f\left(\xi \right)=E_{\xi }{\left[\left(\frac{\partial l(X;\xi )}{\partial \xi }\right)\left(\frac{\partial l(X;\xi )}{\partial \xi }\right)\right]}^{\top }=\frac{1}{{\sigma }^2}\mathbf{C}\left(f\right)=\frac{1}{{\sigma }^2}\left(\begin{array}{cc}{C}_{11}\left(f\right)& C_{12}\left(f\right)\\ C_{21}\left(f\right)& C_{22}\left(f\right)\end{array}\right)<\infty ,$$

where

$$\begin{array}{cc}\hfill C_{11}\left(f\right)& ={\int }_{-\infty }^{\infty }\frac{{\left[f^{\prime }\left(y\right)\right]}^2}{f\left(y\right)}dy\hfill \\ \hfill C_{12}\left(f\right)& =C_{21}\left(f\right)={\int }_{-\infty }^{\infty }y\frac{{\left[f^{\prime }\left(y\right)\right]}^2}{f\left(y\right)}dy\hfill \\ \hfill C_{22}\left(f\right)& ={\int }_{-\infty }^{\infty }{\left(1+y\frac{f^{\prime }\left(y\right)}{f\left(y\right)}\right)}^{2}f\left(y\right)dy={\int }_{-\infty }^{\infty }{y}^2\frac{{\left[f^{\prime }\left(y\right)\right]}^2}{f\left(y\right)}dy-1\hfill \end{array}$$

(3) 

The Fisher information matrix is given by

$$I_f\left(\xi \right)=-E_{\xi }\left[\frac{{\partial }^{2}f(X,\xi )}{\partial \xi \partial {\xi }^{\top }}\right].$$

The proof is given in Appendix.

Proposition 2.

Assume that $0<{\lambda }_0<1$ and $f(·)$ is regular. Then $\left\{g\right(x;\theta ):\theta \in \Omega \}$ given by

$$g(x,\theta )=\frac{{\lambda }_0}{{\sigma }_1}f\left(\frac{x-{\mu }_1}{{\sigma }_1}\right)+\frac{1-{\lambda }_0}{{\sigma }_2}f\left(\frac{x-{\mu }_2}{{\sigma }_2}\right),$$

(2)

where $\theta ={({\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2)}^{\top }$ and $\Omega =R^2\times {R^{+}}^2$, has following properties.

(1) 

$$E_{\theta }\left[\frac{\partial logg(x,\theta )}{\partial \theta }\right]=0;$$

(2) 

$$I\left(\theta \right)=E{\left[\left(\frac{\partial logg(x,\theta )}{\partial \theta }\right)\left(\frac{\partial logg(x,\theta )}{\partial \theta }\right)\right]}^{\top }<\infty ;$$

(3) 

$$I\left(\theta \right)=-E_{\theta }\left[\frac{{\partial }^{2}logg(x;\theta )}{\partial \theta \partial {\theta }^{\top }}\right].$$

The proof is given in Appendix.

We then give the Fisher information matrix of the normal and logistic distribution. For the normal distribution, we have

$$\begin{array}{cc}\hfill C_{11}\left(f\right)& ={\int }_{-\infty }^{\infty }\frac{1}{\sqrt{2\pi }}{e}^{-\frac{y^2}{2}}{y}^{2}dy=1\hfill \\ \hfill C_{12}\left(f\right)& =C_{21}\left(f\right)={\int }_{-\infty }^{\infty }\frac{1}{\sqrt{2\pi }}{e}^{-\frac{y^2}{2}}{y}^{3}dy=0\hfill \\ \hfill C_{22}\left(f\right)& ={\int }_{-\infty }^{\infty }\frac{1}{\sqrt{2\pi }}{e}^{-\frac{y^2}{2}}{y}^{4}dy-1=2\hfill \end{array}$$

Thus the Fisher information matrix of normal distribution is

$$I_f\left(\xi \right)=\frac{1}{{\sigma }^2}\left(\begin{array}{cc}1& 0\\ 0& 2\end{array}\right)$$

Similarly, for the logistic distribution,

$$\begin{array}{cc}\hfill C_{11}\left(f\right)& ={\int }_{-\infty }^{\infty }\frac{e^{-y}{(e^{-y}-1)}^2}{{(1+e^{-y})}^4}dy=1/3\hfill \\ \hfill C_{12}\left(f\right)& =C_{21}\left(f\right)={\int }_{-\infty }^{\infty }\frac{e^{-y}{(e^{-y}-1)}^2}{{(1+e^{-y})}^4}dy=0\hfill \\ \hfill C_{22}\left(f\right)& ={\int }_{-\infty }^{\infty }\frac{e^{-y}{(e^{-y}-1)}^2}{{(1+e^{-y})}^4}dy=\frac{1}{3}+\frac{{\pi }^2}{9}\hfill \end{array}$$

Thus the Fisher information matrix of logistic distribution is

$$I_f\left(\xi \right)=\frac{1}{{\sigma }^2}\left(\begin{array}{cc}\frac{1}{3}& 0\\ 0& \frac{1}{3}+\frac{{\pi }^2}{9}\end{array}\right)$$

2.2. A Posterior p-Value

Now we consider testing the homogeneity of model (1) where $f_1$ and $f_2$ are in ${\mathcal{R}}_f$,

$$f_1=\frac{1}{{\sigma }_1}f\left(\frac{x-{\mu }_1}{{\sigma }_1}\right),f_2=\frac{1}{{\sigma }_2}f\left(\frac{x-{\mu }_2}{{\sigma }_2}\right).$$

This is equivalent to testing equality of the parameters of the two density functions, that is,

$$H_0:{\mu }_1={\mu }_2,{\sigma }_1={\sigma }_{2}v.s.H_1:{\mu }_1\ne {\mu }_2\mathrm{or}{\sigma }_1\ne {\sigma }_2.$$

(3)

Consider the density function

$$g(x;\theta )=\frac{\lambda }{{\sigma }_1}f\left(\frac{x-{\mu }_1}{{\sigma }_1}\right)+\frac{\lambda }{{\sigma }_2}f\left(\frac{x-{\mu }_2}{{\sigma }_2}\right),$$

where $\theta ={({\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2,\lambda )}^{\top }$ is the unknown parameter. $f\left(x\right)$ is the regular density function, then the Fisher information matrix is

$$I_g\left(\theta \right)=E_{\theta }\left[\frac{\partial logg(x;\theta )}{\partial \theta }\right]{\left[\frac{\partial logg(x;\theta )}{\partial \theta }\right]}^{\top },$$

where

$$\frac{\partial logg(x;\theta )}{\partial \theta }=\left(\begin{array}{cc}& \frac{1}{g(x;\theta )}\left[-\frac{\lambda }{{\sigma }_1^2}{f}^{\prime }\left(\frac{x-{\mu }_1}{{\sigma }_1}\right)\right]\\ & \frac{1}{g(x;\theta )}\left[-\frac{1-\lambda }{{\sigma }_2^2}{f}^{\prime }\left(\frac{x-{\mu }_2}{{\sigma }_2}\right)\right]\\ & \frac{1}{g(x;\theta )}\left[-\frac{\lambda }{{\sigma }_1^2}{f}^{\prime }\left(\frac{x-{\mu }_1}{{\sigma }_1}\right)-\frac{\lambda }{{\sigma }_1^3}{f}^{\prime }\left(\frac{x-{\mu }_1}{{\sigma }_1}\right)\right]\\ & \frac{1}{g(x;\theta )}\left[-\frac{1-\lambda }{{\sigma }_2^2}{f}^{\prime }\left(\frac{x-{\mu }_2}{{\sigma }_2}\right)-\frac{1-\lambda }{{\sigma }_2^3}{f}^{\prime }\left(\frac{x-{\mu }_2}{{\sigma }_2}\right)\right]\\ & \frac{1}{g(x;\theta )}\left[\frac{1}{{\sigma }_1}f\left(\frac{x-{\mu }_1}{{\sigma }_1}\right)-\frac{1}{{\sigma }_2}f\left(\frac{x-{\mu }_2}{{\sigma }_2}\right)\right]\end{array}\right)$$

When ${\mu }_1={\mu }_2,{\sigma }_1={\sigma }_2$,

$$\frac{\partial logg(x;\theta )}{\partial \lambda }=0,$$

the last row and column of $I_g\left(\theta \right)$ is zero, which means that $|I_g\left(\theta \right)|=0$ and is non-definite. Thus, we may encounter some difficulties when using some traditional test methods, such as the likelihood ratio test.

We suggest a solution here. First we assume that $\lambda ={\lambda }_0$ is known. Then there are 4 parameters and we still denote them by $\theta ={({\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2)}^{\top }$. We use the estimate of $\lambda$ instead since $\lambda$ is unknown. This is because that when the homogeneity hypothesis holds, the distribution of the population is irrelative to $\lambda$. So the level of the test is irrelative to the estimate of $\lambda$. We then give the inference on $\theta$ below. For the first two samples, the fiducial density of $({\mu }_1,{\sigma }_1)$ and $({\mu }_2,{\sigma }_2)$ are

$$({\mu }_1,{\sigma }_1)\propto \left[\prod _{i=1}^{n_1}\frac{1}{{\sigma }_1}f\left(\frac{x_{1i}-{\mu }_1}{{\sigma }_1}\right)\right]\frac{1}{{\sigma }_1},({\mu }_2,{\sigma }_2)\propto \left[\prod _{j=1}^{n_2}\frac{1}{{\sigma }_2}f\left(\frac{x_{2j}-{\mu }_2}{{\sigma }_2}\right)\right]\frac{1}{{\sigma }_2},$$

(4)

where “∝” denote “proportion to”, see the example 3 of Hannig et al. [19]. Then to combine (4) with the third sample, we regard (4) as the prior distribution. By the Bayes’ theorem

$$\begin{array}{cc}\hfill \theta \propto & \prod _{k=1}^{n_3}\left[\frac{{\lambda }_0}{{\sigma }_1}f\left(\frac{x_{3k}-{\mu }_1}{{\sigma }_1}\right)+\frac{1-{\lambda }_0}{{\sigma }_2}f\left(\frac{x_{3k}-{\mu }_2}{{\sigma }_2}\right)\right]\hfill \\ \hfill ·& \left[\prod _{i=1}^{n_1}\frac{1}{{\sigma }_1}f\left(\frac{x_{1i}-{\mu }_1}{{\sigma }_1}\right)\right]\frac{1}{{\sigma }_1}·\left[\prod _{j=1}^{n_2}\frac{1}{{\sigma }_2}f\left(\frac{x_{2j}-{\mu }_2}{{\sigma }_2}\right)\right]\frac{1}{{\sigma }_2}\hfill \end{array}$$

(5)

Denote the probability measure on the parameter space determined by (5) by $P^{\Theta |x}$, where $x=(x_1^{\top },x_2^{\top },x_3^{\top })$, $x_1={(x_{11},x_{12},\cdots ,x_{1n_1})}^{\top }$, $x_2={(x_{21},x_{22},\cdots ,x_{2n_2})}^{\top }$, $x_3={(x_{31},x_{32},\cdots ,x_{3n_3})}^{\top }$. $\Theta$ denotes the random variable. We can see from the expression (5) that $P^{\Theta |x}$ is the posterior distribution under the prior distribution

$$d\theta =\frac{1}{{\sigma }_1{\sigma }_2}d{\mu }_{1}d{\mu }_{2}d{\sigma }_1{\sigma }_2.$$

Let

$$\mathbf{A}=\left[\begin{array}{cccc}1& -1& 0& 0\\ 0& 0& 1& -1\end{array}\right],\mathbf{b}={[0,0]}^{\top }.$$

Then the hypotheses (3) is equivalent to

$$H_0:\mathbf{A}\theta =\mathbf{b}v.s.H_1:\mathbf{A}\theta \ne \mathbf{b}.$$

(6)

where $\theta ={({\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2)}^{\top }$.

To establish Bernstein-von Mises theorem for multiple samples, we first introduce some necessary assumptions below. Let $l_i\left(\theta \right)$ be the log-likelihood function of the ith sample, where $i=1,2,3$.

Assumption 1.

Given any $\epsilon >0$, there exists $\delta >0$ such that in the expansion

$$l_i\left(\theta \right)=l_i\left({\theta }_0\right)+{(\theta -{\theta }_0)}^{\top }{l}_i^{\prime }\left({\theta }_0\right)-\frac{1}{2}{(\theta -{\theta }_0)}^{\top }[n{I}_i\left({\theta }_0\right)+R_{n_i}\left(\theta \right)](\theta -{\theta }_0),i=1,2,3,$$

where ${\theta }_0$ is the true value of the parameter.$I_i\left({\theta }_0\right)$ is the Fisher Information matrix. The probability of the following event

$$sup\left\{\frac{1}{n}{\lambda }_{max}\left[R_{n_i}\left(\theta \right)\right]:\parallel \theta -{\theta }_0\parallel \le \delta \right\}\ge \epsilon$$

tends to 0 as $n\to \infty$, where $\parallel ·\parallel$ is the Euclidean norm and ${\lambda }_{max}\left(\mathbf{A}\right)$ denotes the largest absolute eigenvalues of a square matrix A.

Assumption 2.

For any $\delta >0$, there exists $\epsilon >0$ such that the probability of the event

$$sup\left\{\frac{1}{n_i}\left[l_i\left(\theta \right)-l_i\left({\theta }_0\right)\right]:\parallel \theta -{\theta }_0\parallel \ge \delta \right\}\le -\epsilon$$

tends to 1 as $n\to \infty$.

Assumption 3.

Under the prior π, there exist $k_0$ such that the integral of $\parallel \theta \parallel$ below exists,

$${\int \parallel \theta \parallel }^2\prod _{i=1}^{k_0}\left[\frac{1}{{\sigma }_1}f\left(\frac{x_{1i}-{\mu }_1}{{\sigma }_1}\right)\right]\prod _{j=1}^{k_0}\left[\frac{1}{{\sigma }_2}f\left(\frac{x_{2i}-{\mu }_2}{{\sigma }_2}\right)\right]\pi \left(\theta \right)d\theta <\infty .$$

Assumption 4.

When $n=n_1+n_2+n_3\to \infty$,

$$\frac{n_i}{n}\to r_i\in (0,1),i=1,2,3.$$

Then we give the Berstein-von Mises theorem for multiple samples as follow.

Theorem 1.

Denote the posterior density of $t=\sqrt{n}(\theta -T_n)$ by ${\pi }^{*}\left(t\right|x)$, where

$$T_n={\theta }_0+\frac{1}{n}{I}^{-1}\left({\theta }_0\right)l^{\prime }\left({\theta }_0\right).$$

If the Assumption 1, 2 and 4 hold, then

$${\int }_{\Omega }\left|{\pi }^{*}\left(t\right|x)-{\left(2\pi \right)}^{-\frac{k}{2}}{\left|I\left({\theta }_0\right)\right|}^{\frac{1}{2}}{e}^{-\frac{1}{2}{t}^{\top }I\left({\theta }_0\right)t}\right|dt\stackrel{P}{\to }0.$$

Furthermore, if Assumption 3 holds, then

$${\int }_{\Omega }\left(1+{\parallel t\parallel }^2\right)\left|{\pi }^{*}\left(t\right|x)-{\left(2\pi \right)}^{-\frac{k}{2}}{\left|I\left({\theta }_0\right)\right|}^{\frac{1}{2}}{e}^{-\frac{1}{2}{t}^{\top }I\left({\theta }_0\right)t}\right|dt\stackrel{P}{\to }0.$$

Then we can define the posterior p-value as follow

Definition 2.

Let

$$\begin{array}{cc}\hfill p\left(x\right)=P^{\Theta |x}(& {(\Theta -{\theta }_B)}^{\top }{A}^{\top }{\left(\mathbf{A}{\Sigma }_B{A}^{\top }\right)}^{-1}A(\Theta -{\mathbf{\theta }}_B)\hfill \\ & ⩾{(b-A{\mathbf{\theta }}_B)}^{\top }{\left(A{\Sigma }_B{\mathbf{A}}^{\top }\right)}^{-1}(b-A{\mathbf{\theta }}_B)),\hfill \end{array}$$

(7)

where $P^{\Theta |x}(·)$ is the probability under the posterior distribution. ${\theta }_B$ is the posterior mean and ${\Sigma }_B$ is the posterior covariance matrix. We call $p\left(x\right)$ as a posterior p-value.

The theorem below gurantees the validity of the posterior p-value.

Theorem 2.

Under the assumption of Theorem 1, the p-value define by (7) satisfies

$$p\left(X\right)\stackrel{d}{⟶}U(0,1).$$

where $“{\stackrel{d}{{}^{\prime }}}^{\prime }$ is the convergence in distribution and $U(0,1)$ is the uniform distribution on the internal $(0,1)$.

The proof is given in Appendix. Then for a given significance level $\alpha$, we may reject the null hypothesis if the p-value is less than $\alpha$.

2.3. Sampling Method

The posterior mean ${\theta }_B$ and the posterior variance ${\Sigma }_B$ in equation (7) can be estimated by the sample mean and variance, respectively. Now the remain problem is how to sample from the posterior distribution. When $\lambda$ is unknown, we first propose an EM algorithm to estimate $\lambda$, then we sample from the posterior distribution where $\lambda$ is fixed to the estimate of $\lambda$. The Markov Chain Monte Carlo (MCMC) method are commonly used. However, as we have mentioned earlier, MCMC method needs to discard a large number of samples in the burn-in period to gurantee the samples accepted sufficiently close to the ones from the real distribution. Fortunately, when $f_1$ and $f_2$ are normal distributions, we find that the posterior distribution can be transformed and sampled by using the approximate Bayesian computation (ABC) method. However, when $f_1$ and $f_2$ are some more common distributions, such as the logistic distributions, the two-block Gibbs sampling proposed by Ren et al. [18] can be an appropriate substitution. We will discuss the details in the following subsection.

2.3.1. Em Algorithm for $\lambda$

In this subsection we propose the EM alogorithm for estimating $\lambda$.

The log-likelihood function of the model is

$$L(x;\theta ,\lambda )=\sum _{i=1}^{n_1}log\left[f_1\left(x_{1i};\theta \right)\right]+\sum _{j=1}^{n_2}log\left[f_2\left(x_{2j};\theta \right)\right]+\sum _{k=1}^{n_3}log\left[p\left(x_{3k};\theta ,\lambda \right)\right],$$

where $f_1$ and $f_2$ are in the same regular location-scale family ${\mathcal{R}}_f$, with parameters $({\mu }_1,{\sigma }_1)$ and $({\mu }_2,{\sigma }_2)$ respectively. $p\left(x_{3k};\theta ,\lambda \right)$ in the log-likelihood function of the third sample is

$$p\left(x_{3k};\theta ,\lambda \right)=\lambda f_1\left(x_{3k};\theta \right)+\left(1-\lambda \right)f_2\left(x_{3k};\theta \right).$$

The EM algorithmis is first proposed by Dempster et al. [20] and broadly applied to a wide variety of parametric models, see McLachlan and Krishnan [21] for a better review.

Assume that we have obtain the estimate of the parameters after m times of iterative, denote them by ${\theta }^{\left(m\right)}={({\mu }_1^{\left(m\right)},{\sigma }_1^{\left(m\right)},{\mu }_2^{\left(m\right)},{\sigma }_2^{\left(m\right)})}^{\top }$ and ${\lambda }^{\left(m\right)}$. We introduce the latent variable $\gamma ={({\gamma }_1,{\gamma }_2,\cdots ,{\gamma }_{n_3})}^{\top }$, the component ${\gamma }_k$ indicates which distribution the sample $x_{3k}$ is drawn from. ${\gamma }_k=1$ when it is drawn from the first distribution $f_1(x_{3k};\theta )$, otherwise, ${\gamma }_k=0$. Then we have

$$P({\gamma }_k=1)=\lambda ,P({\gamma }_k=0)=1-\lambda ,k=1,2,\cdots ,n_3$$

The density of the joint distribution of $X_1,X_2,X_3,\gamma$ is

$$\prod _{i=1}^{n_1}{f}_1(x_{1i};\theta )\prod _{j=1}^{n_2}{f}_2(x_{2j};\theta )\prod _{k=1}^{n_3}{\left[\lambda f_1(x_{3k};\theta )\right]}^{{\gamma }_k}{\left[(1-\lambda )f_2(x_{3k};\theta )\right]}^{1-{\gamma }_k}.$$

Given $X_1=x_1,X_2=x_2,X_3=x_3$, the conditional distribution of ${\gamma }_k$ is

$${\left[\frac{\lambda f_1(x_{3k};\theta )}{\lambda f_1(x_{3k};\theta )+(1-\lambda )f_2(x_{3k};\theta )}\right]}^{{\gamma }_k}{\left[\frac{(1-\lambda )f_2(x_{3k};\theta )}{\lambda f_1(x_{3k};\theta )+(1-\lambda )f_2(x_{3k};\theta )}\right]}^{1-{\gamma }_k},$$

where ${\gamma }_k=0,1$, $k=1,2,\cdots ,n_3$. Thus the conditional expectation of ${\gamma }_k$ is

$$E_{(\theta ,\lambda )}{\gamma }_k=P_{(\theta ,\lambda )}({\gamma }_k=1)=\frac{\lambda f_1(x_{3k};\theta )}{\lambda f_1(x_{3k};\theta )+(1-\lambda )f_2(x_{3k};\theta )}.$$

When $\theta ={\theta }^{\left(m\right)}$ and $\lambda ={\lambda }^{\left(m\right)}$, the conditional expectation of ${\gamma }_k$ can be the estimate of ${\gamma }_k$.

$${\widehat{\gamma }}_k({\theta }^{\left(m\right)},{\lambda }^{\left(m\right)})=\frac{{\lambda }^{\left(m\right)}{f}_1(x_{3k};{\theta }^{\left(m\right)})}{{\lambda }^{\left(m\right)}{f}_1(x_{3k};{\theta }^{\left(m\right)})+(1-{\lambda }^{\left(m\right)})f_2(x_{3k};{\theta }^{\left(m\right)})},$$

$k=1,2,\cdots ,n_3$.

The log likelihood function is

$$\begin{array}{cc}\hfill & \sum _{i=1}^{n_1}log{f}_1\left(x_{1i};\theta \right)+\sum _{j=1}^{n_2}log{f}_2\left(x_{2j};\theta \right)+\sum _{k=1}^{n_3}{\gamma }_{k}log{f}_1\left(x_{3k};\theta ,\right)\hfill \\ \hfill & +\sum _{k=1}^{n_3}(1-{\gamma }_k)log{f}_2\left(x_{3k};\theta ,\right)+\left(\sum _{k=1}^{n_3}{\gamma }_k\right)log\lambda +\left(n_3-\sum _{k=1}^{n_3}{\gamma }_k\right)log(1-\lambda ).\hfill \end{array}$$

Since the latent variable is unknown, we use its conditional expectation. Besides, the MLE of $\lambda$ is

$${\lambda }^{(m+1)}=\frac{{\sum }_{k=1}^{n_3}{\widehat{\gamma }}_k({\theta }^{\left(m\right)},{\lambda }^{\left(m\right)})}{n_3}.$$

Then in the E-step, we calculate the expectation of new parameters conditional on $({\theta }^{\left(m\right)},{\lambda }^{\left(m\right)})$,

$$\begin{array}{cc}\hfill Q(\theta ,\lambda |{\theta }^{\left(m\right)},{\lambda }^{\left(m\right)})=& \sum _{i=1}^{n_1}log{f}_1\left(x_{1i};\theta \right)+\sum _{j=1}^{n_2}log{f}_2\left(x_{2j};\theta \right)\hfill \\ \hfill & +\sum _{k=1}^{n_3}{\widehat{\gamma }}_k\left(({\theta }^{\left(m\right)},{\lambda }^{\left(m\right)})\right)log{f}_1\left(x_{3k};\theta ,\right)\hfill \\ \hfill & +\sum _{k=1}^{n_3}(1-{\widehat{\gamma }}_k\left(({\theta }^{\left(m\right)},{\lambda }^{\left(m\right)})\right))log{f}_2\left(x_{3k};\theta ,\right).\hfill \end{array}$$

Let ${\widehat{\gamma }}_k={\widehat{\gamma }}_k({\theta }^{\left(m\right)},{\lambda }^{\left(m\right)})$, then

$$\begin{array}{cc}\hfill Q(\theta ,\lambda |{\theta }^{\left(m\right)},{\lambda }^{\left(m\right)})=& \sum _{i=1}^{n_1}\left[-log{\sigma }_1+logf\left(\frac{x_{1i}-{\mu }_1}{{\sigma }_1}\right)\right]\hfill \\ \hfill & +\sum _{j=1}^{n_2}\left[-log{\sigma }_2+logf\left(\frac{x_{2j}-{\mu }_2}{{\sigma }_2}\right)\right]\hfill \\ \hfill & +\sum _{k=1}^{n_3}{\widehat{\gamma }}_k\left[-log{\sigma }_1+logf\left(\frac{x_{3k}-{\mu }_1}{{\sigma }_1}\right)\right]\hfill \\ \hfill & +\sum _{k=1}^{n_3}(1-{\widehat{\gamma }}_k)\left[-log{\sigma }_2+logf\left(\frac{x_{3k}-{\mu }_2}{{\sigma }_2}\right)\right].\hfill \end{array}$$

In the M-step we compute the simultaneous equations below to maximize $Q(\theta ,\lambda |{\theta }^{\left(m\right)},{\lambda }^{\left(m\right)})$. The solutions are the new parameters $({\theta }^{(m+1)},{\lambda }^{(m+1)})$. We give the equations of $({\mu }_1,{\sigma }_1)$, similarly can we get $({\mu }_1,{\sigma }_1)$.

$$\left\{\begin{array}{c}0=\sum _{k=1}^{n_3}\frac{{\gamma }_k}{\lambda f_1+\left(1-\lambda \right)f_2}{f}^{\prime }\left(\frac{x_{3k}-{\mu }_1}{{\sigma }_1}\right)+\sum _{i=1}^{n_1}\left[\frac{1}{f_1}{f}^{\prime }\left(\frac{x_{1i}-{\mu }_1}{{\sigma }_1}\right)\right]\hfill \\ 0=\sum _{k=1}^{n_3}\frac{{\gamma }_k}{\lambda f_1+\left(1-\lambda \right)f_2}\left[{\sigma }_1^2{f}_1+f^{\prime }\left(\frac{x_{3k}-{\mu }_1}{{\sigma }_1}\right)\left(x_{3k}-{\mu }_1\right)\right]+\sum _{i=1}^{n_1}\left[{\sigma }_1^2+\left(x_{1i}-{\mu }_1\right)f^{\prime }\left(\frac{x_{1i}-{\mu }_1}{{\sigma }_1}\right)\frac{1}{f_1}\right]\hfill \\ \lambda =\frac{1}{n_3}\sum _{k=1}^{n_3}{\gamma }_k\hfill \end{array}\right.$$

In the simulation study we consider the normal and logistic cases. Then the maximization step of the normal case can be simplified as

$$\begin{array}{cc}\hfill {\mu }_1& =\frac{{\sum }_{k=1}^{n_3}{\gamma }_k{x}_{3k}+{\sum }_{i=1}^{n_1}{x}_{1i}}{{\sum }_{k=1}^{n_3}{\gamma }_k+n_1},\hfill \\ \hfill {\sigma }_1^2& =\frac{{\sum }_{k=1}^{n_3}{\gamma }_k{(x_{3k}-{\mu }_1)}^2+{\sum }_{i=1}^{n_1}{(x_{1i}-{\mu }_1)}^2}{{\sum }_{k=1}^{n_3}{\gamma }_k+n_1},\hfill \end{array}$$

while that of the logistic case is

$$\left\{\begin{array}{c}0=\sum _{k=1}^{n_3}{\gamma }_k\left(\frac{e^{\frac{x_{3k}-{\mu }_1}{{\sigma }_1}}-1}{e^{\frac{x_{3k}-{\mu }_1}{\sigma }}+1}\right)+\sum _{i=1}^{n_1}\frac{e^{\frac{x_{1i}-{\mu }_1}{\sigma }}-1}{e^{\frac{x_{1i}-{\mu }_1}{{\sigma }_1}}+1}\hfill \\ 0=\sum _{k=1}^{n_3}\left[{\gamma }_k\left(\frac{x_{3k}-{\mu }_1}{{\sigma }_1}\frac{e^{\frac{x_{3k}-{\mu }_1}{{\sigma }_1}}-1}{e^{\frac{x_{3k}-{\mu }_1}{{\sigma }_1}}+1}-1\right)\right]+\sum _{i=1}^{n_1}\left[\frac{x_{1i}-{\mu }_1}{{\sigma }_1}\left(\frac{e^{\frac{x_{1i}-{\mu }_1}{{\sigma }_1}}-1}{e^{\frac{x_{1i}-{\mu }_1}{{\sigma }_1}}+1}\right)-1\right].\hfill \end{array}\right.$$

The two steps are repeated sufficiently to gurantee the convergence. Then we can get the MLE of the parameters.

2.3.2. Normal Case

When the estimate of $\lambda$ is obtained, the posterior distribution (5) can be rewritten as

$$\begin{array}{cc}\hfill \pi \left(\theta \right|X_1,X_2,X_3)& =\left(\frac{1}{{\sigma }_1}\prod _{i=1}^{n_1}f(X_{1i};{\mu }_1,{\sigma }_1)\right)\times \left(\frac{1}{{\sigma }_2}\prod _{j=1}^{n_2}f(X_{2j};{\mu }_2,{\sigma }_2)\right)\hfill \\ \hfill & \times \prod _{k=1}^{n_3}\left[\widehat{\lambda }f(X_{3k};{\mu }_1,{\sigma }_1)+(1-\widehat{\lambda })f(X_{3k};{\mu }_2,{\sigma }_2)\right].\hfill \end{array}$$

(8)

This means that the posterior distribution is equivalent to using the first two terms on the right side of the equation as the “prior distribution” and the third term as the likelihood function. For the first term we have

$$\frac{1}{{\sigma }_1}\prod _{i=1}^{n_1}f(X_i;{\mu }_1,{\sigma }_1)=\frac{1}{{\sigma }_1}{\left(\frac{1}{\sqrt{2\pi }{\sigma }_1}\right)}^{n}exp\left[-\frac{{\sum }_{i=1}^{n_1}{(X_i-{\mu }_1)}^2}{2{\sigma }_1^2}\right].$$

Denote the sample mean and variance by $\overline{X}$ and $S_1^2$, respectively, we have

$$\overline{X_1}=\frac{1}{n_1}\sum _{i=1}^{n_1}{X}_{1i},S_1^2=\frac{1}{n-1}\sum _{i=1}^{n_1}{(X_{1i}-\overline{X_1})}^2,$$

which follows a normal and ${\chi }^2\left(n_1-1\right)$ distribution respectively, that is,

$$\overline{X_1}\sim N\left({\mu }_1,\frac{{\sigma }_1^2}{n_1}\right),\frac{(n_1-1)S_1^2}{{\sigma }_1^2}\sim {\chi }^2\left(n_1-1\right).$$

Let $U\sim N(0,1)$ and $V\sim {\chi }^2\left(n_1-1\right)$ be two independent random variables. Then

$$\overline{X_1}={\mu }_1+\frac{{\sigma }_1}{\sqrt{n_1}}U,(n_1-1)S_1^2={\sigma }_1^{2}V.$$

Given $\overline{X_1}=\overline{x_1}$ and $S_1^2=s_1^2$, then ${\mu }_1$ and ${\sigma }_1$ can be regard as the functions of U and V

$${\mu }_1=\overline{x_1}-\frac{{\sigma }_1}{\sqrt{n_1}}U,{\sigma }_1^2=\frac{(n_1-1)s_1^2}{V}.$$

The joint distribution of $(U,V)$ is

$$\frac{1}{\sqrt{2\pi }}{e}^{-\frac{u^2}{2}}\frac{v^{\frac{n_1-1}{2}-1}}{\Gamma \left(\frac{n_1-1}{2}\right)2^{\frac{n_1-1}{2}}}{e}^{-\frac{v}{2}}.$$

Then the joint distribution of $({\mu }_1,{\sigma }_1)$ can be calculated as

$$\pi ({\mu }_1,{\sigma }_1|x_{obs})=\frac{\sqrt{n_1}}{\sqrt{2\pi }{\sigma }_1}{e}^{-\frac{n_1{(\overline{x_1}-{\mu }_1)}^2}{2{\sigma }_1^2}}\frac{{\left(s_1^2\right)}^{\frac{n_1-1}{2}}{(n_1-1)}^{\frac{n_1-1}{2}}}{\Gamma \left(\frac{n_1-1}{2}\right)2^{\frac{n_1-1}{2}}}{\left(\frac{1}{{\sigma }_1^2}\right)}^{\frac{n_1-1}{2}-1+\frac{3}{2}}{e}^{-\frac{(n_1-1)s_1^2}{2{\sigma }_1^2}},$$

(9)

where ${x_1}_{obs}={(x_1,x_2,\cdots ,x_{n_1})}^{\top }$. This coincides with the joint fiducial density proposed by Fisher [22], which means that the fiducial distribution of $({\mu }_1,{\sigma }_1)$ is

$${\mu }_1|{\sigma }_1^2\sim N\left(\overline{x_1},\frac{{\sigma }_1^2}{n_1}\right),\frac{1}{{\sigma }_1^2}\sim \frac{{\chi }^2\left(n_1-1\right)}{(n_1-1)s_1^2}.$$

(10)

Similarly can we get

$${\mu }_2|{\sigma }_2^2\sim N\left(\overline{x_2},\frac{{\sigma }_2^2}{n_1}\right),\frac{1}{{\sigma }_2^2}\sim \frac{{\chi }^2\left(n_2-1\right)}{(n_2-1)s_2^2},$$

(11)

where $\overline{x_2}$ and $s_2^2$ are the sample mean and variance of the second sample and ${x_2}_{obs}={(x_{21},x_{22},\cdots ,x_{2n_2})}^{\top }$.

With the conclusion above, sampling from the posterior distribution (5) can be done by sampling first from the fiducial distribution of the parameters and then combine the information with the likelihood function of the third sample from the mixture model (1). This can be done simply using the approximate Bayesian computation (ABC) method. In this case, we regard the fiducial distributions of ${({\mu }_1,{\sigma }_1,{\mu }_2,{\sigma }_2)}^{\top }$ as the prior distribution. After we drawn samples of parameters from (10) and (11), denote by ${({\mu }_1^{*},{\sigma }_1^{*},{\mu }_2^{*},{\sigma }_2^{*})}^{\top }$, we generate simulations from the model below and denote them by ${x_3}_{sim}={(x_{31},x_{32},\cdots ,x_{3n_3})}^{\top }$,

$$\widehat{\lambda }f(x;{\mu }_1^{*},{\sigma }_1^{*})+(1-\widehat{\lambda })f(x;{\mu }_2^{*},{\sigma }_2^{*})$$

where $\widehat{\lambda }$ is the MLE of $\lambda$ estimated aforehand using the EM algorithm proposed in the last subsection. Then we calculate the distance between the simulations and the observation and accept those whose distance is below a given threshold $\epsilon$. The algorithm is given below.

• Compute the sample mean and variance of the first two samples and denote them by $\overline{x}$, $s_1^2$, $\overline{y}$ and $s_2^2$. Calculate the MLE of $\lambda$ using EM algorithm and denote it by $\widehat{\lambda }$
• Sample $U_1$ and $U_2$ from the standard normal distribution, $V_1$ from the ${\chi }^2\left(n_1-1\right)$ distribution and $V_2$ from ${\chi }^2\left(n_2-1\right)$. respectively. To sample from the fiducial distributions of the parameters, we calculate ${\mu }_1$, ${\sigma }_1$, ${\mu }_2$ and ${\sigma }_2$ using

  $$\begin{array}{cc}\hfill & {\mu }_1=\overline{x_1}-\frac{U_1}{V_1/\sqrt{n_1-1}}\frac{s_1}{\sqrt{n_1}},{\sigma }_1^2=\frac{(n_1-1)s_1^2}{V_1^2},\hfill \\ \hfill & {\mu }_2=\overline{x_2}-\frac{U_2}{V_2/\sqrt{n_2-1}}\frac{s_2}{\sqrt{n_2}},{\sigma }_2^2=\frac{(n_2-1)s_2^2}{V_2^2}.\hfill \end{array}$$

  We denote the samples of the parameters by ${\theta }^{*}={({\mu }_1^{*},{\mu }_2^{*},{\sigma }_1^{*},{\sigma }_2^{*})}^{\top }$.
• Generate a simulation of size $n_3$ from

  $$\widehat{\lambda }f(x;{\mu }_1^{*},{\sigma }_1^{*})+(1-\widehat{\lambda })f(x;{\mu }_2^{*},{\sigma }_2^{*}).$$

  The simulation is represented by ${x_3}_{sim}={(x_{31}^{*},x_{32}^{*},\cdots ,x_{3n_3}^{*})}^{\top }$.
• Calculate the Euclidean distance between the order statistics of the observation $Z_1,Z_2,\cdots ,Z_{n_3}$ and the simulation $z_1^{*},z_2^{*},\cdots ,z_{n_3}^{*}$. We accept the parameters if the distance is below a given threshold $\epsilon$. Otherwise we reject the parameters.
• The procedure is repeated until we accept a certain number of parameters.

A remark that should be noted in this algorithm is that the samples we get from is an approximation to the posterior distribution $\left(\right)$. We actually samples from

$$\pi (\theta ,\epsilon |{x_1}_{obs},{x_2}_{obs},x)\propto \pi \left(\theta \right|{x_1}_{obs},{x_2}_{obs},{x_3}_{obs})\mathbb{I}(\parallel x-{x_3}_{obs}\parallel \le \epsilon ),$$

(12)

where $\mathbb{I}$ is the indicator function. $\epsilon$ controls the proximity of (12) to (5) and can be adjusted to balance the accuracy and computational cost.

2.3.3. General Case

When $f_1$ and $f_2$ are not normal distributions, to sample from the posterior (8), it is natural to use the Markov chain Monte Carlo (MCMC) method. Metropolisi-Hastings (MH) sampling method and Gibbs sampling method are commonly used. An early version of MH algorithm was given by Metropolis et al. [23] in a statistical physics context, with subsequent generalization by Hastings [24], who focused on statistical problems. Some computational problem and solutions can be further see in Owen and Glynn [25].

The initial values of the parameters can be determined by the EM algorithm mentioned above. For the proposal distribution, we choose

$$\begin{array}{cc}\hfill & q\left({\mu }_k\right|{\mu }_k^{\left(\tau \right)})=N(·;{\mu }_k^{\left(\tau \right)},1)\hfill \\ \hfill & q\left({\sigma }_k\right|{\sigma }_k^{\left(\tau \right)})=Ga(·;{\sigma }_k^{\left(\tau \right)},1)\hfill \end{array}$$

(13)

where $Ga(·)$ and $N(·)$ denote the gamma distribution and normal distribution respectively. $k=1,2$ and ${\mu }_k^{\left(\tau \right)}$, ${\sigma }_k^{\left(\tau \right)}$ denotes the parameters accepted in the $\tau$th loop. After we get ${\theta }^{\left(\tau \right)}={({\mu }_1^{\left(\tau \right)},{\sigma }_1^{\left(\tau \right)},{\mu }_2^{\left(\tau \right)},{\sigma }_2^{\left(\tau \right)})}^{\top }$, we can further obtain ${\theta }^{(\tau +1)}$ via the following two-step algorithm.

• Sample ${({\mu }_1^{*},{\sigma }_1^{*},{\mu }_2^{*},{\sigma }_2^{*})}^{\top }$ respectively from the proposal distribution $\left(13\right)$. Compute

  $$\begin{array}{cc}\hfill logQ\left({\theta }^{*}\right|{\theta }^{\tau })& =logq\left({\mu }_1^{*}\right|{\mu }_1^{\left(\tau \right)})+logq\left({\mu }_2^{*}\right|{\mu }_2^{\left(\tau \right)})+logq\left({\sigma }_1^{*}\right|{\sigma }_1^{\left(\tau \right)})+logq\left({\sigma }_2^{*}\right|{\sigma }_2^{\left(\tau \right)})\hfill \\ \hfill & +\sum _{i=1}^{n_1}logf(x_i;{\mu }_1^{*},{\sigma }_1^{*})+\sum _{i=1}^{n_2}logf(y_i;{\mu }_2^{*},{\sigma }_2^{*})\hfill \\ \hfill & +\sum _{i=1}^{n_3}\left[\widehat{\lambda }f(z_i;{\mu }_1^{*},{\sigma }_1^{*})+(1-\widehat{\lambda })f(z_i;{\mu }_2^{*},{\sigma }_2^{*})\right]-log\left({\sigma }_1^{*}{\sigma }_2^{*}\right).\hfill \end{array}$$
• Accept ${\theta }^{*}$ with probability

  $$P({\theta }^{*},{\theta }^{\left(\tau \right)})=exp\left\{min\left[0,logQ\left({\theta }^{*}\right|{\theta }^{\tau })-logQ\left({\theta }^{\left(\tau \right)}\right|{\theta }^{*})\right]\right\}.$$

  and let ${\theta }^{(\tau +1)}={\theta }^{*}$. Otherwise we reject the parameters and return to the first step.

The algorithm should be repeated sufficiently before obtain the samples from the posterior distribution. This cost much more time compare with the ABC algorithm for the normal case. What’s more, in our simulation we found that the MH algorithm may be too conservative. A better subsititution can be the two-block Gibbs sampling proposed by Ren et al. [18]. In this sampling method, $\lambda$ is first estimiated using the EM algorithm, then for each loop, the parameters are updated by the conditional generalized pivotal quantities.

3. Real Data Example

In this section we apply the proposed posterior p-value to the real halibut dataset studied by Hosmer [1], which is given by International Halibut Commission in Seattle, Washington. This dataset consists of the lengths of 208 halibut caught on one of their research cruises, in which 134 are female while the rest 74 are male. The data is summarized by Karunamuni and Wu [26] and given in Table 1. We follow their method and randomly select 14 males and 26 females from the samples and regard them as the first and second sample of the mixture model 1. Then the remain male proportion of 60/168 is approximately indentical to the original male proportion of 74/208, which is 0.3558. 100 replications is generated with the same procedure. Hosmer [1] pointed out that the component for the dataset can be fitted by the normal distribution. A problem of interest is whether the sex effects the length of the halibut.

To test the homogeneity, for each replication we first use the EM algorithm to estimate $\lambda$, then we use the reject-ABC method to generate 8000 samples. We choose a moderate threshold $\epsilon$ to balance the accuracy and the computational cost. For the 100 replications, the mean of estimate of the male proportion $\lambda$ is 0.3381, with the mean squared error 0.0045, which illustrate the accuracy of our EM algorithm The estimates of the location and scale parameter of the male halibut are $\widehat{{\mu }_1}=96.655$ and $\widehat{{\sigma }_1}=12.983$ while that of the female ones are $\widehat{{\mu }_2}=118.806$ and $\widehat{{\sigma }_2}=9.077$. This is close to the estimates of Ren et al. [18]. As with the hypothesis testing $f_1=f_2$, we calculate the posterior p-value of the 100 replications. Given the significance level $\alpha =0.05$, all the p-values are less than $\alpha$. Thus the null hypothesis is rejected, which indicates that there exists association between the sex and length of the halibut.

4. Simulation Study

In this section we give the simulation study of the cases discussed above. We compare the results of the posterior p-value (7) using different sampling methods and the generalized fiducial method proposed by Ren et al. [18]. As we can see from the simulations, the posterior p-value we proposed largely improves the testing of homogeneity.

4.1. Normal Case

When $f_1$ and $f_2$ are normal distributions, we compare the results of three different tests. The first two are the posterior p-value we proposed, but using two-block Gibbs sampling and reject-ABC sampling method, respectively. The last one is the generalized fiducial method proposed by Ren et al. [18]. In the following tables, the first two are denoted by “$T_G$” and “$T_R$”, while the last is denoted by “G”. We fix $f_1$ to $N(0,1)$ while $f_2$ is set to be $N(0,1)$, $N(1,1)$, $N(0,1.5^2)$ and $N(1,1.5^2)$. For each $f_1$ and $f_2$ we consider $\lambda =0.3,0.5,0.7$ and different sample sizes for $n_1,n_2$ and $n_3$. We simulate $N=10000$ repetitions for each case. For the Gibbs sampling, we accept 3000 samples after burning in the first 2000. For the reject-ABC sampling method, we first calculate the estimate of $\lambda$ and accept 4000 parameters with $\epsilon$ set to $\sqrt{n_3}/2$. Then we calculate the posterior p-value using the samples. We set the significance level to $alpha=0.05$ and reject the null hypothesis when the posterior p-value is below $\alpha$. The results are shown in Table 2, Table 3, Table 4 and Table 5. We further give the QQ-plot of $T_R$ in Figure 1, which indicates the correctness of Theorem 2. The first row are the cases of $(n_1,n_2,n_3)$=$(10,10,10)$, $(20,20,20)$ and $(30,30,30)$, while the second row are the cases of $(10,20,30)$, $(30,20,10)$ and $(15,25,150)$.

We can see from the results that the posterior p-value largely improves the testing of homogeneity in normal cases. The Type-I error is controlled as well as the generalized fiducial methods. Moreover, our method significantly improves the power of testing homogeneity, especially when $\sigma$ is different. The reject-ABC sampling method has the advantage of lower computational cost, comparing with the two-block Gibbs sampling method. However, the power of using reject-ABC sampling method is smaller than using two-block Gibbs sampling when $n_3$ is much larger than $n_1$ and $n_2$. Thus, we can use reject-ABC sampling method when the sample size is small or moderate and two-block Gibbs sampling when the sample size is large.

4.2. General Case

For general case we assume that $f_1$ and $f_2$ are logistic distributions. The location and scale parameters of $f_1$ and $f_2$ are set the same as that of the normal case. We simulate 10000 repetitions for each sample size. We propose three methods in this simulation. The first two are the generalized fiducial method proposed by Ren et al. [18] and our posterior p-value using two-blocks Gibbs sampling. They are denoted by “G” and “$T_G$” as last simulation. The last one is the posterior p-value using M-H algorithm, which is denoted by “$T_M$”. First we calculate the MLE of $\lambda$ using the EM algorithm. Then we propose the Metropolis-Hastings algorithm to obtain 12000 samples after first burn-in 8000 ones. To avoid the dependency between the samples, we choose the first one in every three samples, which leaves us 4000 samples. Then we use these samples to calculate the posterior p-value. The algorithm is natural and seems to be feasible. However, from the Table 6 we can see that with this sampling method the results are rather conservative. Given the significance level $\alpha =0.05$, the type-I error of $T_M$ is always much more smaller than 0.05, which makes the power of $T_M$ also smaller than the rest two when $f_1\ne f_2$. However, We find that the two-block Gibbs sampling method can successfully solve the problem. It can be seen that the type-I error of $“T_G^{″}$ can be controlled well while the power is largerly improved, compares with the genralized fiducial method. The results are shown in Table 7, Table 8 and Table 9. We also give the QQ-plot of $“T_G^{″}$ in Figure 2.

5. Conclusions

In this paper, we propose a new posterior p-value for testing the homogeneity of the three-sample problem. We define the regular location-scale family and assume that both $f_1$ and $f_2$ are in the same family. Then testing the homogeneity is equivalent to testing the equality of the location and scale parameters. We use the Bayes’ theorem to obtain the posterior distribution of the parameters and propose the Bernstein-von Mises theorem for multiple samples. Then we propose the posterior p-value for testing the equality of the parameters. To sample from the posterior distribution, we compare different sampling methods. The simulation studies illustrate that reject-ABC sampling method may be a good choice for the normal case while the two-block Gibbs sampling is better for the general ones. It should be noticed that we transform the hypotheses of homogeneity to the hypotheses (6). Then with different matrix $\mathbf{A}$, we can generate our method to a variety of hypotheses.