EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
Vol. 17, No. 1, 2024, 180-200
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Bayesian regression analysis using median rank set
sampling

Inad Nawajah1,∗, Hassan Kanj2,, Yehia Kotb2,, Julian Hoxha2,, Mouhammad
Alakkoumi2,, Kamel Jebreen3,4,5,

1 Department of Mathematics, College of Science and Technology, Hebron University, Pales-
tine
2 College of Engineering and Technology, American University of the Middle East, Egaila
54200, Kuwait
3 Department of Mathematics, Palestine Technical University Kadoorie, Hebron, Palestine
4 Department of Mathematics, An-Najah National University, Nablus, Palestine
5 Unité de Recherche Clinique Saint-Louis Fernand-Widal Lariboisière, APHP, Paris, France

Abstract. Bayesian estimation of the linear regression parameter system is considered by de-
ploying Median Rank Set Sampling (MRSS). The full conditional distributions and the associated
posterior distribution are obtained. Therefore, based on Markov Chain Monte Carlo simulation,
the Bayesian point estimates and credible intervals for the regression parameters are determined.
To measure the efficiency of the obtained Bayesian estimates concerning the frequentist estimates
we compute the asymptotic relative efficiency of the obtained Bayesian estimates using Markov
Chain Monte Carlo simulation.
This study shows that the Bayesian estimation of the simple linear regression parameters under
frequentist MRSS is highly beneficial and much superior to the RSS scheme.

2020 Mathematics Subject Classifications: 62C10, 62J20, 62D05

Key Words and Phrases: Median Ranked Set Sampling; Bayes factor; Regression; Bayesian
approach

1. Introduction

Regression is the process of investigating the relationship between a dependent variable
and some independent variables given a set of historical data. The objective of regression
is to discover a model that accurately captures this relationship and minimizes the dis-
crepancy between predicted values from the model and the actual values observed in the
data.

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v17i1.5015

Email addresses: inadn@hebron.edu (I. Nawaja), hassan.kanj@aum.edu.kw (H. Kanj),
Yehia.Kotb@aum.edu.kw (Y. Kotb), Julian.Hoxha@aum.edu.kw (J. Hoxha), Mouhammad.A@aum.edu.kw
(M. Alakkoumi) K.Jebreen@yahoo.com (K. Jebreen)

https://www.ejpam.com 180 © 2024 EJPAM All rights reserved.



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 181

Regression encompasses various techniques, and one common approach is linear regression,
which takes the form of:

Y = β0 + β1X1 + β2X2 + . . . + βnXn + ϵ (1)

In the Equation 1, the dependent variable Y is related to the independent variables Xi,
where 1 ≤ i ≤ n. The coefficients βi represent the respective coefficients for each inde-
pendent variable, and the term ϵ represents the error. This equation serves as a general
representation of linear regression and can be modified to accommodate polynomial re-
gression, as we will explore later in this section.
As previously mentioned, linear regression can be categorized into different types, and
one such type is simple linear regression. In simple linear regression, the objective is
to establish a model that relates a single dependent variable Y to a single independent
variable X. By substituting n = 1 into Equation 1, we obtain the equation specifically for
simple linear regression:

Y = β0 + β1X+ ϵ (2)

Note that X1 was referred to as X since it is the only independent variable in the data set.
The error in this context is defined as the Euclidean distance between each data point

and the model. The objective of simple linear regression is to determine the optimal values
of β0 and β1 that minimize this error.
Simple linear regression offers several advantages and disadvantages. Some of the ad-
vantages include simplicity, high performance, and interpretability. The simplicity arises
from the model’s relationship between the dependent variable Y and a single indepen-
dent variable X. This simplicity contributes to its high performance, making it suitable
for real-time systems. Additionally, the model is easily interpretable, with the intercept
(β0) representing the expected value of Y when X is zero and the slope coefficient (β1)
indicating the change in Y for each unit increase in X.
Simple linear regression also has some disadvantages. One such disadvantage is the lin-
earity assumption, which assumes that the relationship between the variables is strictly
linear. In reality, observations often exhibit non-linear patterns, making this assumption
restrictive and potentially leading to higher mean square error. Another drawback is the
sensitivity of simple linear regression to outliers. Outliers can disproportionately influence
the model’s estimates, leading to distorted results and increased mean square error.
As mentioned earlier, while Equation 1 represents the multi-variable regression, it can also
be considered as a general equation of regression. As seen before, Equation 2 is a special
case of Equation 1.
Similar to Simple linear regression, polynomial regression is a regression model that estab-
lishes a relationship between a single independent variable X and the dependent variable
Y. In contrast to simple linear regression, polynomial regression offers improved accuracy
by incorporating higher-degree polynomials. Additionally, polynomial regression can be
viewed as a special case of linear regression. The Equation 3 defines polynomial regression
as follows:

Y = β0 + β1X+ β2X2 + . . .+ Xn + ϵ (3)



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 182

Since data is being collected from real life scenarios, there is a possibility of error or
uncertainly about the validity of the gathered data. That is when Bayesian statistics
come into place. Bayesian statistics is a powerful technique that performs predictions
based on history of data. It combines the uncertainty of data with the model which allows
predictions to be realistic and match the real world event. The equation of Bayes rule is
as follows:

P(A|B) =
P(B|A)× P(A)

P(B)
(4)

