# Path: Utility.py
# def l2(x):
# def apply_scale(w, a, b):
# def logabsdiffexp(x, y):
#     def f(x, y):
# def logit_mean(x):
# def digest(x):
# def tree():
# def lift_tree(src, dst):
# def init_latex_rendering():
# def draw_confidence(a, m, S, levels = [0.95]):


import numpy as np
from Utility import logsumexp
        log_Q_X[l] = log_likelihood(X, theta_l)
        for k in range(K):
            log_Q_Y[l,k] = log_likelihood(Y[k], theta_l)
        if verbose:
            print '%.2f: %.2g, %.2g' % \
              (theta_l, np.exp(log_Q_X[l]), np.exp(log_Q_Y[l].max()))
    log_Q_sum_X = logsumexp(log_Q_X)
    log_Q_sum_Y = np.empty(K)
    for k in range(K):
        log_Q_sum_Y[k] = logsumexp(log_Q_Y[:,k])

    # Step over the grid, calculating approximate p-values
    if two_sided:
        log_p_plus = np.empty(L)
    log_p_minus = np.empty(L)
    for l in range(L):
        theta_l = theta_grid[l]
        log_w_l = np.empty(K + 1)

        # X contribution
        if corrected:
            log_w_l[K] = theta_l * (suff * X).sum() - log_Q_sum_X
        else:
            log_w_l[K] = -np.inf

        # Y contribution
        for k in range(K):
            log_w_l[k] = theta_l * (suff * Y[k]).sum() - log_Q_sum_Y[k]

        if two_sided:

# Path: Utility.py
# def l2(x):
# def apply_scale(w, a, b):
# def logabsdiffexp(x, y):
#     def f(x, y):
# def logit_mean(x):
# def digest(x):
# def tree():
# def lift_tree(src, dst):
# def init_latex_rendering():
# def draw_confidence(a, m, S, levels = [0.95]):


# Path: Confidence.py
import numpy as np
from Utility import logsumexp

        log_Q_X[l] = log_likelihood(X, theta_l)
        for k in range(K):
            log_Q_Y[l,k] = log_likelihood(Y[k], theta_l)
        if verbose:
            print '%.2f: %.2g, %.2g' % \
              (theta_l, np.exp(log_Q_X[l]), np.exp(log_Q_Y[l].max()))
    log_Q_sum_X = logsumexp(log_Q_X)
    log_Q_sum_Y = np.empty(K)
    for k in range(K):
        log_Q_sum_Y[k] = logsumexp(log_Q_Y[:,k])

    # Step over the grid, calculating approximate p-values
    if two_sided:
        log_p_plus = np.empty(L)
    log_p_minus = np.empty(L)
    for l in range(L):
        theta_l = theta_grid[l]
        log_w_l = np.empty(K + 1)

        # X contribution
        if corrected:
            log_w_l[K] = theta_l * (suff * X).sum() - log_Q_sum_X
        else:
            log_w_l[K] = -np.inf

        # Y contribution
        for k in range(K):
            log_w_l[k] = theta_l * (suff * Y[k]).sum() - log_Q_sum_Y[k]

        if two_sided:

            log_p_num_plus = logsumexp(log_w_l[I_t_Y_plus[l]])