<jupyter_start><jupyter_text>Comparison with CVX<jupyter_code>class ThresholdTest:
    @staticmethod
    def cvxlb(data, wmin, wmax, theta, rmin=0, rmax=1):
        import cvxpy as cp
        import numpy as np
        from collections import defaultdict
        from scipy.special import xlogy
        from pprint import pformat
        
        assert wmin >= 0
        assert wmin < 1
        assert wmax > 1
        
        crdict = defaultdict(int)
        n = 0
        clogc = 0
        for (ci, wi, ri) in data:
            assert ci >= 0
            assert wi >= wmin and wi <= wmax
            assert ri >= rmin and ri <= rmax
            if ci > 0:
                crdict[(wi, ri)] += ci
            n += ci
            clogc += xlogy(ci,ci)
        del ci, wi, ri
        assert n >= 1
        crdict.default_factory = None
        
        wvec = np.array(list(sorted(set(w for (w, _), _ in crdict.items()))))
        wmaxvec = np.max(wvec)
        rvec = np.array(list(sorted(set(r for (_, r), _ in crdict.items()))))
        
        C = np.array([ [ crdict.get((w, r), 0) for r in rvec ] for w in wvec ])
        W = np.array([ [ w for r in rvec ] for w in wvec ])
        WR = np.array([ [ w*r for r in rvec ] for w in wvec ])
                
        gamma = cp.Variable()
        beta = cp.Variable()
        
        ex = list(set(
            [ (w, r) 
              for (w, r) in [(wmin, rmin), (wmin, rmax), (wmax, rmin), (wmax, rmax) ]
            ]
        ))
        
        constraints = [ gamma + beta * w + (w * r) >= 0 for (w, r) in ex ]
        objective = (
#            cp.exp((theta-clogc)/n + cp.sum(cp.multiply((C/n),cp.log(alpha + beta * W + WR))))
#             cp.exp((theta-clogc)/n) * cp.geo_mean(x=cp.vec(alpha + beta * W + WR), p=C.flatten(order='F')/n)
            n * cp.exp((theta-clogc)/n) * cp.geo_mean(x=cp.vec(gamma + beta * W + WR), 
                                                      p=C.astype('int').flatten(order='F'),
                                                      max_denom=1<<24)
            -gamma
            -beta
        )
        prob = cp.Problem(cp.Maximize(objective), constraints)        
        prob.solve(solver='ECOS')
        assert prob.status == 'optimal', pformat({ 'status': prob.status, 'data': data })
        gammastar = np.asscalar(gamma.value)
        betastar = np.asscalar(beta.value)
        from math import exp, log
        kappastar = (prob.value + gammastar + betastar) / n
        
        return prob.value, {
            'gammastar': gammastar,
            'betastar': betastar,
            'kappastar': kappastar,
        }        
    
    @staticmethod
    def cvxub(data, wmin, wmax, theta, rmin=0, rmax=1):
        import cvxpy as cp
        import numpy as np
        from collections import defaultdict
        from scipy.special import xlogy
        from pprint import pformat
        
        assert wmin >= 0
        assert wmin < 1
        assert wmax > 1
        
        crdict = defaultdict(int)
        n = 0
        clogc = 0
        sumcw = 0
        for (ci, wi, ri) in data:
            assert ci >= 0
            assert wi >= wmin and wi <= wmax
            assert ri >= rmin and ri <= rmax
            if ci > 0:
                crdict[(wi, ri)] += ci
            n += ci
            clogc += xlogy(ci,ci)
            sumcw += ci*wi
        del ci, wi, ri
        assert n >= 1
        crdict.default_factory = None
        
        wvec = np.array(list(sorted(set(w for (w, _), _ in crdict.items()))))
        wmaxvec = np.max(wvec)
        rvec = np.array(list(sorted(set(r for (_, r), _ in crdict.items()))))
        
        C = np.array([ [ crdict.get((w, r), 0) for r in rvec ] for w in wvec ])
        W = np.array([ [ w for r in rvec ] for w in wvec ])
        WR = np.array([ [ w*r for r in rvec ] for w in wvec ])
                
        gamma = cp.Variable()
        beta = cp.Variable()
        
        ex = list(set(
            [ (w, r) 
              for (w, r) in [(wmin, rmin), (wmin, rmax), (wmax, rmin), (wmax, rmax) ]
            ]
        ))
        
        constraints = [ gamma + beta * w - (w * r) >= 0 for (w, r) in ex ]
        objective = (
#            cp.exp((theta-clogc)/n + cp.sum(cp.multiply((C/n),cp.log(alpha + beta * W + WR))))
#             cp.exp((theta-clogc)/n) * cp.geo_mean(x=cp.vec(alpha + beta * W - WR), p=C.flatten(order='F')/n)
            n * cp.exp((theta-clogc)/n) * cp.geo_mean(x=cp.vec(gamma + beta * W - WR), 
                                                      p=C.astype('int').flatten(order='F'),
                                                      max_denom=1<<24)
            -gamma
            -beta
        )
        prob = cp.Problem(cp.Maximize(objective), constraints)        
        prob.solve(solver='ECOS')
        assert prob.status == 'optimal', pformat({ 'status': prob.status, 'data': data })
        gammastar = np.asscalar(gamma.value)
        betastar = np.asscalar(beta.value)
        from math import exp, log
        kappastar = (prob.value + gammastar + betastar) / n

        return -prob.value, {
            'gammastar': gammastar,
            'betastar': betastar,
            'kappastar': kappastar,
        } 

