Marginal distributions given the distribution of range

Question

I'm working with an upper diagonal distribution whose distance from the diagonal is Lomax Pareto (Type II) distribution.

The distance of a point from the diagonal line y = x is $\frac{\sqrt{(x_0-y_0)^2}}{\sqrt{2}} \propto \text{hypotenuse of the right triangle formed by the point and the line } \sqrt{(x_0-y_0)^2}$, i.e. the range between x and y

Now, since the actual distance is distributed Lomax Pareto, and the range of |x - y| is proportional to the distance, the distribution of the range is similarly Lomax

$$ \begin{align} p(r) &= {\alpha \over \lambda} \left[{1 + {r \over \lambda}}\right]^{-(\alpha+1)}, \qquad r \geq 0 \\ &= 2 * \int_0^\infty f(z)f(z+r) dz \text{ (since there are only 2 points)} \end{align} $$

I my first attempt was to try a Lomax marginal for each variable and then work forward.

$$ \begin{align} f_{R_2}(r) &= 2 * \int_0^\infty f(z)f(z+r) dz, \qquad r \geq 0 \\ &= 2 * \int_0^\infty \frac{\alpha}{\lambda} \left\{ 1 + \frac{z}{\lambda} \right\}^{-\alpha-1} \frac{\alpha}{\lambda} \left\{ 1 + \frac{z + r}{\lambda} \right\}^{-\alpha-1} dz \\ &= 2 * \left( \frac{\alpha}{\lambda} \right)^2 \int_0^\infty \left\{ \frac{\lambda + z}{\lambda} \right\}^{-\alpha-1} \left\{ \frac{\lambda + z + r}{\lambda} \right\}^{-\alpha-1} dz \\ &= 2 \alpha^2 \lambda^{\alpha - 1} \int_0^\infty \left\{ \lambda + z \right\}^{-\alpha-1} \left\{ \lambda + z + r \right\}^{-\alpha-1} dz \\ \end{align} $$

...but I cannot push the integral through.

Are there any methods that can be used to identify one or more possible marginal distributions whose range forms a Lomax Pareto (Type II) distribution?

Answer 1

If what you're after is the equivalent of the distribution of $R=\max(X_1,X_2)-\min(X_1,X_2)$ with $X_i\sim \text{LomaxPareto}(\alpha,\lambda)$ and $X_1$ and $X_2$ independent, then the following might be helpful.

Using the functions available in Mathematica one can find the pdf of $R$:

(* Single Lomax Pareto distribution *)
dist = ParetoDistribution[\[Lambda], \[Alpha], 0];

(* Distribution of minimum and maximum for a random sample of size 2 from Lomax Pareto distribution *)
odist = OrderDistribution[{dist, 2}, {1, 2}];

(* Distribution of range *)
tdist = TransformedDistribution[x2 - x1, {x1, x2} \[Distributed] odist];
pdf = FullSimplify[PDF[tdist, z], Assumptions -> {z > 0, \[Lambda] > 0, \[Alpha] > 0}]

Unfortunately, this formula does not work for integer values of $\alpha$. So we must try a few values and hope for a pattern to become evident.

t = Table[dist = ParetoDistribution[\[Lambda], \[Alpha], 0];
  odist = OrderDistribution[{dist, 2}, {1, 2}];
  tdist = TransformedDistribution[x2 - x1, {x1, x2} \[Distributed] odist];
  {\[Alpha],  FullSimplify[PDF[tdist, z], Assumptions -> {z > 0, \[Lambda] > 0, \[Alpha] > 0}]},
  {\[Alpha], 1, 5}];
TableForm[t, TableHeadings -> {None, {"\[Alpha]", "pdf"}}]

From looking at the patterns and consulting oeis.org the following will produce the pdf for any (legal) combination of real values of $z$, $\lambda$, and $\alpha$:

(* oeis.org: A068553 *)
c1[n_] := (LCM @@ Range[2 n])/(n Binomial[2 n, n])
(* oeis.org: A120114 *)
c2[n_] := (LCM @@ Range[2 n + 4])/LCM @@ Range[2 n + 2]
(* oeis.org:  A005430 *)
c3[n_] := n Binomial[2 n, n]

(* Recursive function *)
a[z_, \[Lambda]_, 1] := 1
a[z_, \[Lambda]_, 2] := z^2 - 6 z \[Lambda] - 6 \[Lambda]^2;
a[z_, \[Lambda]_, \[Alpha]_] := c1[\[Alpha]] z^(2 \[Alpha] - 2) - 
  c2[\[Alpha] - 2] \[Lambda] (z + \[Lambda]) a[z, \[Lambda], \[Alpha] - 1] // Expand // Simplify