Equation 4 deals with causalities. The equation evaluates the probability for event A to
happen given that event B has already occurred and been observed. This probability
depends on history. History here is seen as prior knowledge represented by P(A) that is
the the probability of the occurrence of event A. The prior probability is also known to
be the initial belief that even A would occur.
The marginal likelihood P(B), also known as evidence, is the overall probability for event
B to occur. This is particularly important since Bayes studies the causality relationship
between B and A, that is the frequency of B causing A.
The Likelihood of event A causing event B to occur is represented by P(B|A). Finaly the
posterior probability or the most updated probability of event B causing the occurrence
of A is given by P(A|B).
The problem with regression is the fact that it deals with continuous data that could
literally be infinite when observed from a continuously active environment. When data
exceeds certain size, the model could be vulnerable to over-fitting. We use sampling
to make sure regression is performed smoothly and avoid any possible data problems.
Sampling is much more efficient than data sample selection from a given population.
Sampling also ensures data availability since collecting data for entire population might
no be feasible specially for application that involve big data [5].
One of the sampling techniques is the Rank Set Sampling (RSS for short) [22]. RSS is very
useful when data sizes are huge or if the evaluation or measurement of the independent
variables is computationally expensive. Real time system that uses Regression analysis
as an example can make great benefit of RSS sampling since system tasks have deadlines
to meet. Rank Set Sampling performs grouping on observations in the available popula-
tion. Every group has one or more observations. Groups are then ranked and chosen for
sampling.
The structure of this article is as follows: Section 1 provides an introduction to the topic.
Section 2 presents a review of prior research in the field. Section 3 outlines the proposed
model described in this paper. Section 4 introduces the Bayes Estimator and the model
selector proposed in this study. Section 5 discusses the Bayes factor of estimators utilized
in this article. Section 6 describes the simulation methodology and presents the corre-
sponding results. Finally, Section 8 offers a discussion of the findings and concludes the
article [16, 19].



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 183

2. Previous Work

The introduction of Ranked Set Sampling (RSS), a relatively new sampling technique for
estimating population mean, was originally presented in a paper by McIntyre [22]. RSS
involves the following steps: Firstly, a visual inspection is employed to randomly select m2

sample units from the target population. Secondly, these selected units are allocated into
m sets as randomly as possible, with each set comprising m units. Thirdly, an RSS of size
m is constructed for analysis by sequentially selecting the smallest ranked unit from the
first set, the second smallest ranked unit from the second set, and so on, until the largest
ranked unit is chosen from the last set. This process can be repeated r times to obtain a
desired sample size of n = rm [26].
Several authors have evaluated the effectiveness of Bayesian approaches in the context of
RSS, including Al-Saleh et al. [2], Alodat and Al Sagheer [8], Wolfe [27], Kohlschmidt et
al. [20], and Al-saleh and Al-Shrafa [3]. AlOdat et al [6] initially proposed the application
of Bayesian statistics to RSS, and Al-Hamide et al. [4] further studied Bayesian inference
for the linear regression model, assuming a prior distribution for the regression parameters
following the alpha-skew-normal distribution.
Muttlak [23] introduced Median Ranked Set Sampling (MRSS), which involves quantifying
the median of each set from the aforementioned m sets. Within each group, the units are
ranked after randomly distributing the chosen m2 units into m sets of size m. If m is odd,

the m+1
2

th
smallest rank unit (median) is selected from each set, whereas if m is even, the

m
2
th smallest rank unit is chosen from the first m

2 set. This process is repeated r times
until a sample of size n = rm is obtained. AlOdat et al. [6] analyzed parameter estimators
for a simple linear regression model using the MRSS scheme. Additionally, in his another
article AlOdat et al. [7] discussed the large sample properties of the parameter estimators
for simple linear regression based on the MRSS design.
In a study by Al-Hadhrami et al. [1], it was demonstrated that Bayesian estimation of the
mean of a normal distribution using moving extreme ranked set sampling is more efficient
than the frequentist approach of Simple Random Sampling (SRS). Hassan [13] obtained
the maximum likelihood and Bayesian estimators of shape and scale parameters of the
exponentiated exponential distribution based on SRS and RSS.
Li and Balakrshnan [21] developed the best linear unbiased estimators for parameters of
a simple linear regression model using ordered RSS. Haq et al. [12] investigated the best
linear unbiased estimators based on double RSS and ordered double ranked set sampling
(DRSS) for the simple linear regression model with replicated observations.
Yao et al. [28] recently derived the best linear unbiased estimators for simple linear
regression based on moving extremes RSS, which were found to be more efficient than the
estimators obtained under SRS.
In a study by Sazak and Ozel [25], the modified maximum likelihood parameter estimation
of the regression model using bivariate MRSS was investigated. The obtained estimators
were compared with the least squares estimators based on MRSS, as well as with the
modified maximum likelihood and least squares estimators based on RSS.

When estimating the population mean under symmetrical unimodal distributions, me-



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 184

dian ranked set sampling (MRSS), a variation of RSS, outperforms classical RSS [23].
Numerous studies have been conducted on MRSS and its uses [9, 11]. However, in the
case of asymmetric distributions, MRSS can outperform RSS. Moreover, MRSS can im-
prove estimator efficiency since error terms in regression models generally have a normal
distribution. In order to develop effective estimators of regression model parameters, the
Bayesian regression estimators were examined in this study utilizing a median ranked set
sample.

3. Proposed Model

This article focuses on conducting a Bayesian analysis of the simple linear regression model
using Ranked Set Sampling (RSS). We investigate the estimation of regression parameters
by incorporating a prior distribution. Equation 5 represents the modified simple linear
regression model defined by Equation 2, while taking sampling process into consideration.

Yij = β0 + β1Xij + ϵij (5)

where
1 ≤ i ≤ r and 1 ≤ j ≤ n (6)

where r is the iteration number of the optimization process, and n is the group size. The
group size n = 2m− 1 where m is the number of groups. In each group, we have an odd
number of samples, that increases linearly with the number of groups.
ϵij is the error for group i and data sample j. The normal distribution ℵ(µ, σ2) with mean
µ = 0 and standard deviation of σ is used to calculate the error ϵ. β0 is the interception
and β1 is the slope and both are unknown parameters that need to be found through the
regression process.
Assume that we have a performance of the mean ranked set sampling (MRSS) on the
variable of interest and