def testasymptoticci():
    from importlib import reload
    from math import ceil
    import environments.ControlledRangeVariance
    import MLE.MLE
    import numpy as np
    from tqdm import tqdm as tqdm
    import sys
    from pprint import pformat
    from scipy.stats import f
    from scipy.special import xlogy
    
    reload(environments.ControlledRangeVariance)
    reload(MLE.MLE)
    
    wsupport = [ 0, 2, 10 ]
    wmin, wmax = wsupport[0], wsupport[-1]
    env = environments.ControlledRangeVariance.ControlledRangeVariance(seed=45, wsupport=wsupport, expwsq=5)
    
    total = 0
    badcvx = 0

    for ndata in tqdm(map(ceil, np.logspace(1, 7, 14)), file=sys.stderr, total=14):
        for i in range(11):
            (truevalue, data) = env.sample(ndata)
            vhat, qmle = MLE.MLE.estimate(datagen=lambda: data, wmin=wmin, wmax=wmax)
            likelihoodmle = sum(xlogy(c, qmle['qfunc'](c, w, r)) for c, w, r in data if c > 0)
            theta = likelihoodmle - 0.5 * f.isf(q=0.05, dfn=1, dfd=sum(c for c,_,_ in data)-1)

            (lb, ub), (qlb, qub) = MLE.MLE.asymptoticconfidenceinterval(datagen=lambda: data, wmin=wmin, wmax=wmax, 
                                                                        raiseonerr=True)
            
            total += 1
            try:
                cvxlb, cvxqlb = ThresholdTest.cvxlb(data, wmin, wmax, theta)
                cvxub, cvxqub = ThresholdTest.cvxub(data, wmin, wmax, theta)
            except:
                badcvx += 1
                continue

            assert (np.allclose(cvxlb, lb, atol=5e-3) 
                    and np.allclose(cvxub, ub, atol=5e-3)
                   ), pformat(
            {
                'data': [(c, w, r) for c, w, r in data if c > 0],
                'wmin': wmin,
                'wmax': wmax,
                'cvxlb': cvxlb,
                'cvxub': cvxub,
                'cvxqlb': cvxqlb,
                'cvxqub': cvxqub,
                'lb': lb,
                'ub': ub,
                'qlb': qlb,
                'qub': qub,
            })
    
    return { 'badcvx': badcvx, 'total': total }
    
