# Bayesian Inference for Johnson’s SB and Weibull distributions

Mahdi Teimouri  Department of Mathematics and Statistics, Faculty of Science
and Engineering, Gonbad Kavous University, Gonbad Kavous, Iran, Email:
<EMAIL_ADDRESS>
###### Abstract

The four-parameter Johnson’s SB (JSB) and three-parameter Weibull
distributions have received much attention in the field of forestry for
characterizing diameters at breast height (DBH). In this work, we suggest the
Bayesian method for estimating parameters of the JBS distribution. The maximum
likelihood approach uses iterative methods such as Newton-Raphson (NR)
algorithm for maximizing the logarithm of the likelihood function. But there
is no guarantee that the NR method converges. This fact that the NR method for
estimating the parameters of the JSB distribution sometimes fails to converge
was verified through simulation in this study. Further, it was shown that the
Bayesian estimators presented in this work were robust with respect to the
initial values and estimate the parameters of the JSB distribution
efficiently. The performance of the JSB and three-parameter Weibull
distributions was compared in a Bayesian paradigm when these models were
fitted to DBH data of three plots that randomly selected from a study
established in 107 plots of mixed-age ponderosa pine (Pinus ponderosa Dougl.
ex Laws.) with scattered western junipers at the Malheur National Forest in
south end of the Blue Mountains near Burns, Oregon, USA. Bayesian paradigm
demonstrated that JBS was superior model than the three-parameter Weibull for
characterizing the DBH distribution when these models were fitted to the DBH
data of the three plots.

###### keywords:

Bayesian analysis , Diameter distribution , Forest management , Johnson’s SB
distribution , Maximum likelihood method , Weibull distribution

##  1 Introduction

Statistical modelling for the distribution of the diameter at breast height
(DBH) is becoming increasingly popular in order to characterizing the forest
height structure, forest dynamics, and comparing the forest stands (Gorgoso et
al., 2007; Mateus and Tomé, 2011; Özçelik et al., 2016). The statistical
characterization or modelling of the DBH distribution has a long history in
both managed and natural forest stands. Among all statistical models, the
desired is that shows more flexibility, i.e., capturing well the DBH
distribution. This is because different types of the forest stands show
different shape for DBH distribution. For example, two main types of the
forest stands include the even-aged that usually are unimodal (one peak) and
roughly symmetric and uneven-aged that whose DBH distributions often have a
reverse-$J$ shape. Among all statistical distributions, the Johnson’s SB (JSB)
and Weibull have received much attention in the context forest management.
Numerous efforts have been made in the literature for modelling the tree’s DBH
among them are (Bailey and Dell, 1973; Maltamo et al., 1995, 2000; Pretzsch,
2009; Zhang et al., 2010) for two- or three-parameter Weibull distribution,
(Fonseca et al., 2009; Hafley and Buford, 1985; Kiviste et al., 2003; Kudus et
al., 1999; Marto et al., 2009; Mateus and Tomé, 2011; Özçelik et al., 2016;
Parresol, 2003; Rennolls and Wang, 2005; Zhou and McTague, 1996) for JSB
distribution, and (Gorgoso et al., 2012; Hafley and Schreuder, 1977; Palahí et
al., 2007; Zhang et al., 2003) for both of them. In statistical modelling of
DBH, precision of the parameter estimator plays a crucial role in the forest
planning and management (Gorgoso et al., 2007; Mateus and Tomé, 2011; Özçelik
et al., 2016).

The maximum likelihood (ML) approach, as the most popular estimation method,
is obtained with the aid of mathematical optimization tools. These tools
maximize the logarithm of the likelihood (log-likelihood) function using
iterative algorithm such as Newton-Raphson (NR) and so need the initial
values. If the initial values are far away from the true parameter (which is
where the log-likelihood function reaches its global maximum), or when the
log-likelihood function at the initial values becomes large, then there is no
guarantee that the NR method will converge. This means that the ML approach is
sensitive to the initial values and really can be considered as a weakness for
the ML approach. Also, the ML approach may break down when the regularity
conditions fail to exist. The above criticisms may happen when one is
interested in estimating the parameters of JSB and three-parameter Weibull
distributions. Other methods such as moment-based estimators are not as
efficient as the ML approach. For example, moment-based estimators of the
three-parameter Weibull distribution overcome the weaknesses in the ML
approach, but their existence, uniqueness, and consistency are still open
questions (Nagatsuka et al., 2013) or in the case of JBS distribution, the
moment-based estimator are not as efficient as regression-type estimators
(Scolforo et al., 2003). The aim of this study is to derive the Bayesian
estimators for the parameters of the JSB and three-parameter Weibull
distributions, that to the best of our knowledge, the Bayesian estimators of
the parameters of JSB distribution have never been used in the forestry
literature for modeling DBH distributions. This paper is organized as follows.
In what follows we give some preliminaries. The Bayesian paradigm for the JSB
and Weibull distributions are given in Section 2. Section 3 is devoted to the
materials and methods. We give the results and discussion in Sections 4 and 5,
respectively. We conclude the paper in Section 6.

In the following, we give some preliminaries.

###  1.1 The family of JSB and three-parameter Weibull distributions

The probability density function (pdf) and cumulative distribution function
(cdf) of the Johnson SB (JSB) are given, respectively, by (Johnson, 1949;
Norman et al., 1994):

$g_{J}\bigl{(}x\big{|}\Theta\bigr{)}=\frac{\delta\lambda}{\sqrt{2\pi}(x-\xi)(\lambda+\xi-x)}\exp\Biggl{\\{}-\frac{1}{2}\Bigg{[}\gamma+\delta\log\biggl{(}\frac{x-\xi}{\lambda+\xi-x}\biggr{)}\Bigg{]}^{2}\Biggr{\\}},$
(1)

and

$G_{J}\bigl{(}x\big{|}\Theta\bigr{)}=\int_{\xi}^{x}g_{J}\bigl{(}y\big{|}\Theta\bigr{)}dy,$
(2)

where $\Theta=(\delta,\gamma,\lambda,\xi)^{T}$, $\xi<x<\lambda+\xi$,
$\delta>0$, $\lambda>0$, $-\infty<\gamma<\infty$, and $-\infty<\xi<\infty$. As
it is seen, the cdf of the JBS distribution has no closed-form expression. The
pdf and cdf of three-parameter Weibull distribution are given, respectively,
by

$g_{W}\bigl{(}x\big{|}\Theta\bigr{)}=\frac{\alpha}{\beta}\biggl{(}\frac{x-\mu}{\beta}\biggr{)}^{\alpha-1}\exp\Biggl{\\{}-\biggl{(}\frac{x-\mu}{\beta}\biggr{)}^{\alpha}\Biggr{\\}},$
(3)

and

$G_{W}\bigl{(}x\big{|}\Theta\bigr{)}=1-\exp\Biggl{\\{}-\biggl{(}\frac{x-\mu}{\beta}\biggr{)}^{\alpha}\Biggr{\\}},$
(4)

where $\Theta=(\alpha,\beta,\mu)^{T}$, $\mu<x$, $\alpha>0$, and $\beta>0$.
Now, $\alpha$, $\beta$, and $\mu$ are the shape, scale, and location
parameters, respectively.

###  1.2 Bayes theorem

In the Bayesian framework, we assume that the unknown parameter vector
$\Theta$ follows a distribution with pdf $\pi(\Theta)$. Using information
available in random observations $\boldsymbol{x}=(x_{1},\dots,x_{n})^{T}$, a
revision will be made on knowledge about $p(\Theta)$ using the well-known
Bayes’ theorem as $\pi(\Theta|\boldsymbol{x})$. We have:

$\pi(\Theta|\boldsymbol{x})=\frac{g(\boldsymbol{x}|\Theta)\pi(\Theta)}{g(\boldsymbol{x})}$

The expression $\pi(\Theta)$ and $\pi(\Theta|\boldsymbol{x})$ are known in the
literature as prior pdf and posterior pdf of $\Theta$, respectively. Here,
$g(\boldsymbol{x})$ is normalizing constant and so the Bayes’ theorem can be
written as

$\pi(\Theta|\boldsymbol{x})\propto g(\boldsymbol{x}|\Theta)\pi(\Theta)$ (5)

###  1.3 The NR algorithm for JSB distribution

As previously mentioned, the NR algorithm may fail to converge. Unfortunately,
this happens when finding the ML estimators of the JSB distribution is
desired. We performed a simulation study to prove our claim. So, a number of
10,000 samples with different sizes including 20, 50, 100, 250, 500, 1000, and
5,000 were simulated from the JSB distribution with pdf given in (1) and then
we obtained the ML estimators of the parameters using the command optim(.) in
R (team2018r) environment. In each of 10,000 runs, the parameters $\delta$,
$\gamma$, $\lambda$, and $\xi$ were generated from uniform distribution
(0.05,10), (-20,20), (1,100), and (-50,50), respectively. While
$x^{(i)}_{(1)}$ and $x^{(i)}_{(1)}$ denote, respectively, the smallest and
largest values of the $i$-th generated sample, for $i=1,\dots,10000$, the
initial values of $\delta$, $\gamma$, $\lambda$, and $\xi$ were generated from
uniform distribution (0.05,10), (-20,20),
$\Bigl{(}x^{(i)}_{(n)}-x^{(i)}_{(1)},100\Bigr{)}$, and
$\Bigl{(}-50,x^{(i)}_{(1)}\Bigr{)}$, respectively. The results of simulation
are given in Table 1. As it is seen, percentage of failed attempts to reach
convergence through the NR algorithm is considerable (say, on the average
32%).

###  1.4 Gibbs sampler

Due to the complicate nature of the posterior pdf
$\pi(\Theta|\boldsymbol{x})$, we have to sample from the posterior pdf and
then the Bayesian estimators are the ergodic average of the generated sample.
In practice, exploitation of $\pi(\Theta|\boldsymbol{x})$ needs the use of
Markov chain Monte Carlo (MCMC) techniques. The Gibbs sampler is one of the
Markov chain Monte Carlo (MCMC) techniques that enables us to sample from full
conditional pdf, i.e., the pdf of each element of the parameter vector given
the other elements and observed data $\boldsymbol{x}=(x_{1},\dots,x_{n})^{T}$.
Suppose that we have a statistical distribution with unknown parameter vector
$\Theta=(\theta_{1},\theta_{2}\dots,\theta_{k})^{T}$ and observed
$\boldsymbol{x}$. Consider
$\Theta^{(0)}=\bigl{(}\theta^{(0)}_{1},\theta^{(0)}_{2}\dots,\theta^{(0)}_{k}\bigr{)}^{T}$
as the initial values. In order to implement the Gibbs sampler technique, we
generate $\theta^{(1)}_{1}$ from the full conditional pdf
$\pi\bigl{(}\theta_{1}\big{|}\theta^{(0)}_{2},\theta^{(0)}_{3},\dots,\theta^{(0)}_{k},\boldsymbol{x}\bigr{)}$,
$\theta^{(1)}_{2}$ from the full conditional pdf
$\pi\bigl{(}\theta_{2}\big{|}\theta^{(1)}_{1},\theta^{(0)}_{3},\dots,\theta^{(0)}_{k},\boldsymbol{x}\bigr{)}$,
and so on up to $\theta^{(1)}_{k}$ from the full conditional pdf
$\pi\bigl{(}\theta_{k}\big{|}\theta^{(1)}_{1},\theta^{(1)}_{2},\dots,\theta^{(1)}_{k-1},\boldsymbol{x}\bigr{)}$.
In the first iteration we obtain the sample
$\Theta^{(1)}=\bigl{(}\theta^{(1)}_{1},\theta^{(1)}_{2}\dots,\theta^{(1)}_{k}\bigr{)}^{T}$.
Under mild regularity conditions (Roberts and Smith, 1994), after a
sufficiently large number of iterations, say $N$, the ergodic average of the
Markov chain yields a consistent estimator of $\Theta$.

