<jupyter_start><jupyter_text>People always ask: "can you randomize several times and use the proportion of selection, instead ofjust one randomization"?Let's try to figure this out.<jupyter_code>import numpy as np
import regreg.api as rr
import seaborn as sns
%matplotlib inline
%load_ext rpy2.ipython


import matplotlib.pyplot as plt
import scipy.stats
import statsmodels.api as sm

from selection.distributions.discrete_family import discrete_family

ntries, sigma, q = 21, 1, 0.3

def algorithm(Z, ntries=ntries, q=q):
    
    proportion = 0
    for _ in range(ntries):
        proportion += ((Z + sigma * np.random.standard_normal() > 2) +
                       (Z + sigma * np.random.standard_normal() < -2)) > 0

    proportion /= ntries
    return proportion > q

Z = np.linspace(-8, 8, 1001)
def fit_algorithm(algorithm, B=5000, ntries=ntries, q=q, Zval=Z, link='logit'):
    
    Z = np.random.standard_normal(B) * 2
    Z = np.hstack([Z, 
                   np.random.standard_normal(B), 
                   np.random.standard_normal(B) * 3, 
                   np.random.standard_normal(B) * 0.5])

    Y = np.array([algorithm(z, ntries=ntries, q=q) for z in Z])

    %R -i Y,Z,Zval 
    %R Z = as.numeric(Z*1)
    if link == 'probit':
        %R M2 = glm(Y ~ poly(Z, 2), family=binomial(link=probit))
    else:
        %R M2 = glm(Y ~ poly(Z, 2), family=binomial(link=logit))

    %R W = predict(M2, newdata=data.frame(Z=Zval), type='link')
    W = %R W
    if link == 'probit':
        return scipy.stats.norm.cdf(W) 
    else:
        return np.exp(W) / (1 + np.exp(W))

def simulate(ntries=ntries, sigma=sigma, truth=0):
               
    while True:
        Z = np.random.standard_normal() + truth
        if algorithm(Z, ntries, q=q):
            return Z

Z = np.linspace(-8, 8, 1001)
W1 = fit_algorithm(algorithm, ntries=ntries, q=q, Zval=Z)
plt.plot(Z, np.log(W1))
selective_law1 = discrete_family(Z, W1 * scipy.stats.norm.pdf(Z))


def pivot1(z, truth=0):
    return 1 - selective_law1.cdf(truth, z)


P0 = []

for _ in range(1000):
    P0.append((pivot1(simulate()),
              1 - scipy.stats.norm.cdf(simulate())))
P0 = np.array(P0)
U = np.linspace(0, 1, 101)
plt.plot(U, sm.distributions.ECDF(P0[:,0])(U), 'c', label='fit')
plt.plot(U, sm.distributions.ECDF(P0[:,1])(U), 'y', label='naive')

plt.plot([0, 1], [0, 1], 'k--')
plt.legend()
PA = []
truth = 1
for _ in range(1000):
    PA.append((pivot1(simulate(truth=truth), truth=truth),
               1 - scipy.stats.norm.cdf(simulate() - 1)))
PA = np.array(PA)
U = np.linspace(0, 1, 101)
plt.plot(U, sm.distributions.ECDF(PA[:,0])(U), 'c', label='fit')
plt.plot(U, sm.distributions.ECDF(PA[:,1])(U), 'y', label='naive')
plt.legend()
plt.plot([0, 1], [0, 1], 'k--')
Z0 = np.linspace(-2,2,501)

LU1 = []
for z in Z0:
    selective_law = discrete_family(Z, W1 * scipy.stats.norm.pdf(Z))
    LU1.append(selective_law.equal_tailed_interval(z))
LU1 = np.array(LU1)
plt.plot(Z0, LU1[:,0], 'c', label='fit')
plt.plot(Z0, LU1[:,1], 'c')
plt.legend()
coverage, ncover, truth = 0, 500, 0
lengths = []
for _ in range(ncover):
    z = simulate(truth=truth)
    selective_law = discrete_family(Z, W1 * scipy.stats.norm.pdf(Z))
    L, U = selective_law.equal_tailed_interval(z)
    coverage += (L < truth) * (U > truth)
    lengths.append(U-L)
    
coverage / ncover, np.mean(lengths), np.std(lengths)
coverage, ncover, truth = 0, 500, 2.5
lengths = []
for _ in range(ncover):
    z = simulate(truth=truth)
    selective_law = discrete_family(Z, W1 * scipy.stats.norm.pdf(Z))
    L, U = selective_law.equal_tailed_interval(z)
    coverage += (L < truth) * (U > truth)
    lengths.append(U-L)
    