testasymptoticci()<jupyter_output>100%|██████████| 14/14 [00:30<00:00,  2.57s/it]<jupyter_text>Comparison with Baselines<jupyter_code>class CITest:
    class ClopperPearson:
        @staticmethod
        def interval(data, wmax, alpha, **kwargs):    
            from scipy.special import betaincinv

            k = sum((c*w/wmax)*r for c, w, r in data)
            n = sum(c for c, w, _ in data)

            lb = betaincinv(k, n - k + 1, alpha/2) if k > 0 and n - k + 1 > 0 else 0
            ub = betaincinv(k + 1, n - k, 1-alpha/2) if n > k else 1
            
            lb = min(1, max(0, wmax * lb))
            ub = min(1, max(0, wmax * ub))

            return min(lb, ub), max(lb, ub)
        
    class AsymptoticGaussian:
        def interval(data, wmax, alpha, **kwargs):
            from math import sqrt
            from scipy.stats import norm  
            
            scale = norm.ppf(1-alpha/2)
 
            n = sum(c for c, _, _ in data)
            tot = sum(c*w*r for c, w, r in data)
            muhat = tot / n
            totsq = sum(c*(w*r - muhat)**2 for c, w, r in data)
            sigmahat = sqrt(totsq) / n
            
            lb = min(1, max(0, muhat - scale*sigmahat))
            ub = min(1, max(0, muhat + scale*sigmahat))
            
            return min(lb, ub), max(lb, ub)

for data in [
    [(2, 1.0101010101010102, 1), (8, 0, 0)],
    [(9, 1.0101010101010102, 0), (1, 1.0101010101010102, 1)],
    [(1, 0, 0), (99, 1000/999, 1)],
    [(21, 0.0, 0), (247, 1.0101010101010102, 0), (1, 1000.0, 0)],
    [(10, 2, 0)]
]:
    from importlib import reload
    import MLE.MLE

    reload(MLE.MLE)
    
    print("****")
    print(CITest.ClopperPearson.interval(data, wmax=2000, alpha=0.05))
    print(CITest.AsymptoticGaussian.interval(data, wmax=2000, alpha=0.05))
    print(MLE.MLE.asymptoticconfidenceinterval(datagen=lambda: data, wmin=0, wmax=2000, alpha=0.05))
from importlib import reload
import environments.ControlledRangeVariance
import MLE.MLE

reload(environments.ControlledRangeVariance)
reload(MLE.MLE)

def produceresults(env, method, minexp=1, maxexp=6, numpts=15, ndataperpt=100):
    from math import ceil
    import numpy
    
    wmin, wmax = env.range()

    for ndata in map(ceil, numpy.logspace(minexp, maxexp, numpts)):
        bounds=[]
        pc=0
        for i in range(1, ndataperpt+1):   
            (truevalue, data) = env.sample(ndata)

            try:
                ci = method(data=data, wmin=wmin, wmax=wmax, alpha=0.05)
                assert numpy.isfinite(ci[0]) and numpy.isfinite(ci[1])
                assert ci[1] >= ci[0]
                assert ci[1] <= 1 + 1e-4
                assert ci[0] >= -1e-4