Yj(m) = Median(Yj1, . . . ,Yjn) (7)

where Yj(m) is the median of group m in iteration j.
The regression of Y

′s
j(m) on X

′s
j can be written as follows:

Yj(m) = B0 +B1Xj(m) + ϵj(m), (8)

where
ϵj(m) = Median(ϵj1, . . . , ϵjn) (9)

It can be noted that ϵj(m)′s are independently and identically distribution (iid) errors and

have constant variance with the probability density function (pdf):

f(ϵ) =
(2m− 1)!

(m− 1)(m− 1)!
Φ(

ϵ

σ
)m−1×



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 185

Φ(
−ϵ

σ
)m−1Θ(

ϵ

σ
)
1

σ
;

with−∞ < ϵ < ∞ (10)

where Θ is the pdf of a standard normal random variable (i.e ℵ(0, 1)) and Φ a cumulative
standard Normal distribution (CDF).
Utilizing the MRSS scheme, the ordinary least squares estimators of β0 and β1 can be
easily obtained as follows: 

β̂0M = Ȳm − β̂1M x̄,

β̂0M =
∑r

j=1(Yj(m)−Ȳm
)(xj−x̄)

Sxx
,

(11)

where

X̄ =
1

r

r∑
j=1

xj , (12)

Ȳm =
1

r

r∑
j=1

Yj(m), (13)

and

Sxx =
r∑

j=1

(xj − x̄)2 (14)

In this work, the Bayesian estimator of the parameters for the simple linear regression
model is obtained employing MRSS and it is compared to SRS setup.

4. Bayes Estimator and Model selector

In this section, we find the Bayes estimator of each parameter in the model defined by
Equation 5.

4.1. Bayesian Model

Based on the likelihood function and the prior distribution, a Bayesian model is described,
where the likelihood function is the conditional distribution of the response variable, given
all the parameters and covariates. We have the following prior distribution of each pa-
rameter, (β0 ,β1 ,σ2 ):

m1 :


Yj(m) = β0 + β1Xj + ϵj(m), 1 ≤ j ≤ r

β0 ∼ ℵ(a, b2)
β1 ∼ ℵ(δ, g2)
σ2 ∼ IG(α, β)

(15)



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 186

Note that, ℵ(µ, σ2) is the normal distribution with mean µ and variance σ2 and IG(α, β)
is the inverse γ distribution. As far as the prior is concerned, all the hyperparameters
a,b2,δ, g2, α, and β are considered a priori (conditionally) independent. Here, the fixed
effects parameters have weakly informative marginal priors, i.e., in this case, the normal
distribution ℵ is centered at 0, with a large variance. Particularly, we let that: a, δ ∼
ℵ(0, 50),b,g ∼ Uniform(0,10) and α, β ∼ γ(2, 2). Distributions that are flat along the whole
real number line are considered weakly informative priors because they don’t provide any
information.
Let πj(θ), j = 1, 2, 3, be the unknown parameters of the prior distribution and define the
marginal density of the random variable Y, we obtain this equation:

m1y =

∫
f2(y | θ)∂y, (16)

where
θ = (β0, β1, σ

2) (17)

The distribution of Yj(m) is given by:

f(yj(m); θ) =
Cm

σ
Φ(

yj(m) − β0 − β1xj

σ
)m−1

× [1− Φ(
yj(m) − β0 − β1xj

σ
)]m−1

×Θ(
yj(m) − β0 − β1xj

σ
)

(18)

where
θ = (β0, β1, σ

2) (19)

and

Cm =
(2m− 1)!

(m− 1)!2.
(20)

We assume that θ has the following prior distribution:

π(θ) = π(β0, β1, σ
2) = π1(β0)π2(β1)π3(σ

2), (21)

where

π1(β0) =
1√
2πb

Exp{− 1

2b2
(β0 − a)2} (22)

and
−∞ < β0 < ∞ (23)

π2(β1) =
1√
2πb

Exp{− 1

2g2
(β1 − δ)2} (24)



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 187

and
−∞ < β1 < ∞ (25)

π3(σ
2) =

1

Γ(α)(σ2)α+1
Exp{− 1

σ2β
} (26)

and σ2 > 0. Equation 22, 24 and 26 are our prior distributions of the three independent
parameters of the vector θ.
If y

m
(y1(m), y2(m), ...., yr(m), ), then y

m
has the following joint pdf:

fy
m
(y

m
; θ) =

Cr
m

σr

r∏
j=1

[Φ(
yj(m) − β0 − β1xj

σ
)m−1]

× [Φ(
−yj(m) − β0 − β1xj

σ
)m−1

×Θ(
−yj(m) − β0 − β1xj

σ
)].

(27)

The joint pdf is formed as a product of individual probabilities because we have assumed
the observations are independently drawn from the underlying distribution. The term
inside the product notation corresponds to the probability density of each individual ob-
servation yj(m). Each of these is distributed according to a modified normal distribution,
with the parameters of the distribution depending on β0, β1, and xj (value of the predictor

variable for the jth group). The term Φ(
yj(m)−β0−β1xj

σ )m−1 accounts for the cumulative
likelihood of each residual from the model raised to the power of (m − 1). Similarly,

Φ(
−yj(m)−β0−β1xj

σ )m−1 captures the remaining part of cumulative distribution from the

observed value, also raised to the power of (m − 1). Θ(
−yj(m)−β0−β1xj

σ ) is the pdf of the
standardized residuals. This term models the density of the residuals around the model
fit.
The posterior pdf of θ given Ym = y

m
is proportional to the product of the likelihood

function, fYm
(Ym; θ), and the prior distribution π(θ). After some manipulation we can

