def lin_reg_var_organinc_lasso(X, y, pen_val, coef, intercept=None,
                               sample_weight=None):
    """
    Estimates the linear regression variance using the organic Lasso estimate (Yu and Bien, 2019). This requires first fitting a L1 squared penalized model.

    Parameters
    ----------
    X: array-like, shape (n_samples, n_features)
        The training covariate data.

    y: array-like, shape (n_samples, ) or (n_samples, n_responses)
        The training response data.

    pen_val: float
        The squared L1 penalty value.

    coef: array-like, shape (n_features, )
        The Lasso estimated coefficient.

    intercept: None, float
        The (optional) estimated intercept.

    sample_weight: None or array-like,  shape (n_samples,)
        (Optional) Individual weights for each sample.

    Output
    ------
    sigma_sq: float
        An estimate of the noise variance.

    References
    ----------
    Yu, G. and Bien, J., 2019. Estimating the error variance in a high-dimensional linear model. Biometrika, 106(3), pp.533-546.
    """

    if sample_weight is not None:
        raise NotImplementedError

    # See equation (17)
    y_hat = X @ coef
    if intercept is not None:
        y_hat += intercept

    RSS = ((y - y_hat) ** 2).sum()
    L1_sq = abs(coef).sum() ** 2
    sigma_sq = (1 / X.shape[0]) * RSS + 2 * pen_val * L1_sq

    return sigma_sq