根据两个预先计算的直方图报告两个样本的 K-S 统计量

Problem:

在这里，我绘制了存储在文本文件中的 2 个数据集（在列表中）dataset）每个包含 218 亿个数据点。这使得数据太大而无法作为数组保存在内存中。我仍然能够将它们绘制为直方图，但我不确定如何通过2 样本KS测试 http://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.stats.ks_2samp.html。这是因为我无法弄清楚如何访问 plt 对象中的每个直方图。

Example:

这是一些生成虚拟数据的代码：

mu = [100, 120]
sigma = 30
dataset = ['gsl_test_1.txt', 'gsl_test_2.txt']
for idx, file in enumerate(dataset):
    dist = np.random.normal(mu[idx], sigma, 10000)
    with open(file, 'w') as g:
        for s in dist:
            g.write('{}\t{}\t{}\n'.format('stuff', 'stuff', str(s)))

这会生成我的两个直方图（使here https://stackoverflow.com/questions/37082402/how-to-build-or-precompute-a-histogram-from-a-file-too-large-for-memory):

chunksize = 1000
dataset = ['gsl_test_1.txt', 'gsl_test_2.txt']
for fh in dataset:
    # find the min, max, line qty, for bins
    low = np.inf
    high = -np.inf

    loop = 0
    for chunk in pd.read_table(fh, header=None, chunksize=chunksize, delimiter='\t'):
        low = np.minimum(chunk.iloc[:, 2].min(), low)
        high = np.maximum(chunk.iloc[:, 2].max(), high)
        loop += 1
    lines = loop*chunksize

    nbins = math.ceil(math.sqrt(lines))   

    bin_edges = np.linspace(low, high, nbins + 1)
    total = np.zeros(nbins, np.int64)  # np.ndarray filled with np.uint32 zeros, CHANGED TO int64

    for chunk in pd.read_table(fh, header=None, chunksize=chunksize, delimiter='\t'):

        # compute bin counts over the 3rd column
        subtotal, e = np.histogram(chunk.iloc[:, 2], bins=bin_edges)  # np.ndarray filled with np.int64

        # accumulate bin counts over chunks
        total += subtotal


    plt.hist(bin_edges[:-1], bins=bin_edges, weights=total)
    plt.savefig('gsl_test_hist.svg')

问题：

Most KS 统计的示例 https://stackoverflow.com/questions/10884668/two-sample-kolmogorov-smirnov-test-in-python-scipy使用两个原始数据/观察/点/等数组，但我没有足够的内存来使用这种方法。根据上面的示例，我如何访问这些预先计算的垃圾箱（从'gsl_test_1.txt' and 'gsl_test_2.txt'计算两个分布之间的 KS 统计量？

奖励业力：在图表上记录 KS 统计量和 p 值！

我稍微清理了一下你的代码。写信给StringIO所以它比写入文件更精简。设置默认氛围 w/seaborn代替matplotlib使其看起来更现代。这bins如果您希望统计测试保持一致，则两个样本的阈值应该相同。我认为如果你迭代并以这种方式制作垃圾箱，整个事情可能会比需要的时间更长。Counter可能很有用，因为您只需循环一次...而且您将能够制作相同的垃圾箱大小。将浮点数转换为整数，因为您将它们合并在一起。from collections import Counter then C = Counter() and C[value] += 1。你会有一个dict最后你可以从那里制作垃圾箱list(C.keys())。这会很好，因为你的数据是如此粗糙。另外，你应该看看是否有办法做chunksize with numpy代替pandas b/c numpy索引速度更快。尝试一个%timeit for DF.iloc[i,j] and ARRAY[i,j]你就会明白我的意思了。我将其中大部分内容编写为一个函数，以尝试使其更加模块化。

import numpy as np
import pandas as pd
import math
import matplotlib.pyplot as plt
from io import StringIO
from scipy.stats import ks_2samp
import seaborn as sns; sns.set()

%matplotlib inline

#Added seaborn b/c it looks mo betta

mu = [100, 120]
sigma = 30

def write_random(file,mu,sigma=30):
    dist = np.random.normal(mu, sigma, 10000)
    for i,s in enumerate(dist):
        file.write('{}\t{}\t{}\n'.format("label_A-%d" % i, "label_B-%d" % i, str(s)))
    return(file)

#Writing to StringIO instead of an actual file
gs1_test_1 = write_random(StringIO(),mu=100)
gs1_test_2 = write_random(StringIO(),mu=120)

chunksize = 1000

def make_hist(fh,ax):
    # find the min, max, line qty, for bins
    low = np.inf
    high = -np.inf

    loop = 0

    fh.seek(0)
    for chunk in pd.read_table(fh, header=None, chunksize=chunksize, sep='\t'):
        low = np.minimum(chunk.iloc[:, 2].min(), low) #btw, iloc is way slower than numpy array indexing
        high = np.maximum(chunk.iloc[:, 2].max(), high) #you might wanna import and do the chunks with numpy
        loop += 1
    lines = loop*chunksize

    nbins = math.ceil(math.sqrt(lines))   

    bin_edges = np.linspace(low, high, nbins + 1)
    total = np.zeros(nbins, np.int64)  # np.ndarray filled with np.uint32 zeros, CHANGED TO int64

    fh.seek(0)
    for chunk in pd.read_table(fh, header=None, chunksize=chunksize, delimiter='\t'):

        # compute bin counts over the 3rd column
        subtotal, e = np.histogram(chunk.iloc[:, 2], bins=bin_edges)  # np.ndarray filled with np.int64

        # accumulate bin counts over chunks
        total += subtotal

    plt.hist(bin_edges[:-1], bins=bin_edges, weights=total,axes=ax,alpha=0.5)

    return(ax,bin_edges,total)

#Make the plot canvas to write on to give it to the function
fig,ax = plt.subplots()

test_1_data = make_hist(gs1_test_1,ax)
test_2_data = make_hist(gs1_test_2,ax)

#test_1_data[1] == test_2_data[1] The bins should be the same if you're going try and compare them...
ax.set_title("ks: %f, p_in_the_v: %f" % ks_2samp(test_1_data[2], test_2_data[2]))