write the posterior pdf as below:

π(θ | Ym)α fYm
(Ym; θ)π(θ)

α
1

σr
[

r∏
j=1

Ej ]
1

(σ2)α+1
[Exp{− 1

2b2
(β0 − a)2}

×Exp{− 1

2g2
(β1 − δ)2} × Exp{− 1

σ2β
}],

(28)

where

Ej = [Φ(
yj(m) − β0 − β1xj

σ
)Φ×



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 188

(
β0 + β1xj − yj(m)

σ
)]m−1

× θ(
yj(m) − β0 − β1xj

σ
) (29)

The posterior pdf is then compactly represented as the product of three components give
as:

π(β0, β1, σ
2 | y

m
)α

1

(σ2)
r
2
+α+1

exp(
−
∑r

j=1(yj(m) − β0 − β1xj)
2

2σ2
− (β0 − a)2

2b2

− (β1 − δ)

2g2
− 1

σ2β
)x

r∏
j=1

Φ(
yj(m) − β0 − β1xj

σ
)m−1

Φ(
−yj(m) + β0 + β1xj

σ
)m−1. (30)

were 1

(σ2)
r
2+α+1 is a normalization term, the second term is the exponential function which

has four components in the exponent, each representing a component of the prior distri-
butions. Lastly, the third term is the product over all residuals of two standard normal
cumulative distribution functions (CDFs) evaluated at the standardized residuals, each
raised to the power of (m-1).
For facilitating the derivation of posterior distributions we first reformulate the sum of
residuals using some simple algebraic manipulations as:

r∑
j=1

(yj(m) − β0 − β1xj)
2 =

r∑
j=1

(yj(m) − β1xj − ȳm + β1x̄)
2

+ r(ȳ − β1x̄− β0)
2 (31)

then the marginal conditional distribution of β0 given β1, σ
2 and ym is:

π1(β0 | β1, σ2, y
m
)

