Implement logprob derivation for order statistics

[ricardoV94]

Description

The idea would be to derive the logprob for the following type of graphs

import aesara.tensor as at
import pymc as pm

y_rv = at.max(at.random.uniform(0, 1, size=3))   # or min
y_rv = at.sort(at.random.uniform(0, 1, size=3))
y_rv = at.sort(at.random.uniform(0, 1, size=3))[idx]  # max / min correspond to idx==-1 or idx==0

y_vv = y_rv.clone()
pm.logprob.factorized_joint_logprob({y_rv: y_vv})

https://en.wikipedia.org/wiki/Order_statistic#Probabilistic_analysis

This might be a bit far-fetched / difficult to find a good general solution that goes beyond a few simple cases (e.g, order statistics with non-i.i.d RVs):

y_rv = at.max(at.random.uniform([0, 1], [2, 3]))
y_rv = at.max(at.stack([at.random.uniform(0, 1), at.random.normal(0, 1)]))

Probability for order statistics of IID variables are pretty straightforward to obtain, requiring expressions that depend on the CDFs. For non IID, things grow quickly in complexity: https://en.wikipedia.org/wiki/Bapat%E2%80%93Beg_theorem

[ricardoV94]

Superseded by #7121