<commit_before>from __future__ import print_function
import numpy as np, sys

import regreg.api as rr
from selection.tests.instance import gaussian_instance
from scipy.stats import norm as ndist
from selection.randomized.api import randomization
from selection.adjusted_MLE.selective_MLE import M_estimator_map, solve_UMVU
from statsmodels.distributions.empirical_distribution import ECDF

def test_lasso_approx_var(n=100, p=50, s=5, signal=5., B=1000, lam_frac=1., randomization_scale=1., sigma= 1.):

    while True:
        X, y, beta, nonzero, sigma = gaussian_instance(n=n, p=p, s=s, rho=0., signal=signal, sigma=sigma,
                                                       random_signs=True, equicorrelated=False)
        n, p = X.shape
        lam = lam_frac * np.mean(np.fabs(np.dot(X.T, np.random.standard_normal((n, 2000)))).max(0)) * sigma

        loss = rr.glm.gaussian(X, y)
        epsilon = 1./np.sqrt(n)
        W = np.ones(p) * lam
        penalty = rr.group_lasso(np.arange(p),
                                 weights=dict(zip(np.arange(p), W)), lagrange=1.)

        randomizer = randomization.isotropic_gaussian((p,), scale=randomization_scale)
        M_est = M_estimator_map(loss, epsilon, penalty, randomizer, randomization_scale=randomization_scale, sigma=sigma)

        M_est.solve_map()
        active = M_est._overall

        true_target = np.linalg.inv(X[:, active].T.dot(X[:, active])).dot(X[:, active].T).dot(X.dot(beta))
        nactive = np.sum(active)

        if nactive > 0:
            approx_MLE, value, var, mle_map = solve_UMVU(M_est.target_transform,
                                                         M_est.opt_transform,
                                                         M_est.target_observed,
                                                         M_est.feasible_point,
                                                         M_est.target_cov,
                                                         M_est.randomizer_precision)

            boot_sample = np.zeros((B, nactive))
            resid = y - X[:, active].dot(M_est.target_observed)
            for b in range(B):
                boot_indices = np.random.choice(n, n, replace=True)
                boot_vector = (X[boot_indices, :][:, active]).T.dot(resid[boot_indices])
                target_boot = np.linalg.inv(X[:, active].T.dot(X[:, active])).dot(boot_vector) + M_est.target_observed
                boot_sample[b, :] = mle_map(target_boot)[0]

            print("estimated sd", boot_sample.std(0), np.sqrt(np.diag(var)))
            return np.true_divide((approx_MLE - true_target), boot_sample.std(0)), \
                   ((approx_MLE - true_target).sum()) / float(nactive), \
                   np.true_divide((approx_MLE - true_target), np.sqrt(np.diag(var)))
            break

    return np.true_divide((approx_MLE - true_target),np.sqrt(np.diag(var))), (approx_MLE - true_target).sum()/float(nactive)


if __name__ == "__main__":
    import matplotlib.pyplot as plt

    ndraw = 500
    bias = 0.
    pivot_obs_info= []
    pivot_bootstrap = []
    for i in range(ndraw):
        approx = test_lasso_approx_var(n=300, p=50, s=5, signal=3.5)
        if approx is not None:
            pivot_boot = approx[0]
            pivot_approx_info = approx[2]
            bias += approx[1]
            for j in range(pivot_boot.shape[0]):
                pivot_obs_info.append(pivot_approx_info[j])
                pivot_bootstrap.append(pivot_boot[j])

        sys.stderr.write("iteration completed" + str(i) + "\n")
        sys.stderr.write("overall_bias" + str(bias / float(i+1)) + "\n")
        #print("pivots", pivot_approx_info, pivot_boot)

    #if i % 10 == 0:
    plt.clf()
    ecdf_approx = ECDF(ndist.cdf(np.asarray(pivot_obs_info)))
    ecdf_boot = ECDF(ndist.cdf(np.asarray(pivot_bootstrap)))
    grid = np.linspace(0, 1, 101)
    print("ecdf", ecdf_boot(grid))
    plt.plot(grid, ecdf_approx(grid), c='red', marker='^')
    plt.plot(grid, ecdf_boot(grid), c='blue', marker='^')
    plt.plot(grid, grid, 'k--')
    #plt.show()
    plt.savefig("/Users/snigdhapanigrahi/Desktop/Boot_pivot_n300_p50_amp3.5_sigma2.png")<commit_msg>changed map for sigma not equal to 1.<commit_after>from __future__ import print_function
import numpy as np, sys

