def optimal_davies_pvalue(q, mu, var, kur, w, remain_var, df, trho, grid, pmin=None):
    r"""Joint significance of statistics derived from chi2-squared distributions.

    Parameters
    ----------
    q : array_like
        Most significant of the independent test statistics, qₘᵢₙ(𝜌ᵥ) [1].
    mu : float
        Mean of the linear combination of the chi-squared distributions.
    var : float
        Variance of the linear combination of the chi-squared distributions.
    kur : float
        Kurtosis of the linear combination of the chi-squared distributions.
    w : array_like
        Weights of the linear combination.
    remain_var : float
        Remaining variance assigned to a Normal distribution.
    df : float
        Overall degrees of freedom.
    trho : array_like
        Weight between the combination of chi-squared distributions and independent
        chi-squared distributions.
    grid : array_like
        Grid parameters.
    pmin : float
        Boundary of the possible final p-values.

    Returns
    -------
    float
        Estimated p-value.

    References
    ----------
    [1] Lee, Seunggeun, Michael C. Wu, and Xihong Lin. "Optimal tests for rare variant
    effects in sequencing association studies." Biostatistics 13.4 (2012): 762-775.
    """
    q = asarray(q, float)
    mu = float(mu)
    var = float(var)
    kur = float(kur)

    w = asarray(w, float)
    remain_var = float(remain_var)
    df = float(df)
    trho = asarray(trho, float)
    grid = asarray(grid, float)

    lambda_threshold = sum(w) * 10 ** 4
    chi2s = [ChiSquared(i, 0.0, 1) for i in w]
    args = (q, mu, var, kur, lambda_threshold, remain_var, df, trho, grid, chi2s)
    try:
        u = _find_upper_bound(args)
        re = quad(
            _davies_function, 0, u, args, limit=1000, epsabs=_EPSABS, full_output=1
        )
        pvalue = 1 - re[0]
        if re[1] > 1e-6:
            pvalue = _skat_liu_pvalue(
                q, mu, var, kur, w, remain_var, df, trho, grid, pmin
            )

    except RuntimeError:
        pvalue = _skat_liu_pvalue(q, mu, var, kur, w, remain_var, df, trho, grid, pmin)

    if pmin is not None:
        if pmin * len(grid) < abs(pvalue):
            pvalue = pmin * len(grid)
    return pvalue