###  1.5 Metropolis-Hastings algorithm

The Metropolis-Hastings (MH) algorithm, is an MCMC technique and efficient
method for drawing samples from a given posterior distribution. Assume that we
want to simulate sample from pdf given by (5). Firstly, we need to choose a
proposal distribution $q(.|.)$, that changes the location of the chain at each
iteration of the algorithm. The proposal distribution is arbitrary and can be
chosen so that is easy to simulate from. Secondly, we follow the steps given
by the following.

  * 1.

Choose initial value $\theta^{(0)}$ and set i=1;

    1. (a)

Sample $\theta^{*}$ from $q\bigl{(}.\big{|}\theta^{(i-1)}\bigr{)}$;

    2. (b)

Set $\theta^{(i)}=\theta^{*}$ with probability

$\eta=\min\left\\{1,\frac{\pi\bigl{(}\theta^{*}\big{|}\boldsymbol{x}\bigr{)}q\bigl{(}\theta^{(i-1)}\big{|}\theta^{*}\bigr{)}}{\pi\bigl{(}\theta^{(i-1)}\big{|}\boldsymbol{x}\bigr{)}q\bigl{(}\theta_{*}\big{|}\theta^{(i-1)}\bigr{)}}\right\\},$

otherwise $\theta^{(i)}=\theta^{(i-1)}$;

    3. (c)

Set i=i+1 and go to step 1;

  * 2.

Stop the algorithm if $i=N$ ($N$ is a sufficiently large integer value).

Since realization generated at each iteration is used to generate sample at
next step, the chain constitutes a correlated stochastic process, but after
$N$ numbers of generations, we hope that the chain produces uncorrelated
samples and converges to the target distribution $\pi(\theta|\boldsymbol{x})$
as desired.

###  1.6 Adaptive rejection sampling

For a given sample $\boldsymbol{x}=(x_{1},\dots,x_{n})^{T}$, suppose we are
interested in sampling from posterior pdf $\pi(\theta|\boldsymbol{x})$. The
well-known accept-reject sampling needs a suitable upper bound, called $M$
that satisfies in the following inequality.

$\underset{\theta}{\operatorname{sup}}\frac{\pi(\theta|\boldsymbol{x})}{g(\theta)}\leq
M,$

where $g(\theta)$ is an arbitrary pdf that is easy to sample from and its
support includes the support of $\pi(\theta|\boldsymbol{x})$. If there exists
a suitable choice for $M$, then the accept-reject sampling is efficient
otherwise we refer to another Monte Carlo simulation technique, called
adaptive rejection sampling (ARS) algorithm. The ARS algorithm is used to
simulate realization when posterior pdf is log-concave (i. e., the second
derivative of $\pi(\theta|\boldsymbol{x})$ with respect to $\theta$ is
negative). In such a case, the ARS algorithm developed by Gilks and Wild
(1992) is highly efficient.

##  2 Bayesian analysis

Here, we give the Bayesian paradigm for the JSB and Wribull distributions,
respectively.

###  2.1 Gibbs sampler for JSB distribution

For Bayesian inference of the JSB distribution parameters, we assume that all
four priors are statistically independent and so the full Bayesian model
(joint posterior pdf) up to proportionality becomes

$\displaystyle\pi(\Theta|\boldsymbol{x})$ $\displaystyle\propto
g_{J}(\boldsymbol{x}|\Theta)\pi(\delta,\gamma,\lambda,\xi)$
$\displaystyle=\Pi_{i=1}^{n}g_{J}(x_{i}|\Theta)\pi(\delta)\pi(\gamma)\pi(\lambda)\pi(\xi)$
$\displaystyle=\frac{\delta^{n}\lambda^{n}}{(2\pi)^{\frac{n}{2}}\Pi_{i=1}^{n}(x_{i}-\xi)(\lambda+\xi-
x_{i})}\exp\left\\{-\frac{1}{2}\sum_{i=1}^{n}\Bigg{[}\gamma+\delta\log\biggl{(}\frac{x_{i}-\xi}{\lambda+\xi-
x_{i}}\biggr{)}\Bigg{]}^{2}\right\\}\pi(\delta)\pi(\gamma)\pi(\lambda)\pi(\xi).$
(6)

We note that we produce the Gibbs sampler directly from the full conditionals
by choosing improper priors, i.e., bypassing the propriety of the posterior.
This is due to the fact that the full conditionals are well-defined and the
Bayesian model under study is enough complex (Robert et al., 2010). Details
for generating from full conditionals of $\delta$, $\gamma$, $\lambda$, and
$\xi$ are given A, B, C, and D, respectively.

###  2.2 Gibbs sampler for three-parameter Weibull distribution

The Bayesian paradigm for the three-parameter Weibull distribution originally
was developed by Smith and Naylor (1987) and Green et al. (1994). Here, we
give a slightly different version of the Bayesian paradigm developed by Green
et al. (1994). We mention that the only difference occurs in updating the
location parameter at each iteration if the chain. Assume that
$\boldsymbol{x}=(x_{1},\dots,x_{n})^{T}$ denotes the vector of $n$ independent
observations each follows distribution with pdf given in (3). We consider the
Jeffreys’ prior (Jeffreys, 1961) for $\alpha$ and $\beta$, i.e.,
$\pi(\alpha)\propto 1/\alpha$ and $\pi(\beta)\propto 1/\beta$. Also, we allow
the prior for $\mu$ to be uniformly over ${{\mathbb{R}}}$. So, the full
Bayesian model is given by

$\displaystyle\pi(\Theta|\boldsymbol{x})$ $\displaystyle\propto
g_{W}(\boldsymbol{x}|\Theta)\pi(\alpha,\beta,\mu)$
$\displaystyle=\Pi_{i=1}^{n}g_{W}(x_{i}|\Theta)\pi(\alpha)\pi(\beta)\pi(\mu)$
$\displaystyle=\frac{\alpha^{n-1}}{\beta^{n+1}}\Pi_{i=1}^{n}\biggl{(}\frac{x_{i}-\mu}{\beta}\biggr{)}^{\alpha-1}\exp\left\\{-\sum_{i=1}^{n}\biggl{(}\frac{x_{i}-\mu}{\beta}\biggr{)}^{\alpha}\right\\}.$
(7)

The full conditionals of $\alpha$, $\beta$, and $\mu$ are (up to
proportionality) given by

$\displaystyle\pi(\alpha|\beta,\mu,\boldsymbol{x})$
$\displaystyle\propto\alpha^{n-1}\Pi_{i=1}^{n}\biggl{(}\frac{x_{i}-\mu}{\beta}\biggr{)}^{\alpha-1}\exp\left\\{-\sum_{i=1}^{n}\biggl{(}\frac{x_{i}-\mu}{\beta}\biggr{)}^{\alpha}\right\\},$
(8) $\displaystyle\pi(\beta|\alpha,\mu,\boldsymbol{x})$
$\displaystyle\propto\beta^{-n\alpha-1}\exp\left\\{-\sum_{i=1}^{n}\biggl{(}\frac{x_{i}-\mu}{\beta}\biggr{)}^{\alpha}\right\\},$
(9) $\displaystyle\pi(\mu|\alpha,\beta,\boldsymbol{x})$
$\displaystyle\propto\Pi_{i=1}^{n}\bigl{(}x_{i}-\mu\bigr{)}^{\alpha-1}\exp\left\\{-\sum_{i=1}^{n}\biggl{(}\frac{x_{i}-\mu}{\beta}\biggr{)}^{\alpha}\right\\}.$
(10)

As pointed out by Green et al. (1994), since
$\pi(\alpha|\beta,\mu,\boldsymbol{x})$ given in (8) is log-concave, the ARS
algorithm is highly efficient technique for simulating from this full
conditional. Assuming that we are currently at the $t$-th iteration of the
chain and we have just obtained $\alpha^{(t+1)}$ by simulating from full
conditional given in (8). In order to generate from the full conditional
$\pi\bigl{(}\beta\big{|}\alpha^{(t+1)},\mu^{(t)},\boldsymbol{x}\bigr{)}$ given
in (9), it suffices to simulate a realization from gamma distribution with
shape parameter $n$, say $z$, and then update $\beta^{(t)}$ as $\beta^{(t+1)}$
using the relation

$\displaystyle\beta^{(t+1)}=\Biggl{(}\frac{\sum_{i=1}^{n}\bigl{(}x_{i}-\mu^{(t)}\bigr{)}^{\alpha^{(t+1)}}}{z}\Biggr{)}^{\frac{1}{\alpha^{(t+1)}}}.$

For the location parameter with full conditional given in (10), we do not
follow the accept-reject sampling method proposed by Green et al. (1994). Our
study revealed that when $\alpha$ is small (say $\alpha<2$) the accept-reject
sampling method does not work efficiently (the chin takes too much time for
updating). Instead, we use the MH algorithm for generating from full
conditional of $\mu$. For this purpose, we use the uniform distribution on
$\bigl{(}x_{(1)}-\beta,x_{(1)}\bigr{)}$ as the proposal. The steps of the MH
algorithm are given by the following.

  1. 1.

Suppose we are currently at the $t$-th iteration of the chain. Set $i=1$ and
choose the initial value as $\mu^{(0)}=x_{(1)}-1/n$;

  2. 2.

Sample $\mu^{*}$ from proposal distribution with pdf
$q\bigl{(}\mu)=1/\beta^{(t+1)}$, for $x_{(1)}-\beta^{(t+1)}<\mu<x_{(1)}$;

  3. 3.

Compute $\eta$ as

$\eta=\min\left\\{1,\frac{\pi\bigl{(}\mu^{*}\big{|}\alpha^{(t+1)},\beta^{(t+1)},\boldsymbol{x}\bigr{)}}{\pi\bigl{(}\mu_{(i-1)}\big{|}\alpha^{(t+1)},\beta^{(t+1)},\boldsymbol{x}\bigr{)}}\right\\}.$

  4. 4.

Generate an uniform random variable on $(0,1)$, say $u$. If $u<\eta$, then
$\xi^{(i)}=\xi^{*}$, otherwise $\xi^{(i)}=\xi^{(i-1)}$;

  5. 5.

If $i=N$, then go to the next step. Otherwise set $i=i+1$ and go to step 2;

  6. 6.

Accept $\mu_{(N)}$ as a generation form pdf
$\pi\bigl{(}\mu\big{|}\alpha^{(t+1)},\beta^{(t+1)},\boldsymbol{x}\bigr{)}$,
i.e., $\mu^{(t+1)}=\mu_{(N)}$ and stop the MH algorithm.

We note that the MH algorithm adopted for simulating from full conditional of
$\mu$ is faster than the accept-reject sampling method proposed by Green et
al. (1994) when $\alpha<2$. Also, by choosing a uniform proposal, we allow the
location parameter to vary over real line. This feature of our proposed prior
for location parameter will appeal to a wide range of study fields in which
modelling data through three-parameter Weibull distribution with negative
location parameter occurs frequently.

##  3 Materials and methods

###  3.1 Materials

Since our work motivated by the widespread use and application of the
statistical distributions in the forest management, we conducted a study in
the context of the forestry. A study was established in 107 plots of mixed-age
ponderosa pine (Pinus ponderosa Dougl. ex Laws.) with scattered western
junipers (_Juniperus occidentalis_ Hook.) that located in Malheur National
Forest in south end of the Blue Mountains near Burns, Oregon, USA Kerns et al.
(2017). These data include tree height, diameter, and growth for a prescribed
burning study with unburned controls. Of these variables, we only used the DBH
(measured at a height of 1.3 m) of all live trees in three randomly selected
plots (plots 9, 73, and 81) each of size 0.08 ha for statistical validation of
the Bayesian approach. The plots summary statistics are given in Table 2.