import regreg.api as rr
from selection.tests.instance import gaussian_instance
from scipy.stats import norm as ndist
from selection.randomized.api import randomization
from selection.adjusted_MLE.selective_MLE import M_estimator_map, solve_UMVU
from statsmodels.distributions.empirical_distribution import ECDF

def test_lasso_approx_var(n=100, p=50, s=5, signal=5., B=1000, lam_frac=1., randomization_scale=1., sigma= 1.5):

    while True:
        X, y, beta, nonzero, sigma = gaussian_instance(n=n, p=p, s=s, rho=0., signal=signal, sigma=sigma,
                                                       random_signs=True, equicorrelated=False)
        n, p = X.shape
        lam = lam_frac * np.mean(np.fabs(np.dot(X.T, np.random.standard_normal((n, 2000)))).max(0)) * sigma

        loss = rr.glm.gaussian(X, y)
        epsilon = 1./np.sqrt(n)
        W = np.ones(p) * lam
        penalty = rr.group_lasso(np.arange(p),
                                 weights=dict(zip(np.arange(p), W)), lagrange=1.)

        randomizer = randomization.isotropic_gaussian((p,), scale=randomization_scale)
        M_est = M_estimator_map(loss, epsilon, penalty, randomizer, randomization_scale=randomization_scale, sigma=sigma)

        M_est.solve_map()
        active = M_est._overall

        true_target = np.linalg.inv(X[:, active].T.dot(X[:, active])).dot(X[:, active].T).dot(X.dot(beta))
        nactive = np.sum(active)

        if nactive > 0:
            approx_MLE, value, var, mle_map = solve_UMVU(M_est.target_transform,
                                                         M_est.opt_transform,
                                                         M_est.target_observed,
                                                         M_est.feasible_point,
                                                         M_est.target_cov,
                                                         M_est.randomizer_precision)

            boot_sample = np.zeros((B, nactive))
            resid = y - X[:, active].dot(M_est.target_observed)
            for b in range(B):
                boot_indices = np.random.choice(n, n, replace=True)
                boot_vector = (X[boot_indices, :][:, active]).T.dot(resid[boot_indices])
                target_boot = np.linalg.inv(X[:, active].T.dot(X[:, active])).dot(boot_vector) + M_est.target_observed
                boot_sample[b, :] = mle_map(target_boot)[0]

            print("estimated sd", boot_sample.std(0), np.sqrt(np.diag(var)))
            return np.true_divide((approx_MLE - true_target), boot_sample.std(0)), \
                   ((approx_MLE - true_target).sum()) / float(nactive), \
                   np.true_divide((approx_MLE - true_target), np.sqrt(np.diag(var)))
            break

    return np.true_divide((approx_MLE - true_target),np.sqrt(np.diag(var))), (approx_MLE - true_target).sum()/float(nactive)


if __name__ == "__main__":
    import matplotlib.pyplot as plt

    ndraw = 500
    bias = 0.
    pivot_obs_info= []
    pivot_bootstrap = []
    for i in range(ndraw):
        approx = test_lasso_approx_var(n=300, p=50, s=5, signal=3.5)
        if approx is not None:
            pivot_boot = approx[0]
            pivot_approx_info = approx[2]
            bias += approx[1]
            for j in range(pivot_boot.shape[0]):
                pivot_obs_info.append(pivot_approx_info[j])
                pivot_bootstrap.append(pivot_boot[j])

        sys.stderr.write("iteration completed" + str(i) + "\n")
        sys.stderr.write("overall_bias" + str(bias / float(i+1)) + "\n")
        #print("pivots", pivot_approx_info, pivot_boot)

    #if i % 10 == 0:
    plt.clf()
    ecdf_approx = ECDF(ndist.cdf(np.asarray(pivot_obs_info)))
    ecdf_boot = ECDF(ndist.cdf(np.asarray(pivot_bootstrap)))
    grid = np.linspace(0, 1, 101)
    print("ecdf", ecdf_boot(grid))
    plt.plot(grid, ecdf_approx(grid), c='red', marker='^')
    plt.plot(grid, ecdf_boot(grid), c='blue', marker='^')
    plt.plot(grid, grid, 'k--')
    #plt.show()
    plt.savefig("/Users/snigdhapanigrahi/Desktop/Boot_pivot_n300_p50_amp3.5_sigma1.5.png")
{'fraction_non_alphanumeric': 0.12961855793375582, 'fraction_numerical': 0.01446849216788234, 'mean_word_length': 2.2169230769230768, 'pattern_counts': {'":': 2, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2264796', 'n_tokens_mistral': 2869, 'n_tokens_neox': 2746, 'n_words': 484}