<jupyter_start><jupyter_text>Likelihood of different distributions<jupyter_code>def loglikelihood_ratio(loglikelihoods1, loglikelihoods2,
        nested=False, normalized_ratio=False):
    """
    Calculates a loglikelihood ratio and the p-value for testing which of two
    probability distributions is more likely to have created a set of
    observations.

    Parameters
    ----------
    loglikelihoods1 : list or array
        The logarithms of the likelihoods of each observation, calculated from
        a particular probability distribution.
    loglikelihoods2 : list or array
        The logarithms of the likelihoods of each observation, calculated from
        a particular probability distribution.
    nested : bool, optional
        Whether one of the two probability distributions that generated the
        likelihoods is a nested version of the other. False by default.
    normalized_ratio : bool, optional
        Whether to return the loglikelihood ratio, R, or the normalized
        ratio R/sqrt(n*variance)

    Returns
    -------
    R : float
        The loglikelihood ratio of the two sets of likelihoods. If positive,
        the first set of likelihoods is more likely (and so the probability
        distribution that produced them is a better fit to the data). If
        negative, the reverse is true.
    p : float
        The significance of the sign of R. If below a critical value
        (typically .05) the sign of R is taken to be significant. If above the
        critical value the sign of R is taken to be due to statistical
        fluctuations.
    """
    from numpy import sqrt
    from scipy.special import erfc

    n = float(len(loglikelihoods1))

    if n==0:
        R = 0
        p = 1
        return R, p
    from numpy import asarray
    loglikelihoods1 = asarray(loglikelihoods1)
    loglikelihoods2 = asarray(loglikelihoods2)

    #Clean for extreme values, if any
    from numpy import inf, log
    from sys import float_info
    min_val = log(10**float_info.min_10_exp)
    loglikelihoods1[loglikelihoods1==-inf] = min_val
    loglikelihoods2[loglikelihoods2==-inf] = min_val

    R = sum(loglikelihoods1-loglikelihoods2)

    from numpy import mean
    mean_diff = mean(loglikelihoods1)-mean(loglikelihoods2)
    variance = sum(
            ( (loglikelihoods1-loglikelihoods2) - mean_diff)**2
            )/n

    if nested:
        from scipy.stats import chi2
        p = 1 - chi2.cdf(abs(2*R), 1)
    else:
        p = erfc( abs(R) / sqrt(2*n*variance))

    if normalized_ratio:
        R = R/sqrt(n*variance)

    return R, p

def find_best_fit(func,v):
    '''
    Find the best initial guess for fitting distribution in scipy.
    '''
    print(func)
    if func.shapes: 
        n=len(func.shapes)
    else:
        n = 0
    
    
    def fit_once(n):
        in_guesses = [np.random.rand()*(5*np.log10(max(v))) for i in range(2+n)]
        par = func.fit(v,*in_guesses[:-2], loc = in_guesses[-1],scale=in_guesses[-2])
        l = sum(func.logpdf(v,*par))
        return l,par
        
        
    
    res = Parallel(n_jobs=10)(delayed(fit_once)(n) for _ in tqdm.tqdm(range(50)))
    l,par = sorted(res,key=lambda x:x[0])[-1]

    return par

from scipy import stats
##LISTS OF FUNCTIONS 
names = ['lognormal','truncated power-law','powerlaw',"log-logistic","log gamma","log laplace","log-weibull"]
funcs = [stats.lognorm, 
         stats.gamma,
         stats.powerlaw, 
         stats.fisk,
         stats.loggamma,
         stats.loglaplace,
         stats.gumbel_r]
def significance(p):
    if p<0.001:
        return "***"
    elif p<0.01:
        return "**"
    elif p<0.05:
        return "*"
    else:
        return "(ns)"

def compute_loglikelihood_ratios(values,
                                 funcs,
                                 names,
                                 compare_to=0):
    parameters = [find_best_fit(func, v) for func in funcs]
    likelihoods = [func.logpdf(v, *parameters[n3]) for n3,func in enumerate(funcs) ]
    ratios = [loglikelihood_ratio(likelihoods[compare_to],likelihoods[n]) for n in range(len(funcs))]
    ratios = [str(int(R))+significance(p) for R,p in ratios]
    return dict(zip(names,ratios)),parameters
