<commit_before>#!/usr/bin/env python

import json
import signal
from time import time

import numpy as np
from scipy.stats import chi2
import matplotlib.pyplot as plt

from BinaryMatrix import approximate_conditional_nll as cond_a_nll_b
from BinaryMatrix import approximate_from_margins_weights as cond_a_sample_b
from Confidence import invert_test, ci_conservative_generic
from Utility import logsumexp, logabsdiffexp
from Experiment import Seed

# Parameters
params = { 'fixed_example': None,#'data/rasch_covariates.json',
           'M': 10,
           'N': 5,
           'theta': 2.0,
           'kappa': -1.628,
           'alpha_min': -0.4,
           'beta_min': -0.86,
           'v_min': -0.6,
           'alpha_level': 0.05,
           'n_MC_levels': [],#[10, 50, 100, 500],
           'wopt_sort': False,
           'is_T': 50,
           'n_rep': 3,
           'L': 601,
           'theta_l': -6.0,
           'theta_u': 6.0,
           'do_prune': False,
           'random_seed': 137,
           'verbose': True }

terminated = False
def sigint_handler(signum, frame):
    print 'Terminating after current trial completes.'
    global terminated
    terminated = True
signal.signal(signal.SIGINT, sigint_handler)

def cond_a_nll(X, w):
    return cond_a_nll_b(X, w, sort_by_wopt_var = params['wopt_sort'])

def cond_a_sample(r, c, w, T = 0):
    return cond_a_sample_b(r, c, w, T, sort_by_wopt_var = params['wopt_sort'])

# Generates (or loads) a particular realization of P_{ij} and then
# repeatedly samples independent Bernoulli random variables according
# to these cell probabilities.
def generate_data(params, seed):
    # Advance random seed for parameter and covariate construction
    seed.next()

    if not params['fixed_example']:
        # Generate parameters and covariates
        M, N = params['M'], params['N']
        alpha = np.random.uniform(size = (M,1)) + params['alpha_min']
        beta = np.random.uniform(size = (1,N)) + params['beta_min']
        kappa = params['kappa']
        v = np.random.uniform(size = (M,N)) + params['v_min']
    else:
        # Load parameters and covariates
        with open(params['fixed_example'], 'r') as example_file:
            example = json.load(example_file)

            v = np.array(example['nu'])
            M, N = v.shape

            alpha = np.array(example['alpha']).reshape((M,1))
            beta = np.array(example['beta']).reshape((1,N))
            kappa = example['kappa']

    # Generate Bernoulli probabilities from logistic regression model
    logit_P = np.zeros((M,N)) + kappa
    logit_P += alpha
    logit_P += beta
    logit_P += params['theta'] * v
    P = 1.0 / (1.0 + np.exp(-logit_P))

    while True:
        # Advance random seed for data generation
        seed.next()

        # Generate data for this trial
        X = np.random.random((M,N)) < P

        # Pruning rows and columns of 0's and 1's; this may improve
        # the quality of the approximation for certain versions of the
        # sampler
        if params['do_prune']:
            X_p = X.copy()
            v_p = v.copy()
            while True:
                r, c = X_p.sum(1), X_p.sum(0)
                r_p = (r == 0) + (r == N)
                c_p = (c == 0) + (c == M)
                pruning = np.any(r_p) or np.any(c_p)

                X_p = X_p[-r_p][:,-c_p]
                v_p = v_p[-r_p][:,-c_p]

                if not pruning:
                    break

            yield X_p.copy(), v_p.copy()

        yield X, v

def timing(func):
    def inner(*args, **kwargs):
        start_time = time()
        val = func(*args, **kwargs)
        elapsed = time() - start_time
        return val, elapsed

    return inner

def plot_statistics(ax, theta_grid, test_val, crit):
    # Compute confidence interval from test statistics
    ci_l, ci_u = invert_test(theta_grid, test_val, crit)
    ci_l = max(ci_l, params['theta_l'])
    ci_u = min(ci_u, params['theta_u'])

    ax.plot(theta_grid, test_val, color = 'b')
    ax.hlines(crit, theta_grid[0], theta_grid[-1], linestyle = 'dotted')
    ax.hlines(crit, ci_l, ci_u, color = 'r')
    ax.hlines(2.0 * crit, theta_grid[0], theta_grid[-1],
              color = 'w', linewidth = 9, zorder = 99)
    ax.hlines(2.0 * crit, ci_l, ci_u, color = 'r', linewidth = 9, zorder = 100,
              alpha = 1.0 / ((1 - params['alpha_level']) * params['n_rep']))
    ax.vlines(ci_l, 2.0 * crit, crit, color = 'r', linestyle = 'dotted')
    ax.vlines(ci_u, 2.0 * crit, crit, color = 'r', linestyle = 'dotted')
    ax.set_ylim(2.0 * crit, 0)

# Set up plots
fig_cmle_a, ax_cmle_a = plt.subplots()
fig_cmle_is, ax_cmle_is = plt.subplots()

@timing
def ci_cmle_a(X, v, theta_grid, alpha_level):
    crit = -0.5 * chi2.ppf(1 - alpha_level, 1)

    cmle_a = np.empty_like(theta_grid)
    for l, theta_l in enumerate(theta_grid):
        logit_P_l = theta_l * v
        cmle_a[l] = -cond_a_nll(X, np.exp(logit_P_l))

    plot_statistics(ax_cmle_a, theta_grid, cmle_a - cmle_a.max(), crit)<commit_msg>Final (hopefully) bug fix for coverage alpha...<commit_after>#!/usr/bin/env python

import json
import signal
from time import time

import numpy as np
from scipy.stats import chi2
import matplotlib.pyplot as plt

