<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Create dataset
Select city
Step1:  Read data
Use the following to remove an existing store file, if needed. (Use C-M y to make the cell runnable).
rm -f /u/gciampag/data/tbdc15/trip_accidents_store.hdf
Step2:  <hr style="height
Step3:  Scatter plot of trips vs trip starts
Distributions of accidents
Load the data
Step4:  Histogram grouping data $>0$ in bins of size $9$
Step5:  Fit the $>0$ data with an exponential law (with binning), and geometric distribution
Step6:  Zero-inflated geometric distribution
Step7:  Simulate from a zero-inflated geometric
Step8:  Compare the trip distribution for cells with $0$ accidents with a normal cell (with $\ge 0$ accidents)
<ASSISTANT_TASK:>
Python Code:
cities = ['bari', 'milano', 'napoli', 'palermo', 'roma', 'torino', 'venezia']
location = 'gfs'
boxes_location = 'data'
store_location = 'data'
scale = 1000.0
Explanation: Create dataset
Select city
End of explanation
ll -h ~/data/tbdc15/
for city in cities:
    logging.info("City: {}".format(city))
    make_city_frame(city, 
                    loc=location, 
                    scale=scale, 
                    boxesloc=boxes_location,
                    storeloc=store_location, 
                    break_at=None)
Explanation: Read data
Use the following to remove an existing store file, if needed. (Use C-M y to make the cell runnable).
rm -f /u/gciampag/data/tbdc15/trip_accidents_store.hdf
End of explanation
%matplotlib inline
import numpy
from pylab import *
# scatter plot
def scatter_trips_accidents(df, city, use_trips_starts=False):
    fig = figure()
    if use_trips_starts:
        xcol = 'trips_start'
    else:
        xcol = 'trips'
            
    df.plot(x=xcol, y='accidents', kind='scatter', marker='x', alpha=.2, color='k', fig=fig)
    # trend line
    emin = numpy.log10(df[xcol].min())
    emax = numpy.log10(df[xcol].max())
    bins = numpy.logspace(max(emin, 0), emax, 20)
    print bins
    df.groupby(numpy.digitize(df[xcol], bins=bins)).mean().plot(x=xcol, y='accidents',
                                                                color='r', linestyle='solid',
                                                                marker='o', ax=gca(), alpha=.5, 
                                                                linewidth=2, fig=fig)
    grid('off')
    title(city)
    if use_trips_starts:
        xlabel('Traffic (start of trip)')
    else:
        xlabel('Traffic')
    ylabel('Accidents')
    xscale('log')
    xlim(1, xlim()[1])
    tight_layout()
    legend()
    savefig('trips_accidents_scatter_{}.pdf'.format(city))
    savefig('trips_accidents_scatter_{}.png'.format(city))
    show()
def hist_accidents(df, city):
    fig = figure()
    ax = gca()
    ax.hist(df['accidents'].values, log=True, bins=60, color='white')
    ylim(.1, ylim()[1])
    xlabel('Accidents')
    ylabel('Frequency')
    title(city)
    tight_layout()
    legend()
    savefig('accidents_histogram_{}.pdf'.format(city))
    savefig('accidents_histogram_{}.png'.format(city))
    show()
    
    
def plot_all(city):
    paths = getpaths(city, 
                     loc=location, 
                     boxesloc=boxes_location,
                     storeloc=store_location)
    df = pandas.read_hdf(paths['store'], city)
    df = df.groupby(['i', 'j']).sum().filter(['trips', 'trips_start', 'accidents'])
    scatter_trips_accidents(df, city)
    scatter_trips_accidents(df, city, use_trips_starts=True)
    hist_accidents(df, city)