#                 assert ci[1] >= truevalue
#                 assert ci[0] <= truevalue
            except:
                print('truevalue was {}'.format(truevalue))
                print('data was {}'.format(data))
                print('wmin, wmax was {} {}'.format(wmin, wmax))
                print('ci was {}'.format(ci))
                raise
                                                
            bounds.append(
                ( truevalue,
                  data,
                  [ 
                      ci[0], 'n/a', ci[1]
                  ]
                )  
            )
            
        uppers = [ d[2][2] for d in bounds ]
        lowers = [ d[2][0] for d in bounds ]
        upperbounded = [ 1 if d[0] <= d[2][2] else 0 for d in bounds ]
        lowerbounded = [ 1 if d[2][0] <= d[0] else 0 for d in bounds ]
        covered = [ 1 if u*l > 0 else 0 for (u, l) in zip(upperbounded, lowerbounded) ]
        width = [ d[2][2] - d[2][0] for d in bounds ]    
        yield (ndata,
                { 
                    'cov': numpy.mean(covered),
                    'covstd': numpy.std(covered, ddof=1) / numpy.sqrt(len(covered) - 1),
                    'ubcov': numpy.mean(upperbounded),
                    'lbcov': numpy.mean(lowerbounded),
                    'width': numpy.mean(width),
                    'widthstd': numpy.std(width, ddof=1) / numpy.sqrt(len(width) - 1),
                    'widthlo': numpy.quantile(width, q=[0.05])[0],
                    'widthhi': numpy.quantile(width, q=[0.95])[0],
                    'ub': numpy.mean(uppers),
                    'lb': numpy.mean(lowers),
                },
              )

def getenv():
    wsupport = [ 0, 2, 1000 ]
    env = environments.ControlledRangeVariance.ControlledRangeVariance(seed=45, wsupport=wsupport, expwsq=100)
    return env, env.getpw(), env.range(), env.expectedwsq()

allres = []
for (name, method) in [ 
        ('MLE CI', lambda data, **kwargs: MLE.MLE.asymptoticconfidenceinterval(datagen=lambda: data, **kwargs)[0]),
        ('Asymptotic Gaussian', lambda data, wmax, alpha, **kwargs: CITest.AsymptoticGaussian.interval(data, wmax, alpha)),
        ('Binomial CI', lambda data, wmax, alpha, **kwargs: CITest.ClopperPearson.interval(data, wmax, alpha)),
    ]:
    print('****** {} ******'.format(name))
    res = []
    for zzz in produceresults(getenv()[0], method, ndataperpt=10000, numpts=20):
        wmax = getenv()[2][1]
        res.append( (zzz[0] / wmax, zzz[1]) )
        print('{}'.format(( zzz[0], { y: round(z, 4) for y, z in zzz[1].items() })), flush=True)
        del wmax
    allres.append((name, res))
import pickle
pickle.dump( allres, open( "epsilongreedy_ci_allres_cvxopt.p", "wb" ) )
%matplotlib inline

