BLBVar_methods.py

"""
    BLBVar:
Input: data D of size n, private estimator theta_tilde, nonprivate plug-in theta,
       subset size m, privacy eps_var, smoothing rho, variance bound sigma2_max,
       MC resamples Nmc, confidence level 1-alpha

1. Partition D into s = floor(n/m) disjoint subsets D_1,...,D_s of size m
2. For each subset D_i:
      a. Compute theta_hat_i = theta(D_i)
      b. For b = 1,...,Nmc:
            - Draw bootstrap resample D_i^{*(b)} of size n from empirical(D_i)
            - Compute u_{i,b} = sqrt(n) * (theta_tilde(D_i^{*(b)}) - theta_hat_i)
      c. Set v_i = average_b ||u_{i,b}||^2
3. Compute sigma2_hat = PrivMedian(v_1,...,v_s; eps_var, rho, [0, sigma2_max])
4. Compute theta_priv = theta_tilde(D)
5. Return CI = theta_priv ± z_{1-alpha/2} * sqrt(sigma2_hat / n)"""
import numpy as np
from scipy.stats import norm

from Moshes_experiments.implementations.optimized_quantile_inverse_sensitivity import vectorized_quantile_combined
from experiment_utils import build_estimator_from_enum, create_logreg_est
from scripts.clean_one import runs_per_boot_sens_unstud

# Experiment‑level parameters
MC_EXP = 1.5
MIN_BOOTS = 100
MAX_BOOTS = 1000
CONST_INV_SENS = 10
STD_UP = 10



def BLBVar(X, eps, estimator_to_build, dist, params):
    """Compute a private estimate of the variance of a private estimator using the BLB approach.
    todo: handle multiple experiments.
    """
    L, U, significance = params.L, params.U, params.significance
    if np.ndim(X)==1:
        X = X.reshape(1,-1)

    n = X.shape[1]
    true_dens = dist.pdf(dist.ppf(0.5))
    true_std = 1 / (2 * true_dens) # todo: it is tailored to the median, need to be generalized
    s, b, run_per_boots = runs_per_boot_sens_unstud(eps, n)

    eps_var, eps_est = eps * 0.5, eps * 0.5
    priv_estimator = build_estimator_from_enum(estimator_to_build, eps_est, L, U, wrap_3d=False, axis=-1, is_sorted=True)
    non_priv_estimator = build_estimator_from_enum(estimator_to_build, np.inf, L, U, wrap_3d=False, axis=-1,is_sorted=True)

    Lb, Ub = [], []
    for i in range(X.shape[0]):
        m, p = blb_var_helper(X[i], eps_var, priv_estimator, non_priv_estimator, n, run_per_boots, s, true_std, significance)
        Lb.append(m)
        Ub.append(p)
    return np.column_stack([Lb, Ub])


def blb_var_helper(Xi, eps_var, priv_est, non_priv_est, n, run_per_boots, s, true_std, significance):
    # Partition D into s disjoint subsets
    # shuffle the data to ensure randomness in the subsets
    subsets = np.array_split(np.random.permutation(Xi), s)

    bag_vars = np.zeros(int(s))
    for i, subset in enumerate(subsets):
        theta_hat = non_priv_est(np.sort(subset))

        # Bootstrap resampling
        # todo: handle the other bootstrap sizes (n-> m)
        resamples = np.random.choice(subset, size=(run_per_boots, n), replace=True)
        resamples.sort(axis=-1)
        theta_tildes = priv_est(resamples)

        bag_vars[i] = np.mean((theta_tildes - theta_hat) ** 2)

    # Compute sigma2_hat using a private median mechanism
    sigma_for_bounds = (STD_UP * true_std / np.sqrt(n)) ** 2
    sigma2_hat = np.maximum(vectorized_quantile_combined(bag_vars, L=0, U=sigma_for_bounds, sigma= (2.0 / max(eps_var, 1e-12)), quantile=0.5, is_sorted=False), 1e-10)

    theta_priv_cent = priv_est(np.sort(Xi))
    z = norm.ppf(1 - significance / 2)
    return theta_priv_cent - z * np.sqrt(sigma2_hat), theta_priv_cent + z * np.sqrt(sigma2_hat)

def BLBVar_logreg(X, y, eps, center_mech_enum, settings, log_exp_setting, params):
    B, n = X.shape
    s, _, run_per_boots = runs_per_boot_sens_unstud(eps, n)
    eps_var, eps_est = 0.5 * eps, 0.5 * eps
    true_sd_ = log_exp_setting.fisher_n[n]

    priv_estimator, sigma_center = create_logreg_est(center_mech_enum, params, n, n, 1, eps_est, log_exp_setting)
    print("sigma_center: *", sigma_center)
    non_priv_estimator, _ = create_logreg_est(None, params, n, n, 1, np.inf, log_exp_setting)
    Lb, Ub = [], []
    for b in range(B):
        m, p = BLBvar_logreg_helper(X[b],y[b], eps_var, priv_estimator, non_priv_estimator, n, run_per_boots, s, true_sd_, params.significance)
        Lb.append(m)
        Ub.append(p)
    return np.column_stack([Lb, Ub])

def BLBvar_logreg_helper(Xi, yi, eps_var, priv_est, non_priv_est, n, run_per_boots, s, sigma_for_bounds, significance):
    # Partition D into s disjoint subsets
    perm = np.random.permutation(n)
    # subsets_X = np.array_split(Xi[perm], s)
    # subsets_y = np.array_split(yi[perm], s)
    idx_subsets = np.array_split(perm, s)


    bag_vars = np.zeros(int(s))
    for i, idx in enumerate(idx_subsets):
        X_sub = Xi[idx]
        y_sub = yi[idx]

        theta_hat = non_priv_est(X_sub, y_sub)

        # Bootstrap resampling
        boot_idx_local = np.random.choice(len(idx), size=(run_per_boots, n), replace=True)
        theta_tildes = np.empty(run_per_boots, dtype=float)
        for r in range(run_per_boots):
            j = boot_idx_local[r]
            theta_tildes[r] = priv_est(X_sub[j], y_sub[j])

        bag_vars[i] = np.mean((theta_tildes - theta_hat) ** 2)

    # Private aggregation of bag variances via private median
    sigma2_hat = np.maximum(vectorized_quantile_combined(bag_vars, L=0, U=sigma_for_bounds, sigma= (2.0 / max(eps_var, 1e-12)), quantile=0.5, is_sorted=False), 1e-10)

    # Private center on full sample
    theta_priv_cent = priv_est(Xi, yi)

    z = norm.ppf(1 - significance / 2.0)
    radius = z * np.sqrt(sigma2_hat)
    return theta_priv_cent - radius, theta_priv_cent + radius