Assessing Multinomial Distributions with a Bayesian Approach

Abstract

This paper introduces a unified Bayesian approach for testing various hypotheses related to multinomial distributions. The method calculates the Kullback–Leibler divergence between two specified multinomial distributions, followed by comparing the change in distance from the prior to the posterior through the relative belief ratio. A prior elicitation algorithm is used to specify the prior distributions. To demonstrate the effectiveness and practical application of this approach, it has been applied to several examples.

1. Introduction

Multinomial distribution tests are a crucial statistical tool in many fields, especially when data can be categorized into multiple groups. These tests were first proposed by Karl Pearson in 1890 and have since been widely used to analyze and make inferences about the probabilities or proportions associated with each category in the multinomial distribution [1].

Let the sample space $\mathcal{A}$ of a random experiment be the union of a finite number k of mutually disjoint sets (categories) $A_1,\dots ,A_k$. Assume that $P\left(A_j\right)={\theta }_j$, $j=1,\dots ,k$, where ${\sum }_{j=1}^k{\theta }_j=1$. Here ${\theta }_j$ represents the probability that the outcome is an element of the set $A_j$. The random experiment is to be repeated n independent times. Define the random variables $Y_j$ to be the number of times the outcome is an element of set $A_j$, $j=1,\dots ,k$. That is, $Y_1,\dots ,Y_k=n-Y_1-Y_2-\cdots -Y_{k-1}$ denote the frequencies with which the outcome belongs to $A_1,\dots ,A_k$, respectively. Then the joint probability mass function (pmf) of $Y_1,\dots ,Y_k$ is the multinomial with parameters $n,{\theta }_1,\dots ,{\theta }_k$ [2]. It is desired to test the null hypothesis:

$$\begin{array}{c}\hfill H_0^1:{\theta }_j={\theta }_{j0}, \mathrm{for} j=1,\dots ,k\end{array}$$

(1)

against all alternatives, where ${\theta }_{j0}$ are known constants. Within the classical frequentist framework, to test $H_0^1$, it is common to use the test statistic [3]:

$$\begin{array}{c}\hfill {\chi }^2=\sum _{j=1}^k\frac{{(Y_j-n{\theta }_{j0})}^2}{n{\theta }_{j0}}.\end{array}$$

(2)

It is known that, under $H_0^1$, the limiting distribution of ${\chi }^2$ is chi-squared with $k-1$ degrees of freedom. When $H_0^1$ is true, $n{\theta }_{j0}$ represents the expected value of $Y_j$. This implies the observed value ${\chi }^2$ should not be too large if $H_0^1$ is true. For a given significance level $\alpha$, an approximate test of size $\alpha$ is to reject $H_0^1$ if the observed ${\chi }^2>{\chi }_{k-1}^2\left(\alpha \right)$, where the ${\chi }_{k-1}^2\left(\alpha \right)$ is the $1-\alpha$ quantile of the chi-squared distribution with $k-1$ degrees of freedom; otherwise, fail to reject $H_0^1$. Other possible tests for $H_0^1$ include Fisher’s exact test and likelihood ratio tests [4].

If there are r independent samples, then the interest is to test whether the r samples come from the same multinomial population or that r multinomial populations are different. Let $A_1,A_2,\dots ,A_k$ denote k possible types of categories in the ith sample, $i=1,\dots ,k$. Let the probability that an outcome of category $A_j$ will occur for the ith population (or ith sample) be denoted by ${\theta }_{j|i}$. Note that, ${\sum }_{j=1}^k{\theta }_{j|i}=1$ for each $i=1,\dots ,r$. Moreover, let $Y_{j|i}$ be the number of times the outcome is an element of $A_j$ in sample i. Consider the completely specified hypothesis:

$$\begin{array}{c}\hfill H_0^2:{\theta }_{j|i}={\theta }_{j0|i}, \mathrm{for} j=1,2,\dots ,k.\end{array}$$

(3)

Under $H_0^2$, the test statistic in (2) can be extended to

$$\begin{array}{c}\hfill {\chi }^2=\sum _{i=1}^r{\chi }_i^2=\sum _{i=1}^r\sum _{j=1}^k\frac{{(Y_{j|i}-n_i{\theta }_{j0|i})}^2}{n_i{\theta }_{j0|i}}.\end{array}$$

(4)

If $H_0^2$ is true, then ${\chi }^2$ in (4) has an approximately chi-squared distribution with $r(k-1)$ degrees of freedom. Likewise, for a given significance level $\alpha$, an approximate test of size $\alpha$ is to reject $H_0^2$ if the observed ${\chi }^2$ is bigger than ${\chi }_{r(k-1)}^2\left(\alpha \right)$; otherwise, fail to reject $H_0^2$ [5].

A third and more common hypothesis is to test whether the r multinomial populations are the same without specifying the values of the ${\theta }_{j|i}$. That is, we consider the null hypothesis:

$$\begin{array}{c}\hfill H_0^3:{\theta }_{j|1}={\theta }_{j|2}=\cdots ={\theta }_{j|r}={\theta }_j, \mathrm{for} j=1,2,\dots ,k.\end{array}$$

(5)

The test statistics to test $H_0^3$ are given by

$$\begin{array}{c}\hfill {\chi }^2=\sum _{i=1}^r{\chi }_i^2=\sum _{i=1}^r\sum _{j=1}^k\frac{{(Y_{j|i}-n_i{\widehat{\theta }}_{j|i})}^2}{n_i{\widehat{\theta }}_{j|i}},\end{array}$$

(6)

where ${\widehat{\theta }}_{j|i}=\frac{{\sum }_{i=1}^r{Y}_{j|i}}{{\sum }_{i=1}^r{n}_i}=\frac{{\sum }_{i=1}^r{Y}_{j|i}}{n}$. Here, $n_i$ denotes the sample size of sample i and ${\sum }_{i=1}^r{Y}_{j|i}$ represents the total in category $A_j$. Note that ${\widehat{\theta }}_{j|i}$ represents the pooled maximum likelihood estimator (MLE) of ${\theta }_j$ under $H_0^3$. It is known that the limiting distribution of ${\chi }^2$ in (6) is a chi-squared distribution with $(r-1)(k-1)$ degrees of freedom. So, for a given significance level $\alpha$, an approximate test of size $\alpha$ is to reject $H_0^3$ if the observed ${\chi }^2>{\chi }_{(r-1)(k-1)}^2\left(\alpha \right)$; otherwise, fail to reject $H_0^3$ [5]. It is worth mentioning that several other frequentist methods for testing the multinomial distribution have been proposed, utilizing different distance measures. These methods include the Euclidean distance proposed by [6], the smooth total variation distance introduced by [7], and $\varphi$-divergences discussed by [8]. These approaches provide alternative ways to assess the goodness-of-fit of the multinomial distribution using distance metrics.