class FlassPlot:
    from matplotlib import scale as mscale
    
    class ExpScale(mscale.ScaleBase):
        from matplotlib import transforms as mtransforms

        """
        Scales data in range -Infinity, 1 using log(1 - x)
        """

        # The scale class must have a member ``name`` that defines the
        # string used to select the scale.  For example,
        # ``gca().set_yscale("mercator")`` would be used to select this
        # scale.
        name = 'exp'

        def __init__(self, axis, **kwargs):
            """
            Any keyword arguments passed to ``set_xscale`` and
            ``set_yscale`` will be passed along to the scale's
            constructor.
            """
            from matplotlib import scale as mscale
            myscale = kwargs.pop("expscale", 1.0)
            mscale.ScaleBase.__init__(self, axis, **kwargs)
            self.myscale = myscale

        def get_transform(self):
            """
            Override this method to return a new instance that does the
            actual transformation of the data.
            """
            return FlassPlot.ExpScale.ExpTransform(self.myscale)

        def set_default_locators_and_formatters(self, axis):
            """
            Override to set up the locators and formatters to use with the
            scale.  This is only required if the scale requires custom
            locators and formatters.  Writing custom locators and
            formatters is rather outside the scope of this example, but
            there are many helpful examples in ``ticker.py``.

            In our case, the Mercator example uses a fixed locator from
            -90 to 90 degrees and a custom formatter class to put convert
            the radians to degrees and put a degree symbol after the
            value::
            """

            return axis

        def limit_range_for_scale(self, vmin, vmax, minpos):
            """
            Override to limit the bounds of the axis to the domain of the
            transform.
            """
            return vmin, vmax

        class ExpTransform(mtransforms.Transform):
            # There are two value members that must be defined.
            # ``input_dims`` and ``output_dims`` specify number of input
            # dimensions and output dimensions to the transformation.
            # These are used by the transformation framework to do some
            # error checking and prevent incompatible transformations from
            # being connected together.  When defining transforms for a
            # scale, which are, by definition, separable and have only one
            # dimension, these members should always be set to 1.
            input_dims = 1
            output_dims = 1
            is_separable = True

            def __init__(self, scale):
                from matplotlib import transforms as mtransforms
                
                mtransforms.Transform.__init__(self)
                self.scale = scale

            def transform_non_affine(self, a):
                """
                This transform takes an Nx1 ``numpy`` array and returns a
                transformed copy.  Since the range of the Mercator scale
                is limited by the user-specified threshold, the input
                array must be masked to contain only valid values.
                ``matplotlib`` will handle masked arrays and remove the
                out-of-range data from the plot.  Importantly, the
                ``transform`` method *must* return an array that is the
                same shape as the input array, since these values need to
                remain synchronized with values in the other dimension.
                """
                import numpy as np
                
                return np.exp(self.scale * a)

            def inverted(self):
                """
                Override this method so matplotlib knows how to get the
                inverse transform for this transform.
                """
                return FlassPlot.ExpScale.InvertedExpTransform(self.scale)

        class InvertedExpTransform(mtransforms.Transform):
            input_dims = 1
            output_dims = 1
            is_separable = True

            def __init__(self, scale):
                from matplotlib import transforms as mtransforms
                
                mtransforms.Transform.__init__(self)
                self.scale = scale

            def transform_non_affine(self, a):
                import numpy as np
                
                return np.log(a / self.scale)

            def inverted(self):
                return FlassPlot.ExpScale.ExpTransform(self.scale)

    @staticmethod
    def yo():
        import matplotlib.pyplot as plt
        plt.figure(figsize=(4,2))
        
    @staticmethod
    def forpaper():
        import matplotlib.pyplot as plt
        from matplotlib import scale as mscale

        mscale.register_scale(FlassPlot.ExpScale)

        SMALL_SIZE = 16
        MEDIUM_SIZE = 22
        BIGGER_SIZE = 24

        plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
        plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
        plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
        plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
        plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
        plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
        plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title
        
    @staticmethod
    def exppic(x, y, label, expscale, **kwargs):
        import matplotlib.pyplot as plt
        import matplotlib.ticker as ticker
        plt.ylim(0.0, 1.01)
        plt.plot(x, y, label=label, **kwargs)
        plt.gca().set_xscale('log')
        plt.gca().set_yscale('exp', expscale=expscale)
        plt.ylim(0.0, 1.01)
        plt.gca().yaxis.set_major_locator(ticker.FixedLocator([0.38,0.68,0.86,0.95,1]))

        plt.legend()
        
    @staticmethod
    def pic(x, y, label): 
        import matplotlib.pyplot as plt
        plt.semilogx(x, y, label=label)

        #plt.loglog(x, y, label=label)
        plt.legend()
        
    @staticmethod
    def axeslabel(xlabel, ylabel):  
        import matplotlib.pyplot as plt
        plt.xlabel(xlabel)
        plt.ylabel(ylabel)
        
    @staticmethod
    def title(title):
        import matplotlib.pyplot as plt
        plt.title(title)
        
    @staticmethod
    def legend(**kwargs):
        import matplotlib.pyplot as plt
        plt.legend(**kwargs)
        
    @staticmethod
    def savefig(filename):
        import matplotlib.pyplot as plt
        plt.savefig(filename, bbox_inches='tight')
        
    @staticmethod
    def plt():
        import matplotlib.pyplot as plt
        return plt