from BinaryMatrix import approximate_conditional_nll as cond_a_nll_b
from BinaryMatrix import approximate_from_margins_weights as cond_a_sample_b
from Confidence import invert_test, ci_conservative_generic
from Utility import logsumexp, logabsdiffexp
from Experiment import Seed

# Parameters
params = { 'fixed_example': 'data/rasch_covariates.json',
           'M': 20,
           'N': 10,
           'theta': 2.0,
           'kappa': -1.628,
           'alpha_min': -0.4,
           'beta_min': -0.86,
           'v_min': -0.6,
           'alpha_level': 0.05,
           'n_MC_levels': [10, 50, 100, 500],
           'wopt_sort': False,
           'is_T': 50,
           'n_rep': 100,
           'L': 601,
           'theta_l': -6.0,
           'theta_u': 6.0,
           'do_prune': False,
           'random_seed': 137,
           'verbose': True }

terminated = False
def sigint_handler(signum, frame):
    print 'Terminating after current trial completes.'
    global terminated
    terminated = True
signal.signal(signal.SIGINT, sigint_handler)

def cond_a_nll(X, w):
    return cond_a_nll_b(X, w, sort_by_wopt_var = params['wopt_sort'])

def cond_a_sample(r, c, w, T = 0):
    return cond_a_sample_b(r, c, w, T, sort_by_wopt_var = params['wopt_sort'])

# Generates (or loads) a particular realization of P_{ij} and then
# repeatedly samples independent Bernoulli random variables according
# to these cell probabilities.
def generate_data(params, seed):
    # Advance random seed for parameter and covariate construction
    seed.next()

    if not params['fixed_example']:
        # Generate parameters and covariates
        M, N = params['M'], params['N']
        alpha = np.random.uniform(size = (M,1)) + params['alpha_min']
        beta = np.random.uniform(size = (1,N)) + params['beta_min']
        kappa = params['kappa']
        v = np.random.uniform(size = (M,N)) + params['v_min']
    else:
        # Load parameters and covariates
        with open(params['fixed_example'], 'r') as example_file:
            example = json.load(example_file)

            v = np.array(example['nu'])
            M, N = v.shape

            alpha = np.array(example['alpha']).reshape((M,1))
            beta = np.array(example['beta']).reshape((1,N))
            kappa = example['kappa']

    # Generate Bernoulli probabilities from logistic regression model
    logit_P = np.zeros((M,N)) + kappa
    logit_P += alpha
    logit_P += beta
    logit_P += params['theta'] * v
    P = 1.0 / (1.0 + np.exp(-logit_P))

    while True:
        # Advance random seed for data generation
        seed.next()

        # Generate data for this trial
        X = np.random.random((M,N)) < P

        # Pruning rows and columns of 0's and 1's; this may improve
        # the quality of the approximation for certain versions of the
        # sampler
        if params['do_prune']:
            X_p = X.copy()
            v_p = v.copy()
            while True:
                r, c = X_p.sum(1), X_p.sum(0)
                r_p = (r == 0) + (r == N)
                c_p = (c == 0) + (c == M)
                pruning = np.any(r_p) or np.any(c_p)

                X_p = X_p[-r_p][:,-c_p]
                v_p = v_p[-r_p][:,-c_p]

                if not pruning:
                    break

            yield X_p.copy(), v_p.copy()

        yield X, v

def timing(func):
    def inner(*args, **kwargs):
        start_time = time()
        val = func(*args, **kwargs)
        elapsed = time() - start_time
        return val, elapsed

    return inner

def plot_statistics(ax, theta_grid, test_val, crit):
    # Compute confidence interval from test statistics
    ci_l, ci_u = invert_test(theta_grid, test_val, crit)
    ci_l = max(ci_l, params['theta_l'])
    ci_u = min(ci_u, params['theta_u'])
    cov_alpha = min(1.0, 1.0 / ((1 - params['alpha_level']) * params['n_rep']))

    ax.plot(theta_grid, test_val, color = 'b')
    ax.hlines(crit, theta_grid[0], theta_grid[-1], linestyle = 'dotted')
    ax.hlines(crit, ci_l, ci_u, color = 'r')
    ax.hlines(2.0 * crit, theta_grid[0], theta_grid[-1],
              color = 'w', linewidth = 9, zorder = 99)
    ax.hlines(2.0 * crit, ci_l, ci_u, color = 'r', linewidth = 9, zorder = 100,
              alpha = cov_alpha)
    ax.vlines(ci_l, 2.0 * crit, crit, color = 'r', linestyle = 'dotted')
    ax.vlines(ci_u, 2.0 * crit, crit, color = 'r', linestyle = 'dotted')
    ax.set_ylim(2.0 * crit, 0)

# Set up plots
fig_cmle_a, ax_cmle_a = plt.subplots()
fig_cmle_is, ax_cmle_is = plt.subplots()

@timing
def ci_cmle_a(X, v, theta_grid, alpha_level):
    crit = -0.5 * chi2.ppf(1 - alpha_level, 1)

    cmle_a = np.empty_like(theta_grid)
    for l, theta_l in enumerate(theta_grid):
        logit_P_l = theta_l * v
        cmle_a[l] = -cond_a_nll(X, np.exp(logit_P_l))

    plot_statistics(ax_cmle_a, theta_grid, cmle_a - cmle_a.max(), crit)
{'fraction_non_alphanumeric': 0.131722293750619, 'fraction_numerical': 0.01792611666831732, 'mean_word_length': 2.428862478777589, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, '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}}
{'dir': 'git-commits-cleaned', 'id': '600231', 'n_tokens_mistral': 3709, 'n_tokens_neox': 3483, 'n_words': 900}