####  3.1.1 Methods

The Bayesian approach were applied to the DBH data addressed in the previous
subsection. We compared the performance of the JBS and three-parameter Weibull
distributions for modelling DBH data when the parameters of both models were
estimated using the Bayesian approach. Figure 1 and Figure 2 display
histograms of the samples drawn from the full conditionals and the pairwise
scatterplots of the sampler output for the JSB and three-parameter Weibull
distributions, respectively, when these distributions fitted to DBH data of
the plot 9. Accordingly, Figure 3 and Figure 4 show histograms of the samples
drawn from the full conditionals and pairwise scatterplots of the sampler
output for the JSB and three-parameter Weibull distributions, respectively,
when these distributions fitted to the DBH data in the plot 81. Finally, the
graphical visualizations when the JSB and three-parameter Weibull
distributions fitted to the DBH data in the plot 73 are shown in Figure 5 and
Figure 6, respectively. We assumed that the sampler’s convergence has been
attained before 5,000 iterations in all three plots for both the JSB and
three-parameter Weibull distributions. Therefore, we removed the first 5000
samples from the sampler output when the sampler were repeated for 10,000
times. Based on the average of the final 5,000 samples, we obtained the
Bayesian estimators of the models parameters. The estimated parameters are
given in Table 4. The goodness-of-fit statistics including Anderson-Darling
(AD), Cramér–von Mises (CM), Kolmogorov-Smirnov (KS), and log-likelihood (LL)
statistics were used as criterion for choosing the better model. These
criteria are defined as

$\displaystyle AD$ $\displaystyle=-n-\sum_{i=1}^{n}\frac{2i-1}{n}\biggl{[}\log
G\bigl{(}x_{(i)}\big{|}\hat{\Theta}\bigr{)}+\log\Bigl{(}1-G\bigl{(}x_{(n+i-1)}\big{|}\hat{\Theta}\bigr{)}\Bigr{)}\biggr{]},$
$\displaystyle CM$
$\displaystyle=\frac{1}{12n}+\sum_{i=1}^{n}\biggl{[}G\bigl{(}x_{(i)}\big{|}\hat{\Theta}\bigr{)}-\frac{2i-1}{2n}\biggr{]}^{2},$
$\displaystyle KS$
$\displaystyle=\underset{x_{i}\in{\boldsymbol{x}}}{\operatorname{sup}}\Big{|}G\bigl{(}x_{i}\big{|}\hat{\Theta}\bigr{)}-G_{n}(x_{i})\Big{|},$
$\displaystyle LL$ $\displaystyle=\sum_{i=1}^{n}\log
g\bigl{(}x_{i}\big{|}\hat{\Theta}\bigr{)},$

where $G_{n}(.)$ denotes the empirical distribution function, $x_{(i)}$ is th
$i$-th order statistic in the random sample of size $n$, and $\hat{\Theta}$ is
the estimated vector parameter. The computed goodness-of-fit statistics are
given in Table 5. For analyzing the DBH data given in Table 2, we used the R
package `ForestFit` (Teimouri et al., 2020) developed for R (R Core Team,
2018) and uploaded to CRAN (Comprehensive R Archive Network) at
`https://cran.r-project.org/web/packages/Fores`
`tFit/index.html`.

##  4 Results

The DBH values ranged from 9.1 to 88.6 cm and the mean of DBH ranged between
24.27 and 33.35 cm in DBH (Table 2). The DBH distributions were usually
continuous with peaks at the lower, middle, and higher bins (Fig. 7). Mixed
forests with ponderosa pine and scattered western junipers are usually a
mosaic composed of different forest patches. The analysis of three randomly
selected plots represented the following results.

Plot 9 was characterized by DBH distribution that had trees in the size
between 10.4 and 55.9 cm in DBH (Fig. 2). The understory cohort was
characterized by a distribution with an understory cohort that has one peak
around 15 cm in DBH, and truncated at 60 cm in DBH (Figure 7). Overall, the
JSB model was the superior model than the Weibull model for DBH distribution
(Table 4). The DBH distribution was best captured better by JSB distribution
than the Weibull. The JSB was the superior model in the sense of AD, CM, KS,
and LL measures (Table 4). The shapes of the fitted distributions to the DBH
data of plot 9 were different in the right tail (Fig. 7). In fact, the JSB was
characterized the tail of the DBH distribution better than the Weibull
distribution.

Plot 44 was characterized by DBH distribution that had trees in the size
between 11.9 and 83.8 cm in DBH (Table 2). This plot represented a one-storied
forest that was characterized by a DBH distribution with three modes around
15, 35, and 60 cm in DBH. Ignoring these modes, the understory cohort was
characterized by a broadly uniform distribution with an understory cohort
between 9-60 cm in DBH (Figure 7). Overall, the JSB model was the superior
model than the Weibull model for DBH distribution (Table 4). The shape of the
fitted densities corresponding to the JSB and Weibull distributions have
differences in the middle and the right tail (Figure 7). This is largely
because the JSB has two threshold parameters that makes its pdf bounded from
the left and right.

Plot 73 was characterized by DBH distribution that had trees in the size
between 9.1 and 88.6 cm in DBH (Table 2). This plot was characterized by a
generally negative exponential or reverse-_J_ DBH distribution. Overall, this
type of DBH distribution occurs in uneven-aged, complex structure, and old-
growth forests where the number of trees declines sharply with increasing tree
size (Figure 7). There are two peaks around 15 and 35 cm in DBH. In contrast
to the plot 9 and 44, the DBH distribution of plot 73 was much wider with a
long and skewed to the right tail. The JSB model was characterized the DBH
distribution better than the Weibull model (Table 4). Evidently, the JSB was
the superior model than the Weibull model for DBH observations at the first
peak, between two peaks, and right tail (Figure 7). In all three plots, the
JSB outperformed the Weibull model in terms of all AD, CM, KS, and LL measures
(Table 4).

##  5 Discussion

For implementing the Bayesian paradigm, we carried out a study for choosing
the initial values. In the case of the JSB distribution, it is known that the
first order statistic, i.e., $x_{(1)}$ is a sufficient statistic for $\xi$.
So, choosing $\xi^{(0)}=x_{(1)}-1/n$ in which $n$ is the sample size would be
quite reasonable as the initial value for $\xi$. In the same fashion as for
$\xi$, a good initial value for $\lambda$, since $\xi<x<\xi+\lambda$, is given
by $\lambda^{(0)}=x_{(n)}-x_{(1)}+2/n$ in which $x_{(n)}$ is the maximum value
of DBH in the sample, i.e., $x_{(n)}$. Here, the constants $1/n$ and $2/n$ for
$\xi^{(0)}$ and $\lambda^{(0)}$ has been used to avoid the possible
singularity problems. A suitable initial value for $\gamma$ is obtained by
considering the relation
$\gamma=\delta\log\bigl{(}1/{\boldsymbol{y}}_{0.5}-1\bigr{)}$ where
${\boldsymbol{y}}_{0.5}$ is median of the transformation
${\boldsymbol{y}}=({\boldsymbol{x}}-\xi)/\lambda$ with
${\boldsymbol{x}}=\bigl{(}x_{1},\dots,x_{n}\bigr{)}^{T}$ (Özçelik et al.,
2016). Therefore,
$\gamma^{(0)}=\delta^{(0)}\log\bigl{(}1/{\boldsymbol{y}}_{0.5}-1\bigr{)}$
where ${\boldsymbol{y}}_{0.5}$ is median of the transformation
${\boldsymbol{y}}=\bigl{(}{\boldsymbol{x}}-\xi^{(0)}\bigr{)}/\lambda^{(0)}$.
It should be noted that we took $\delta^{(0)}=1$ as the initial value for
$\delta$.

Additionally, our study revealed that the Bayesian paradigm presented in this
work is robust with respect to $\delta^{(0)}$, so that it can be started well
away from the true value of $\delta$. In the case of the three-parameter
Weibull distribution, similar to the JSB distribution, we used
$\xi^{(0)}=x_{(1)}-1/n$ as the initial value for $\mu$. The initial values of
the shape and scale parameters obtained by using the method of moments (Norman
et al., 1994). We also performed a simulation study to check the robustness of
the Bayesian paradigm with respect to the initial values for estimating the
parameters of the JSB distribution. For this purpose, we confine ourselves to
the case in which we have simulated 300 samples each of size 100 from JSB
distribution with parameter vector $\Theta=(2,2,20,0)^{T}$, i.e., $\delta=2$,
$\gamma=2$, $\lambda=20$, and $\xi=0$. The results of simulation are given in
Table 3. We note that in each of 200 runs, the initial values were not chosen
by the method suggested above. Instead, the initial values were generated
randomly from uniform distribution. We used this scenario in order to check
the robustness of the Gibbs sampler. The initial values for $\delta$,
$\gamma$, $\lambda$ ,and $\xi$ were generated from uniform distribution
(0.1,15), (-15,15), (20.1,60), and (-10,10), respectively. For example, the
general motion of the Gibbs sampler has been shown in Figure 8, when the
initial values were chosen as $\delta^{(0)}=15$, $\gamma^{(0)}=-15$,
$\lambda^{(0)}=60$, and $\xi^{(0)}=-10$ to show the robustness of the Bayesian
paradigm.

##  6 Conclusion

We have derived the Bayesian estimators for the four-parameter Johnson’s SB
(JSB) distribution. The maximum likelihood (ML) approach is the most commonly
used method for estimating the model parameters, but it has been shown using
simulation study that percentage of failed attempts to reach the convergence
through this method is on the average 32%. So, we suggest to use the Bayesian
paradigm to estimate the parameters of the JSB distribution. We have shown
that the proposed Bayesian approach works efficiently and is robust with
respect to the initial values. Furthermore, we have considered the Bayesian
estimators for parameters of the three-parameter Weibull distribution that
proposed by Green et al. (1994). Our algorithm for sampling from full
conditional of the location parameter is faster to that of Green et al.
(1994). The We have fitted both of the JSB and three-parameter Weibull
distributions to the diameters at breast height (DBH) obtained from 3 plots
out of 107 plots of size 0.08 ha established in mixed-age ponderosa pine
(Pinus ponderosa Dougl. ex Laws.) forests with scattered western junipers
located in the Malheur National Forest on the south end of the Blue Mountains
near Burns, Oregon, USA. The estimation results indicated that the JBS model
outperformed the three-parameter Weibull distribution and so characterized
more accurately the DBH distribution. As a possible future work, we are
interested in estimating the parameters of the bivariate JSB distribution
using the Bayesian method. The users can access the `R` package `ForestFit`
that is available at address
`https://cran.r-project.org/web/packages/ForestFit/index.html`.

## Appendix A

We have

$\displaystyle\pi(\delta|\gamma,\lambda,\xi,\boldsymbol{x})$
$\displaystyle\propto\frac{\delta^{n}\lambda^{n}}{(2\pi)^{\frac{n}{2}}\Pi_{i=1}^{n}(x_{i}-\xi)(\lambda+\xi-
x_{i})}\exp\left\\{-\frac{1}{2}\sum_{i=1}^{n}\Bigg{[}\gamma+\delta\log\biggl{(}\frac{x_{i}-\xi}{\lambda+\xi-
x_{i}}\biggr{)}\Bigg{]}^{2}\right\\}$
$\displaystyle\propto\delta^{n}\exp\left\\{-\frac{k_{2}}{2}\biggl{[}\delta+\frac{\gamma
k_{1}}{k_{2}}\biggr{]}^{2}\right\\},$

