<jupyter_start><jupyter_text>k-NN, Function Expectation, Density Estimation<jupyter_code>from experiment_framework.helpers import build_convergence_curve_pipeline
from empirical_privacy.one_bit_sum import GenSampleOneBitSum

# from empirical_privacy import one_bit_sum_joblib as one_bit_sum
# from empirical_privacy import lsdd
# reload(one_bit_sum)

def B_pmf(k, n, p):
    return binom(n, p).pmf(k)
def B0_pmf(k, n, p):
    return B_pmf(k, n-1, p)
def B1_pmf(k, n, p):
    return B_pmf(k-1, n-1, p)
def sd(N, P):
    return 0.5*np.sum(abs(B0_pmf(i, N, P) - B1_pmf(i, N, P)) for i in range(N+1))
def optimal_correctness(n, p):
    return 0.5 + 0.5*sd(n, p)
n_max = 2**10
ntri=30
n=7
p=0.5
sd(n,p)
B0 = [B0_pmf(i, n, p) for i in range(n+1)]
B1 = [B1_pmf(i, n, p) for i in range(n+1)]
dif = np.abs(np.array(B0)-np.array(B1))
sdv = 0.5*np.sum(dif)
pc = 0.5+0.5*sdv
print(f'n={n} coin flips p={p} probability of heads'\
      '\nB0 has first outcome=0, B1 has first outcome=1')
print(f'Statistic is the total number of heads sum')
print(f'N_heads=\t{"     ".join(np.arange(n+1).astype(str))}')
print(f'PMF of B0=\t{B0}\nPMF of B1=\t{B1}')
print(f'|B0-B1|=\t{dif}')
print(f'sd = 0.5 * sum(|B0-B1|) = {sdv}')
print(f'P(Correct) = 0.5 + 0.5*sd = {pc}')
ccc_kwargs = {
    'confidence_interval_width':10,
    'n_max':2**13,
    'dataset_settings' : {
        'n_trials':n,
        'prob_success':p,
        'gen_distr_type':'binom'
    },
    'validation_set_size' : 2000
}

CCCs = []
Fits = ['knn', 'density', 'expectation']
for fit in Fits:
    CCCs.append(build_convergence_curve_pipeline(
                    GenSampleOneBitSum,
        gensample_kwargs = {'generate_in_batch':True},
        fitter=fit,
        fitter_kwargs={} if fit=='knn' else {'statistic_column':0}
                )(**ccc_kwargs)
   )
luigi.build(CCCs, local_scheduler=True, workers=4, log_level='ERROR')

colors = cm.Accent(np.linspace(0,1,len(CCCs)+1))

ax = plt.figure(figsize=(10,5))
ax = plt.gca()
leg_handles = []
for (i, CC) in enumerate(CCCs):
    with CC.output().open() as f:
        res = pickle.load(f)
    handle=sns.tsplot(res['sd_matrix'], ci='sd', color=colors[i], ax=ax, legend=False, time=res['training_set_sizes'])

j=0
for i in range(len(CCCs), 2*len(CCCs)):
    handle.get_children()[i].set_label('{}'.format(Fits[j]))
    j+=1
plt.semilogx()
plt.axhline(optimal_correctness(n, p), linestyle='--', color='r', label='_nolegend_')
plt.axhline(0.5, linestyle='-', color='b', label='_nolegend_')
plt.title('n={n} p={p} $\delta$={d:.3f}'.format(n=n, p=p, d=sd(n,p)), fontsize=20)
plt.xlabel('num samples')
plt.ylabel('Correctness Rate')
plt.legend(loc=(0,1.1))<jupyter_output>/opt/conda/lib/python3.6/site-packages/seaborn/timeseries.py:183: UserWarning: The `tsplot` function is deprecated and will be removed in a future release. Please update your code to use the new `lineplot` function.
  warnings.warn(msg, UserWarning)