plot_all('bari')
plot_all('milano')
plot_all('napoli')
plot_all('palermo')
plot_all('roma')
plot_all('torino')
plot_all('venezia')
Explanation: <hr style="height: .1em"/>
Plot the data
End of explanation
city = 'palermo'
df = pandas.read_hdf('/u/gciampag/data/tbdc15/trip_accidents_store.hdf', city)
df = df.groupby(['i', 'j']).sum().filter(['trips', 'trips_start', 'accidents'])
Explanation: Scatter plot of trips vs trip starts
Distributions of accidents
Load the data
End of explanation
bin_size = 9  # lower bound on bin size
max_accidents = df['accidents'].max()
start = 1
stop = 1 + ceil((max_accidents - 1) / bin_size) * bin_size
num = (stop - start) / bin_size + 1
bins = numpy.linspace(start, stop, num, endpoint=True)
bins = numpy.hstack([[0,], bins])
nh, __, ___ = hist(df['accidents'].values, bins=bins, color='lightgray', log=True, normed=True, histtype='bar')
xlim(-5, stop)
ylim(1e-7, 1)
xlabel('Accidents')
ylabel('Frequency')
title(city.title())
tick_params(axis='both', direction='out', which='both')
tick_params(axis='x', which='minor', bottom='on', top='off')
tick_params(axis='y', which='both', right='off')
tick_params(axis='x', which='major', top='off')
tick_params(axis='x', which='minor', bottom='on')
Explanation: Histogram grouping data $>0$ in bins of size $9$
End of explanation
from scipy.stats import expon, geom, poisson
group_size = 9
df_nza = df.query('accidents > 0')
a = df_nza.groupby(ceil(df_nza['accidents'] / group_size)).count()['accidents']
x = a.index.values
p = a / a.sum()
vlines(x, 0, p, color='LightGray')
plot(x, p, 'wo ', label='Data', ms=8)
# expected number of accidents (computed as a weighted average of the frequencies)
exp_accidents = np.sum(p.values * a.index.values)
#x = np.hstack([[0]])
rv = expon(loc=0, scale=exp_accidents ** -1)
plot(x, rv.cdf(x + 1) - rv.cdf(x), 'x ', color='k', mew=2, label='Exponential')
rv = geom(exp_accidents ** -1, loc=0)
plot(x, rv.pmf(x), '+ ', color='gray', mew=1.5, ms=10, label='Geometric')
rv = poisson(exp_accidents ** -1, loc=0)
plot(x, rv.pmf(x), marker=(6, 2, 0), ls=' ', color='LightGray', mew=1, ms=10, label='Poisson')
xlim(0, xlim()[1] + 1)
xlabel(r'$\left\lceil\rm{Accidents} \,/\, %d\right\rceil$' % group_size, fontsize='large')
ylabel('Probability')
yscale('log')
ylim(ylim()[0], 2)
title(city.title())
legend(loc='best', frameon=False)
tick_params(axis='both', direction='out', which='both')
tick_params(axis='x', which='minor', bottom='on', top='off')
tick_params(axis='y', which='both', right='off')
tick_params(axis='x', which='major', top='off')
tick_params(axis='x', which='minor', bottom='on')
savefig("accidents_grouped_fit_{}.png".format(city))
Explanation: Fit the $>0$ data with an exponential law (with binning), and geometric distribution
End of explanation
import numpy
from scipy.stats import rv_discrete
class zigeom_gen(rv_discrete):
    def _pmf(self, k, pi, p):
        s = numpy.sign(k)
        return (1 - s) * (pi + (1.0 - pi) * p) + s * (1.0 - pi) * (1.0 - p) ** k * p
            
zigeom = zigeom_gen()
Explanation: Zero-inflated geometric distribution
End of explanation
from scipy.optimize import minimize
def fit(data):
    def _llk(args):
        return - zigeom(*args).logpmf(data).sum()
    N = float(len(data))
    pi0 = (data == 0).sum() / N
    x0 = (pi0, .5)
    ret = minimize(_llk, x0, method='Nelder-Mead')
    if ret['success']:
        return ret['x']
    else:
        raise RuntimeError(ret['message'])
pi = .2
p = .3
data = zigeom(pi, p).rvs(size=1000)
print fit(data)
p = 0.1
pi = 0.5
data = zigeom(pi, p).rvs(size=2000)
data_max = data.max()
hist(data, bins=data_max, color='white', log=True, normed=True)
pih, ph = fit(data)
x = np.arange(data_max)
px = zigeom(pih, ph).pmf(x)
plot(x + .5, px, 'r-')
title('$\pi = {}, p = {}$'.format(pi, p))
Explanation: Simulate from a zero-inflated geometric
End of explanation
import numpy
from truthy_measure.plotting import plot_pdf_log2, plot_cdf
trips_all = df['trips'].values
trips_zac = df.query('accidents == 0')['trips'].values
num_points = 20
bins = numpy.logspace(0, numpy.log2(trips_all.max()), num=num_points, base=2)
hist_all, __ = numpy.histogram(trips_all, bins=bins, normed=True)
hist_zac, __ = numpy.histogram(trips_zac, bins=bins, normed=True)
plot(bins[1:], numpy.log(hist_zac) - numpy.log(hist_all), 'ko:', mfc='LightGray')
axhline(0, ls='--', color='gray', lw=2)
xscale('log')
xlabel('Trips $t$')
ylabel('$\log\Pr\{T = t | A = 0\} - \log\Pr\{T\}$')
yb = max(numpy.abs(ylim()))
ylim(-yb, yb)
title(city.title())
tight_layout()
savefig("logratio_{}.png".format(city))
Explanation: Compare the trip distribution for cells with $0$ accidents with a normal cell (with $\ge 0$ accidents)
End of explanation
<END_TASK>