5  Analyze

With datasets in hand, the next step is to turn raw counts into distributions and summary statistics. The analyze function adds two items to a dataset, modifying it in place: a 'distribution' item, holding the relative frequencies at each checkpoint, and a 'statistics' item, holding the mean, second moment, variance, standard deviation, median, and mode. It, and the two ancillary functions it is built on, live in src/primes_in_intervals/statistics.py; from the shell, pii analyze (see docs/cli.md).

Every chapter of this book runs in its own session, so the folded cell below rebuilds the running datasets of Chapter 4 from the database.

import numpy as np

# NumPy 2 displays scalars inside containers as np.float64(...); restore the
# plain numeric display so statistics read cleanly.
np.set_printoptions(legacy="1.25")

from primes_in_intervals import nest, partition, retrieve


def pick(found, lower_bound):
    """Return the retrieved dataset whose data begins at the given lower bound.

    Introduced in the Datasets chapter: accepts either return shape of
    ``retrieve`` (one dataset, or a list of them).
    """
    if isinstance(found, dict):
        found = [found]
    return next(d for d in found if d['header']['lower_bound'] == lower_bound)


# Example 1: interval length 100 over (2 * 10**6, 3 * 10**6].
A1, B1, H1 = 2*10**6, 3*10**6, 100
retrieve_data1disj = retrieve(H1, 'disjoint')
retrieve_data1olap = pick(retrieve(H1, 'overlap'), A1)
partition(retrieve_data1disj)
partition(retrieve_data1olap)
nest_data1disj = nest(retrieve(H1, 'disjoint'))
nest_data1olap = nest(pick(retrieve(H1, 'overlap'), A1))

# Example 2: interval length 85 about exp(18), and the companion length-90
# datasets used for nesting.
N2, H2 = int(np.exp(18)), 85
retrieve_data2disj = retrieve(H2, 'disjoint')
retrieve_data2olap = retrieve(H2, 'overlap')
partition(retrieve_data2disj)
partition(retrieve_data2olap)
nest_data2disj = nest(retrieve(90, 'disjoint'))
nest_data2olap = nest(pick(retrieve(90, 'overlap'), 65559979))

from primes_in_intervals import analyze

The first ancillary function, dictionary_sort, returns a copy of a dictionary sorted by its keys (the checkpoints’ frequency dictionaries carry their keys in order of first appearance, not in numerical order).

The second, dictionary_statistics, takes a frequency dictionary, whose keys are numbers and each key’s value is the number of times (frequency of) the key occurs in some data, and computes everything the analysis needs without ever expanding the data into a list. Its core is a single pass over the sorted frequencies:

    for s in frequencies.keys():
        number_of_objects_counted += frequencies[s]
        mean += s*frequencies[s]
        second_moment += (s**2)*frequencies[s]
        if frequencies[s] == M:
            mode.append(s)
    mean = mean/number_of_objects_counted
    second_moment = second_moment/number_of_objects_counted
    variance = second_moment - mean**2
    standard_deviation = np.sqrt(variance)

Here M is the maximum frequency, so the mode collects every key attaining it, and is therefore a list; and the variance is the second moment minus the square of the mean, over the whole population of intervals. The median comes from cumulative frequencies, in a small subroutine: with \(n\) values counted in all, walk the keys in increasing order, accumulating their frequencies; if \(n\) is odd, the median is the first key at which the cumulative count exceeds \(n/2\), and if \(n\) is even, it is the average of the first keys at which the cumulative count reaches \(n/2\) and \(n/2 + 1\). Finally, dividing each frequency by \(n\) gives the relative frequencies, and the function returns a dictionary with the distribution under 'dist' and the summary statistics under 'mean', '2ndmom', 'var', 'sdv', 'med', and 'mode'.

analyze applies dictionary_statistics across a whole dataset. It handles both configurations. For a flat dataset, it runs over the checkpoints of the 'data' item, skipping the first: the first checkpoint’s counts are all zero, so its distribution and statistics are stored as empty dictionaries (no meaningful statistics for the trivial item). For a nested dataset, it runs over every pair in 'nested_interval_data', none of which is trivial. Either way, the results land in the new 'distribution' and 'statistics' items, keyed exactly as the source data is, the header’s 'contents' grows accordingly, and the modified dataset is returned. Guards cover the remaining cases: a dataset that has already been analyzed is left alone with a message, as is one with no data at all.

Example 1

analyze returns the whole dataset, which by now runs to many hundreds of lines, so here we display the header, in which 'contents' records the two new items.