Refs. [9,10,11,12,13] made early advances in Bayesian methods for analyzing categorical data, focusing on smoothing proportions in contingency tables and inference about odds ratios, respectively. These methods typically employed conjugate beta and Dirichlet priors. Ref. [14,15] extended these methods to develop Bayesian analogs of small-sample frequentist tests for $2\times 2$ tables, also using such priors. Ref. [16] recommended the use of the uniform prior for predictive inference, but other priors were also suggested by discussants of his paper. The Jeffreys prior is the most commonly used prior for binomial inference, partially due to its invariance to the scale of measurement for the parameter. Reference priors (see [17]), such as the Jeffreys prior for the binomial parameter (see [18]), are viable options, but their specification can be computationally complex. Ref. [10] may have been the first to utilize an empirical Bayesian approach with contingency tables, estimating parameters in gamma and log-normal priors for association factors. Empirical Bayes involves estimating the prior distribution from the observed data itself and is particularly useful when dealing with large amounts of data. Refs. [19,20] derived integral expressions for the posterior distributions of the difference, ratio, and odds ratio under independent beta priors. Ref. [19] introduced Bayesian highest posterior density (HPD) confidence intervals for these measures. The HPD approach ensures that the posterior probability matches the desired confidence level, and the posterior density is higher inside the interval than outside. Ref. [21] discussed Bayesian confidence intervals for association parameters in $2\times 2$ tables. They argued that to achieve good coverage performance in the frequentist sense across the entire parameter space, it is advisable to use relatively diffuse priors. Even uniform priors are often too informative, and they recommended the use of the Jeffreys prior. Bayesian methods for analyzing categorical data have been extensively surveyed in the literature, including comprehensive reviews by [22,23] with a focus on contingency table analysis. Refs. [24,25,26] proposed tests based on Bayesian nonparametric methods using Dirichlet process priors.

We build on the recent work of [27] by extending their Bayesian approach for hypothesis testing on one-sample proportions based on Kullback–Leibler divergence and relative belief ratio, using a uniform (0, 1) prior on binomial proportions, to multinomial distributions. Our goal is to provide a comprehensive Bayesian approach for testing hypotheses $H_0^1$, $H_0^2$, and $H_0^3$. We derive distance formulas and use the Dirichlet distribution as a prior on probabilities. To ensure proper values of the prior’s hyperparameters, we employ the elicitation algorithm developed by [28]. The proposed approach offers several advantages, including computational simplicity, ease of interpretation, evidence in favor of the null hypothesis, and no requirement to specify a significance level.

The paper is structured as follows. Section 2 provides an overview of the relative belief ratio inference and KL divergence. Section 3 details the proposed approach, including the formulas and computational algorithms. In Section 4, several examples are presented to illustrate the approach. Finally, Section 5 contains concluding remarks and discussions.

2. Relevant Background

2.1. Inferences Using Relative Belief

Ref. [29] introduced the relative belief ratio, which has become a popular tool in statistical hypothesis testing theory. Several works have employed this approach, including [30,31,32,33,34,35].

Suppose we have a statistical model with a density function $\{f_{\theta }\left(\mathit{y}\right):\theta \in \mathsf{\Theta }\}$ with respect to the Lebesgue measure on the parameter space $\mathsf{\Theta }$. Let $\pi \left(\theta \right)$ be a prior on $\mathsf{\Theta }$. After observing the data $\mathit{y}$, the posterior distribution of $\theta$ can be expressed as

$$\pi \left(\theta \right|\mathit{y})=\frac{f_{\theta }\left(\mathit{y}\right)\pi \left(\theta \right)}{m\left(\mathit{y}\right)},$$

where $m\left(\mathit{y}\right)={\int }_{\mathsf{\Theta }}{f}_{\theta }\left(\mathit{y}\right)\pi \left(\theta \right)d\theta$.

Assume that the goal is to draw inferences about the parameter $\theta$. If the prior $\pi (·)$ and the posterior $\pi (·|\mathit{y})$ are continuous at $\theta$, then the relative belief ratio for a hypothesized value ${\theta }_0$ of $\theta$ can be expressed as follows:

$$RB\left({\theta }_0\right|\mathit{y})=\frac{\pi \left({\theta }_0\right|\mathit{y})}{\pi \left({\theta }_0\right)},$$

the ratio of the posterior density to the prior density at ${\theta }_0$. In other words, $RB\left({\theta }_0\right|\mathit{y})$ quantifies how the belief in ${\theta }_0$ being the true value has changed from prior to posterior. It is worth noting that when $\pi (·)$ and $\pi (·|\mathit{y})$ are discrete, the relative belief ratio is defined through limits, and further details can be found in [29].

The relative belief ratio $RB\left({\theta }_0\right|\mathit{y})$ provides a measure of evidence for ${\theta }_0$ being the true value. A value of $RB\left({\theta }_0\right|\mathit{y})>1$ indicates evidence in favor of ${\theta }_0$ being the true value, whereas $RB\left({\theta }_0\right|\mathit{y})<1$ indicates evidence against ${\theta }_0$ being the true value. If $RB\left({\theta }_0\right|\mathit{y})=1$, there is no evidence in either direction.

Once the relative belief ratio is calculated, it is important to determine the strength of the evidence in favor of or against $H_0:\theta ={\theta }_0$. A common way to quantify this is by computing the tail probability [29]:

$$Str\left({\theta }_0\right|\mathit{y})=\Pi (RB\left(\theta \right|\mathit{y})\le RB\left({\theta }_0\right|\mathit{y})|\mathit{y})={\int }_{\{\theta \in \mathsf{\Theta }:RB\left(\theta \right|\mathit{y})\le RB\left({\theta }_0\right|\mathit{y})\}}\pi (·|\mathit{y})d\theta ,$$

(7)

where $\Pi (·|\mathit{y})$ in (7) is the posterior cumulative distribution function with posterior density $\pi (·|\mathit{y})$. Therefore, equation (7) represents the posterior probability that the true value of $\theta$ has a relative belief ratio no greater than that of the hypothesized value ${\theta }_0$. When $RB\left({\theta }_0\right|\mathit{y})<1$, there is evidence against ${\theta }_0$. A small value of $Str\left({\theta }_0\right|\mathit{y})$ indicates a high posterior probability that the true value has a relative belief ratio greater than $RB\left({\theta }_0\right|\mathit{y})$, indicating strong evidence against ${\theta }_0$. Conversely, when $RB\left({\theta }_0\right|\mathit{y})>1$, there is evidence in favor of ${\theta }_0$. A large value of $Str\left({\theta }_0\right|\mathit{y})$ indicates a low posterior probability that the true value has a relative belief ratio greater than $RB\left({\theta }_0\right|\mathit{y})$, indicating strong evidence in favor of ${\theta }_0$. A small value of $Str\left({\theta }_0\right|\mathit{y})$ indicates weak evidence in favor of ${\theta }_0$.