coverage / ncover, np.mean(lengths), np.std(lengths)
coverage, ncover, truth = 0, 500, -1.
lengths = []
for _ in range(ncover):
    z = simulate(truth=truth)
    selective_law = discrete_family(Z, W1 * scipy.stats.norm.pdf(Z))
    L, U = selective_law.equal_tailed_interval(z)
    coverage += (L < truth) * (U > truth)
    lengths.append(U-L)
    
coverage / ncover, np.mean(lengths), np.std(lengths)<jupyter_output><empty_output><jupyter_text>Increasing number of tries<jupyter_code>ntries, sigma, q = 41, 0.5, 0.65
Z = np.linspace(-8, 8, 1001)

def pivot(z, truth=0):
    return 1 - selective_law.cdf(truth, z)

def pivot0(z, truth=0):
    return 1 - selective_law0.cdf(truth, z)

def algorithm(Z, ntries=ntries, q=q):
    
    proportion = 0
    for _ in range(ntries):
        proportion += ((Z + sigma * np.random.standard_normal() > 2) + 
                       (Z + sigma * np.random.standard_normal() < -2)) > 0

    proportion /= ntries
    return proportion > q


W1 = fit_algorithm(algorithm, ntries=ntries, q=q, Zval=Z)
selective_law1 = discrete_family(Z, W1 * scipy.stats.norm.pdf(Z))

def pivot1(z, truth=0):
    return 1 - selective_law1.cdf(truth, z)

pivot(simulate())

P0 = []
truth = 0.05
for _ in range(1000):
    P0.append((pivot1(simulate(ntries=ntries, sigma=sigma, truth=truth)),
              1-scipy.stats.norm.cdf(simulate(ntries=ntries, sigma=sigma, truth=truth) - truth)))
P0 = np.array(P0)
U = np.linspace(0, 1, 101)
plt.plot(U, sm.distributions.ECDF(P0[:,0])(U), 'c', label='fit')
plt.plot(U, sm.distributions.ECDF(P0[:,1])(U), 'y', label='naive')
plt.plot([0, 1], [0, 1], 'k--')
plt.legend()
truth = -1
PA = []

for _ in range(1000):
    PA.append((pivot1(simulate(ntries=ntries, sigma=sigma, truth=truth), truth=truth),
              1-scipy.stats.norm.cdf(simulate(ntries=ntries, sigma=sigma, truth=truth) - truth)))
PA = np.array(PA)
U = np.linspace(0, 1, 101)
plt.plot(U, sm.distributions.ECDF(PA[:,0])(U), 'c', label='fit', linewidth=2)
plt.plot(U, sm.distributions.ECDF(PA[:,1])(U), 'y', label='naive')
plt.plot([0, 1], [0, 1], 'k--')
plt.legend()

truth = -2
PA = []

for _ in range(1000):
    PA.append((pivot1(simulate(ntries=ntries, sigma=sigma, truth=truth), truth=truth),
              1-scipy.stats.norm.cdf(simulate(ntries=ntries, sigma=sigma, truth=truth) - truth)))
PA = np.array(PA)
U = np.linspace(0, 1, 101)
plt.plot(U, sm.distributions.ECDF(PA[:,0])(U), 'c', label='fit', linewidth=2)
plt.plot(U, sm.distributions.ECDF(PA[:,1])(U), 'y', label='naive')
plt.plot([0, 1], [0, 1], 'k--')
plt.legend()


truth = 1
PA = []

for _ in range(1000):
    PA.append((pivot1(simulate(ntries=ntries, sigma=sigma, truth=truth), truth=truth),
              1-scipy.stats.norm.cdf(simulate(ntries=ntries, sigma=sigma, truth=truth) - truth)))
PA = np.array(PA)
U = np.linspace(0, 1, 101)
plt.plot(U, sm.distributions.ECDF(PA[:,0])(U), 'c', label='fit', linewidth=2)
plt.plot(U, sm.distributions.ECDF(PA[:,1])(U), 'y', label='naive')
plt.plot([0, 1], [0, 1], 'k--')
plt.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18532876099597437, 'fraction_numerical': 0.03548531385120024, 'mean_word_length': 3.585099111414901, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '17278', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2938, 'n_tokens_neox': 2693, 'n_words': 560, 'path': 'doc/notebooks/learning/Non convex region.ipynb', 'repo_name': 'yiling-h/selective-inference', 'size': '6707'}