/opt/conda/lib/python3.6/site-packages/seaborn/timeseries.py:183: UserWarning: The `tsplot` function is deprecated and will be removed in a future release. Please update your code to use the new `lineplot` function.
  warnings.warn(msg, UserWarning)
/opt/conda/lib/python3.6/site-packages/seaborn/timeseries.py:183: UserWarning: The `tsplot` function is deprecated and will be removed in a future release. Please update your code to use the new `lineplot` function.
  warnings.warn(msg, UserWarning)<jupyter_text>Repeat the above using joblib to make sure the luigi implementation is correct<jupyter_code>from math import ceil, log
one_bit_sum.n_jobs=1
N = int(ceil(log(n_max) / log(2)))
N_samples = np.logspace(4,N,num=N-3, base=2).astype(np.int)

ax = plt.figure(figsize=(10,5))
ax = plt.gca()
AlgArg = namedtuple('AlgArg', field_names=['f_handle', 'f_kwargs'])
algs = [
    AlgArg(one_bit_sum.get_knn_correctness_rate_cached, {'neighbor_method':'sqrt'}),
    AlgArg(one_bit_sum.get_knn_correctness_rate_cached, {'neighbor_method':'sqrt_random_tiebreak'}),

    AlgArg(one_bit_sum.get_density_est_correctness_rate_cached, {'bandwidth_method':None}),
    AlgArg(one_bit_sum.get_expectation_correctness_rate_cached, {'bandwidth_method':None}),
    AlgArg(one_bit_sum.get_lsdd_correctness_rate_cached, {})
    #AlgArg(one_bit_sum.get_knn_correctness_rate_cached, {'neighbor_method':'cv'})
]
colors = cm.Accent(np.linspace(0,1,len(algs)+1))
leg_handles = []
for (i,alg) in enumerate(algs):
    res = one_bit_sum.get_res(n,p,ntri, alg.f_handle, alg.f_kwargs, n_max=n_max)
    handle=sns.tsplot(res, ci='sd', color=colors[i], ax=ax, legend=False, time=N_samples)
#     f, coef = get_fit(res, N_samples)
#     print alg, coef
#     lim = coef[0]
#     plt.plot(N_samples, f(N_samples), linewidth=3)
#    plt.text(N_samples[-1], lim, '{:.3f}'.format(lim),fontsize=16)
j=0
for i in range(len(algs), 2*len(algs)):
    #print i, i/2-1 if i%2==0 else (i)/2
    handle.get_children()[i].set_label('{} {}'.format(algs[j].f_handle.func.__name__, algs[j].f_kwargs))
    j+=1
plt.semilogx()
plt.axhline(optimal_correctness(n, p), linestyle='--', color='r', label='_nolegend_')
plt.axhline(0.5, linestyle='-', color='b', label='_nolegend_')
plt.title('n={n} p={p} $\delta$={d:.3f}'.format(n=n, p=p, d=sd(n,p)), fontsize=20)
plt.xlabel('num samples')
plt.ylabel('Correctness Rate')
plt.legend(loc=(0,1.1))
#print ax.get_legend_handles_labels()<jupyter_output>[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.5s remaining:    0.0s
[Parallel(n_jobs=1)]: Done   2 out of   2 | elapsed:    1.3s remaining:    0.0s
[Parallel(n_jobs=1)]: Done   3 out of   3 | elapsed:    2.0s remaining:    0.0s
[Parallel(n_jobs=1)]: Done   3 out of   3 | elapsed:    2.0s finished
/opt/conda/lib/python3.6/site-packages/seaborn/timeseries.py:183: UserWarning: The `tsplot` function is deprecated and will be removed in a future release. Please update your code to use the new `lineplot` function.
  warnings.warn(msg, UserWarning)