α× Exp{−r(β0 − ȳ + β1x̄)
2

2σ2
− (β0 − a)2

2b2

×
r∏

j=1

Φ(
yj(m) − β0 − β1xj

σ
)m−1

× Φ(
−yj(m) + β0 + β1xj

σ
)m−1,

αExp{−β2
0

2
(
r

σ2
+

1

b2
)+



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 189

β0(r
β1x̄− ȳm

σ2
+

a

b2
)}

×
r∏

j=1

Φ(
yj(m) − β0 − β1xj

σ
)m−1

× Φ(
−yj(m) + β0 + β1xj

σ
)m−1,

(32)

This marginal posterior distribution represents our updated beliefs about the parameter
β0, given the observed data and the specified values for the other parameters. After some
manipulations we compactly write it as:
π1(β0 | β1, σ2, y

m
)

αExp{−
( r
σ2 + 1

¯
2)

2

× (β0 −
r β1x̄−ȳm

σ2 + a
b2

r
σ2 + 1

b2

)2}

×
r∏

j=1

Φ(
yj(m) − β0 − β1xj

σ
)m−1

× Φ(
−yj(m) + β0 + β1xj

σ
)m−1. (33)

We can easily identify two terms; the first one is the Gaussian likelihood function with
the required adjustment for the prior belief about the parameter β0; second term used to
capture the cumulative distribution of the residuals. Overall, this expression models how
the residuals are distributed around the model fit and combines this with prior information
and likelihood derived from the data to give a posterior belief about β0.



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 190

The conditional distribution for β1 given β0, σ
2 and ym is:

π2(β1 | β0, σ2, y
m
)

α exp

{
−
∑r

j=1(yj(m) − β0 − β1xj)
2

2σ2

}

× exp

{
−(β1 − δ)2

2g2

}
×

r∏
j=1

[
Φ

(
yj(m) − β0 − β1xj

σ

)m−1

× Φ

(−yj(m) + β0 + β1xj

σ

)m−1
]

α exp

{
−2
∑r

j=1 xj(yj(m) − β0)β1

2σ2

}

× exp

{
−
β2
1

∑r
j=1 x

2
j

2σ2

}

× exp

{
− β2

1

2g2
+

δβ1
g2

}
×

r∏
j=1

[
Φ

(
yj(m) − β0 − β1xj

σ

)m−1

× Φ

(−yj(m) + β0 + β1xj

σ

)m−1
]
.

(34)



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 191

By some simplification,

π2(β1 | β0, σ2, y
m
) α exp

{
−β2

2

2

(∑r
j=1 x

2
j

σ2
+

1

g2

)

+

(
β2
1

∑r
j=1 x

2
j

σ2
+

δ

g2

)
β1

}

×
r∏

j=1

[
Φ

(
yj(m) − β0 − β1xj

σ

)m−1

×Φ

(−yj(m) + β0 + β1xj

σ

)m−1
]

α exp

−

∑r
j=1 x

2
j

σ2 + 1
g2

2

+

β1 −
∑r

j=1 xj(yj(m)−β0)

σ2 + δ
g2∑r

j=1 x
2
j

σ2 + 1
g2

2


×
r∏

j=1

[
Φ

(
yj(m) − β0 − β1xj

σ

)m−1

×Φ

(−yj(m) + β0 + β1xj

σ

)m−1
]
.

(35)

and finally we have that:

π2(β1 | β0, σ2, y
m
) α exp

−

∑r
j=1 x

2
j

σ2 + 1
g2

2

+

β1 −
∑r

j=1 xj(yj(m)−β0)

σ2 + σ2δ
g2∑r

j=1 x
2
j +

σ2

g2

2


×
r∏

j=1

[
Φ

(
yj(m) − β0 − β1xj

σ

)m−1

×Φ

(−yj(m) + β0 + β1xj

σ

)m−1
]
.

(36)

The conditional distribution of σ2 given β0, β1 and ym is given in the simplified compact



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 192

form as:

π3(σ
2 | β0, β1, ym) α

1

(σ2)
r
2
+α+1

× exp


−

∑r
j=1(yj(m)−β0−β1xj)

2

2 + 1
β

σ2


×

r∏
j=1

[
Φ

(
yj(m) − β0 − β1xj

σ

)m−1

×Φ

(−yj(m) + β0 + β1xj

σ

)m−1
]
.

(37)

5. Bayes Factor of Estimators in MRSS

The Bayes factor is the ratio of the marginal densities of the two cases. The Bayes
factor is given by:

Bji =
mj(y)

mi(y)
=

∫
fj(y | θ)πj(θ)dθ∫
fi(y | θ)πi(θ)dθ

(38)

Here, this ratio evaluates the modification of the odds of mj(y) against mi(y) due to the
observation and can naturally be compared to 1, although an exact comparison scale can
only be based upon a loss function. And the Bayes factor depends on prior information.
It measures the strength of evidence for a model in a way that takes into account both the
goodness of fit and the complexity of the model. However, it requires the computation of
multi-dimensional integrals, which can be quite difficult in practice, especially for complex
models with many parameters. For this reason, it’s often approximated using techniques
like Markov chain Monte Carlo (MCMC) methods.
We can simplify the Bays factor for the estimator depending on their posterior density as:

β̂0b =

∫∞
−∞

∫∞
−∞

∫∞
−∞

β0

σr [
∏r

j=1Em]× 1
(σ2)α+1×∫∞

−∞
∫∞
−∞

∫∞
−∞[

∏r
j=1Em]dβ0dβ1dσ2

×
[
exp

(
− 1

2b2
(β0 − a)2

)
× exp

(
− 1

2g2
(β1 − δ)2

)
× exp

(
− 1

σ2β

)]
dβ0dβ1dσ

2

(39)

In this particular equation, the expectation of β0 (β̂0b) under the posterior distribution
is being calculated, which serves as the Bayesian point estimate for this parameter. This
estimate takes into account both the likelihood of the data given the parameters and the



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 193

prior beliefs about the parameters.

β̂1b =

∫∞
−∞

∫∞
−∞

∫∞
−∞

β1

σr [
∏r

j=1Em]× 1
(σ2)α+1×∫∞

−∞
∫∞
−∞

∫∞
−∞[

∏r
j=1Em]dβ0dβ1dσ2

×
[
exp

(
− 1

2b2
(β0 − a)2

)
× exp

(
− 1

2g2
(β1 − δ)2

)
× exp

(
− 1

σ2β

)]
dβ0dβ1dσ

2

(40)

σ̂2
b =

∫∞
−∞

∫∞
−∞

∫∞
−∞

σ2

σr [
∏r

j=1Em]× 1
(σ2)α+1×∫∞

−∞
∫∞
−∞

∫∞
−∞[

∏r
j=1Em]dβ0dβ1dσ2

×
[
exp

(
− 1

2b2
(β0 − a)2

)
× exp

(
− 1

2g2
(β1 − δ)2

)
× exp

(
− 1

σ2β

)]
dβ0dβ1dσ

2

(41)

All the above three equations involve a triple integral over all possible values of β0,
β1, and σ2

b . In the Bayesian context, the integral goes from (−∞,+∞) because we’re
integrating over all possible values that the parameters could take. In practice, the values
are constrained to make sense in the context of the problem in hand.
The equations are written in such a way to clearly distinguish three common terms in-
volved; in the numerator, the terms for each equation β0

σr ,
β1

σr and σ2

σr are part of the
function we want to find the expected value of, where r is the number of data points. The
product

∏r
j=1Em represents the likelihood function for the model. The terms involving

exponential functions are the prior distributions of β0, respectively. We multiply these
three parts together and integrate over the parameters to get the expected value of β0, β1,
and σ2

b under the posterior distribution. In the denominator, the integral of the likelihood
over all possible parameter values, which serves as a normalization factor to ensure that
the result is a valid probability distribution. This is often called the evidence or marginal
likelihood in Bayesian statistics. The ratio of the numerator to the denominator gives the
expected value of β0, β1, and σ2

b which is an estimate for this parameter under the pos-
terior distribution. The integrals above are typically not solvable analytically, especially
for complex models. Thus, numerical methods like MCMC are often used to approximate
these integrals.

6. Simulation Study

In this study, the MCMC numerical method is used to approximate the calculation
of the posterior expectation of the parameters β0, β1, and σ2

b , derived in section 5. This
serves as the Bayesian point estimate for this parameters. We run the MCMC simulation
for 200,000 iterations. An iteration here corresponds to one cycle through the algorithm
– proposing new values for the parameters, checking whether these new values are likely



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 194

given the data (using the likelihood function), and deciding whether to accept these new
values or keep the old ones [15]. The first 5,000 iterations are discarded. This is known as
the ”burn-in” period. The idea is that early iterations are based on our initial (potentially
poor) guess, and so we have to give the algorithm some time to converge towards the
true values before we start collecting data. After the burn-in period, we start collecting
data, but not at every iteration. Instead, we only keep the parameter values every 50th
iteration. This is known as ”thinning”. The purpose is to reduce the correlation between
successive samples, which can bias the results. After the burn-in and thinning, you are
left with 5,000 samples. These are the values of the parameters at different points in the
MCMC simulation. Each sample represents a possible set of parameters that could have
generated the data. Finally, we need to check whether the algorithm has converged –
that is, whether it has found the true posterior distribution. This is typically done using
various statistical tests, such as the Geweke diagnostic, which compares the mean of the
early and late portions of the chain, and the Heidelberger-Welch diagnostics, which involve
a stationarity test (to see if the chain has reached equilibrium) and a half-width test (to
see if the chain has run long enough to achieve a desired precision) [17, 24].
Through this process, we have generated many different samples of parameter values. The
collection of these samples approximates the posterior distribution of the parameters and
the triple integral in the formula for the posterior mean of β0 can then be approximated
by taking the average of β0 over these samples.
In order to assess the performance of the simulation and the convergence of the chain, in
Figure 1 we show the trace plots for the three parameters of interest β0,β1 and σ2

b . The
trace plot of each parameter displays the sequence of sampled values (on the y-axis) at
each step in the MCMC chain (on the x-axis). It enables us to visualize the path taken by
the Markov Chain over the iterations. By looking at the trace plot, we can assess whether
the chain has converged to the target distribution. Ideally, the plot should look like a
’hairy caterpillar’ – it should oscillate around a constant mean without any trend, and it
should cover the entire range of plausible values for the parameter, meaning it mixes well.
Figure 1 show the values of β0, β1 and σ2

b within a stable range after a certain number of
iterations. This indicates that the MCMC simulation has ”converged” on an estimate for
β0, β1 and σ2

b . It also show that these three parameters values are ”mixing” well, meaning
they move freely and explore the entire range of plausible values rather than getting stuck
in particular regions [18]. This indicates that the MCMC simulation is effectively exploring
the full range of possible values for β0, β1 and σ2

b . We conclude from Figure 1 that the
MCMC simulation is likely to provide a very good estimate for this three parameter.

In this study, the benefits of Bayesian estimates based on the Median Ranked Set
Samples (MRSS) have been explored. This exploration is carried out via MCMC numerical
simulations following the steps outlined below:

• Generate Ranked Set Samples (RSS) of size ’n’ from the full posterior distribution
for the case when m = 1. Here, ’m’ refers to the number of cycles used in the
sampling process. In many practical applications, one cycle (m = 1) is commonly
used.



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 195

β0 β1

σ2

Figure 1: Trace plot of the parameters β0, β1, and σ2

• Use these RSS samples to calculate the Bayesian estimates as described in Section
5. These estimates rely on the MRSS samples.

The estimates are then analyzed in two ways:

• Point Estimation: The posterior mean and standard error (SE) for each parameter
is computed for different values of ’m’. These values are presented in Table 1. This
table includes the mean and SE for parameters β0, β1 and σ2

b , under different ’m’
values. The ’m’ values are increased from 3 to 7. A key observation from this table
is that as the value of ’m’ increases, the Bayesian standard error for all parameters
decreases. This indicates that as the cycle of sampling increases, the accuracy of the
estimates also increases.

• Interval Estimation: 95% credible intervals for the mean posterior of the parameters
are also calculated. These intervals, representing a range of values within which the
true parameter value lies with 95% certainty, are presented in Table 2. Here, for each
parameter and for different ’m’ values, the 2.5th percentile and the 97.5th percentile
values are reported. The range between these two percentile values forms the 95%
credible interval.

The tables illustrate that as the number of cycles ’m’ increases, the posterior means
become more accurate (SE decreases) and the credible intervals become narrower. These
outcomes suggest that by using more cycles in the MRSS sampling, we can achieve more
precise estimates and more confidence about where the true parameter values lie.
From the tables provided, we can also infer that as the number of cycles (m) increases,
not only does the precision of our estimates improve (seen by decreasing standard error),
but our credible intervals become narrower. This indicates that our certainty about the
true values of the parameters increases. The findings support the use of larger ’m’ values
in applications where precision is critical. However, larger ’m’ values may involve more



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 196

complex computations or more extensive data collection. The fact that standard errors
decrease as ’m’ increases suggests that the estimator becomes more efficient with more
cycles, which is advantageous in providing more reliable estimates.

Table 1: The Posterior mean and the standard error (SE) for the parameters at different m.

β̂0 β̂1 σ̂2

Mean SE Mean SE Mean SE

m=3 1.5 0.65 2.24 0.54 0.17 0.057

m=5 1.24 0.55 2.13 0.41 0.09 0.0096

m=7 1.12 0.34 2.03 0.30 0.07 0.0067

The 95% credible intervals for the parameters at different m.

β̂0 β̂1 σ̂2

2.5% 97.5% 2.5% 97.5% 2.5% 97.5%

m=3 1.22 2.87 2.04 3.55 0.05 0.18

m=5 1.12 2.06 1.85 2.34 0.06 0.11

m=7 0.97 1.76 1.67 2.65 0.05 0.10

Table 3 compare the efficiency of Bayesian and frequentist Mean Rank Set Sampling
(MRSS) methods for estimating the parameters of a simple regression model. The effi-
ciency is calculated as the ratio of the Mean Squared Errors (MSEs) of the parameter
estimates from the two methods. Specifically, the MSE of a parameter estimate using the
frequentist MRSS method (denoted by θ̂2) is divided by the MSE of the same parameter
estimate using the Bayesian MRSS method (denoted by θ̂1).
The table lists the efficiency ratios for the three parameters of interest: β0, β1 and σ2

b .
Each row in the table corresponds to a different value of ’m’, the size of the cycle in MRSS.
Each efficiency value in the table (like 1.29, 1.76, etc.) is calculated as the MSE of the
frequentist MRSS estimate divided by the MSE of the Bayesian MRSS estimate for a given
parameter and ’m’ value. This efficiency ratio measures how well the Bayesian method
performs in comparison to the frequentist method. The MSE is a common measure of an
estimator’s quality: the smaller the MSE, the better the estimator. The efficiencies are
all greater than one, indicating that the Bayesian MRSS method is more efficient than
the frequentist MRSS method for all parameter estimates and all values of ’m’. This
is because a smaller MSE (which implies better performance) for the Bayesian method
results in a larger efficiency ratio. As ’m’ increases from 3 to 7, the efficiencies also in-
crease. This shows that the advantage of the Bayesian MRSS method over the frequentist
MRSS method becomes more pronounced as ’m’ increases. In summary from Table 3, we
conclude that in this particular application, the Bayesian MRSS method provides more
efficient (lower MSE) estimates of the regression parameters than the frequentist MRSS
method, and this advantage increases as the size of the MRSS cycle (’m’) increases.



I. Nawajah, H. Kanj, Y. Kotb, J. Hoxha, M. Alakkoumi, K. Jebreen / Eur. J. Pure Appl. Math, 17 (1) (2024), 180-200 197

Table 3: The efficiency of Bayesian regression parameters using MRSS concerning MRSS with different m.

efficiency (θ̂1, θ̂2)= MSE(θ̂1)/MSE(θ̂2)

efficiency (β̂0, ˆβ0M ) efficiency (β̂1, ˆβ1M ) efficiency (σ̂2, σ̂2
M )

m=3 1.29 1.76 1.01

m=5 1.78 2.88 1.98

m=7 2.21 3.07 2.34

7. Discussion

In the present study, we developed a Bayesian model for estimating the parameters of
a simple linear regression model through median ranked set sampling (MRSS). Our results
demonstrate that Bayesian estimators derived from this approach are more efficient than
those from the frequentist MRSS method detailed in [6].
This study marks the first exploration of a Bayesian model for estimating coefficients of
a simple linear regression model via MRSS. The primary goal was to investigate the po-
tential of Bayesian statistical analysis for a straightforward and efficient estimation of the
simple linear regression model via MRSS, especially when compared with the traditional
frequentist MRSS approach.
The paper aligns with the results of Alodat et al., [6] and AlSaleh and AlSharafat [3],
emphasizing that Bayesian estimators prove to be more efficient than those derived from
simple random sampling. We also noticed a decrease in the standard error (SE) of Bayesian
estimates based on MRSS as the sample size m increases, reinforcing findings by De Iorio
et al. [10] and Helo et al., [14], both of whom demonstrated a similar trend in ranked set
sampling and the Bayesian estimation of Weibull parameters respectively.
This study assumes a symmetric distribution for random errors. When compared to the
frequentist model in Alodat et al. [6], our Bayesian model’s estimates exhibit a noticeable
advantage: they are simpler both in mathematical formulation and calculation. This ease
and efficiency in handling underline the potential advantages of Bayesian estimation when
dealing with MRSS in simple linear regression models.
The findings of this study bear significant implications for statistical modeling and infer-
ential analysis. They establish the value of Bayesian approaches in providing robust and
efficient estimates in the context of a simple linear regression model, using median ranked
set sampling (MRSS).
One key implication of the study is that it demonstrates the efficiency of Bayesian estima-
tors in comparison to their frequentist MRSS counterparts. This points to the potential
of Bayesian approaches to be more widely applied in statistical modeling and inferential
analysis, providing researchers with more robust and precise tools to analyze their data.
In addition, this study suggests that as the sample size (m) increases, the standard error
(SE) of Bayesian estimates based on MRSS decreases. This trend could inform future
data collection strategies, by encouraging the collection of larger sample sizes to improve
the precision of Bayesian estimates.



REFERENCES 198

8. Conclusion

Related to our knowledge, in the literature, there are no previous studies on Bayesian
estimating of the linear regression parameters using MRSS. Thus, the emphasis of this pa-
per is to estimate the simple linear regression parameters via MRSS in Bayesian approach.

This study investigated the Bayesian estimation of parameters in a simple linear re-
gression model, employing MRSS. Using Markov Chain Monte Carlo (MCMC) numerical
simulations, Bayesian estimates were obtained and analyzed. The findings underscored
the efficiency of Bayesian estimators obtained through MRSS, surpassing their frequentist
counterparts using the same MRSS methodology, at least for one MRSS design. Addi-
tionally, it was observed that the Bayesian standard error for all parameters diminished
as the values of m increased.
However, this study is not without its limitations. The efficiency of the Bayesian estima-
tors was established in the context of one MRSS specific design and this may not hold
true for all MRSS design. MRSS is a specific type of RSS where the median item from
each subset is chosen. The way these subsets are chosen and ranked can vary, and these
specifics constitute the MRSS design. Further research should investigate the performance
of Bayesian estimators in a more diverse range of MRSS design. Additionally, this study
has taken a symmetric distribution of random errors as an assumption. Future investi-
gations could benefit from considering asymmetric or heavy-tailed error distributions to
expand the scope of this research.
Moreover, while this study highlights the increased efficiency and simplicity of the Bayesian
approach, it is crucial to further explore the computational aspects and feasibility of these
methods, especially in the context of large-scale data. Future research should also delve
into the practical implications and applications of the methodology proposed in this study,
in various fields and for a variety of research questions.
In conclusion, this study has paved the way for further research in the field of Bayesian
estimation using MRSS, providing a springboard for more extensive studies, wider appli-
cations, and refined methodologies.

Acknowledgements

The authors are grateful to the Hebron University, Palestine American University of
the Middle East, Egaila 54200, Kuwait, and Palestine Technical University – Kadoorie,
Palestine, for their support in completing this research.

References

[1] Said Ali Al-Hadhrami and Amer Ibrahim Al-Omari. Bayesian inference on the variance of
normal distribution using moving extremes ranked set sampling. Journal of Modern Applied
Statistical Methods, 8(1):25, 2009.

[2] M Fraiwan Al-Saleh, Khalaf Al-Shrafat, and H Muttlak. Bayesian estimation using ranked set



REFERENCES 199

sampling. Biometrical Journal: Journal of Mathematical Methods in Biosciences, 42(4):489–
500, 2000.

[3] Mohammad Fraiwan Al-Saleh and Khalaf Al-Shrafat. Estimation of average milk yield using
ranked set sampling. Environmetrics: The official journal of the International Environmetrics
Society, 12(4):395–399, 2001.

[4] AA Alhamide, Kamarulzaman Ibrahim, MT Alodat, and Wan Zawiah Wan Zin. Bayesian
inference for linear regression under alpha-skew-normal prior. Sains Malaysiana, 48(1):227–
235, 2019.

[5] Iftikhar Ali, Nabaa Muhammad Diaa, Muhammad Haroon Aftab, Muhammad Waheed Ra-
heed, Kamel Jebreen, and Hassan Kanj. Topological effects of chiral pamam dendrimer for
the treatment of cancer. 31(2).

[6] MT Alodat, MT AL-Rawwash, and IM Nawajah. Analysis of simple linear regression via
median ranked set sampling. Metron, 67(1):57–74, 2009.

[7] MT Alodat, MY Al-Rawwash, and IM Nawajah. Inference about the regression parame-
ters using median-ranked set sampling. Communications in Statistics—Theory and Methods,
39(14):2604–2616, 2010.

[8] MT Alodat and OA Al-Sagheer. Estimation the location and scale parameters using ranked
set sampling. J. Appl. Statis. Sci, 15:245–252, 2007.

[9] Kara Gülay Begüm and Demirel Neslihan. Two-layer median ranked set sampling. Hacettepe
Journal of Mathematics and Statistics, page 1–10.

[10] Maria De Iorio, Peter Müller, Gary L Rosner, and Steven N MacEachern. An anova model for
dependent random measures. Journal of the American Statistical Association, 99(465):205–
215, 2004.

[11] Begüm Kara Gülay and Neslihan Demirel. Two-layer median ranked set sampling. Hacettepe
Journal of Mathematics and Statistics, 48:1–10, 2019.

[12] Abdul Haq, Jennifer Brown, and Elena Moltchanova. Improved best linear unbiased esti-
mators for the simple linear regression model using double ranked set sampling schemes.
Communications in Statistics-Theory and Methods, 45(12):3541–3561, 2016.

[13] Amal S Hassan. Maximum likelihood and bayes estimators of the unknown parameters for
exponentiated exponential distribution using ranked set sampling. International Journal of
Engineering Research and Applications, 3(1):720–725, 2013.

[14] Amal Helu, Muhammad Abu-Salih, and Osama Alkam. Bayes estimation of weibull distri-
bution parameters using ranked set sampling. Communications in Statistics—Theory and
Methods, 39(14):2533–2551, 2010.

[15] Kamel Jebreen. Modèles graphiques pour la classification et les séries temporelles. Aix-
Marseille.

[16] Kamel Jebreen, Muhammad Aftab, Iftikhar Ali, Mohammad Sowaity, and Hassan Kanj. Topo-
logical aspects investigated from m-polynomial of -sheet of boron clusters. International Jour-
nal of Chemical and Biochemical Sciences, 24(4):469–477, 2023.



REFERENCES 200

[17] Kamel Jebreen and Badih Ghattas. Bayesian network classification: Application to epilepsy
type prediction using PET scan data. In 2016 15th IEEE International Conference on Machine
Learning and Applications (ICMLA), pages 965–970.

[18] Kamel Jebreen, Mohamad Motasem Nawaf, Amjad Barham, and Badih Ghattas. Inferring
linear and nonlinear interaction networks using neighborhood support vector machines. In
2021 International Conference on Engineering and Emerging Technologies (ICEET), pages
1–6. ISSN: 2409-2983.

[19] Hassan Kanj, Hifza Iqbal, Muhammad Haroon Aftab, Hasnain Raza, Kamel Jebreen, and
Mohammed Issa Sowaity. Topological characterization of hexagonal network and non-kekulean
benzenoid hydrocarbon. 16(4):2187–2197.

[20] Jessica K Kohlschmidt, Elizabeth A Stasny, and Douglas A Wolfe. Ranked set sampling for a
population proportion: Allocation of sample units to each judgment order statistic. Pakistan
Journal of Statistics and Operation Research, pages 511–530, 2012.

[21] Tao Li and Narayanaswamy Balakrishnan. Best linear unbiased estimators of parameters of
a simple linear regression model based on ordered ranked set samples. Journal of Statistical
Computation and Simulation, 78(12):1267–1278, 2008.

[22] GA McIntyre. A method for unbiased selective sampling, using ranked sets. Australian journal
of agricultural research, 3(4):385–390, 1952.

[23] HA Muttlak. Median ranked set sampling with concomitant variables and a comparison with
ranked set sampling and regression estimators. Environmetrics: The official journal of the
International Environmetrics Society, 9(3):255–267, 1998.

[24] Martyn Plummer et al. Jags: A program for analysis of bayesian graphical models using
gibbs sampling. In Proceedings of the 3rd international workshop on distributed statistical
computing, volume 124, pages 1–10. Vienna, Austria., 2003.

[25] Hakan Savaş Sazak and Melis Zeybek. The modified maximum likelihood estimators for the
parameters of the regression model under bivariate median ranked set sampling. Computa-
tional Statistics, 37(3):1069–1109, 2022.

[26] Arun Kumar Sinha. Ranked set sampling: As a cost-effective and more efficient data collection
method. Statistical Journal of the IAOS, 32(4):607–611, 2016.

[27] Douglas A Wolfe. Ranked set sampling: its relevance and impact on statistical inference.
International scholarly research notices, 2012, 2012.

[28] Dong-sen Yao, Wang-xue Chen, and Chun-xian Long. Parametric estimation for the sim-
ple linear regression model under moving extremes ranked set sampling design. Applied
Mathematics-A Journal of Chinese Universities, 36(2):269–277, 2021.