##sizes

ratios = dict()
parameters = dict()
for level in range(2,6):
    v = data_size[data_size['level']==level]['size'].values
    res = compute_loglikelihood_ratios(values,
                                                 funcs,
                                                 ["R lognormal/{}".format(i) for i in names])
    ratios[("D1","level={}".format(level))] = res[0]
    parameters[level] = res[1]
###PLOT
fig, ax = plt.subplots(1,4, figsize = (3.60236*2,2))
plt.subplots_adjust(bottom = 0.25, wspace=0.4,left=0.05,right=0.99)

y1 = ax[0].get_position().get_points()[-1][1]+0.02

for level in range(0,4):
    
    ##VALUES
    v = data_size[data_size['level']==level+2]['size'].values

    
    ##WRITE TITLE
    x1 =  ax[level].get_position().get_points()[0][0]
    x2 = ax[level].get_position().get_points()[1][0]
    ax[level].set_title("D1, level={}".format(level+2),y = y1,x = (x2+x1)/2, transform=fig.transFigure)
    
    
    
    #PLOT HISTOGRAM
    X = np.logspace(np.log10(min(v)*2),np.log10(max(v)),25)
    hist,edges = np.histogram(v,bins = X, density = True)
    x, y = zip(*[(i,j) for (i,j) in zip((edges[:-1]+edges[1:])/2, hist) if j>0])
    ax[level].loglog(x,y, color = 'k',label = 'data')
    
    
    ##FIND PARAMETERS
    aa = [func.pdf(X,*parameters[level][n]) for n, func in enumerate(funcs[:2])]
    for n, i in enumerate(aa):
        ax[level].loglog(X,i,label = names[n])
    ax[level].set_xlabel('container size (m)')
    ax[level].set_ylabel('pdf')
    ymin, ymax = ax[level].get_ylim()
    ax[level].set_ylim(min(y)*0.1,max(aa[1]*2))

        

        
##WRITE LETTERS


y1 = ax[0].get_position().get_points()[-1][1]+0.02

letters = ['a','b','c','d']
positions = [ax[i].get_position().get_points()[0][0] for i in range(4)] 
positions2 = [y1]*(4)
for n,a in enumerate(ax):
    a.text(positions[n], positions2[n], 
           letters[n], 
           transform=fig.transFigure, 
           fontsize = 8, 
           weight="bold", 
           fontdict={ 'family':'Helvetica 55 Roman'} )
    
largest_x = ax[-1].get_position().get_points()[1][0]



plt.savefig("../outputs/Extended_data/Extended_data_3_a_d.pdf",dpi=600)
plt.savefig("../outputs/Extended_data/Extended_data_3_a_d.svg",dpi=600)

#Print to file
table = pd.DataFrame(ratios).loc[["R lognormal/{}".format(i) for i in names[1:]]]


doc = docx.Document()
t = doc.add_table(table.shape[0]+1,table.shape[1])
for j in range(table.shape[-1]):
    t.cell(0,j).text = table.columns[j]

for i in range(table.shape[0]):
    for j in range(table.shape[-1]):
        t.cell(i+1,j).text = str(table.values[i,j])
for row in t.rows:
    for cell in row.cells:
        paragraphs = cell.paragraphs
        for paragraph in paragraphs:
            for run in paragraph.runs:
                font = run.font 
                font.name = 'Helvetica 55 Roman'
                font.size = Pt(7)
doc.save("../outputs/Extended_data/Extended_data_2_D1.docx")
##TIME
import warnings 
warnings.filterwarnings("ignore")
from collections import defaultdict
ratios_time = dict()
for level in range(2,6):
    v = data_time[data_time['trip_scale']==level]['inter_time'].values
    ratios_time[("D1","level={}".format(level))] = compute_loglikelihood_ratios(values,
                                                 funcs[:2],
                                                 ["R lognormal/{}".format(i) for i in names[:2]])