where

$\displaystyle k_{1}=\sum_{i=1}^{n}\log\Bigl{(}\frac{x_{i}-\xi}{\lambda+\xi-
x_{i}}\Bigr{)},$ (11)

and

$\displaystyle
k_{2}=\sum_{i=1}^{n}\bigg{[}\log\Bigl{(}\frac{x_{i}-\xi}{\lambda+\xi-
x_{i}}\Bigr{)}\biggr{]}^{2}.$ (12)

The full conditional pdf of $\delta$ is given by

$\displaystyle\pi(\delta|\gamma,\lambda,\xi,\boldsymbol{x})=\text{C}\delta^{n}\exp\left\\{-\frac{k_{2}}{2}\biggl{[}\delta+\frac{\gamma
k_{1}}{k_{2}}\biggr{]}^{2}\right\\},$

where C is a normalizing constant independent of $\delta$. The First and
second derivatives of $\log\pi(\delta|\gamma,\lambda,\xi,\boldsymbol{x})$ with
respect to $\delta$ are

$\displaystyle\frac{\partial\pi(\delta|\gamma,\lambda,\xi,\boldsymbol{x})}{\partial\delta}=\frac{n}{\delta}-k_{2}\biggl{[}\delta+\frac{\gamma
k_{1}}{k_{2}}\biggr{]},$

and

$\displaystyle\frac{\partial^{2}\pi(\delta|\gamma,\lambda,\xi,\boldsymbol{x})}{\partial\delta^{2}}=-\frac{n}{\delta^{2}}-k_{2}.$
(13)

Since $k_{2}>0$, the right-hand side of (13) is always negative and so
$\pi(\delta|\gamma,\lambda,\xi,\boldsymbol{x})$ is log-concave. Assume that we
are currently at $t$-th iteration of the sampler, sampling from
$\pi\bigl{(}\delta\big{|}\gamma^{(t)},\lambda^{(t)},\xi^{(t)},\boldsymbol{x}\bigr{)}$
is carried out through the ARS algorithm. Here, $\gamma^{(t)}$,
$\lambda^{(t)}$, and $\xi^{(t)}$ denote the generated values, respectively,
from $\gamma$, $\lambda$, and $\xi$, at $t$-th iteration.

## Appendix B

We have

$\displaystyle\pi(\gamma|\delta,\lambda,\xi,\boldsymbol{x})\propto\exp\left\\{-\frac{1}{2}\sum_{i=1}^{n}\Bigg{[}\gamma+\delta\log\biggl{(}\frac{x_{i}-\xi}{\lambda+\xi-
x_{i}}\biggr{)}\Bigg{]}^{2}\right\\}\propto\exp\left\\{-\frac{n}{2}\biggl{[}\gamma+\frac{\delta
k_{1}}{n}\biggr{]}^{2}\right\\},$

where $k_{1}$ is defined in (11). Assume that we are currently at $t$-th
iteration of the sampler, for sampling from full conditional pdf of $\gamma$,
it is enough to sample from Gaussian distribution with mean

$\displaystyle-\frac{\delta^{(t+1)}}{n}\sum_{i=1}^{n}\log\biggl{(}\frac{x_{i}-\xi^{(t)}}{\lambda^{(t)}+\xi^{(t)}-x_{i}}\biggr{)},$

and variance $1/n$.

## Appendix C

It is easy to see that the full conditional of $\lambda$ is (up to
proportionality)

$\displaystyle\pi(\lambda|\delta,\gamma,\xi,\boldsymbol{x})$
$\displaystyle\propto\Pi_{i=1}^{n}\biggl{(}\frac{\lambda}{\lambda+\xi-
x_{i}}\biggr{)}\exp\left\\{-\frac{1}{2}\sum_{i=1}^{n}\Bigg{[}\gamma+\delta\log\biggl{(}\frac{x_{i}-\xi}{\lambda+\xi-
x_{i}}\biggr{)}\Bigg{]}^{2}\right\\}$
$\displaystyle=\lambda^{n}\Pi_{i=1}^{n}\Bigl{(}\lambda+\xi-
x_{i}\Bigr{)}^{\delta\gamma-1}\exp\left\\{-\frac{\delta^{2}}{2}\sum_{i=1}^{n}\Bigg{[}\log\biggl{(}\frac{x_{i}-\xi}{\lambda+\xi-
x_{i}}\biggr{)}\Bigg{]}^{2}\right\\},$

As it is seen, the structure of the full conditional pdf of $\lambda$ is
complicated. Additionally, $\pi(\lambda|\delta,\gamma,\xi,\boldsymbol{x})$ is
not always log-concave. So, we use the MH algorithm for sampling from
$\pi(\lambda|\delta,\gamma,\xi,\boldsymbol{x})$ by choosing
$\exp\Bigl{\\{}-\bigl{(}\lambda-x_{(n)}+\xi\bigr{)}\Bigr{\\}}$ as the proposal
pdf for $\lambda>x_{(n)}-\xi$ wherein $x_{(n)}=\max\\{x_{1},\dots,x_{n}\\}$.
The six-step MH algorithm is given a follows.

  1. 1.

Suppose we are currently at $t$-th iteration of the sampler. Choose the
initial value as $\lambda_{(0)}=x_{(n)}-\xi^{(t)}+1/n$ and set i=1;

  2. 2.

Sample $\lambda^{*}$ from proposal distribution $q\bigl{(}\lambda)$ with pdf
$\exp\Bigl{\\{}-\bigl{(}\lambda-x_{(n)}+\xi^{(t)}\bigr{)}\Bigr{\\}}$, for
$\lambda>x_{(n)}-\xi^{(t)}$;

  3. 3.

Compute $\eta$ as

$\eta=\min\left\\{1,\frac{\pi\bigl{(}\lambda_{*}\big{|}\delta^{(t+1)},\gamma^{(t+1)},\xi^{(t)},\boldsymbol{x}\bigr{)}\exp\bigl{\\{}-\lambda_{(i-1)}\bigr{\\}}}{\pi\bigl{(}\lambda_{(i-1)}\big{|}\delta^{(t+1)},\gamma^{(t+1)},\xi^{(t)},\boldsymbol{x}\bigr{)}\exp\bigl{\\{}-\lambda_{*}\bigr{\\}}}\right\\}.$

  4. 4.

Generate an uniform random variable on $(0,1)$, say $u$. If $u<\eta$, then
$\lambda_{(i)}=\lambda_{*}$, otherwise $\lambda_{(i)}=\lambda_{(i-1)}$;

  5. 5.

If $i=N$, then go to the next step. Otherwise set $i=i+1$ and go to step 2;

  6. 6.

Accept $\lambda_{(N)}$ as a generation form
$\pi\bigl{(}\lambda\big{|}\delta^{(t)},\gamma^{(t)},\xi^{(t)},\boldsymbol{x}\bigr{)}$,
i.e., $\lambda^{(t+1)}=\lambda_{(N)}$ and stop the MH algorithm.

## Appendix D

Similar to the $\lambda$, the full conditional pdf of $\xi$ has complicated
structure. We have

$\displaystyle\pi(\xi|\delta,\gamma,\lambda,\boldsymbol{x})$
$\displaystyle\propto\Pi_{i=1}^{n}\Biggl{(}\frac{(\lambda+\xi-
x_{i})^{\delta\gamma-1}}{(x_{i}-\xi)^{\delta\gamma+1}}\Biggr{)}\exp\left\\{-\frac{\delta^{2}}{2}\sum_{i=1}^{n}\Bigg{[}\log\biggl{(}\frac{x_{i}-\xi}{\lambda+\xi-
x_{i}}\biggr{)}\Bigg{]}^{2}\right\\},$

where $x_{(n)}-\lambda<\xi<x_{(1)}$ in which
$x_{(1)}=\min\\{x_{1},\dots,x_{n}\\}$. Since the full conditional pdf of $\xi$
is not always log-concave, so we use the MH algorithm to sample from it. For
this aim, we choose the uniform distribution on
$\bigl{(}x_{(n)}-\lambda,x_{(1)}\bigr{)}$ as the proposal. The six-step MH
algorithm is given by the following.

  1. 1.

Suppose we are currently at $t$-th iteration of the sampler. Set $i=1$ and
generate a random variable from uniform distribution on
$\bigl{(}x_{(n)}-\lambda^{(t+1)},x_{(1)}\bigr{)}$, say $u$. Set $u=\xi^{(0)}$;

  2. 2.

Sample $\xi^{*}$ from proposal distribution with pdf
$q\bigl{(}\xi)=1/\bigl{(}x_{(1)}-x_{(n)}+\lambda^{(t+1)}\bigr{)}$;

  3. 3.

Compute $\eta$ as

$\eta=\min\left\\{1,\frac{\pi\bigl{(}\xi^{*}\big{|}\delta^{(t+1)},\gamma^{(t+1)},\lambda^{(t+1)},\boldsymbol{x}\bigr{)}}{\pi\bigl{(}\xi^{(i-1)}\big{|}\delta^{(t+1)},\gamma^{(t+1)},\lambda^{(t+1)},\boldsymbol{x}\bigr{)}}\right\\}.$

  4. 4.

Generate an uniform random variable, say $u$, on $(0,1)$. If $u<\eta$, then
$\xi^{(i)}=\xi^{*}$, otherwise $\xi^{(i)}=\xi^{(i-1)}$;

  5. 5.

If $i=N$, then go to the next step. Otherwise set $i=i+1$ and go to step 2;

  6. 6.

Accept $\xi_{(N)}$ as a generation form pdf
$\pi\Bigl{(}\xi\big{|}\delta^{(t+1)},\gamma^{(t+1)},\lambda^{(t+1)},\boldsymbol{x}\Bigr{)}$,
i.e., $\xi^{(t+1)}=\xi_{(N)}$ and stop the MH algorithm.