analyze(retrieve_data1disj)['header']

{'interval_type': 'disjoint',
 'lower_bound': 2000000,
 'upper_bound': 3000000,
 'interval_length': 100,
 'no_of_checkpoints': 101,
 'contents': ['data', 'partition', 'distribution', 'statistics']}

retrieve_data1disj['distribution'][2500000]

{0: 0.0,
 1: 0.0024,
 2: 0.0094,
 3: 0.0316,
 4: 0.0752,
 5: 0.1366,
 6: 0.186,
 7: 0.1986,
 8: 0.1602,
 9: 0.1104,
 10: 0.0538,
 11: 0.0246,
 12: 0.0082,
 13: 0.0024,
 14: 0.0006,
 15: 0.0,
 17: 0.0}

retrieve_data1disj['statistics'][2500000]

{'mean': 6.8278,
 '2ndmom': 50.6258,
 'var': 4.006947160000003,
 'sdv': 2.00173603654428,
 'med': 7.0,
 'mode': [7]}

analyze(retrieve_data1olap)['header']

{'interval_type': 'overlap',
 'lower_bound': 2000000,
 'upper_bound': 3000000,
 'interval_length': 100,
 'no_of_checkpoints': 101,
 'contents': ['data', 'partition', 'distribution', 'statistics']}

retrieve_data1olap['distribution'][2500000], retrieve_data1olap['statistics'][2500000]

({0: 0.000156,
  1: 0.001936,
  2: 0.009184,
  3: 0.031292,
  4: 0.076508,
  5: 0.135348,
  6: 0.186364,
  7: 0.19746,
  8: 0.162852,
  9: 0.108816,
  10: 0.05598,
  11: 0.023608,
  12: 0.008028,
  13: 0.002004,
  14: 0.000428,
  15: 3.6e-05,
  16: 0.0,
  17: 0.0,
  18: 0.0},
 {'mean': 6.827924,
  '2ndmom': 50.59066,
  'var': 3.9701138502239957,
  'sdv': 1.9925144542070443,
  'med': 7.0,
  'mode': [7]})

analyze(nest_data1disj)
nest_data1disj['distribution'][(2490000, 2510000)]

{0: 0.0,
 1: 0.0,
 2: 0.02,
 3: 0.03,
 4: 0.055,
 5: 0.135,
 6: 0.22,
 7: 0.21,
 8: 0.165,
 9: 0.105,
 10: 0.025,
 11: 0.01,
 12: 0.015,
 13: 0.01,
 14: 0.0,
 15: 0.0,
 17: 0.0}

analyze(nest_data1olap)
nest_data1olap['statistics'][(2490000, 2510000)]

{'mean': 6.7601,
 '2ndmom': 49.3741,
 'var': 3.675147989999992,
 'sdv': 1.917067549670588,
 'med': 7.0,
 'mode': [6]}

Example 2

analyze(retrieve_data2disj)['header']

{'interval_type': 'disjoint',
 'lower_bound': 65557969,
 'upper_bound': 65761969,
 'interval_length': 85,
 'no_of_checkpoints': 201,
 'contents': ['data', 'partition', 'distribution', 'statistics']}

analyze(retrieve_data2olap)['header']

{'interval_type': 'overlap',
 'lower_bound': 65557969,
 'upper_bound': 65761969,
 'interval_length': 85,
 'no_of_checkpoints': 201,
 'contents': ['data', 'partition', 'distribution', 'statistics']}

analyze(nest_data2disj)
Cnest = list(nest_data2disj['nested_interval_data'].keys())
nest_data2disj['statistics'][Cnest[0]]

{'mean': 5.045454545454546,
 '2ndmom': 28.40909090909091,
 'var': 2.9524793388429735,
 'sdv': 1.718278015585072,
 'med': 5.0,
 'mode': [7]}

analyze(nest_data2olap)
Cnest = list(nest_data2olap['nested_interval_data'].keys())
Cnest[20], nest_data2olap['distribution'][Cnest[20]]

((65639179, 65680759),
 {0: 0.0012987012987012987,
  1: 0.016378066378066377,
  2: 0.06955266955266955,
  3: 0.1303030303030303,
  4: 0.197017797017797,
  5: 0.2342953342953343,
  6: 0.16998556998557,
  7: 0.10134680134680135,
  8: 0.047546897546897546,
  9: 0.023136123136123135,
  10: 0.0075998075998075995,
  11: 0.001443001443001443,
  12: 9.62000962000962e-05,
  13: 0.0})