pd.DataFrame.from_dict(ratios_time,orient='columns')
#Print to file
table = pd.DataFrame(ratios_time).iloc[1:]


doc = docx.Document()
t = doc.add_table(table.shape[0]+1,table.shape[1])
for j in range(table.shape[-1]):
    t.cell(0,j).text = table.columns[j]

for i in range(table.shape[0]):
    for j in range(table.shape[-1]):
        t.cell(i+1,j).text = str(table.values[i,j])
for row in t.rows:
    for cell in row.cells:
        paragraphs = cell.paragraphs
        for paragraph in paragraphs:
            for run in paragraph.runs:
                font = run.font 
                font.name = 'Helvetica 55 Roman'
                font.size = Pt(7)
doc.save("../outputs/Extended_data/Extended_data_3_D1.docx")<jupyter_output><empty_output><jupyter_text>charachteristics of lognormal<jupyter_code>!pip install statsmodels
import scipy
from uncertainties import ufloat
from uncertainties.umath import exp,pow,log,sqrt


from statsmodels.base.model import GenericLikelihoodModel


class Lognorm(GenericLikelihoodModel):

    def loglike(self, params):
        shape=params[0]
        scale=params[-1]
        return scipy.stats.lognorm.logpdf(self.endog, shape,0,scale).sum()

def median(mu,sigma):
    return exp(mu)
    
def mode(mu,sigma):
    return exp(mu-pow(sigma,2))

def mean(mu,sigma):
    return exp(mu+pow(sigma,2)/2)

def variance(mu,sigma):
    a = (exp(2*mu+pow(sigma,2)))
    b = (exp(pow(sigma,2)) - 1)
    return a*b

def coeff(mu,sigma):
    return sqrt(exp(pow(sigma,2))-1)


values = []

for level in range(0,4):
    
    ##VALUES
    data = data_size[data_size['level']==level+2]['size'].values

    params = scipy.stats.lognorm.fit(data, floc=0)
    shape=params[0]
    scale=params[-1]
    model_res = Lognorm(data).fit(start_params=[shape,scale])
    
    
    expmu = ufloat(model_res.params[1], model_res.bse[1])                          
    mu = log(expmu)
    sigma = ufloat(model_res.params[0], model_res.bse[0])
    
    values.append(['{:%.3f}'.format(mu),
                  '{:%.3f}'.format(sigma),
                  '{:%.3f}'.format(median(mu,sigma)),
                  '{:%.3f}'.format(mode(mu,sigma)),
                  '{:%.3f}'.format(mean(mu,sigma)),
                  '{:%.3f}'.format(sqrt(variance(mu,sigma))),
                  '{:%.3f}'.format(coeff(mu,sigma))]
                 )
    print()
d = pd.DataFrame(values,columns = ['\mu','\sigma','median','mode','mean','std','coeff_variation'])
d.index = range(2,6)
d[['\mu','\sigma','median','coeff_variation','mode']]
#Print to file
table = pd.DataFrame(d)

doc = docx.Document()
t = doc.add_table(table.shape[0]+1,table.shape[1])
for j in range(table.shape[-1]):
    t.cell(0,j).text = table.columns[j]

for i in range(table.shape[0]):
    for j in range(table.shape[-1]):
        t.cell(i+1,j).text = str(table.values[i,j])
for row in t.rows:
    for cell in row.cells:
        paragraphs = cell.paragraphs
        for paragraph in paragraphs:
            for run in paragraph.runs:
                font = run.font 
                font.name = 'Helvetica 55 Roman'
                font.size = Pt(7)
doc.save("../outputs/Extended_data/Extended_data_5_D1.docx")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14342518956907713, 'fraction_numerical': 0.021361198446458295, 'mean_word_length': 2.6648593697051846, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '340564', 'license': 'MIT', 'n_tokens_mistral': 3826, 'n_tokens_neox': 3681, 'n_words': 798, 'path': 'Paper_figures/Figure1_and_Supplementary_Note_1.ipynb', 'repo_name': 'janbenetka/scales_human_mobility', 'size': '10814'}