2.2. KL Divergence

The KL divergence, also referred to as relative entropy, is a measure of dissimilarity between two probability distributions that quantifies how far apart they are from each other. It was introduced by Solomon Kullback and Richard Leibler in 1951. Let P and Q be two discrete cumulative distribution functions (cdf’s) on the same probability space $\phi$, with corresponding probability mass functions (pmf’s) p and q (with respect to the counting measure). The KL divergence between p and q is given by:

$$\begin{array}{c}\hfill d(p,q)=\sum _{x\in \phi }p\left(x\right)\log \left(\frac{p\left(x\right)}{q\left(x\right)}\right).\end{array}$$

The KL divergence is always non-negative, and it attains its minimum value when $p=q$ almost surely. This property makes it a useful tool in many areas of machine learning and information theory, such as hypothesis testing, model selection, and clustering. One interpretation of the KL divergence is that it measures how much information is lost when using Q to approximate P. It is worth noting that the KL divergence is not symmetric: $d(p,q)$ and $d(q,p)$ are generally not equal. Therefore, it is important to specify which distribution is the “true” or “target” distribution and which is the “approximating” or “predicted” distribution for some applications when using KL divergence in practice, as noted by [36].

The following lemma is essential to the proposed approach.

Lemma 1. 

Let $p(y_1,y_2,\dots ,y_r)=p_1\left(y_1\right)p_2\left(y_2\right)\cdots p_r\left(y_r\right)$ and $q(y_1,y_2,\dots ,y_r)=$$q_1\left(y_1\right)q_2\left(y_2\right)\cdots q_r\left(y_r\right)$, where $p_i\left(y_i\right)$ and $q_i\left(y_i\right)$ are probability mass functions with supports $y_i=1,\dots ,n_i,i=1,\dots ,r$. Then

$$d(p,q)=\sum _{i=1}^{r}d(p_i,q_i).$$

Proof. 

We have

$$\begin{array}{ccc}\hfill d(p,q)& =& \sum _{y_1=0}^{n_1}\cdots \sum _{y_r=0}^{n_r}p(y_1,y_2,\dots ,y_r)\log \frac{p(y_1,y_2,\dots ,y_r)}{q(y_1,y_2,\dots ,y_r)}\hfill \\ & =& \sum _{y_1=0}^{n_1}\cdots \sum _{y_r=0}^{n_r}{p}_1\left(y_1\right)\cdots p_r\left(y_r\right)\log \frac{p_1\left(y_1\right)\cdots p_r\left(y_r\right)}{q_1\left(y_1\right)\cdots q_r\left(y_r\right)}\hfill \\ & =& \sum _{y_1=0}^{n_1}\cdots \sum _{y_r=0}^{n_r}{p}_1\left(y_1\right)\cdots p_r\left(y_r\right)[\log p_1\left(y_1\right)+\cdots +\log p_r\left(y_r\right)-\log q_1\left(y_1\right)-\cdots -\log q_r\left(y_r\right)]\hfill \\ & =& \sum _{y_1=0}^{n_1}\cdots \sum _{y_r=0}^{n_r}{p}_1\left(y_1\right)\cdots p_r\left(y_r\right)\log p_1\left(y_1\right)+\cdots +\sum _{y_1=0}^{n_1}\cdots \sum _{y_r=0}^{n_r}{p}_1\left(y_1\right)\cdots p_r\left(y_r\right)\log p_r\left(y_r\right)\hfill \\ & & -\sum _{y_1=0}^{n_1}\cdots \sum _{y_r=0}^{n_r}{p}_1\left(y_1\right)\cdots p_r\left(y_r\right)\log q_1\left(y_1\right)-\cdots -\sum _{y_1=0}^{n_1}\cdots \sum _{y_r=0}^{n_r}{p}_1\left(y_1\right)\cdots p_r\left(y_r\right)\log q_r\left(y_r\right).\hfill \end{array}$$

Since, for $i=1,\dots ,r$, ${\sum }_{y_i=1}^{n_i}{p}_i\left(y_i\right)={\sum }_{y_i=1}^{n_i}{q}_i\left(y_i\right)=1$, we have

$$\begin{array}{ccc}\hfill d(p,q)& =& \sum _{y_1=0}^{n_1}{p}_1\left(y_1\right)\log p_1\left(y_1\right)+\cdots +\sum _{y_r=0}^{n_r}{p}_r\left(y_r\right)\log p_r\left(y_r\right)\hfill \\ & & -\sum _{y_1=0}^{n_1}{p}_1\left(y_1\right)\log q_1\left(y_1\right)-\cdots -\sum _{y_r=0}^{n_r}{p}_r\left(y_r\right)\log q_r\left(y_r\right)\hfill \\ & =& \sum _{y_1=0}^{n_1}{p}_1\left(y_1\right)\log \frac{p_1\left(y_1\right)}{q_1\left(y_1\right)}+\cdots +\sum _{y_r=0}^{n_r}{p}_r\left(y_r\right)\log \frac{p_r\left(y_r\right)}{q_r\left(y_r\right)}\hfill \\ & =& d(p_1,q_1)+\cdots +d(p_r,q_r).\hfill \end{array}$$

□

3. The Approach

3.1. Bayesian One-Sample Multinomial

Let $\mathit{Y}=(Y_1,\dots ,Y_k)\sim multinomial(n,{\theta }_1,\dots ,{\theta }_k)$. The joint pmf of $Y_1,\dots ,Y_k$ is given by

$$\begin{array}{c}\hfill p(y_1,\dots ,y_k)=\left(\genfrac{}{}{0pt}{}{n}{y_1,y_2,\dots ,y_k}\right)\prod _{j=1}^k{\theta }_j^{y_j},\end{array}$$

(8)

where $\left(\genfrac{}{}{0pt}{}{n}{y_1,y_2,\dots ,y_k}\right)=\frac{n!}{y_1!\cdots y_k!}$, ${\sum }_{j=1}^k{\theta }_j=1$, and ${\sum }_{j=1}^k{y}_j=n.$

To test the null hypothesis $H_0^1$ as defined in (1), we first compute the Kullback–Leibler (KL) divergence between $p(y_1,\dots ,y_k)$ and the pmf under $H_0^1$ represented by

$$\begin{array}{ccc}\hfill q(y_1,\dots ,y_k)& =& \left(\genfrac{}{}{0pt}{}{n}{y_1,y_2,\dots ,y_k}\right)\prod _{j=1}^k{\theta }_{j0}^{y_i}.\hfill \end{array}$$

(9)

Here, ${\theta }_{j0}$ denotes the null hypothesis value for ${\theta }_j$. The following proposition provides the formula for the KL divergence between p and q.

Proposition 1. 