[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.7s remaining:    0.0s
[Parallel(n_jobs=1)]: Done   2 out of   2 | elapsed:    1.3s remaining:    0.0s
[Parallel(n_jobs=1)]: Done   3 out of   3 | elapsed:    2.1s remaining:    0.0s
[Parallel(n_jobs=1)]: Done   3 out of   3 | elapsed:    2.1s finished
/opt/conda/lib/python3.6/site-packages/seaborn/timeseries.py:183: UserWarning: The `tsplot` function is deprecated and will be rem[...]<jupyter_text>Timing GenSamples Without halving: 7.5secWith halving: 8.1sec (i.e. not much overhead)<jupyter_code>from luigi_utils.sampling_framework import GenSamples
import time

class GS(GenSamples(GenSampleOneBitSum, generate_in_batch=True)):
    pass
GSi = GS(dataset_settings = ccc_kwargs['dataset_settings'],
  random_seed='0',
  generate_positive_samples=True,
  num_samples=2**15)

start = time.time()
luigi.build([GSi], local_scheduler=True, workers=8, log_level='ERROR')
cputime = time.time() - start
print(cputime)
res['training_set_sizes'].shape
np.concatenate((np.array([]), np.array([1,2,3])))<jupyter_output><empty_output><jupyter_text>More exp<jupyter_code>def get_fit(res, N_samples):
    ntri, nsamp = res.shape
    sqrt2 = np.sqrt(2)
    Xlsq = np.hstack((np.ones((nsamp,1)),
                      sqrt2/(N_samples.astype(np.float)**0.25)[:, np.newaxis]))
    y = 1.0 - res.reshape((nsamp*ntri, 1))
    
    Xlsq = reduce(lambda x,y: np.vstack((x,y)), [Xlsq]*ntri)
    coef = np.linalg.lstsq(Xlsq, y)[0].ravel()
    f = lambda n: 1.0 - coef[0] - coef[1]*sqrt2/n.astype(np.float)**0.25, coef
    return f

trial=0
num_samples=2**11
bandwidth_method=None
from scipy.stats import gaussian_kde

X0, X1, y0, y1 = one_bit_sum.gen_data(n, p, num_samples, trial)
X0 = X0.ravel()
X1 = X1.ravel()
bw = None
if hasattr(bandwidth_method, '__call__'):
    bw = float(bandwidth_method(num_samples)) / num_samples  # eg log
if type(bandwidth_method) == float:
    bw = num_samples**(1-bandwidth_method)
f0 = gaussian_kde(X0, bw_method = bw)
f1 = gaussian_kde(X1, bw_method = bw)
#Omega = np.unique(np.concatenate((X0, X1)))
_min = 0
_max = n
x = np.linspace(_min, _max, num=10*num_samples)
print('difference of densities=',0.5 + 0.5 * 0.5 * np.mean(np.abs(f0(x)-f1(x))))
denom = f0(x)+f1(x)
numer = np.abs(f0(x)-f1(x))
print('expectation = ',0.5 + 0.5*np.mean(numer/denom))<jupyter_output><empty_output><jupyter_text>Uniforml distributed random variables $$g_0 = U[0,0.5]+\sum_{i=1}^{n-1} U[0,1]$$$$g_1 = U[0.5,1.0]+\sum_{i=1}^{n-1} U[0,1]$$ Let $\mu_n = \frac{n-1}{2}$ and $\sigma_n = \sqrt{\frac{n-0.75}{12}}$By the CLT $g_0\sim N(\mu_n+0.25, \sigma_n)$ and $g_1\sim N(\mu_n+0.75, \sigma_n)$.<jupyter_code>from math import sqrt
n=3
x = np.linspace(n/2.0-sqrt(n), n/2.0+sqrt(n))
sigma = sqrt((n-0.75)/12.0)
sqrt2 = sqrt(2)
mu = (n-1.0)/2

def g0_pdf(x):
    return norm.pdf(x, loc=mu+0.25, scale=sigma)
def g1_pdf(x):
    return norm.pdf(x, loc=mu+0.75, scale=sigma)
def d_pdf(x):
    return norm.pdf(x, loc=-0.5, scale=sigma*sqrt2)
def g_int(n):
    sigma = sqrt((n-0.75)/12.0)
    mu = (n-1.0)/2
    N0 = norm(loc=mu+0.25, scale=sigma)
    N1 = norm(loc=mu+0.75, scale=sigma)
    I0 = N0.cdf(n*0.5)-N0.cdf(0)
    I1 = N1.cdf(n*0.5)-N1.cdf(0)
    return 2*(I0-I1)
def g_stat_dist(n):
    return 0.5 * g_int(n)
def g_optimal_correctness(n):
    return 0.5 + 0.5*g_stat_dist(n)

plt.plot(x, g0_pdf(x), label='$g_0$')
plt.plot(x, g1_pdf(x), label='$g_1$')
#plt.plot(x, d_pdf(x), label='$d$')
plt.axvline(x=n/2.0, color='r')
assert g0_pdf(n/2.0)==g1_pdf(n/2.0)
plt.legend()
print(g_optimal_correctness(n))
from math import ceil, log
if n_max >= 2**13:
    one_bit_sum.n_jobs=1
else:
    one_bit_sum.n_jobs=-1
    
N = int(ceil(log(n_max) / log(2)))
N_samples = np.logspace(4,N,num=N-3, base=2).astype(np.int)

ax = plt.figure(figsize=(10,5))
ax = plt.gca()
AlgArg = namedtuple('AlgArg', field_names=['f_handle', 'f_kwargs'])
algs = [
    AlgArg(one_bit_sum.get_knn_correctness_rate_cached, {'neighbor_method':'sqrt'}),
    AlgArg(one_bit_sum.get_knn_correctness_rate_cached, {'neighbor_method':'sqrt_random_tiebreak'}),

    AlgArg(one_bit_sum.get_density_est_correctness_rate_cached, {'bandwidth_method':None}),
    AlgArg(one_bit_sum.get_expectation_correctness_rate_cached, {'bandwidth_method':None}),
    AlgArg(one_bit_sum.get_lsdd_correctness_rate_cached, {})
    #AlgArg(one_bit_sum.get_knn_correctness_rate_cached, {'neighbor_method':'cv'})
]
for A in algs:
    A.f_kwargs['type']='norm'
colors = cm.Accent(np.linspace(0,1,len(algs)+1))
leg_handles = []
for (i,alg) in enumerate(algs):
    res = one_bit_sum.get_res(n,p,ntri, alg.f_handle, alg.f_kwargs, n_max=n_max)
    handle=sns.tsplot(res, ci='sd', color=colors[i], ax=ax, legend=False, time=N_samples)
#     f, coef = get_fit(res, N_samples)
#     print alg, coef
#     lim = coef[0]
#     plt.plot(N_samples, f(N_samples), linewidth=3)
#     plt.text(N_samples[-1], lim, '{:.3f}'.format(lim),fontsize=16)
j=0
for i in range(len(algs), 2*len(algs)):
    #print i, i/2-1 if i%2==0 else (i)/2
    handle.get_children()[i].set_label(algs[j].f_handle.func.__name__)
    j+=1
    #print handle.get_children()[i].get_label()
plt.semilogx()
plt.axhline(g_optimal_correctness(n), linestyle='--', color='r', label='_nolegend_')
plt.axhline(0.5, linestyle='-', color='b', label='_nolegend_')
plt.title('n={n} $\delta$={d:.3f}'.format(n=n, d=g_stat_dist(n)), fontsize=20)
plt.xlabel('num samples')
plt.ylabel('Correctness Rate')
plt.legend(loc=(1.1,0))
#print ax.get_legend_handles_labels()
true_value = g_optimal_correctness(n)
print(true_value)
trial=0
num_samples=2**15
bandwidth_method=None
from scipy.stats import gaussian_kde

X0, X1, y0, y1 = one_bit_sum.gen_data(n, p, num_samples, trial, type='norm')
X0 = X0.ravel()
X1 = X1.ravel()
bw = None
if hasattr(bandwidth_method, '__call__'):
    bw = float(bandwidth_method(num_samples)) / num_samples  # eg log
if type(bandwidth_method) == float:
    bw = num_samples**(1-bandwidth_method)
f0 = gaussian_kde(X0, bw_method = bw)
f1 = gaussian_kde(X1, bw_method = bw)
#Omega = np.unique(np.concatenate((X0, X1)))
_min = 0
_max = n
x = np.linspace(_min, _max, num=num_samples)

print('difference of densities=',0.5 + 0.5 * 0.5 * integrate.quad(lambda x: np.abs(f0(x)-f1(x)), -np.inf, np.inf)[0])
X = np.concatenate((X0,X1))
f0x = f0(X)
f1x = f1(X)
denom = (f0x+f1x+np.spacing(1))
numer = np.abs(f0x-f1x)
print('expectation = ',0.5 + 0.5*np.mean(numer/denom))
print('exact=',g_optimal_correctness(n))
plt.plot(x, f0(x),label='$\hat g_0$', linestyle='--')
plt.plot(x, f1(x),label='$\hat g_1$', linestyle='--')
plt.plot(x, g0_pdf(x), label='$g_0$')
plt.plot(x, g1_pdf(x), label='$g_1$')
plt.legend(loc=(1.05,0))<jupyter_output><empty_output><jupyter_text>Comparing different numerical integration techniques<jupyter_code>to_int = [f0,f1]

print 'Quad'
# for (i,f) in enumerate(to_int):
#     intr = integrate.quad(f, -np.inf, np.inf)
#     print 'func={0} err={1:.3e}'.format(i, abs(1-intr[0]))
g_int(n)-integrate.quad(lambda x: np.abs(f0(x)-f1(x)), -np.inf, np.inf)[0]
to_int = [f0,f1]

print 'Quad'
g_int(n)-integrate.quad(lambda x: np.abs(f0(x)-f1(x)), -np.inf, np.inf)[0]
g_int(n)
print 'Simps'
def delta(x):
    return np.abs(f0(x)-f1(x))
X = np.unique(np.concatenate((X0,X1)))
y = delta(X)
g_int(n)-integrate.simps(y,X)
import empirical_privacy.lsdd
rtv = lsdd.lsdd(X0[np.newaxis, :], X1[np.newaxis, :])
plt.hist(rtv[1])
np.mean(rtv[1])<jupyter_output><empty_output><jupyter_text>Sympy-based analysis<jupyter_code>import sympy as sy
n,k = sy.symbols('n k', integer=True)
#k = sy.Integer(k)
p = sy.symbols('p', real=True)
q=1-p

def binom_pmf(k, n, p):
    return sy.binomial(n,k)*(p**k)*(q**(n-k))
def binom_cdf(x, n, p):
    return sy.Sum([binom_pmf(j, n, p) for j in sy.Range(x+1)])

B0 = binom_pmf(k, n-1, p)
B1 = binom_pmf(k-1, n-1, p)
def stat_dist(N,P):
    return 0.5*sum([sy.Abs(B0.subs([(n,N),(p,P), (k,i)])-B1.subs([(n,N),(p,P), (k,i)])) for i in range(N+1)])
def sd(N, P):
    return 0.5*np.sum(abs(B0(i, N, P) - B1(i, N, P)) for i in range(N+1))
stat_dist(50,0.5)
sd(5000,0.5)
N=2
terms =[(B0.subs([(n,N), (k,i)]).simplify(),B1.subs([(n,N), (k,i)]).simplify()) for i in range(N+1)]
print terms
0.5*sum(map(lambda t: sy.Abs(t[0]-t[1]), terms)).subs([(p,0.5)])
stat_dist(4,0.5)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18936830231246474, 'fraction_numerical': 0.03715454032712916, 'mean_word_length': 5.018243529910904, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 28, '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': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '668401', 'license': 'MIT', 'n_tokens_mistral': 6454, 'n_tokens_neox': 6220, 'n_words': 1107, 'path': 'Notebooks/development/1-bit sum.ipynb', 'repo_name': 'maksimt/empirical_privacy', 'size': '14184'}