## References

  * Bailey and Dell (1973) Bailey, R.L., Dell, T., 1973\.  Quantifying diameter distributions with the Weibull function.  Forest Science 19, 97–104.
  * Fonseca et al. (2009) Fonseca, T.F., Marques, C.P., Parresol, B.R., 2009.  Describing maritime pine diameter distributions with Johnson’s SB distribution using a new all-parameter recovery approach.  Forest Science 55, 367–373.
  * Gilks and Wild (1992) Gilks, W.R., Wild, P., 1992\.  Adaptive rejection sampling for Gibbs sampling.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 41, 337–348.
  * Gorgoso et al. (2007) Gorgoso, J., González, J.Á., Rojo, A., Grandas-Arias, J., 2007\.  Modelling diameter distributions of betula alba l. stands in northwest Spain with the two-parameter Weibull function.  Forest Systems 16, 113–123.
  * Gorgoso et al. (2012) Gorgoso, J., Rojo, A., Cámara-Obregón, A., Diéguez-Aranda, U., 2012.  A comparison of estimation methods for fitting Weibull, Johnson’s SB and beta functions to Pinus pinaster, Pinus radiate and Pinus sylvestris stands in northwest Spain.  Forest Systems 21, 446–459.
  * Gove et al. (2008) Gove, J.H., Ducey, M.J., Leak, W.B., Zhang, L., 2008\.  Rotated sigmoid structures in managed uneven-aged northern hardwood stands: a look at the Burr type III distribution.  Forestry 81, 161–176.
  * Green et al. (1994) Green, E.J., Roesch Jr, F.A., Smith, A.F., Strawderman, W.E., 1994\.  Bayesian estimation for the three-parameter Weibull distribution with tree diameter data.  Biometrics , 254–269.
  * Hafley and Buford (1985) Hafley, W., Buford, M., 1985\.  A bivariate model for growth and yield prediction.  Forest Science 31, 237–247.
  * Hafley and Schreuder (1977) Hafley, W., Schreuder, H., 1977\.  Statistical distributions for fitting diameter and height data in even-aged stands.  Canadian Journal of Forest Research 7, 481–487.
  * Jeffreys (1961) Jeffreys, H., 1961.  Theory of Probability.  Oxford University Press London.
  * Johnson (1949) Johnson, N.L., 1949.  Systems of frequency curves generated by methods of translation.  Biometrika 36, 149–176.
  * Kerns et al. (2017) Kerns, B.K., Westlind, D.J., Day, M.A., 2017.  Season and interval of burning and cattle exclusion in the southern Blue Mountains, Oregon: Overstory tree height, diameter and growth.  USDA Forest Service Research Data Archive, https://doi.org/10.2737/RDS-2017-0041 .
  * Kiviste et al. (2003) Kiviste, A., Nilson, A., Hordo, M., Merenäkk, M., 2003\.  Diameter distribution models and height-diameter equations for Estonian forests.  Modelling Forest Systems , 169–179.
  * Kudus et al. (1999) Kudus, K.A., Ahmad, M., Lapongan, J., 1999.  Nonlinear regression approach to estimating Johnson SB parameters for diameter data.  Canadian Journal of Forest Research 29, 310–314.
  * Maltamo et al. (2000) Maltamo, M., Kangas, A., Uuttera, J., Torniainen, T., Saramäki, J., 2000.  Comparison of percentile based prediction methods and the Weibull distribution in describing the diameter distribution of heterogeneous Scots pine stands.  Forest Ecology and Management 133, 263–274.
  * Maltamo et al. (1995) Maltamo, M., Puumalainen, J., Päivinen, R., 1995.  Comparison of beta and Weibull functions for modelling basal area diameter distribution in stands of pinus sylvestris and picea abies.  Scandinavian Journal of Forest Research 10, 284–295.
  * Marto et al. (2009) Marto, M., Palma, J., Mateus, A., Tomé, M., 2009\.  Computer program for estimation of Johnson’s SB parameters using a parameter recovery method.  Publicações Científicas Forchange PC-X/2009. Centro de Estudos Florestais, Instituto Superior de Agronomia, Universidade Técnica de Lisboa, Lisboa .
  * Mateus and Tomé (2011) Mateus, A., Tomé, M., 2011\.  Modelling the diameter distribution of eucalyptus plantations with Johnson’s SB probability density function: parameters recovery from a compatible system of equations to predict stand variables.  Annals of Forest Science 68, 325–335.
  * Nagatsuka et al. (2013) Nagatsuka, H., Kamakura, T., Balakrishnan, N., 2013.  A consistent method of estimation for the three-parameter Weibull distribution.  Computational Statistics & Data Analysis 58, 210–226.
  * Norman et al. (1994) Norman, L.J., Kotz, S., Balakrishnan, N., 1994.  Continuous Univariate Distributions. volume I.  John Wiley & Sons.
  * Özçelik et al. (2016) Özçelik, R., Fidalgo Fonseca, T.J., Parresol, B.R., Eler, Ü., 2016\.  Modeling the diameter distributions of Brutian Pine stands using Johnson’s SB distribution.  Forest Science 62, 587–593.
  * Palahí et al. (2007) Palahí, M., Pukkala, T., Blasco, E., Trasobares, A., 2007\.  Comparison of beta, Johnson’s SB, Weibull and truncated Weibull functions for modeling the diameter distribution of forest stands in Catalonia (north-east of Spain).  European Journal of Forest Research 126, 563–571.
  * Parresol (2003) Parresol, B.R., 2003.  Recovering parameters of Johnson’s SB distribution.  Res. Pap. SRS-31. Asheville, NC: US Department of Agriculture, Forest Service, Southern Research Station. 9 p. 31\.
  * Pretzsch (2009) Pretzsch, H., 2009.  Forest dynamics, growth, and yield, in: Forest dynamics, growth and yield. Springer, pp. 1–39.
  * R Core Team (2018) R Core Team, 2018.  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing. Vienna, Austria.  URL: https://www.R-project.org/.
  * Rennolls and Wang (2005) Rennolls, K., Wang, M., 2005\.  A new parameterization of Johnson’s SB distribution with application to fitting forest tree diameter data.  Canadian Journal of Forest Research 35, 575–579.
  * Robert et al. (2010) Robert, C.P., Casella, G., Casella, G., 2010.  Introducing Monte Carlo methods with R. volume 18.  Springer.
  * Roberts and Smith (1994) Roberts, G.O., Smith, A.F., 1994\.  Simple conditions for the convergence of the Gibbs sampler and Metropolis-Hastings algorithms.  Stochastic Processes and Their Applications 49, 207–216.
  * Scolforo et al. (2003) Scolforo, J.R.S., Tabai, F.C.V., de Macedo, R.L.G., Acerbi Jr, F.W., de Assis, A.L., 2003.  SB distribution’s accuracy to represent the diameter distribution of Pinus taeda, through five fitting methods.  Forest Ecology and Management 175, 489–496.
  * Smith and Naylor (1987) Smith, R.L., Naylor, J., 1987\.  A comparison of maximum likelihood and Bayesian estimators for the three-parameter weibull distribution.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 36, 358–369.
  * Teimouri et al. (2020) Teimouri, M., Doser, J.W., Finley, A.O., 2020.  Forestfit: An R package for modeling plant size distributions.  Environmental Modelling & Software , 104668\.
  * Zhang et al. (2003) Zhang, L., Packard, K.C., Liu, C., 2003.  A comparison of estimation methods for fitting Weibull and Johnson’s SB distributions to mixed spruce fir stands in northeastern North America.  Canadian Journal of Forest Research 33, 1340–1347.
  * Zhang et al. (2010) Zhang, X., Lei, Y., Cao, Q.V., 2010.  Compatibility of stand basal area predictions based on forecast combination.  Forest Science 56, 552–557.
  * Zhou and McTague (1996) Zhou, B., McTague, J.P., 1996\.  Comparison and evaluation of five methods of estimation of the Johnson system parameters.  Canadian Journal of Forest Research 26, 928–935.

Table 1: Percentage of runs that NR method truly converged for the JSB distribution. Sample size | 20 | 50 | 100 | 250 | 500 | 1000 | 5000
---|---|---|---|---|---|---|---
Percentage | 68.5% | 68.7% | 68.3% | 68.3% | 68.1% | 68.4% | 67.9%
Table 2: Summary statistics for DBH data. Plot | Plot size | Min | 1st Quar. | Median | Mean | 3rd Quar. | Max. | St.Dev. | Skewness
---|---|---|---|---|---|---|---|---|---
9 | 52 | 10.4 | 17.32 | 27.55 | 28.01 | 36.42 | 55.9 | 12.37 | 0.48
44 | 42 | 11.9 | 14.78 | 17.40 | 24.27 | 26.50 | 83.8 | 15.06 | 2.13
73 | 35 | 9.1 | 13.85 | 21.10 | 30.35 | 39.75 | 88.6 | 20.94 | 1.15
Table 3: Descriptive statistics for the sampler output. Parameter | Min | 1st Quar. | Median | Mean | 3rd Quar. | Max. | St.Dev. | Skewness
---|---|---|---|---|---|---|---|---
$\delta$ | 1.242 | 1.724 | 1.850 | 1.874 | 1.996 | 15.000 | 0.3406 | 20.123
$\gamma$ | -15.000 | 1.407 | 1.615 | 1.609 | 1.804 | 15.982 | 0.450 | 7.225
$\lambda$ | 12.399 | 15.941 | 17.215 | 17.155 | 18.410 | 60.000 | 1.719 | 1.469
$\xi$ | -10.000 | -0.1270 | 0.245 | 0.160 | 0.549 | 1.399 | 0.542 | -1.460
Table 4: Bayesian estimators for the parameters of the JBS and three-parameter Weibull distributions. | JSB distribution |  | Weibull distribution
---|---|---|---
Plot | ${\delta}$ | ${\gamma}$ | ${\lambda}$ | $\xi$ |  | $\alpha$ | $\beta$ | $\mu$
9 | 0.772 | 0.545 | 52.311 | 8.719 |  | 1.682 | 23.120 | 7.436
44 | 0.875 | 0.203 | 64.162 | 3.642 |  | 2.145 | 34.496 | 2.171
73 | 0.641 | 0.978 | 88.592 | 8.182 |  | 1.005 | 22.746 | 8.278
Table 5: Computed goodness-of-fit statistics for modelling DBH data using JBS and three-parameter Weibull distributions. | JSB distribution |  | Weibull distribution
---|---|---|---
Plot | AD | CM | KS | LL |  | AD | CM | KS | LL
9 | 0.167 | 0.023 | 0.059 | -196.808 |  | 0.310 | 0.048 | 0.082 | -199.579
44 | 0.216 | 0.028 | 0.100 | -120.072 |  | 0.302 | 0.040 | 0.115 | -121.846
73 | 0.105 | 0.031 | 0.075 | -141.751 |  | 0.257 | 0.045 | 0.084 | -143.299

Figure 1: Left-hand side: pairwise scatterplots of trimmed output of the Gibbs
sampler for estimation parameters of the JSB distribution fitted to the DBH
observations in plot 9. These outputs suggest that there is little dependence
between $\delta$ and $\lambda$. Right-hand side: histograms of the full
conditionals (a) $\pi(\delta|\gamma,\lambda,\xi,\boldsymbol{x})$, (b)
$\pi(\gamma|\delta,\lambda,\xi,\boldsymbol{x})$, (c)
$\pi(\lambda|\delta,\gamma,\xi,\boldsymbol{x})$, and (d)
$\pi(\xi|\delta,\gamma,\lambda,\boldsymbol{x})$ produced by the Gibbs sampler
for 10,000 runs.

Figure 2: Left-hand side: pairwise scatterplots of trimmed output of the Gibbs
sampler for estimation parameters of the three-parameter Weibull distribution
fitted to the DBH observations in plot 9. These outputs suggest that there is
considerable dependence between all pairs. Right-hand side: histograms of the
full conditionals (a) $\pi(\alpha|\beta,\mu,\boldsymbol{x})$, (b)
$\pi(\beta|\alpha,\mu,\boldsymbol{x})$, and (c)
$\pi(\mu|\alpha,\beta,\boldsymbol{x})$ produced by the Gibbs sampler for
10,000 runs.

Figure 3: Left-hand side: pairwise scatterplots of trimmed output of the Gibbs
sampler for estimation parameters of the JSB distribution fitted to the DBH
observations in plot 44. These outputs suggest that there is little dependence
between between pairs $(\delta,\lambda)$ and $(\lambda,\xi)$. Right-hand side:
histograms of the full conditionals (a)
$\pi(\delta|\gamma,\lambda,\xi,\boldsymbol{x})$, (b)
$\pi(\gamma|\delta,\lambda,\xi,\boldsymbol{x})$, (c)
$\pi(\lambda|\delta,\gamma,\xi,\boldsymbol{x})$, and (d)
$\pi(\xi|\delta,\gamma,\lambda,\boldsymbol{x})$ produced by the Gibbs sampler
for 10,000 runs.

Figure 4: Left-hand side: pairwise scatterplots of trimmed output of the Gibbs
sampler for estimation parameters of the three-parameter Weibull distribution
fitted to the DBH observations in plot 44. These outputs suggest that there is
considerable dependence between all pairs. Right-hand side: histograms of the
full conditionals (a) $\pi(\alpha|\beta,\mu,\boldsymbol{x})$, (b)
$\pi(\beta|\alpha,\mu,\boldsymbol{x})$, and (c)
$\pi(\mu|\alpha,\beta,\boldsymbol{x})$ produced by the Gibbs sampler for
10,000 runs.