Let $p(y_1,\dots ,y_k)$ and $q(y_1,\dots ,y_k)$ be two joint probability mass functions as defined in (8) and (9), respectively. We have,

$$\begin{array}{ccc}\hfill d(p,q)& =& n\sum _{j=1}^k\left[{\theta }_j\log \right(\frac{{\theta }_j}{{\theta }_{j0}}\left)\right].\hfill \end{array}$$

Proof. 

Let the support of $Y_j$ be $1,2,\dots ,n_j$, $j=1,\dots ,k$. We have

$$\begin{array}{ccc}\hfill d(p,q)& =& \sum _{y_1=0}^{n_1}\cdots \sum _{y_k=0}^{n_k}p(y_1,\dots ,y_k)\log \frac{p(y_1,\dots ,y_k)}{q(y_1,\dots ,y_k)}\hfill \\ \hfill & =& \sum _{y_1=0}^{n_1}\cdots \sum _{y_k=0}^{n_k}p(y_1,\dots ,y_k)\log \frac{\left(\genfrac{}{}{0pt}{}{n}{y_1,y_2,\cdots ,y_k}\right){\prod }_{j=1}^k{\theta }_j^{y_j}}{\left(\genfrac{}{}{0pt}{}{n}{y_1,y_2,\cdots ,y_k}\right){\prod }_{j=1}^k{\theta }_{j0}^{y_j}}\hfill \\ \hfill & =& \sum _{y_1=0}^{n_1}\cdots \sum _{y_k=0}^{n_k}p(y_1,\dots ,y_k)\log \prod _{j=1}^k{\left[\frac{{\theta }_j}{{\theta }_{j0}}\right]}^{y_j}.\hfill \end{array}$$

Using the properties of logarithmic function, we get

$$\begin{array}{ccc}\hfill d(p,q)& =& \sum _{y_1=0}^{n_1}\cdots \sum _{y_k=0}^{n_k}p(y_1,\dots ,y_k)\sum _{j=1}^k{y}_j\log \left[\frac{{\theta }_j}{{\theta }_{j0}}\right]\hfill \\ \hfill & =& \sum _{y_1=0}^{n_1}\cdots \sum _{y_k=0}^{n_k}p(y_1,\dots ,y_k)\times y_1\times \log \left[\frac{{\theta }_1}{{\theta }_{01}}\right]\hfill \\ \hfill & & \cdots +\sum _{y_1=0}^{n_1}\cdots \sum _{y_k=0}^{n_k}p(y_1,\dots ,y_k)\times y_k\times \left[\frac{{\theta }_k}{{\theta }_{k0}}\right]\hfill \\ \hfill & =& \mathbb{E}\left[Y_1\right]\times \log \left[\frac{{\theta }_1}{{\theta }_{01}}\right]+\cdots +\mathbb{E}\left[Y_k\right]\times \log \left[\frac{{\theta }_k}{{\theta }_{k0}}\right]\hfill \\ \hfill & =& \sum _{j=1}^k\mathbb{E}\left[Y_j\right]\log \left[\frac{{\theta }_j}{{\theta }_{j0}}\right].\hfill \end{array}$$

Since the marginal probability mass function of $Y_j,j=1,\dots ,k$, is the binomial with parameters n and ${\theta }_j$, we get

$$\begin{array}{ccc}\hfill d(p,q)& =& \sum _{j=1}^{k}n{\theta }_j\log \left[\frac{{\theta }_j}{{\theta }_{j0}}\right]=n\sum _{j=1}^k\left[{\theta }_j\log \right(\frac{{\theta }_j}{{\theta }_{j0}}\left)\right].\hfill \end{array}$$

□

To connect the distance formula presented in Proposition 1 with the test statistic ${\chi }^2$ in (2), we use the Taylor series expansion of the function $f\left(x\right)=x\log \frac{x}{x_0}$ about $x_0$. This gives us

$$\begin{array}{c}\hfill f\left(x\right)=(x-x_0)+0.5{(x-x_0)}^2\frac{1}{x_0}+\cdots \end{array}$$

If $H_0^1$ is true and n is large, then we can approximate the distance $d(p,q)$ as

$$\begin{array}{ccc}\hfill d(p,q)& \approx & n\sum _{j=1}^k({\theta }_j-{\theta }_{j0})+0.5n\sum _{j=1}^k\frac{{({\theta }_j-{\theta }_{j0})}^2}{{\theta }_{j0}}.\hfill \end{array}$$

(10)

Since the probabilities sum to 1, the first term in (10) equals 0. The second term in (10) can be expressed as

$$\begin{array}{c}\hfill 0.5\sum _{j=1}^k\frac{{(n{\theta }_j-n{\theta }_{j0})}^2}{n{\theta }_{j0}}=0.5\sum _{j=1}^k\frac{{(E\left(Y_j\right)-n{\theta }_{j0})}^2}{n{\theta }_{j0}}.\end{array}$$

This shows a direct connection between the KL divergence and ${\chi }^2$.

For the proposed Bayesian test, the probabilities ${\theta }_1,\dots ,{\theta }_k$ are now random. The suggested prior for the joint probabilities $({\theta }_1,\dots ,{\theta }_k)$ is the Dirichlet distribution with parameters ${\alpha }_1,\dots {\alpha }_k$. That is,

