def METHOD_NAME(shape, value, axis):
    """Test whether the logprob for ```pt.max``` produces the corrected
    The fact that order statistics of i.i.d. uniform RVs ~ Beta is used here:
        U_1, \\dots, U_n \\stackrel{\text{i.i.d.}}{\\sim} \text{Uniform}(0, 1) \\Rightarrow U_{(k)} \\sim \text{Beta}(k, n + 1- k)
    for all 1<=k<=n
    """
    x = pt.random.uniform(0, 1, size=shape)
    x.name = "x"
    x_max = pt.max(x, axis=axis)
    x_max_value = pt.scalar("x_max_value")
    x_max_logprob = logp(x_max, x_max_value)
    assert_no_rvs(x_max_logprob)
    test_value = value
    n = np.prod(shape)
    beta_rv = pt.random.beta(n, 1, name="beta")
    beta_vv = beta_rv.clone()
    beta_rv_logprob = logp(beta_rv, beta_vv)
    np.testing.assert_allclose(
        beta_rv_logprob.eval({beta_vv: test_value}),
        (x_max_logprob.eval({x_max_value: test_value})),
        rtol=1e-06,
    )