import pickle
allres = pickle.load(open( "epsilongreedy_ci_allres_cvxopt.p", "rb" ) )

renameit = { "Binomial CI": "Binomial", "Asymptotic Gaussian": "Gaussian", "MLE CI": "EL" }
FlassPlot.forpaper()
for name, res in allres:
    x = [ x[0] for x in res ]
    y = [ x[1]['width'] for x in res ]
    ylo = [ x[1]['widthlo'] for x in res ]
    yhi = [ x[1]['widthhi'] for x in res ]
    FlassPlot.plt().semilogx(x, y, label=renameit.get(name, name))
    FlassPlot.plt().fill_between(x, ylo, yhi, alpha=0.15)

FlassPlot.legend(loc="center", bbox_to_anchor=[0.8,0.825])
FlassPlot.axeslabel('n / $w_{max}$', 'CI width')
FlassPlot.plt().savefig("epsilongreedy_ci.pdf", bbox_inches='tight')
%matplotlib inline

import pickle
allres = pickle.load(open( "epsilongreedy_ci_allres_cvxopt.p", "rb" ) )

renameit = { "Binomial CI": "Binomial", "Asymptotic Gaussian": "Gaussian", "MLE CI": "EL" }
FlassPlot.forpaper()
for name, res in allres:
    x = [ x[0] for x in res ]
    y = [ x[1]['cov'] for x in res ]
    ylo = [ x[1]['cov'] - 1.96 * x[1]['covstd'] for x in res ]
    yhi = [ x[1]['cov'] + 1.96 * x[1]['covstd'] for x in res ]
    FlassPlot.exppic(x, y, label=renameit.get(name, name), expscale=7.0)
    FlassPlot.plt().fill_between(x, ylo, yhi, alpha=0.25)
FlassPlot.exppic([ x[0] for x in res ], [ 0.95 for x in res ], label='_nolegend_', expscale=7.0, ls='dashed')
FlassPlot.legend(loc="center", bbox_to_anchor=[0.5,0.425])
FlassPlot.axeslabel('n / $w_{max}$', 'empirical coverage')
FlassPlot.plt().savefig("epsilongreedy_coverage_ci.pdf", bbox_inches='tight')
from importlib import reload
import environments.ControlledRangeVariance
import MLE.MLE

reload(environments.ControlledRangeVariance)
reload(MLE.MLE)

def megaproduceresults(env, method):
    from math import ceil
    import numpy
    
    wmin, wmax = env.range()

    for ndata in [100000]:
        bounds=[]
        for i in range(1, 101):
            (truevalue, data) = env.sample(ndata)

            try:
                ci = method(data=data, wmin=wmin, wmax=wmax, alpha=0.05)
                assert numpy.isfinite(ci[0]) and numpy.isfinite(ci[1])
#                 assert ci[1] >= ci[0]
#                 assert ci[1] >= truevalue
#                 assert ci[0] <= truevalue
            except:
                print('truevalue was {}'.format(truevalue))
                print('data was {}'.format(data))
                print('ci was {}'.format(ci))
                raise
                                                
            bounds.append(
                ( truevalue,
                  data,
                  [ 
                      ci[0], 'n/a', ci[1]
                  ]
                )  
            )
            
        uppers = [ d[2][2] for d in bounds ]
        lowers = [ d[2][0] for d in bounds ]
        upperbounded = [ 1 if d[0] <= d[2][2] else 0 for d in bounds ]
        lowerbounded = [ 1 if d[2][0] <= d[0] else 0 for d in bounds ]
        covered = [ 1 if u*l > 0 else 0 for (u, l) in zip(upperbounded, lowerbounded) ]
        width = [ d[2][2] - d[2][0] for d in bounds ]    
        yield (ndata,
                { 
                    'coverage': round(sum(covered) / len(covered), 4),
                    'ubcoverage': round(sum(upperbounded) / len(upperbounded), 4),
                    'lbcoverage': round(sum(lowerbounded) / len(lowerbounded), 4),
                    'width': round(sum(width) / len(width), 4),
                    'avub': round(sum(uppers) / len(uppers), 4),
                    'avlb': round(sum(lowers) / len(lowers), 4),
                },
              )