Figure 5: Left-hand side: pairwise scatterplots of trimmed output of the Gibbs
sampler for estimation parameters of the JSB distribution fitted to the DBH
observations in plot 73. These outputs suggest that there is little dependence
between $\delta$ and $\gamma$. Right-hand side: histograms of the full
conditionals (a) $\pi(\delta|\gamma,\lambda,\xi,\boldsymbol{x})$, (b)
$\pi(\gamma|\delta,\lambda,\xi,\boldsymbol{x})$, (c)
$\pi(\lambda|\delta,\gamma,\xi,\boldsymbol{x})$, and (d)
$\pi(\xi|\delta,\gamma,\lambda,\boldsymbol{x})$ produced by the Gibbs sampler
for 10,000 runs.

Figure 6: Left-hand side: pairwise scatterplots of trimmed output of the Gibbs
sampler for estimation parameters of the three-parameter Weibull distribution
fitted to the DBH observations in plot 73. These outputs suggest that there is
little dependence between $\alpha$ and $\mu$. Right-hand side: histograms of
the full conditionals (a) $\pi(\alpha|\beta,\mu,\boldsymbol{x})$, (b)
$\pi(\beta|\alpha,\mu,\boldsymbol{x})$, and (c)
$\pi(\mu|\alpha,\beta,\boldsymbol{x})$ produced by the Gibbs sampler for
10,000 runs.

Figure 7: Histograms of DBH data in plots 9, 44, and 73. Superimposed in each
subfigure are estimated probability density functions of the JSB (blue solid
line) and Weibull (red dashed line) distributions. Figure 8: Output of (a)
$\delta$, (b) $\gamma$, (c) $\lambda$, and (d) $\xi$ from Gibbs sampler for
the robust analysis. Each subfigure shows the motion of the Gibbs sampler
across 10,000 iterations.


multicitedelim, compcitedelim, postnotedelim mydate mydate 210mm5.5in* [] []
*26pc1.618**1*** 2.5 *2* 2ex10mm0pt [nearest] 0pt manaartplain manaart0.5
manaartPorta Manalog-likelihood and cross-validation manaart0 plain0

This document is designed for screen reading and two-up printing on A4 or
Letter paper

plain figure * #1

# A relation between log-likelihood
and cross-validation log-scores

P.G.L. Porta Mana Kavli Institute, Trondheim, Norway <ntnu.no> piero.mana

(18 August 2019; updated 20 August 2019)

It is shown that the log-likelihood of a hypothesis or model given some data
is equal to an average of all leave-one-out cross-validation log-scores that
can be calculated from all subsets of the data. This relation can be
generalized to any $k$-fold cross-validation log-scores.
_Note: Dear Reader & Peer, this manuscript is being peer-reviewed by you.
Thank you._

##  1 Log-likelihoods and cross-validation log-scores

The probability calculus unequivocally tells us how our degree of belief in a
hypothesis $H_{h}$ given data $D$ and background information or assumptions
$I$, that is, $\mathrm{P}(H_{h}\nonscript\,|\nonscript\;\mathopen{}D\,I)$, is
related to our degree of belief in observing those data when we entertain that
hypothesis as true, that is,
$\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H_{h}\,I)$:

$\displaystyle\mathrm{P}(H_{h}\nonscript\,|\nonscript\;\mathopen{}D\,I)$
$\displaystyle=\frac{\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H_{h}\,I)\;\mathrm{P}(H_{h}\nonscript\,|\nonscript\;\mathopen{}I)}{\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}I)}$
(1a)
$\displaystyle=\frac{\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H_{h}\,I)\;\mathrm{P}(H_{h}\nonscript\,|\nonscript\;\mathopen{}I)}{\sumop\slimits@_{h^{\prime}}\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H_{h^{\prime}}\,I)\;\mathrm{P}(H_{h^{\prime}}\nonscript\,|\nonscript\;\mathopen{}I)}.$
(1b)

$D$, $H_{h}$, $I$ denote propositions, which are usually about numeric
quantities. I use the terms ‘degree of belief’, ‘belief’, and ‘probability’ as
synonyms. By ‘hypothesis’ I mean either a scientific (physical, biological,
etc.) hypothesis – a state or development of things capable of experimental
verification, at least in a thought experiment – or more generally some
proposition, often not precisely specified, which leads to quantitatively
specific distributions of beliefs for any contemplated data set. In the latter
case we often call $H_{h}$ a ‘(probabilistic) model’ rather than a
‘hypothesis’.

Expression (1b) assumes that we have a set $\\{H_{h}\\}$ of mutually exclusive
and exhaustive hypotheses under consideration, which is implicit in our
knowledge $I$. In fact it’s only valid if

$\mathrm{P}\bigl{(}\mathop{\textstyle\bigvee}\nolimits_{h}H_{h}\nonscript\,|\nonscript\;\mathopen{}I\bigr{)}=1,\qquad\mathrm{P}(H_{h}\land
H_{h^{\prime}}\nonscript\,|\nonscript\;\mathopen{}I)=0\quad\text{if $h\neq
h^{\prime}$}.$ (2)

Only rarely does the set of hypotheses $\\{H_{h}\\}$ encompass and reflect the
extremely complex and fuzzy hypotheses lying in the backs of our minds.
They’re simplified pictures. That’s also why they’re called ‘models’.

Expression (1a) is universally valid instead, but it’s rarely possible to
quantify its denominator $\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}I)$
unless we simplify our inferential problem by introducing a possibly
unrealistic exhaustive set of hypotheses, thus falling back to (1b). We can
bypass this problem if we are content with comparing our beliefs about any two
hypotheses through their ratio, so that the term
$\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}I)$ cancels out. See Jaynes’s
[§§ 4.3–4.4]jaynes1994_r2003 insightful remarks about such binary comparisons,
and also Good’s [§ 6.3–6.6]good1950.

The term $\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H_{h}\,I)$ in eq.
(1) is called the _likelihood_ of the hypothesis given the data [good1950, §
6.1 p. 62]. Its logarithm is surprisingly called log-likelihood:

$\log\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H_{h}\,I),$ (3)

where the logarithm can be taken in an arbitrary basis (Turing, Good
[good1985, good1950, good1969, e.g.], Jaynes [jaynes1994_r2003, § 4.2]
recommend base $\cramped{1}0^{1/10}$, leading to a measurement in decibels;
see the cited works for the practical advantages of such choice).

The ratio of the likelihoods of two hypotheses, called _relative Bayes factor_
, or its logarithm, the _relative weight of evidence_ ,111[ch.
6]good1950good1975,good1981,good1985, and many other works in [good1983]; [§
1.4]osteyeeetal1974mackay1992,kassetal1995; see also [jeffreys1939_r1983, chs
V, VI, A]. are often used to quantify how much the data favour our belief in
one versus the other hypothesis (that is, assuming at least momentarily that
they be exhaustive). ‘It is historically interesting that the expression
“weight of evidence”, in its technical sense, anticipated the term
“likelihood” by over forty years’ [osteyeeetal1974, § 1.4.2 p. 12].

######

Recent literature [kassetal1995, for example] seems to exclusively deal with
_relative_ Bayes factors. I’d like to recall, lest it fades from the memory,
the definition of the non-relative Bayes factor for a hypothesis $H_{h}$
provided by data $D$:[good1981, § 2]

$\frac{\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H_{h}\;I)}{\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}\lnot
H_{h}\;I)}\equiv\frac{\mathrm{O}(H_{h}\nonscript\,|\nonscript\;\mathopen{}D\;I)}{\mathrm{O}(H_{h}\nonscript\,|\nonscript\;\mathopen{}I)}=\frac{\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H_{h}\;I)\;[1-\mathrm{P}(H_{h}\nonscript\,|\nonscript\;\mathopen{}I)]}{\sumop\slimits@_{h^{\prime}}^{h^{\prime}\neq
h}\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H_{h^{\prime}}\;I)\;\mathrm{P}(H_{h^{\prime}}\nonscript\,|\nonscript\;\mathopen{}I)},$
(4)

where the _odds_ $\mathrm{O}$ is defined as
$\mathrm{O}\coloneqq\mathrm{P}/(1-\mathrm{P})$. Looking at the expression on
the right, which can be derived from the probability rules, it’s clear that
the Bayes factor for a hypothesis involves the likelihoods of _all_ other
hypotheses as well as their pre-data probabilities. This quantity and its
logarithm, the (non-relative) weight of evidence, have important properties
which relative Bayes factors and relative weights of evidence don’t enjoy. For
example, the expected weight of evidence for a correct hypothesis is always
positive, and for a wrong hypotheses always negative[good1950, § 6.7]. See
Jaynes[§§ 4.3–4.4]jaynes1994_r2003 for further discussion and a numeric
example.

The literature in probability and statistics has also employed and debated
other ad-hoc measures to quantify how the data relate to the hypotheses – or
even to select one hypothesis for further use, discarding the
others[bernardoetal1994_r2000, §§ 3.4, 6.1.6 gives the clearest motivation and
explanation][see
also][]stone1977,geisseretal1979,vehtarietal2012,vehtarietal2002,krnjajicetal2011,krnjajicetal2014,gelmanetal2014,gronauetal2019,chandramoulietal2019.
Here I consider one measure in particular: the _leave-one-out cross-validation
log-score_ 11footnotemark: 1, which I’ll just call ‘log-score’ for brevity:

$\frac{1}{d}\sumop\slimits@_{i=1}^{d}\log\mathrm{P}(D_{i}\nonscript\,|\nonscript\;\mathopen{}D_{-i}\,H_{h}\,I)$
(5)

where every $D_{i}$ is one datum in the data
$D\equiv\bigwedge_{i=1}^{d}D_{i}$, and $D_{-i}$ denotes the data with datum
$D_{i}$ excluded. The intuition behind this score can be colloquially
expressed thus: ‘let’s see what my belief in one datum would be, on average,
once I’ve observed the other data, if I consider $H_{h}$ as true’. ‘On
average’ means considering such belief for every single datum in turn, and
then taking the geometric mean of the resulting beliefs. Other variants of
this score use more general partitions of the data into two disjoint
subsets11footnotemark: 1.

If you find this you can claim a postcard from me.

My purpose is to show an exact relation between the log-likelihood (3) and the
leave-one-out cross-validation log-score (5). This relation doesn’t seem to
appear in the literature, and I find it very intriguing because it portrays
the log-likelihood as a sort of full-scale use of the log-score: it says that
_the log-likelihood is the sum of all averaged log-scores that can be formed
from all data subsets_. The relation can be extended to more general cross-
validation log-scores, and it can be of interest for the debate about the
soundness of log-scores in deciding among hypotheses.

##  2 A relation between log-likelihood and log-score

We can obviously write the likelihood as the $d$th root of its $d$th power:

$\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H\,I)\equiv\bigl{[}\,\underbracket{\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H\,I)\times\dotsm\times\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H\,I)}_{\text{$d$
times}}\,\bigr{]}^{1/d}$ (6)

where we have dropped the subscript h for simplicity. By the rules of
probability we have

$\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H\,I)=\mathrm{P}(D_{i}\nonscript\,|\nonscript\;\mathopen{}D_{-i}\,H_{h}\,I)\times\mathrm{P}(D_{-i}\nonscript\,|\nonscript\;\mathopen{}H_{h}\,I)$
(7)

no matter which specific $i\in\\{1,\dotsc,d\\}$ we choose (temporal ordering
and similar matters are completely irrelevant in the formula above: it’s a
logical relation between propositions). So let’s expand each of the $d$
factors in the identity (6) using the product rule (7), using a different $i$
for each of them. The result can be thus displayed:

$\displaystyle\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H\,I)\equiv{}\bigl{[}$
$\displaystyle\mathrm{P}(D_{1}\nonscript\,|\nonscript\;\mathopen{}D_{-1}\,H\,I)\times\mathrm{P}(D_{-1}\nonscript\,|\nonscript\;\mathopen{}H\,I)\times{}$
(8)
$\displaystyle\mathrm{P}(D_{2}\nonscript\,|\nonscript\;\mathopen{}D_{-2}\,H\,I)\times\mathrm{P}(D_{-2}\nonscript\,|\nonscript\;\mathopen{}H\,I)\times{}$
$\displaystyle\hphantom{\mathrm{P}(D_{2}\nonscript\,|\nonscript\;\mathopen{}D_{-2}\,H\,I)}\mathrlap{\dotso}\hphantom{{}\times\mathrm{P}(D_{-2}\nonscript\,|\nonscript\;\mathopen{}H\,I)}\times{}$
$\displaystyle\mathrm{P}(D_{d}\nonscript\,|\nonscript\;\mathopen{}\underbracket{D_{-d}}_{\mathclap{\begin{subarray}{c}\Big{\uparrow}\\\
\text{this column leads to the log-
score}\end{subarray}}}H\,I)\times\mathrm{P}(D_{-d}\nonscript\,|\nonscript\;\mathopen{}H\,I)\,\bigr{]}^{1/d}.$

Upon taking the logarithm of this expression, the $d$ factors vertically
aligned on the left add up to the log-score (5), as indicated. But the
mathematical reshaping we just did for
$\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H\,I)$ – that is, the root-
product identity (6) and the expansion (8) – can be done for each of the
remaining factors $\mathrm{P}(D_{-i}\nonscript\,|\nonscript\;\mathopen{}H\,I)$
vertically aligned on the right in the expression above; and so on
recursively. Here is an explicit example for $d=3$:

$\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H\,I)\equiv{}\\\
\begin{aligned}
b]{2}\Bigl{\\{}\,&\mathrm{P}(D_{1}\nonscript\,|\nonscript\;\mathopen{}D_{2}\,D_{3}\,H\,I)\times\bigl{[}&&\mathrm{P}(D_{2}\nonscript\,|\nonscript\;\mathopen{}D_{3}\,H\,I)\times\mathrm{P}(D_{3}\nonscript\,|\nonscript\;\mathopen{}H\,I)\times{}\\\\[-3.0pt]
&&&\mathrm{P}(D_{3}\nonscript\,|\nonscript\;\mathopen{}D_{2}\,H\,I)\times\mathrm{P}(D_{2}\nonscript\,|\nonscript\;\mathopen{}H\,I)\bigr{]}^{1/2}\times{}\\\\[6.0pt]
&\mathrm{P}(D_{2}\nonscript\,|\nonscript\;\mathopen{}D_{1}\,D_{3}\,H\,I)\times\bigl{[}&&\mathrm{P}(D_{1}\nonscript\,|\nonscript\;\mathopen{}D_{3}\,H\,I)\times\mathrm{P}(D_{3}\nonscript\,|\nonscript\;\mathopen{}H\,I)\times{}\\\\[-3.0pt]
&&&\mathrm{P}(D_{3}\nonscript\,|\nonscript\;\mathopen{}D_{1}\,H\,I)\times\mathrm{P}(D_{1}\nonscript\,|\nonscript\;\mathopen{}H\,I)\bigr{]}^{1/2}\times{}\\\\[6.0pt]
&\mathrm{P}(D_{3}\nonscript\,|\nonscript\;\mathopen{}D_{1}\,D_{2}\,H\,I)\times\bigl{[}&&\mathrm{P}(D_{1}\nonscript\,|\nonscript\;\mathopen{}D_{2}\,H\,I)\times\mathrm{P}(D_{2}\nonscript\,|\nonscript\;\mathopen{}H\,I)\times{}\\\\[-3.0pt]
&&&\mathrm{P}(D_{2}\nonscript\,|\nonscript\;\mathopen{}D_{1}\,H\,I)\times\mathrm{P}(D_{1}\nonscript\,|\nonscript\;\mathopen{}H\,I)\bigr{]}^{1/2}\Bigr{\\}}^{1/3}.\end{aligned}$
(9)

In this example the logarithm of the three vertically aligned factors in the
left column is, as already noted, the log-score (5). The logarithm of the six
vertically aligned factors in the central column is an average of the log-
scores calculated for the three distinct subsets of pairs of data
$\\{D_{1}\,D_{2}\\}$, $\\{D_{1}\,D_{3}\\}$, $\\{D_{2}\,D_{3}\\}$. Likewise,
the logarithm of the six factors vertically aligned on the right is the
average of the log-scores for the three subsets of data singletons
$\\{D_{1}\\}$, $\\{D_{2}\\}$, $\\{D_{3}\\}$.

In the general case with $d$ data there are $\binom{d}{k}$ subsets with $k$
data points. We therefore obtain

$\log\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H\,I)\equiv\frac{1}{d}\sumop\slimits@_{i=1}^{d}\log\mathrm{P}(D_{i}\nonscript\,|\nonscript\;\mathopen{}D_{-i}\,H\,I)+{}\\\
\frac{1}{d}\sumop\slimits@_{i\in\\{1,\dotsc,d\\}}\frac{1}{d-1}\sumop\slimits@_{j\in\\{1,\dotsc,d\\}}^{j\neq
i}\\!\\!\\!\log\mathrm{P}(D_{-i,j}\nonscript\,|\nonscript\;\mathopen{}D_{-i,-j}\,H\,I)+{}\\\\[3.0pt]
\binom{d}{d-2}^{-1}\\!\\!\\!\sumop\slimits@_{i,j\in\\{1,\dotsc,d\\}}^{i<j}\frac{1}{d-2}\sumop\slimits@_{k\in\\{1,\dotsc,d\\}}^{k\neq
i,j}\\!\\!\\!\log\mathrm{P}(D_{-i,-j,k}\nonscript\,|\nonscript\;\mathopen{}D_{-i,-j,-k}\,H\,I)+{}\\\\[6.0pt]
{}\dotsb+{}\\\
\binom{d}{2}^{-1}\\!\\!\\!\sumop\slimits@_{i,j\in\\{1,\dotsc,d\\}}^{i<j}\frac{1}{2}\bigl{[}\log\mathrm{P}(D_{i}\nonscript\,|\nonscript\;\mathopen{}D_{j}\,H\,I)+\log\mathrm{P}(D_{j}\nonscript\,|\nonscript\;\mathopen{}D_{i}\,H\,I)\bigr{]}+{}\\\\[3.0pt]
\frac{1}{d}\sumop\slimits@_{i=1}^{d}\log\mathrm{P}(D_{i}\nonscript\,|\nonscript\;\mathopen{}H\,I),$
(10)

which can be compactly written

$\log\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H\,I)\equiv\sumop\slimits@_{k=1}^{d}\binom{d}{k}^{-1}\\!\smashoperator[]{\sumop\slimits@_{\begin{subarray}{c}\text{ordered}\\\
\text{$k$-tuples}\end{subarray}}^{}}\;\;\;\frac{1}{k}\;\;\smashoperator[]{\sumop\slimits@_{\begin{subarray}{c}\text{cyclic}\\\
\text{permutations}\end{subarray}}^{}}\log\mathrm{P}(D_{i_{1}}\nonscript\,|\nonscript\;\mathopen{}D_{i_{2}}\dotsm
D_{i_{k}}\,H\,I).$ (11)

That is, _the log-likelihood is the sum of all averaged log-scores that can be
formed from all (non-empty) data subsets with $k$ elements_, the average for
log-scores over $k$ data being taken over the $\binom{d}{k}$ subsets having
the same cardinality $k$.

There’s also an equivalent form with a slightly different cross-validating
interpretation: We take each datum $D_{j}$ in turn and calculate our log-
belief in it conditional on all possible subsets of remaining data, from the
empty subset with no data (term $k=0$), to the only subset $D_{-j}$ with all
data except $D_{j}$ (term $k=d-1$). These log-beliefs are averaged over the
$\binom{d-1}{k}$ subsets having the same cardinality $k$. The result can be
expressed as

$\log\mathrm{P}(D\nonscript\,|\nonscript\;\mathopen{}H\,I)\equiv\frac{1}{d}\sumop\slimits@_{j=1}^{d}\sumop\slimits@_{k=0}^{d-1}\binom{d-1}{k}^{-1}\smashoperator[]{\sumop\slimits@_{\begin{subarray}{c}\text{ordered}\\\
\text{$k$-tuples,}\\\ \text{$j$
excluded}\end{subarray}}^{}}\log\mathrm{P}(D_{j}\nonscript\,|\nonscript\;\mathopen{}D_{i_{1}}\dotsm
D_{i_{k}}\,H\,I).$ (12)

##  3 Brief discussion

It’s remarkable that the individual log-scores in expressions (11) and (12)
above are computationally expensive, but their sum results in the log-
likelihood, which is less expensive.

The relation (11) invites us to see the log-likelihood as a refinement and
improvement of the log-score. The log-likelihood takes into account not only
the log-score for the whole data, but also the log-scores for all possible
subsets of data. Figuratively speaking it examines the relationship between
data and hypothesis locally, globally, and on all intermediate scales. To me
this property makes the log-likelihood preferable to any single log-score
(besides the fact that the log-likelihood is directly obtained from the
principles of the probability calculus), because our interest is usually in
how the hypothesis $H$ relates to single data points as well as to any
collection of them. I hope to discuss this point, which also involves the
distinction between simple and composite hypotheses[bernardoetal1994_r2000, §
6.1.4], more in detail elsewhere[portamana9999b].

By applying the identity (6) and generalizing the expansion (7) to different
divisions of the data – leave-two-out, leave-three-out, and so on – we see
that the relation (11) can be generalized to any $k$-fold cross-validation
log-scores. Thus the log-likelihood is also equivalent to an average of _all
conceivable_ cross-validation log-scores for all subsets of data, though I
haven’t calculated the weights of such average.

###### Acknowledgements.

…to the Kavli Foundation and the Centre of Excellence scheme of the Research
Council of Norway (Yasser Roudi group) for financial support. To Aki Vehtari
for some references. To the staff of the NTNU library for their always prompt
assistance. To Mari, Miri, Emma for continuous encouragement and affection,
and to Buster Keaton and Saitama for filling life with awe and inspiration. To
the developers and maintainers of Open Science Framework, LaTeX, Emacs,
AUCTeX, Python, Inkscape, Sci-Hub for making a free and impartial scientific
exchange possible.

prenote(‘de $X$’ is listed under D, ‘van $X$’ under V, and so on, regardless
of national conventions.)


# Distribution of Square-Prime Numbers

Raghavendra N Bhat
Raghavendra N Bhat
University of Illinois, Urbana Champaign
Department of Mathematics
1409 West Green Street
Urbana, IL 61801<EMAIL_ADDRESS>(Corresponding Author)

###### Abstract.

For $a\neq 1$ and $p$ prime, we define numbers of the form $pa^{2}$ to be
Square-Prime (SP) Numbers. For example, 75 = 3 $\cdot$ 25; 108 = 3 $\cdot$ 36;
45 = 5 $\cdot$ 9\. These numbers are listed in the OEIS as A228056. We study
the properties of these numbers, their distribution/density and also develop a
few claims on their distribution/density. We rely on computer programs to
verify some conjectures up to large numbers.

##  1\. Introduction

