#!/usr/bin/env python3 # 04 Fabian Ying 2016-11-13 Added discrete example # 03 Keith Briggs 2016-10-25 # 02: Added ccdf option to plot function # Switched order for spans, the variable for CI (now: [Lower bound, Best estimate, Upper bound]) # Added a new plot function with shading # 01: Added best estimate point for each quantile (3rd column of spans). Added plot cdf tool # Keith Briggs 2016-08-08 was quantile_estimation_02.py # see slide 25 in ~/KTP_TVWS_Project/tvws_ktn_london_slides-0.2/Spectrum_sharing_deliverable_2013-10-11.pdf import numpy as np import matplotlib.pyplot as plt from sys import exit,stderr,stdout from math import sqrt from random import random from scipy.stats import binom class Quantile_estimator: def __init__(self,f,**f_params): ''' Instantiate a Quantile_estimator: f_params are extra keywords arguments passed to the sampling function f, which must also take a size keyword argument ''' self.f=f self.f_params=f_params def _predict_indices_for_confidence_interval(self,sample_size,quantile,conf=0.95): # internal use only # cf. ~/KTP_TVWS_software/quantsim/src/binom.cc rv=binom(sample_size,quantile) r,s=rv.interval(conf) return max(0,int(r)),min(sample_size-1,int(s)) def estimate(self,quantiles=(0.5,),conf=0.95,eps_abs=1e-3,eps_rel=1e-2,batch_size=1000,max_n_batches=1000): ''' Take samples from callable object f until all quantiles are sufficiently accurately determined, and return the confidence intervals as rows of an np.array, and also the number of samples needed. The default quantiles=(0.5,) finds the median. ''' if type(quantiles) is float: quantiles=[quantiles,] for quantile in quantiles: if not 0.0x\$]',img_fn_base,ccdf) print(spans,'\nsample size needed was %d \n'%sample_size) # Continuous distributions for name,fun in continuous_tests: print(fmt%(100.0*conf,str(quantiles),name,)) qe=Quantile_estimator(fun,scale=scale) spans,sample_size=qe.estimate(quantiles,conf,eps_rel=eps_rel) print(spans,'\nsample size needed was %d \n'%sample_size) img_fn_base='%s/qe_%s_%s'%(save_path,name,('cdf','ccdf')[ccdf],) qe.plot_cdf(spans,quantiles,False,'\$x\$','Prob[\$X>x\$]',img_fn_base,ccdf) # Discrete distributions name='Uniform (discrete)' print(fmt%(100.0*conf,str(quantiles),name,)) qe=Quantile_estimator(np.random.choice, a=a) spans,sample_size=qe.estimate(quantiles,conf,eps_rel=eps_rel) print(spans,'\nsample size needed was %d \n'%sample_size) img_fn_base='%s/qe_uniform_discrete_%s'%(save_path,('cdf','ccdf')[ccdf],) qe.plot_cdf(spans,quantiles,False,'\$x\$','Prob[\$X>x\$]',img_fn_base,ccdf) name = 'Geometric' print(fmt % (100.0*conf, str(quantiles), name,)) qe = Quantile_estimator(np.random.geometric, p=p) spans, sample_size = qe.estimate(quantiles, conf, eps_rel=eps_rel) print(spans, '\nsample size needed was %d \n' % sample_size) img_fn_base = '%s/qe_geometric_%s_%s' % (save_path, 'p=%d' % p, ('cdf', 'ccdf')[ccdf],) qe.plot_cdf(spans, quantiles, False, '\$x\$', 'Prob[\$X>x\$]', img_fn_base, ccdf) name = 'Poisson' for lam in [lam1, lam2]: print(fmt % (100.0*conf, str(quantiles), name,)) qe = Quantile_estimator(np.random.poisson, lam=lam) # spans,sample_size=qe.estimate(quantiles,conf,eps_abs=2.0) # absolute tolerance spans, sample_size = qe.estimate(quantiles, conf, eps_rel=eps_rel) print(spans, '\nsample size needed was %d \n' % sample_size) img_fn_base = '%s/qe_poisson_%s_%s' % (save_path, 'lam=%d' % lam, ('cdf', 'ccdf')[ccdf],) qe.plot_cdf(spans, quantiles, False, '\$x\$', 'Prob[\$X>x\$]', img_fn_base, ccdf) #qe.plot_cdf_alt(spans,quantiles,exponential_quantiles=exponential_quantiles,savefig=True,img_filename=img_filename_alt,ccdf=ccdf) # non-uniformly-spaces quantiles... # quantiles_exponential=1.0-np.power(10.0,np.arange(-30.0,-5.0,5.0)/10.0) #quantiles_exponential=np.power(10.0,np.arange(-30.0,-5.0,2.0)/10.0) test_02()