pdf[z_, \[Lambda]_, \[Alpha]_] := If[z <= 0 || \[Lambda] <= 0 || \[Alpha] <= 0, 0,
  If[IntegerQ[\[Alpha]],
   2*\[Alpha] \[Lambda]^\[Alpha] (z (z + 2 \[Lambda]) a[z, \[Lambda], \[Alpha]]/(z + \[Lambda])^\[Alpha] - 
    (-1)^\[Alpha] c1[\[Alpha]] c3[\[Alpha]] \[Lambda]^\[Alpha]*
    Log[\Lambda]/(z + \[Lambda])])/(c1[\[Alpha]] z^(2 \[Alpha] + 1)),
    (2 \[Alpha] \[Lambda]^\[Alpha] (-((4^\[Alpha] Sqrt[\[Pi]]* 
    z^(-2 \[Alpha]) \[Lambda]^\[Alpha] Csc[\[Pi] \[Alpha]] Gamma[1/2 + \[Alpha]])/
    Gamma[\[Alpha]]) + (z + \[Lambda])^-\[Alpha] *
    (-1 + 2 Hypergeometric2F1[1, \[Alpha], 1 - \[Alpha], \[Lambda]/(z + \[Lambda])])))/z]]

Here is an example:

Plot[pdf[z, 2, 3/2], {z, 0, 20}, WorkingPrecision -> 50, PlotRange -> {All, {0, All}}]

I used a high WorkingPrecision as the formulas for the pdf's are very numerically sensitive (or maybe "unstable with machine precision" is a better way to phrase it). For example, using 32/10 (a rational number in Mathematica) for either $\lambda$ or $\alpha$ is recommended over using 3.2.

Answer 2

Edit: In the meantime I've tried a variety of sets of parameters and found the result is highly numerically unstable. I'll investigate more and update.

Here is a more direct answer that does not rely on as many special functions in Mathematica.

First, the initial integral you were solving:

Integrate[(2 \[Alpha]^2 ((z + \[Lambda])/\[Lambda])^(-1 - \[Alpha]) ((z + r + \[Lambda])/\[Lambda])^(-1 - \[Alpha]))/\[Lambda]^2,
  {z, 0, \[Infinity]}, Assumptions -> {\[Alpha] > 0, \[Lambda] > 0, r > 0}]

(* (r^(-1 - 2 \[Alpha]) ((-1)^(1 + \[Alpha]) Sqrt[\[Pi]] Beta[-(\[Lambda]/r), -\[Alpha], -\[Alpha]] + 
    4^\[Alpha] Gamma[-\[Alpha]] Gamma[1/2 + \[Alpha]]))/Sqrt[\[Pi]] *)

$$\int_0^{\infty } \frac{2 \alpha ^2 \left(\frac{\lambda +z}{\lambda }\right)^{-\alpha -1} \left(\frac{\lambda +r+z}{\lambda }\right)^{-\alpha -1}}{\lambda ^2} \, dz = \frac{r^{-2 \alpha -1} \left(4^{\alpha } \Gamma (-\alpha ) \Gamma \left(\alpha +\frac{1}{2}\right)+\sqrt{\pi } (-1)^{\alpha +1} B_{-\frac{\lambda }{r}}(-\alpha ,-\alpha )\right)}{\sqrt{\pi }}$$

where $\Gamma[.]$ is the gamma function and $B_z(a,b)$ is the incomplete beta function.

Note that in my previous answer the special Mathematica functions for random variables results in an equivalent but more complicated representation of the result. That happens. And will many times require more than machine precision numbers in part depending on how close $\alpha$ is to an integer.

So this function also does not work for integer values of $\alpha$. However, each desired positive integer value of $\alpha$ can be found by a separate integration:

$\alpha=1$

Integrate[(2 \[Alpha]^2 ((z + \[Lambda])/\[Lambda])^(-1 - \[Alpha])*
  ((z + r + \[Lambda])/\[Lambda])^(-1 - \[Alpha]))/\[Lambda]^2 /. \[Alpha] -> 1,
  {z, 0, \[Infinity]}, Assumptions -> {\[Lambda] > 0, z > 0, r > 0}]

$$\frac{2 \lambda \left(\frac{\lambda r}{\lambda +r}+2 \lambda \log \left(\frac{\lambda }{\lambda +r}\right)+r\right)}{r^3}$$

$\alpha=5$

Integrate[(2 \[Alpha]^2 ((z + \[Lambda])/\[Lambda])^(-1 - \[Alpha])*
  ((z + r + \[Lambda])/\[Lambda])^(-1 - \[Alpha]))/\[Lambda]^2 /. \[Alpha] -> 5,
  {z, 0, \[Infinity]}, Assumptions -> {\[Lambda] > 0, z > 0, r > 0}]

$$\frac{5 \lambda ^{10} \left(\frac{r (2 \lambda +r) \left(1260 \lambda ^8+33 \lambda ^2 r^6-126 \lambda ^3 r^5+672 \lambda ^4 r^4+4410 \lambda ^5 r^3+7350 \lambda ^6 r^2-9 \lambda r^7+2 r^8+5040 \lambda ^7 r\right)}{\lambda ^5 (\lambda +r)^5}+2520 \log \left(\frac{\lambda }{\lambda +r}\right)\right)}{r^{11}}$$

The general form for integer values of $\alpha$ can be found in my previous answer.