def getenv(expwsq):
    import numpy
    wsupport = [ 0, 2, 1000 ]
    env = environments.ControlledRangeVariance.ControlledRangeVariance(seed=45, wsupport=wsupport, expwsq=expwsq)
    return env, env.getpw(), env.range(), env.expectedwsq()
    
allres = []
for (name, method) in [ ('Binomial', lambda data, wmax, alpha, **kwargs: CITest.ClopperPearson.interval(data, wmax, alpha)),
                        ('MLE', lambda data, **kwargs: MLE.MLE.asymptoticconfidenceinterval(datagen=lambda: data, **kwargs)[0]),
                      ]:
    import numpy
    print('****** {} ******'.format(name))
    res = []
    for expwsq in numpy.linspace(100, 999, 10):
        for zzz in megaproduceresults(getenv(expwsq)[0], method):
            res.append((expwsq, zzz))
            print('{}'.format((expwsq, zzz)), flush=True)
    allres.append((name, res))
import pickle
pickle.dump( allres, open( "epsilongreedy_ci_expwsqsweep.p", "wb" ) )
import pickle
allres = pickle.load(open( "epsilongreedy_ci_expwsqsweep.p", "rb" ) )

%matplotlib inline
        
class MegaFlassPlot:
    @staticmethod
    def forpaper():
        import matplotlib.pyplot as plt

        SMALL_SIZE = 16
        MEDIUM_SIZE = 22
        BIGGER_SIZE = 24

        plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
        plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
        plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
        plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
        plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
        plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
        plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title
        
    @staticmethod
    def pic(x, y, label):
        import matplotlib.pyplot as plt
        
        plt.plot(x, y, label=label)
        plt.legend()
        
    @staticmethod
    def axeslabel(xlabel, ylabel):
        import matplotlib.pyplot as plt
        
        plt.xlabel(xlabel)
        plt.ylabel(ylabel)
        
    @staticmethod
    def title(title):
        import matplotlib.pyplot as plt
        
        plt.title(title)
        
    @staticmethod
    def savefig(filename):
        import matplotlib.pyplot as plt
        
        plt.savefig(filename, bbox_inches='tight')       

MegaFlassPlot.forpaper()
for name, res in allres:
    MegaFlassPlot.pic([ x[0] for x in res], [ x[1][1]['width'] for x in res], label=name)
MegaFlassPlot.axeslabel('$E[w^2]$', 'CI width')
from importlib import reload
import environments.Chirpy
import MLE.MLE

reload(environments.Chirpy)
reload(MLE.MLE)

def megaproduceresults(env, method):
    from math import ceil
    import numpy
    
    wmin, wmax = env.range()

    for ndata in [10000]:
        bounds=[]
        for i in range(1, 101):
            (truevalue, data) = env.sample(ndata)

            try:
                ci = method(data=data, wmin=wmin, wmax=wmax, alpha=0.05)
                expwsq = sum(c*w*w for c, w, _ in data) / sum(c for c, _, _ in data)
                assert numpy.isfinite(ci[0]) and numpy.isfinite(ci[1])