$$\begin{array}{c}\hfill ({\theta }_1,{\theta }_2,\dots ,{\theta }_k)\sim \mathrm{Dirichlet}({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_k).\end{array}$$

(11)

To elicit the prior, we use the elicitation algorithm developed by [28], which requires some domain knowledge to provide a lower bound for each ${\theta }_i$. For convenience, we have made this algorithm available on Shiny at the following link: https://bayesian-chi-square-test.shinyapps.io/DirichletprocessKyusonlim/ (accessed on 23 May 2023). For comparison purposes, we also considered the non-informative (uniform) prior, and Jeffreys prior; see Section 4. For the proposed Bayesian approach, when $({\theta }_1,\dots ,{\theta }_k)$ has the prior defined in (11), we put

$$\begin{array}{c}\hfill D=n\sum _{j=1}^k\left[{\theta }_j\log \right(\frac{{\theta }_j}{{\theta }_{j0}}\left)\right].\end{array}$$

(12)

We also have that the posterior distribution of $({\theta }_1,\dots ,{\theta }_k)$ given the observed data $\mathit{y}=(y_1,\dots ,y_k)$ is $\mathrm{Dirichlet}({\alpha }_1+y_1,{\alpha }_2+y_2,\dots ,{\alpha }_k+y_k)$. We write

$$\begin{array}{c}\hfill D_y=n\sum _{j=1}^k\left[{\theta }_j\log \right(\frac{{\theta }_j}{{\theta }_{j0}}\left)\right].\end{array}$$

(13)

Note that,

$$E\left({\theta }_j\right|\mathit{y})=\frac{{\alpha }_j+y_j}{{\sum }_{j=1}^k({\alpha }_j+y_j)}=\frac{{\alpha }_j+y_j}{n+{\sum }_{j=1}^k{\alpha }_j}=\frac{n}{n+{\sum }_{j=1}^k{\alpha }_j}\frac{y_j}{n}+(1-\frac{n}{n+{\sum }_{j=1}^k{\alpha }_j})\frac{{\alpha }_j}{{\sum }_{j=1}^k{\alpha }_j},$$

which is a convex combination between the sample proportion (MLE) and the prior mean. As $n\to \infty$, the weak law of large numbers ensures that $E\left({\theta }_j\right|\mathit{y})$ converges in probability to the true value of ${\theta }_j$. Hence, if $H_0^1$ is true, then $D_y\stackrel{a.s.}{\to }0$. Conversely, if $H_0^1$ is false, then $D_y\stackrel{a.s.}{\to }c$, where $c>0$. Proposition 1 establishes that $d(p,q)=0$ if and only if ${\theta }_j={\theta }_{j0}$. Therefore, testing $H_0^1:{\theta }_j={\theta }_{j0}$ is equivalent to testing $d(p,q)=0$. It follows that when $H_0^1$ is true, the distribution of $D_y$ should be more concentrated around 0 than that of D. So, the proposed test involves comparing the distributions of D and $D_y$ around 0 using the relative belief ratio:

$$\begin{array}{c}\hfill R{B}_D\left(0\right|\mathit{y})=\frac{{\pi }_D\left(0\right|\mathit{y})}{{\pi }_D\left(0\right)},\end{array}$$

(14)

where ${\pi }_D\left(0\right)$ and ${\pi }_D\left(0\right|\mathit{y})$ represent the probability density functions of D and $D_y$, respectively. If $RBD\left(0\right|\mathit{y})>1$, it provides evidence in favor of $H_0$ (since the distribution of $D_y$ is more concentrated around 0 than that of D). If $R{B}_D\left(0\right|\mathit{y})<1$, there is evidence against $H_0^1$ (as the distribution of $D_y$ is less concentrated around 0 than that of D). Additionally, we compute the strength of evidence $St{r}_D\left(0\right|\mathit{y})={\Pi }_D(RB\left(d\right|\mathit{y})\le RB\left(0\right|\mathit{y})|\mathit{y})$, where ${\Pi }_D(·|\mathit{y})$ is the cumulative distribution function of $D_y$. As ${\pi }_{{}_D}(·|\mathit{y})$ and ${\pi }_{{}_D}(·)$ in (14) have no closed forms, $R{B}_D\left(0\right|\mathit{y})$ and $St{r}_D\left(0\right|\mathit{y})$ need to be approximated. The following Algorithm 1 summarizes the steps required to test $H_0^1$.

Algorithm 1 RB test for $H_0^1$

(i) Generate $({\theta }_1,{\theta }_2,\dots ,{\theta }_k)$ from $\mathrm{Dirichlet}({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_k)$ based on the algorithm of [28] and compute D as defined in (12). (ii) Repeat step (ii) to obtain a sample of $r_1$ values of D. (iii) Generate $({\theta }_1,{\theta }_2,\dots ,{\theta }_k)$ given the observed data $\mathit{y}=(y_1,\dots ,y_k)$ from $\mathrm{Dirichlet}({\alpha }_1+y_1,{\alpha }_2+y_2,\dots ,{\alpha }_k+y_k)$ and compute $D_y$ as defined in (13). (iv) Repeat step (iii) to obtain a sample of $r_2$ values of $D_y$. (v) Compute the relative belief ratio and the strength as follows: (a) Let L be a positive number. Let ${\widehat{F}}_D$ denote the empirical cdf of D based on the prior sample in (3) and for $i=0,\dots ,L,$ let ${\widehat{d}}_{i/L}$ be the estimate of $d_{i/L},$ the $(i/L)$-the prior quantile of D. Here ${\widehat{d}}_0=0$, and ${\widehat{d}}_1$ is the largest value of D. Let ${\widehat{F}}_D(·|\mathit{y})$ denote the empirical cdf based on $D_y$. For $d\in [{\widehat{d}}_{i/L},{\widehat{d}}_{(i+1)/L})$, estimate $R{B}_D\left(d\right|\mathit{y})={\pi }_D\left(d\right|\mathit{y})/{\pi }_D\left(d\right)$ by $${\widehat{RB}}_D\left(d\right|\mathit{y})=L\{{\widehat{F}}_D\left({\widehat{d}}_{(i+1)/L}\right|\mathit{y})-{\widehat{F}}_D\left({\widehat{d}}_{i/L}\right|\mathit{y})\},$$ (15) the ratio of the estimates of the posterior and prior contents of $[{\widehat{d}}_{i/L},{\widehat{d}}_{(i+1)/L}).$ Thus, we estimate $R{B}_D\left(0\right|\mathit{y})={\pi }_D\left(0\right|\mathit{y})/{\pi }_D\left(0\right)$ by ${\widehat{RB}}_D\left(0\right|\mathit{y})=$$L{\widehat{F}}_D\left({\widehat{d}}_{p_0}\right|\mathit{y})$ where $p_0=i_0/L$ and $i_0$ are chosen so that $i_0/L$ is not too small (typically $i_0/L\approx 0.05)$. (b) Estimate the strength ${\Pi }_D(R{B}_D\left(d\right|\mathit{y})\le R{B}_D\left(0\right|\mathit{y})|\mathit{y})$ by the finite sum $$\sum _{\{i\ge i_0:{\widehat{RB}}_D\left({\widehat{d}}_{i/L}\right|\mathit{y})\le {\widehat{RB}}_D\left(0\right|\mathit{y})\}}({\widehat{F}}_D\left({\widehat{d}}_{(i+1)/L}\right|\mathit{y})-{\widehat{F}}_D\left({\widehat{d}}_{i/L}\right|\mathit{y})).$$ (16) For fixed $L,$ as $r_1\to \infty ,r_2\to \infty ,$ then ${\widehat{d}}_{i/L}$ converges almost surely to $d_{i/L}$ and (15) and (16) converge almost surely to $R{B}_D\left(d\right|\mathit{y})$ and ${\Pi }_D(R{B}_D\left(d\right|\mathit{y})\le R{B}_D\left(0\right|\mathit{y})|\mathit{y})$, respectively. See [34] for the details.

3.2. Bayesian r-Sample Multinomial Test with a Completely Specified Null Hypothesis

Consider r independent samples ${\mathit{Y}}_1,{\mathit{Y}}_2,\dots ,{\mathit{Y}}_r$ where each ${\mathit{Y}}_i=(Y_{1|i},\dots ,Y_{k|i})$ follows a multinomial distribution with parameters $n_i$ and ${\theta }_i=({\theta }_{1|i},\dots ,{\theta }_{k|i})$, where ${\sum }_{j=1}^k{\theta }_{j|i}=1$ for $i=1,\dots ,r$. Here, ${\theta }_{j|i}$ denotes the probability of an outcome falling in category j for the ith sample, and $Y_{j|i}$ represents the number of times the outcome falls in category j in the ith sample. The null hypothesis to be tested is $H_0^2:{\theta }_{j|i}={\theta }_{j0|i}$ for $j=1,2,\dots ,k$, where ${\theta }_{j0|i}$ are known constants.

Let the joint distribution of ${\mathit{Y}}_1,{\mathit{Y}}_2,\dots ,{\mathit{Y}}_r$ be

$$\begin{array}{ccc}\hfill p({\mathit{y}}_1,{\mathit{y}}_2,\dots ,{\mathit{y}}_r)& =& \prod _{i=1}^{r}p\left({\mathit{y}}_i\right)=\prod _{i=1}^r\left\{\right(\genfrac{}{}{0pt}{}{n_i}{y_{1|i},y_{2|i},\dots ,y_{k|i}}\left)\prod _{j=1}^k{\theta }_{j|i}^{y_{j|i}}\right\}.\hfill \end{array}$$

(17)

The proposed test is based on measuring the KL divergence between p and

$$\begin{array}{ccc}\hfill q({\mathit{y}}_1,{\mathit{y}}_2,\dots ,{\mathit{y}}_r)& =& \prod _{i=1}^{r}q\left({\mathit{y}}_i\right)=\prod _{i=1}^r\left\{\right(\genfrac{}{}{0pt}{}{n_i}{y_{1|i},y_{2|i},\dots ,y_{k|i}}\left)\prod _{j=1}^k{\theta }_{j0|i}^{y_{j0|i}}\right\}.\hfill \end{array}$$

(18)

The following proposition provides the expression for the KL divergence between p and q. The proof follows directly from Lemma 1 and Proposition 1.

Proposition 2. 

Let $p({\mathit{y}}_1,{\mathit{y}}_2,\dots ,{\mathit{y}}_r)$ and $q({\mathit{y}}_1,{\mathit{y}}_2,\dots ,{\mathit{y}}_r)$ be the two joint probability mass functions as defined in (17) and (18), respectively. Then

$$\begin{array}{ccc}\hfill d(p,q)& =& \sum _{i=1}^r\left\{n_i\sum _{j=1}^k\right[{\theta }_{j|i}\log \left(\frac{{\theta }_{j|i}}{{\theta }_{j0|i}}\right)\left]\right\}.\hfill \end{array}$$

Proposition 2 provides a connection between the KL divergence formula and the test statistic ${\chi }^2$ in (6). Using a Taylor series expansion, we can approximate the distance $d(p,q)$ as follows:

$$\begin{array}{ccc}\hfill d(p,q)& \approx & 0.5\sum _{i=1}^r\sum _{j=1}^k\frac{{(n_i{\theta }_{j|i}-n_i{\theta }_{j0|i})}^2}{n_i{\theta }_{j0|i}}=0.5\sum _{i=1}^r\sum _{j=1}^k\frac{{(E\left(Y_{j|i}\right)-n_i{\theta }_{j0|i})}^2}{n_i{\theta }_{j0|i}},\hfill \end{array}$$

which reveals a close connection to ${\chi }^2$.

For the proposed Bayesian test of $H_0^2$, we adopt the prior $({\theta }_{1|i},\dots ,{\theta }_{k|i})\sim \mathrm{Dirichlet}({\alpha }_{1|i},{\alpha }_{2|i},\dots ,$${\alpha }_{k|i})$, and use the algorithm developed in [28] to elicit the hyperparameters ${\alpha }_{j|i}$ for $j=1,\dots ,k$ and $i=1,\dots ,r$. In this case, we define the divergence measure as

$$\begin{array}{c}\hfill D=\sum _{i=1}^r\left\{n_i\sum _{j=1}^k\right[{\theta }_{j|i}\log \left(\frac{{\theta }_{j|i}}{{\theta }_{j0|i}}\right)\left]\right\},\end{array}$$

(19)

where ${\theta }_{j0|i}$ is the hypothesized value of ${\theta }_{j|i}$ under the null hypothesis.

The posterior distribution of $({\theta }_{1|i},\dots ,{\theta }_{k|i})$ given the observed data ${\mathit{y}}_i=(y_{1|i},\dots ,y_{k|i})$ is then $\mathrm{Dirichlet}({\alpha }_{1|i}+y_{1|i},{\alpha }_{2|i}+y_{2|i},\dots ,{\alpha }_{k|i}+y_{k|i})$. In this case, the divergence measure becomes

$$\begin{array}{c}\hfill D_y=\sum _{i=1}^r\left\{n_i\sum _{j=1}^k\right[{\theta }_{j|i}\log \left(\frac{{\theta }_{j|i}}{{\theta }_{j0|i}}\right)\left]\right\}.\end{array}$$

(20)

The following algorithm outlines the steps required to test $H_0^2$ using the proposed Bayesian test (Algorithm 2):

Algorithm 2 RB test for $H_0^2$

(i) For $i=1,\dots ,r$, generate $({\theta }_{1|i},{\theta }_{2|i},\dots ,{\theta }_{k|i})$ from $\mathrm{Dirichlet}({\alpha }_{1|i},{\alpha }_{2|i},\dots ,{\alpha }_{k|i})$ based on the algorithm of [28] and compute D as defined in (19). (ii) Repeat step (i) to obtain a sample of $r_1$ values of D. (iii) For $i=1,\dots ,r$, generate $({\theta }_{1|i},{\theta }_{2|i},\dots ,{\theta }_{k|i})$ given the observed data ${\mathit{y}}_i=(y_{1|i},\dots ,y_{k|i})$ from $\mathrm{Dirichlet}({\alpha }_{1|i}+y_{1|i},{\alpha }_{2|i}+y_{2|i},\dots ,{\alpha }_{k|i}+y_{k|i})$ and compute $D_y$ as defined in (20). (iv) Repeat step (iii) to obtain a sample of $r_2$ values of $D_y$. (v) Compute the relative belief ratio and strength as described in Algorithm 1.

3.3. Bayesian Test for Homogeneity in r-Sample Multinomial Data

Consider r independent samples ${\mathit{Y}}_1,{\mathit{Y}}_2,\dots ,{\mathit{Y}}_r$, where for $i=1,\dots ,r$, ${\mathit{Y}}_i=(Y_{1|i},\dots ,Y_{k|i})$$\sim \mathrm{multinomial}(n_i,{\theta }_{1|i},\dots ,{\theta }_{k|i})$ with ${\sum }_{j=1}^k{\theta }_{j|i}=1$. To test the null hypothesis $H_0^3$ as defined in (5), it is required to measure the KL divergence between p and q as defined in (17) and (18) with ${\theta }_{j0|i}$ is replaced by ${\theta }_j$. This requirement is offered in the following proposition.

Proposition 3. 

Consider the probability mass functions p and q as defined in (17) and (18) with ${\theta }_{j0|i}$ is replaced by ${\theta }_j$. Then

$$\begin{array}{ccc}\hfill d(p,q)& =& \sum _{i=1}^r\left\{n_i\sum _{j=1}^k\right[{\theta }_{j|i}\log \left(\frac{{\theta }_{j|i}}{{\theta }_j}\right)\left]\right\}\hfill \end{array}$$

(21)

and

$$\begin{array}{c}\hfill {\theta }_j^{★}=\arg \underset{{\theta }_j}{\min }d(p,q)=\frac{{\sum }_{i=1}^r{n}_i{\theta }_{j|i}}{{\sum }_{i=1}^r{n}_i}=\frac{{\sum }_{i=1}^r{n}_i{\theta }_{j|i}}{n}.\end{array}$$

(22)

Proof. 

(21) follows directly from Lemma 2 by setting ${\theta }_{j|i}={\theta }_j$. To prove (22), we use we use Lagrange multiplier with the constraint ${\sum }_{j=1}^r{\theta }_j=1$:

$$\begin{array}{c}\hfill L=L({\theta }_j,\lambda )=\sum _{i=1}^r\left\{n_i\sum _{j=1}^k\right[{\theta }_{j|i}\log \left(\frac{{\theta }_{j|i}}{{\theta }_j}\right)\left]\right\}+\lambda (\sum _{j=1}^r{\theta }_j-1).\end{array}$$

Now,

$$\begin{array}{c}\hfill \frac{\partial L}{\partial {\theta }_j}=-\sum _{i=1}^r{n}_i\frac{{\theta }_{j|i}}{{\theta }_j}+\lambda ,j=1,\dots ,k.\end{array}$$

Setting $\frac{\partial L}{\partial {\theta }_j}=0$ gives

$$\begin{array}{c}\hfill {\theta }_j=\frac{{\sum }_{i=1}^r{n}_i{\theta }_{j|i}}{\lambda },j=1,\dots ,k.\end{array}$$

Summing over both sides and applying the constraint gives

$$\lambda =\sum _{j=1}^k\sum _{i=1}^r{n}_i{\theta }_{j|i}=\sum _{i=1}^r{n}_i\sum _{j=1}^k{\theta }_{j|i}=\sum _{i=1}^r{n}_i=n.$$

Hence,

$$\begin{array}{c}\hfill {\theta }_j^{★}=\frac{{\sum }_{i=1}^r{n}_i{\theta }_{j|i}}{n}.\end{array}$$

□

Note that ${\theta }_j^{★}$ represents the weighted average of ${\theta }_{j|i}$. Substituting ${\theta }_j^{★}$ into (21), we get

$$\begin{array}{ccc}\hfill \widehat{d}(p,q)& =& \sum _{i=1}^r\left\{n_i\sum _{j=1}^k\right[{\theta }_{j|i}\log \left(\frac{n{\theta }_{j|i}}{{\sum }_{i=1}^r{n}_i{\theta }_{j|i}}\right)\left]\right\},\hfill \end{array}$$

(23)

which is equal to 0 under $H_0^3$. We can also establish a connection between (23) and the test statistic ${\chi }^2$ in (6). By Taylor series expansion, when n is large and under $H_0^3$, we have

$$\begin{array}{ccc}\hfill \widehat{d}(p,q)& \approx & 0.5\sum _{i=1}^r{n}_i\sum _{j=1}^k\frac{{({\theta }_{j|i}-\frac{{\sum }_{i=1}^r{n}_i{\theta }_{j|i}}{n})}^2}{\frac{{\sum }_{i=1}^r{n}_i{\theta }_{j|i}}{n}}\hfill \\ & =& 0.5\sum _{i=1}^r\sum _{j=1}^k\frac{{(n_i{\theta }_{j|i}-n_i\frac{{\sum }_{i=1}^r{n}_i{\theta }_{j|i}}{n})}^2}{n_i\frac{{\sum }_{i=1}^r{n}_i{\theta }_{j|i}}{n}}\hfill \\ & =& 0.5\sum _{i=1}^r\sum _{j=1}^k\frac{{(E\left(Y_{j|i}\right)-n_i\frac{{\sum }_{i=1}^{r}E\left(Y_{j|i}\right)}{n})}^2}{n_i\frac{{\sum }_{i=1}^{r}E\left(Y_{j|i}\right)}{n}},\hfill \end{array}$$

which is closely linked to ${\chi }^2$. The proposed Bayesian test for $H_0^3$ uses the prior described in Section 3.1. We write

$$\begin{array}{c}\hfill \widehat{D}=\sum _{i=1}^r\left\{n_i\sum _{j=1}^k\right[{\theta }_{j|i}\log \left(\frac{n{\theta }_{j|i}}{{\sum }_{i=1}^r{n}_i{\theta }_{j|i}}\right)\left]\right\}.\end{array}$$

(24)

Moreover, for the posterior distribution of $({\theta }_{1|i},\dots ,{\theta }_{k|i})$ given the observed data ${\mathit{y}}_i=(y_{1|i},\dots ,y_{k|i})$, we write

$$\begin{array}{c}\hfill {\widehat{D}}_y=\sum _{i=1}^r\left\{n_i\sum _{j=1}^k\right[{\theta }_{j|i}\log \left(\frac{n{\theta }_{j|i}}{{\sum }_{i=1}^r{n}_i{\theta }_{j|i}}\right)\left]\right\}.\end{array}$$

(25)

The following algorithm is used to test $H_0^3$ (Algorithm 3).

Algorithm 3 RB test for $H_0^3$

(i) For $i=1,\dots ,r$, generate $({\theta }_{1|i},{\theta }_{2|i},\dots ,{\theta }_{k|i})$ from $\mathrm{Dirichlet}({\alpha }_{1|i},{\alpha }_{2|i},\dots ,{\alpha }_{k|i})$ based on the algorithm of [28] and compute $\widehat{D}$ as defined in (24). (ii) Repeat step (i) to obtain a sample of $r_1$ values of $\widehat{D}$. (iii) For $i=1,\dots ,r$, generate $({\theta }_{1|i},{\theta }_{2|i},\dots ,{\theta }_{k|i})$ given the observed data ${\mathit{y}}_i=(y_{1|i},\dots ,y_{k|i})$ from $\mathrm{Dirichlet}({\alpha }_{1|i}+y_{1|i},{\alpha }_{2|i}+y_{2|i},\dots ,{\alpha }_{k|i}+y_{k|i})$ and compute ${\widehat{D}}_y$ as defined in (25). (iv) Repeat step (iii) to obtain a sample of $r_2$ values of ${\widehat{D}}_y$. (v) Compute the relative belief ratio and strength using Algorithm 1, but replace D and $D_y$ with $\widehat{D}$ and ${\widehat{D}}_y$, respectively.

4. Examples

This section presents three examples that demonstrate the effectiveness of our approach in evaluating $H_0^1,H_0^2,$ and $H_0^3$. We use Algorithms 1–3, with fixed values of $L=20$, $i_0=1$, and $r_1=r_2={10}^4$. To further investigate the efficacy of our approach, we consider three different prior distributions: uniform prior, Jeffreys prior, and an elicited prior based [28]. Additionally, we compute the p-values using the test statistics discussed in Section 1 of this paper. The approach was implemented using R (version 4.2.1), and the code is available upon request from the corresponding author.

Example 1 

(Rolling Die; [5]). We roll a die 60 times and seek to test whether it is unbiased, that is, whether $H_0^1:{\theta }_j=1/6$ for $j=1,\dots ,k$. The

Table 1. Data of Example 1.

	1	2	3	4	5	6	Total
Observed	8	11	5	12	15	6	60

below presents the recorded data:

We will use a Bayesian approach to address this problem. We employ three priors: the uniform prior represented by $\mathrm{Dirichlet}(1,1,1,1,1,1)$, Jeffreys prior represented by $\mathrm{Dirichlet}(0.5,0.5,$$0.5,0.5,0.5,0.5)$, and the elicited prior Dirichlet (5.83, 5.83, 5.83, 5.83, 5.83, 5.83) obtained using the algorithm proposed by [28], with a lower bound of 0.05 applied to all probabilities. It is worth noting that setting the lower bound in [28] to 0 yields the uniform prior. Additionally, we will include the p-value for the corresponding frequentist test as a reference. The results of our analysis are presented in

Table 2. The RB and its strength (Str) for Example 1.

Prior	RB (Strength)	Decision
Uniform	15.125 (1)	Strong evidence in favor of $H_0^1$
Jeffreys	19.824 (1)	Strong evidence in favor of $H_0^1$
Evan et al.	1.900 (1)	Strong evidence in favor of $H_0^1$
p-value	0.3027	Fail to reject $H_0^1$ at $\alpha =0.05$

. Clearly, both the proposed Bayesian approach, considering the three priors, and the frequentist approach lead to the same conclusion. It should be noted that the uniform prior and the Jefferey prior have a wider spread around zero compared to the elicited prior. As a result, they have higher relative belief ratios in this example. However, this is not practically significant in our case as we calibrate the relative belief ratio through the strength. See Figure 1.

Example 2 

(Operation Trial [5]). In a system consisting of four independent components, let ${\theta }_{j|i}$ denote the probability of successful operation of the ith component, $i=1,2,3,4$. We will test the null hypothesis $H_0^2:{\theta }_{1|1}=0.9,{\theta }_{2|1}=0.1,{\theta }_{1|2}=0.9,{\theta }_{2|2}=0.1,{\theta }_{1|3}=0.8$,${\theta }_{2|3}=0.2$, ${\theta }_{1|4}=0.8$, ${\theta }_{2|4}=0.2$, given that in 50 trials, the components operated as follows (

Table 3. Data of Example 2.

Component	Successful	Failure
1	40	10
2	48	2
3	45	5
4	40	10

):

We use the priors: $\mathrm{Dirichlet}(1,1)$, $\mathrm{Dirichlet}(0.5,0.5)$, and $\mathrm{Dirichlet}(33.38,5.62)$. We obtain the latter prior using algorithm of [28], with lower bounds of ${\theta }_{1|i}=0.7$ and ${\theta }_{2|i}=0.1$ for all $i=1,2,3,4$.

Table 4. The RB and its strength (Str) for Example 2.

Prior	RB (Strength)	Decision
Uniform	20(1)	Strong evidence in favor of $H_0^2$
Jeffreys	19.998 (0.000)	Weak evidence in favor of $H_0^2$
Evan et al.	0.592 (0.047)	Strong evidence against $H_0^2$
p-value	0.030	Fail to reject $H_0^1$ at $\alpha =0.05$

displays the results of our analysis. As in Example 1, both the uniform prior and the Jefferey prior exhibit less concentration around zero when compared to the elicited prior. This, in turn, leads to a notably different conclusion than that of the elicited prior and the p-value calculated using the chi-square test. See also Figure 2.

Example 3 

(Clinical Trial; [37]). A study was performed to determine whether the type of cancer differed between blue-collar, white-collar, and unemployed workers. A sample of 100 of each type of worker diagnosed as having cancer was categorized into one of three types of cancer. The results are shown in

Table 5. Data of Example 3.

Occupation	Type of Cancer
	Lung	Stomach	Other	Total
Blue collar	53	17	30	100
White Collar	10	67	23	100
Unemployed	30	30	40	100

. See also Table 12.6 of [37]. The hypothesis to be tested is that the proportions of the three cancer types are the same for all three occupation groups. That is, $H_0^3:{\theta }_{j|1}={\theta }_{j|2}={\theta }_{j|3}$ for all j (types of cancer), where ${\theta }_{j|i}$ is the probability of occupation i having cancer type j.

Similar to the previous two examples, we utilize the uniform prior $\mathrm{Dirichlet}(1,1,1)$, Jeffreys prior $\mathrm{Dirichlet}(0.5,0.5,0.5)$, and the elicited prior $\mathrm{Dirichlet}(3,3,3)$. We obtained the elicited prior using the algorithm of [28] by setting a lower bound of $0.05$ for all probabilities.

Table 6. The RB and its strength (Str) for Example 3.

Prior	RB (Strength)	Decision
Uniform	0.102 (0.010)	Strong evidence against $H_0$
Jeffreys	1.874 (0.129)	Weak evidence in favor of $H_0$
Evan et al.	0.012 (0.000)	Strong evidence against $H_0$
p-value	0.000	Reject $H_0^1$ at $\alpha =0.05$

summarizes the results of our analysis. Similar to the previous examples, Jeffreys prior is not sufficiently concentrated around zero, which makes it inefficient when there is evidence against $H_0$. See Figure 3.

5. Concluding Remarks

This study presents a Bayesian method for testing hypotheses related to multinomial distributions. Our approach involves calculating the Kullback–Leibler divergence between two multinomial distributions and comparing the change in distance from the prior to the posterior through the relative belief ratio. To specify the prior distributions, we employ a prior elicitation algorithm. We recommend avoiding the use of Jeffreys prior or the uniform prior unless there is a valid reason to use them. Through several examples, we demonstrated the effectiveness of our approach. Future research may expand our approach to include testing for independence and other related cases.