SP Numbers can be generated in linear time allowing us to study some
interesting properties and generalize a few conjectures and claims that would
be more difficult to prove for the set of prime numbers. For example, there
are many unsolved problems along the lines of Goldbach and twin prime
conjectures. The intention of this paper is to study SP numbers, their density
with respect to primes and to also study their distribution on the number
line.

##  2\. Terminologies, Axioms and Trivial Results

  1. (1)

There are infinitely many SP Numbers. This is trivially proved owing to the
existence of infinite primes and infinite squares.

  2. (2)

There can be consecutive natural numbers who are both SP. An example is 27 and
28 (27 is $3\cdot 9$ and 28 is $7\cdot 4$)

  3. (3)

The product of two SP Numbers is always non-SP.

  4. (4)

An SP gap is defined as the difference between two consecutive SP Numbers.

  5. (5)

Consecutive natural numbers that are both SP are called SP Twins. These SP
pairs have gap 1.

DOI: 10.35834/YYYY/VVNNPPP

MSC2020: 11N25, 11P99

Key words and phrases: generalization of primes, distribution of prime-like
numbers

##  3\. First Few SP Numbers

Here are the first 100 SP numbers:

$\displaystyle
8,12,18,20,27,28,32,44,45,48,50,52,63,68,72,75,76,80,92,98,99,108,$
$\displaystyle
112,116,117,124,125,128,147,148,153,162,164,171,172,175,176,180,$
$\displaystyle
188,192,200,207,208,212,236,242,243,244,245,252,261,268,272,275,$
$\displaystyle
279,284,288,292,300,304,316,320,325,332,333,338,343,356,363,368,$
$\displaystyle
369,387,388,392,396,404,405,412,423,425,428,432,436,448,450,452,$
$\displaystyle 464,468,475,477,496,500,507,508,512,524,531,539,548,549.$

##  4\. Theorems on SP Numbers

######  Theorem 4.1.

For all natural numbers $n$ large enough, the number of SP Numbers smaller
than $n$ is asymptotic to $(\zeta(2)-1)\frac{n}{\log n}$.

###### Proof.

From the prime number theorem [2, pp. 111–115], we know that the number of
prime numbers $\leq$ $n$ equals $\frac{n}{\log n}+O(\frac{n}{\log^{2}n})$. Let
$SP(n)$ denote the number of SP numbers smaller than $n$. $SP(n)$ is equal to
the number of pairs $(a,p)$ where $a^{2}p$ $\leq$ $n$ where $p$ is a prime
number and $a$ is a natural number $\geq$ 2\. $SP(n)$ is

$\sum_{a=2}^{\sqrt{n/2}}\pi\left(\frac{n}{a^{2}}\right).$

Here, $\pi$ is the prime counting function. For arbitrary $A$ smaller than
$n$, the head term is

$\sum_{a=2}^{A}\pi\left(\frac{n}{a^{2}}\right).$

The tail is

$\sum_{a=A}^{\sqrt{n/2}}\pi\left(\frac{n}{a^{2}}\right).$

Set $A$ to be $\log^{2}n$. Thus, the head term is

$\sum_{a=2}^{\log^{2}n}\pi\left(\frac{n}{a^{2}}\right).$

We simplify this using PNT as follows:

$n\sum_{a=2}^{\log^{2}n}\frac{1}{a^{2}(\log n-\log(a^{2}))}$ $=\frac{n}{\log
n}\sum_{a=2}^{\log^{2}n}\frac{1}{a^{2}(1-\frac{2\log a}{\log n})}.$

Using the properties of geometric series, this simplifies to

$\frac{n}{\log
n}\sum_{a=2}^{\log^{2}n}\frac{1}{a^{2}}\left(1+O\left(\frac{\log a}{\log
n}\right)\right)$

$=\frac{n}{\log
n}\sum_{a=2}^{\log^{2}n}\frac{1}{a^{2}}+O\left(\frac{n}{\log^{2}n}\right).$

Simplifying the main term, yields

$\frac{n}{\log n}\sum_{a=2}^{\infty}\frac{1}{a^{2}}-\frac{n}{\log
n}\sum_{a=\log^{2}n}^{\infty}\frac{1}{a^{2}}$ $=\frac{n}{\log
n}\sum_{a=2}^{\infty}\frac{1}{a^{2}}+O\left(\frac{n}{\log^{2}n}\right)$
$=(\zeta(2)-1)\frac{n}{\log n}+O\left(\frac{n}{\log^{2}n}\right).$

Hence we have the head of the original asymptotic to be

$=(\frac{\pi^{2}}{6}-1)\frac{n}{\log n}+O\left(\frac{n}{\log^{2}n}\right).$

Simplifying the tail yields

$\sum_{a=\log^{2}n}^{\sqrt{n/2}}\pi(\frac{n}{a^{2}})\leq
n\sum_{a=\log^{2}n}^{\infty}\frac{1}{a^{2}(\log
n-\log(a^{2}))}=O\left(\frac{n}{\log^{2}n}\right).$

Finally, we have

$SP(n)=(\zeta(2)-1)\frac{n}{\log n}+O\left(\frac{n}{\log^{2}n}\right).$

This completes the proof. ∎

######  Corollary 4.2.

For a sufficiently large $n$, the number of prime numbers smaller than $n$ is
more than the number of SP numbers smaller than $n$.

###### Proof.

Since $\frac{\pi^{2}}{6}-1$ is less than 1, we have $SP(n)$ $\leq$ $\pi(n)$. ∎

######  Theorem 4.3.

There are infinitely many pairs of SP numbers for any existing SP gap.

###### Proof.

We use Gauss’s result on the general Pell equation [3, Chapter 8, Theorem 2 on
p. 57] that states that $x^{2}-dy^{2}=m$ has infinitely many solutions, if it
has one solution. Thus, if a gap $g$ occurs, we have $P_{1},P_{2},a,b$ such
that

$P_{1}a^{2}-P_{2}b^{2}=g.$

Multiplying by $P_{1}$ yields

$P_{1}a^{2}-(P_{1}P_{2})b^{2}=P_{1}g.$

We now have a general Pell equation

$x^{2}-(P_{1}P_{2})y^{2}=P_{1}g$ (4.1)

that has a solution. Thus, it has infinitely many solutions. In each solution,
$P_{1}$ divides $(P_{1}P_{2})y^{2}$ and $P_{1}g$. Hence, it also divides $x$.
Let $x=P_{1}k$. Thus, dividing every solution to (4.1) by $P_{1}$ yields

$P_{1}k^{2}-P_{2}y^{2}=g.$

Since $P_{1}k^{2}$ and $P_{2}y^{2}$ are SP Numbers which have gap $g$, the
proof is complete. ∎

As an immediate corollary, we have the existence of infinite SP twins, our
existing solution being 27 and 28.

##  5\. Conjectures on SP Numbers

These conjectures are all verified computationally beyond numbers in excess of
$10^{9}$. We look at these claims briefly to get some ideas for future
research.

######  Conjecture 5.1.

Every large natural number is a sum of two SP Numbers.

The Goldbach conjecture is probably the most famous conjecture about primes.
It claims that every even number greater than or equal to 4 can be written as
a sum of two primes. However, for SP Numbers, as there are no even-odd parity
restrictions, we claim that every sufficiently large natural number, whether
even or odd, can be written as a sum of two SP Numbers. This has been verified
beyond $10^{9}$ and is valid for all numbers starting $3931(27+3904)$. E.g.,
4041=116 + 3925, 10216=12 + 10204. [4]

######  Conjecture 5.2.

Every SP Number greater than 27 can be written as a sum of two SP Numbers in
at least one way.

This follows from conjecture 1 but takes effect with a lower starting point
(although specifically restricted to only SPs and not all naturals). E.g., 153
= 28 + 125.

######  Conjecture 5.3.

There exists at least one SP number between consecutive square numbers greater
than 500.

This is already conjectured for primes [5]. Computational verification shows
that it works for SP Numbers as well. E.g., Between 625 and 676 we have 637,
which is 49 $\cdot$ 13.

##  6\. On the Last Digit of SP Numbers

Using the fact that primes are asymptotically equidistributed mod 10 (i.e. we
have approximately the same number of primes ending in 1,3,7 and 9), we can
come up with similar asymptotic estimates for SP numbers ending in 1,3,7 and
9. For example, if we wish to estimate the count of SP numbers ending in 1
(i.e. congruent to 1 mod 10), we wish to look at pairs $(p,a)$ such that
$pa^{2}$ $\equiv$ 1 mod 10. $p$ can be $\equiv$ 1, 3, 7 or 9 mod 10. However,
if $pa^{2}$ has to end in 1, $p$ cannot end in 3 or 7, because square numbers
cannot end in 7 or 3 respectively. Thus we have two cases, $p$ ending in 1 or
9. When $p\equiv 1$ mod 10, we require $a^{2}$ to end in 1 (so that $pa^{2}$
ends in 1 too). Thus, $a$ can be $\equiv$ 1 or 9 mod 10. Thus, the number of
pairs $(p,a)$ such that $pa^{2}\equiv$ 1 mod 10 with $p$ ending in 1 is

$\frac{1}{4}\sum_{k=1}^{(n-1)/10}\pi\left(\frac{n}{(10k+1)^{2}}\right)+\pi\left(\frac{n}{(10k+9)^{2}}\right).$

Using the fact that prime numbers are asymptotically distributed mod 10, this
can be approximated by:

$\frac{1}{4}\frac{n}{\log
n}\sum_{k=1}^{\infty}\frac{1}{(10k+1)^{2}}+\frac{1}{(10k+9)^{2}}$
$=\frac{1}{400}\frac{n}{\log
n}\sum_{k=1}^{\infty}\frac{1}{(k+1/10)^{2}}+\frac{1}{(k+9/10)^{2}}$
$=\frac{1}{400}\frac{n}{\log n}(\zeta(2,1/10)+\zeta(2,9/10)-2).$

Here we use the Hurwitz zeta function which denotes
$\sum_{k=0}^{\infty}\frac{1}{(k+c)^{s}}$ as $\zeta(s,c).$ Similarly, if $p$
ends in 9, we require $a^{2}$ to end in 9. Thus, $a$ can be $\equiv$ 3 or 7
mod 10. Therefore, the asymptotic estimate for SP numbers ending in 1 with $p$
ending in 9 is the following:

$\frac{1}{400}\frac{n}{\log n}(\zeta(2,3/10)+\zeta(2,7/10)-2).$

Combining the above results, we have Theorem 6.1.

######  Theorem 6.1.

The asymptotic estimate of SP numbers ending in 1 is

$\frac{1}{400}\frac{n}{\log
n}(\zeta(2,1/10)+\zeta(2,9/10)+\zeta(2,3/10)+\zeta(2,7/10)-4).$

Using the process employed above, we find the same asymptotic estimate for SP
numbers ending in 3, 7, and 9.

Acknowledgements. I would like to acknowledge the professors, at the
University of Illinois for their valuable advice and support. I also
appreciate the comments from the editor, especially the guidance for proving
the results about the last digits of SP numbers.

## References

  * [1] The On-Line Encyclopedia of Integer Sequences, published electronically at https://oeis.org. Sequence A228056. https://oeis.org/A228057
  * [2] Davenport, H., _Multiplicative Number Theory_ , Third Edition. Springer, 2000.
  * [3] Mordell, L. J., _Diophantine Equations_ , Academic Press, 1969.
  * [4] Helfgott, H., _The ternary Goldbach conjecture is true_ , arXiv:1312.7748
  * [5] Bazzanella, D., _Some conditional results on primes between consecutive squares_ , Functiones et Approximatio Commentarii Mathematici, Funct. Approx. Comment. Math. 45(2), 255-263, (December 2011)