#                 assert ci[1] >= ci[0]
#                 assert ci[1] >= truevalue
#                 assert ci[0] <= truevalue
            except:
                print('truevalue was {}'.format(truevalue))
                print('data was {}'.format(data))
                print('ci was {}'.format(ci))
                raise
                                                
            bounds.append(
                ( truevalue,
                  data,
                  [ 
                      ci[0], 'n/a', ci[1]
                  ],
                  expwsq
                )  
            )
            
        uppers = [ d[2][2] for d in bounds ]
        lowers = [ d[2][0] for d in bounds ]
        upperbounded = [ 1 if d[0] <= d[2][2] else 0 for d in bounds ]
        lowerbounded = [ 1 if d[2][0] <= d[0] else 0 for d in bounds ]
        covered = [ 1 if u*l > 0 else 0 for (u, l) in zip(upperbounded, lowerbounded) ]
        width = [ d[2][2] - d[2][0] for d in bounds ]
        expwsqds = [d[3] for d in bounds]
        yield (ndata,
                { 
                    'coverage': round(sum(covered) / len(covered), 4),
                    'ubcoverage': round(sum(upperbounded) / len(upperbounded), 4),
                    'lbcoverage': round(sum(lowerbounded) / len(lowerbounded), 4),
                    'width': round(sum(width) / len(width), 4),
                    'avub': round(sum(uppers) / len(uppers), 4),
                    'avlb': round(sum(lowers) / len(lowers), 4),
                    'avexpwsq': round(sum(expwsqds) / len(expwsqds), 4),
                },
              )

def getenv(numimp):
    import numpy
    wsupport = [ 0, 2, 1000 ]
    env = environments.Chirpy.Chirpy(seed=45, wsupport=wsupport, expwsq=200, nchirp=numimp // 3)
    return env, env.getpw(), env.range(), env.expectedwsq()
    
allres = []
for (name, method) in [ ('Binomial', lambda data, wmax, alpha, **kwargs: CITest.ClopperPearson.interval(data, wmax, alpha)),
                        ('MLE', lambda data, **kwargs: MLE.MLE.asymptoticconfidenceinterval(datagen=lambda: data, **kwargs)[0]),
                      ]:
    import numpy
    print('****** {} ******'.format(name))
    res = []
    for numimp in numpy.ceil(numpy.geomspace(10, 10000, 10)):
        for zzz in megaproduceresults(getenv(numimp)[0], method):
            res.append((numimp, zzz))
            print('{}'.format((numimp, zzz)), flush=True)
    allres.append((name, res))
import pickle
pickle.dump( allres, open( "epsilongreedy_ci_numimpsweep.p", "wb" ) )
import pickle
allres = pickle.load(open( "epsilongreedy_ci_numimpsweep.p", "rb" ) )

%matplotlib inline
        
class MegaFlassPlot:
    @staticmethod
    def forpaper():
        import matplotlib.pyplot as plt

        SMALL_SIZE = 16
        MEDIUM_SIZE = 22
        BIGGER_SIZE = 24

        plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
        plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
        plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
        plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
        plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
        plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
        plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title
        
    @staticmethod
    def pic(x, y, label):
        import matplotlib.pyplot as plt
        
        plt.semilogx(x, y, label=label)
        plt.legend()
        
    @staticmethod
    def axeslabel(xlabel, ylabel):
        import matplotlib.pyplot as plt
        
        plt.xlabel(xlabel)
        plt.ylabel(ylabel)
        
    @staticmethod
    def title(title):
        import matplotlib.pyplot as plt
        
        plt.title(title)
        
    @staticmethod
    def savefig(filename):
        import matplotlib.pyplot as plt
        
        plt.savefig(filename, bbox_inches='tight')

MegaFlassPlot.forpaper()
for name, res in allres:
    MegaFlassPlot.pic([ x[0] for x in res], [ x[1][1]['width'] for x in res], label=name)
MegaFlassPlot.axeslabel('# unique importance weights', 'CI width')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12128817758075557, 'fraction_numerical': 0.0185142931647055, 'mean_word_length': 1.7949348638061402, 'pattern_counts': {'":': 6, '<': 29, '<?xml version=': 0, '>': 39, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '365828', 'license': 'Unlicense', 'n_tokens_mistral': 10307, 'n_tokens_neox': 9653, 'n_words': 2558, 'path': 'ci.ipynb', 'repo_name': 'pmineiro/elfcb', 'size': '30679'}