## 6 Bayesian Inference on \(R\)  
In this section, the Bayes estimate of \(R\) can be obtained under assumption that the shape parameter \(\alpha\) and scale parameters \(\theta_{1}\) and \(\theta_{2}\) are random variables. According to the likelihood function in Section 3, we note that \(\alpha\) has a positive exponent and \(\theta_{1}\) and \(\theta_{2}\) have negative exponents, respectively, so we can assume that \(\theta_{1}\) and \(\theta_{2}\) have independent Inverse Gamma pdf and \(\alpha\) follows Gamma distribution. We choose this family such that prior-to-posterior updating yields a posterior that is also in the family:  
\[\pi(\theta_{1})=\frac{b_{1}^{a_{1}}}{\gamma(a_{1})}\theta_{1}^{-(1+a_{1})}e^{ -\frac{b_{1}}{b_{1}}},\ \theta_{1}>0, \tag{31}\]  
\[\pi(\theta_{2})=\frac{b_{2}^{a_{2}}}{\gamma(a_{2})}\theta_{2}^{-(1+a_{2})}e^{ -\frac{b_{2}}{b_{2}}},\ \theta_{2}>0, \tag{32}\]  
where all the hyper-parameters \(a_{i}\) and \(b_{i}\ (i=1,2)\) are assumed to be known and non-negative. The prior density function of \(\alpha\) is denoted as \(\pi(\alpha)\), and we assume that it has a Gamma \((0,1)\) distribution. We have the likelihood function based on the above assumptions as  
\[L(data|\alpha,\theta_{1},\theta_{2})=\alpha^{m+n}\theta_{1}^{-n}\theta_{2}^{-m}( \prod_{i=1}^{n}x_{i}^{-\alpha-1})(\prod_{j=1}^{m}y_{j}^{-\alpha-1})e^{-\frac{1} {\theta_{1}}\sum_{i=1}^{m}x_{i}^{-\alpha}}e^{-\frac{1}{\theta_{2}}\sum_{j=1}^{m }y_{j}^{-\alpha}}. \tag{33}\]  
The joint density of the data, \(\alpha\), \(\theta_{1}\) and \(\theta_{2}\) becomes  
\[P(data,\alpha,\theta_{1},\theta_{2})=L(data|\alpha,\theta_{1},\theta_{2})\pi( \alpha)\pi(\theta_{1})\pi(\theta_{2}). \tag{34}\]  
Therefore, the joint posterior density of \(\alpha\), \(\theta_{1}\) and \(\theta_{2}\) given the data is  
\[P(\alpha,\theta_{1},\theta_{2}|data)=\frac{P(data,\alpha,\theta_{1},\theta_{2 })}{\int_{0}^{\infty}\int_{0}^{\infty}L(data,\alpha,\theta_{1},\theta_{2}) \mathrm{d}\alpha\mathrm{d}\theta_{1}\mathrm{d}\theta_{2}}. \tag{35}\]  
Since the expression (35) cannot be written in a closed form, the Bayes estimate of \(R\) and the corresponding credible interval of \(R\) are derived adopting the Gibbs sampling technique:  
\[P(\alpha,\theta_{1},\theta_{2}|data)\propto \alpha^{m+n-1}\theta_{1}^{-(n+1+a_{1})}\theta_{2}^{-(m+1+a_{2})}( \prod_{i=1}^{m}x_{i}^{-\alpha-1})(\prod_{j=1}^{m}y_{j}^{-\alpha-1})\] \[\exp\left\{-\frac{1}{\theta_{1}}(\sum_{i=1}^{n}x_{i}^{-\alpha}+b_ {1})-\frac{1}{\theta_{2}}(\sum_{j=1}^{m}y_{j}^{-\alpha}+b_{2})-\alpha\right\}.\]  
The posterior pdfs of \(\alpha\), \(\theta_{1}\) and \(\theta_{2}\) can be obtained based on the expression \(P(\alpha,\theta_{1},\theta_{2}|data)\) as the following:  
\[\theta_{1}|\alpha,\theta_{2},data\sim IG(n+a_{1},b_{1}+\sum_{i=1}^{n}x_{i}^{- \alpha}),\]  
\[\theta_{2}|\alpha,\theta_{1},data\sim IG(m+a_{2},b_{2}+\sum_{j=1}^{m}y_{j}^{- \alpha}),\]  
\[f_{\alpha}(\alpha|\theta_{1},\theta_{2},data)\propto\alpha^{m+n-1}\prod_{i=1} ^{n}x_{i}^{-\alpha-1}\prod_{j=1}^{m}y_{j}^{-\alpha-1}\exp\left\{-\frac{1}{ \theta_{1}}(\sum_{i=1}^{n}x_{i}^{-\alpha}+b_{1})-\frac{1}{\theta_{2}}(\sum_{j =1}^{m}y_{j}^{-\alpha}+b_{2})-\alpha\right\}.\]  
The posterior pdf of \(\alpha\) is not known. We use the Metropolis-Hastings method with normal proposal distribution to generate random numbers from this distribution. The algorithm of Gibbs sampling is described as follows:  
Step 1: Start with an initial guess \((\alpha^{(0)},\theta_{1}^{(0)},\theta_{2}^{(0)})\). Step 2: Set \(t=1\). Step 3: Generate \(\theta_{1}^{(t)}\) form \(IG(n+a_{1},b_{1}+\sum_{i=1}^{n}x_{i}^{-\alpha(t-1)})\). Step 4: Generate \(\theta_{2}^{(t)}\) form \(IG(m+a_{2},b_{2}+\sum_{j=1}^{m}y_{j}^{-\alpha(t-1)})\). Step 5: Using the Metropolis-Hastings method, generate \(\alpha^{(t)}\) from \(f_{\alpha}(\alpha^{(t-1)}|\theta_{1}^{(t)},\theta_{2}^{(t)},data)\). * Generate a new candidate value \(\delta\) from \(N(\ln\alpha^{(t-1)},1)\). * Set \(\alpha^{\prime}=\exp(\delta)\).