概率分布导致“进程以退出代码 137 完成（被信号 9：SIGKILL 中断）”答案

【问题标题】：Probability distribution leads to the 'Process finished with exit code 137 (interrupted by signal 9: SIGKILL)'概率分布导致“进程以退出代码 137 完成（被信号 9：SIGKILL 中断）”
【发布时间】：2020-02-17 09:21:13
【问题描述】：

我正在尝试为我的地质研究创建某种简化的 Oracle Crystal Ball 应用程序，它将使用 P90（90% 置信度）和 P10（10% 置信度）值作为不同概率场景的输入和返回分布。听起来像蒙特卡洛分布。我是 Python 新手，最近才开始，顺便说一句 :)

本主题将分为四个关键部分：

作品范围的一般描述。
伪编码（不过以前从未尝试过）。
实际 Python 代码。
我在这里的原因或逻辑/代码问题。

PART 1. 作品范围概述。

为简单起见，假设我们只有三个类别，每个类别都有 P90 和 P10 参数，它们之间没有任何步骤：
- cat_1：[1, 2]
- cat_2: [2, 4]
- cat_3：[3, 6]
利用笛卡尔积，我们得到以下 8 个可能场景的列表：
- [1, 2, 3], [1, 2, 6], [1, 4, 3], [1, 4, 6], [2, 2, 3], [2, 2, 6] , [2, 4, 3], [2, 4, 6]
将每个列表中的参数相乘得到以下产品：
- [6、12、12、24、12、24、24、48]
测量每个产品的频率会导致：
- {6: 1, 12: 3, 24: 3, 48: 1}，或考虑百分比：
- {6: 12.5%, 12: 37.5%, 24: 37.5%, 48: 12:5%,}表示出现12或24的概率高于6或48。
这就是我想要得到的结果：知道产品能够获得均值、中值和众数值的概率。
对于我的硬件来说，困难的部分是真实案例中可能出现的大量场景。共有六个类别，在 P90 和 P10 值之间有小步长。考虑公制，P90和P10的取值范围可能如下：
- 正方形面积：0.01 - 100.00 km2，步长0.01；
- 层厚：0.10 - 100.00 m，步长0.1；
- 孔隙率：0.01 - 1.00 p.u.，步长 0.01；
- 饱和度：0.01 - 1.00 p.u.，步长 0.01；
- 压力：1 - 2000 atm，步骤 1 atm；
- 表面：0.01 - 1.00 p.u.，步长 0.01。
通常，实际案例研究将使用更窄的范围，例如 0.1 - 2.0 km2 的面积，1 - 10 m 的厚度，8 - 15 的孔隙度等。尽管如此，即使在这种情况下，它听起来也像考虑到上述步骤，“谷歌”可能的场景数量。结果，我收到以下通知，这是关键问题：

进程以退出代码 137 结束（被信号 9：SIGKILL 中断）。

当总计算量超过 ~10MM 和 ~1 分钟时会发生这种情况（经过实验检查，因此数字是粗略的）。

第 2 部分。伪编码。

良好的做法是在伪编码时不应该进行抽象，但是我在这个领域的经验为零，因此会尽力而为。

User inputs minimum possible values (P90) for total 6 categories
User inputs maximum possible values (P10) for total 6 categories

Total 6 list are created (square area, layer thickness, porosity etc.), 1 per each category that contain a range of possible values and indicated step (P90_category1, P10_category1, step1)

Use a Cartesian product to create a list_of_tuples with possible scenarios

Convert list_of_tuples to the list_of_lists

Create empty_list
for each element in the list_of_lists:
    calculate its product
    append to the empty_list

Round values in the empty_list

Create a dictionary that counts similar values in the empty_list

Calculate a probability of each value according to its repetition frequency in the dictionary

就是这样。还应用了一些基本的统计数据和绘图，但这不是关键时刻。

第 3 部分。实际 Python 代码。

最初的 P90 值（90% 置信度）：

P90_area = float(input('P90 area: '))
P90_thickness = float(input('P90 thickness: '))
P90_porosity = float(input('P90 porosity: '))
P90_saturation = float(input('P90 saturation: '))
P90_pressure = float(input('P90 pressure: '))
P90_surface = float(input('P90 surface: '))

然后是 P10 值（10% 置信度）：

P10_area = float(input('P10 area: '))
P10_thickness = float(input('P10 thickness: '))
P10_porosity = float(input('P10 porosity: '))
P10_saturation = float(input('P10 saturation: '))
P10_pressure = float(input('P10 pressure: '))
P10_surface = float(input('P10 surface: '))

使用特定步骤创建从 P90 到 P10 的值范围

area_values = np.arange(P90_area, P10_area + 0.01, 0.01)
thickness_values = np.arange(P90_thickness, P10_thickness + 0.1, 0.1)
porosity_values = np.arange(P90_porosity, P10_porosity + 0.01, 0.01)
saturation_range = np.arange(P90_saturation, P10_saturation + 0.01, 0.01)
pressure_range = np.arange(P90_pressure, P10_pressure + 1, 1)
surface_range = np.arange(P90_surface, P10_surface + 0.01, 0.01)

将所有列表组合成笛卡尔积（即[（area1，thickness1，porosity1），（area1，thickness1，porosity2）等]）：

list_of_tuples = list(itertools.product(area_values, thickness_values, porosity_values, saturation_range, pressure_range, surface_range)

将元组列表转换为列表列表：

list_of_lists = [list(elem) for elem in list_of_tuples]

创建一个包含相乘值的列表并对它们进行排序（'np.prod' 为每个列表返回一个产品）：

multiplied_values = []
for i in list_of_lists:
    i = np.prod(np.array(i))
    multiplied_values.append(i)
multiplied_values = sorted(multiplied_values)

四舍五入：

rounded_values = [float(Decimal('%.2f' % elem)) for elem in multiplied_values]

创建一个统计所有相似/独特对象的字典：

counts = Counter(rounded_values)

通过将值除以列表中元素的总数来计算概率：

probability_mass = {k: v/total for k, v in counts.items()}

它有效，这里有简单的统计数据和特定案例的图表：

总计算量：4899510
P90 为：5.60
P10 是：43.41
P50（概率最大的值）为：15.24
平均值为：23.80

Figure. Probability distribution diagram

第一个问题很关键，因为它阻塞了大堆数据的计算：

第 4 部分。关键问题。

第一季度。关键问题：

因此，我收到以下通知，这是关键问题：

进程以退出代码 137 结束（被信号 9：SIGKILL 中断）。

根据类似的主题，我的脚本很可能由于 CPU 使用率过高而被操作系统杀死。我在运行代码时使用 'top' 命令检查了 CPU 负载，当它可以处理输入参数时，CPU 负载高达 100%，在某些时候被中断时高达 110%。

问题：无论如何，我怎样才能摆脱这种中断，让脚本运行到需要的时间？我很乐意等待为大型数据堆栈获得正确的分布。为每个参数增加一个步骤是一个硬核选项，我不想这样做。

第二季度。概率分布图看起来不像经典的正态分布，而最大概率值和平均值之间的差异很大。你怎么看，代码的逻辑会不会有问题？

附：我知道这个脚本看起来很颠簸，希望你的眼睛不会流血）

【问题讨论】：

听起来您正在实施一种需要计算所有可能性的方法。没关系，这是一个很好的开始方式。在这一点上，我的建议是首先通过增加每个变量的步长来减少需要计算的值的数量。目前该步长似乎是 0.01，也许可以尝试 0.1。这个想法是让它使用更少数量的值，然后尝试减小步长（增加步数）。此外，一旦你对问题有了感觉，就寻找一种更具分析性或象征性的方法。祝你好运，玩得开心。
@RobertDodier 谢谢你，罗伯特）正确，关键是计算每一个可能的场景。当我们确实有一系列可能的地质条件值（提到的面积、厚度等）时，这在地质学中被广泛使用，但是不知道精确值，因为目前还没有直接测量。所以我做了你写的，从小开始，不断增加参数和减少步骤。这就是我在实际情况中将面临的问题，这就是为什么操作系统的中断非常令人沮丧，因为脚本可以工作，但范围相对较小）
我会查看完成计算所需的时间如何随步数的增加而变化，并尝试估计使用我最初想要的步数完成计算所需的时间。如果结果太多，有一些策略可以尝试保持在时间和空间的限制范围内。例如。用大步搜索到附近，然后小步细化。或应用启发式搜索方法，如模拟退火，或无梯度搜索，如多面体（变形虫）算法，或在存在梯度时使用梯度的方法（LBFGS 等）。
@RobertDodier 这么多术语我不熟悉，但这看起来是一个进一步了解它的好机会）谢谢你的提示和指导！如果我有任何新的正面或负面结果，我会在这里发表评论。

标签： python statistics cpu probability montecarlo

【解决方案1】：

由于您尝试计算每种可能的情况，因此此处所需的计算量会随着您的每个范围内的元素数量呈指数增长。我很想尝试为您调试一个完整的代码，但我需要输入，所以您能否发布完整的代码以及已经指定的输入，以便我们知道可以使用哪些合理的值。

稍微有点不同，与其尝试修复您的代码，我们可以先尝试解决您的原始问题吗？当您说“简化的概率分布计算器”时，您是什么意思？在我们尝试了解如何在 Python 中实现该过程之前，您能否用伪代码编写步骤让我们了解该过程。

根据您对上述问题的回答，我可能会建议您采用抽样方法，而不是评估每一种可能性。查找蒙特卡罗模拟。如果您有一个使用新数据更新的先前分布，并且您想知道后验（最终）分布，那么请考虑使用贝叶斯方法，特别是 Winbugs（不是 Python 但非常适合贝叶斯东西的独立程序）。

PS。我知道我的回答可能更适合写成评论，但显然你需要 +50 声望，而我还没有:(

【讨论】：

感谢您的全面回答，马修！我将使用建议的规格编辑我的帖子，并在此处发表评论，以便您收到通知。不过，这几天会这样做，很可能是在周末。
Matthew，刚刚编辑了最初的帖子并添加了伪代码。实际上，我确实想做某种蒙特卡罗模拟，你说得对。
根据您的伪编码，我有一些后续问题：1）由于您给出了 p10 和 p90，因此您暗示您的参数遵循概率分布。我需要知道是哪一个。如果您说它们是正态分布的，那么我们会遇到负（不可能）值的问题，因为 p10 非常接近于零，以至于分布的下尾包含大量负数。他们更有可能遵循 Weibull 或对数正态分布，或截断正态分布。或者，如果我们绘制一个负数，那么我们可以重新采样。你要哪个？
2) 我仍然不确定您想要实现什么作为查找输出。你能像考试题一样说它吗？例如。 “给定 X 和 Y，Z 是什么”
我必须在消息中指出这一点，但只能使用正参数。我还没有编写任何规则，只是为了在这一点上更容易。

【解决方案2】：

所以我在输入参数的均匀分布、随机抽样和笛卡尔积方面完成了您需要的工作。结果有点像指数分布。它最好由 weibull 分布建模。

我做了一些进一步的分析，因为任何模拟的结果都应该进一步调查，以检查模拟是否足够。为此，我做了一个 10,100,1000,10000,100000,10000000 个样本的蒙特卡罗样本来生成直方图。我们从拟合的 weibull 的 alpha 和 beta 收敛中看出，100 万个样本就足够了。

我相信您对此会有疑问，因此请在下方提问。请注意，直方图是对数刻度，因此在可视化分布时需要牢记这一点（或注释掉 xscale 和 yscale 线）。

结果如下： https://i.stack.imgur.com/viQ9i.png https://i.stack.imgur.com/0kc4n.png

这是生成输出的代码：

import numpy as np
from tqdm import tqdm
import random
import matplotlib.pyplot as plt
import scipy.stats as ss

#these should be user inputs
area_min = 0.01
area_max = 100
thickness_min = 0.1
thickness_max = 100
porosity_min = 0.01
porosity_max = 1
saturation_min = 0.01
saturation_max = 1
pressure_min = 1
pressure_max = 2000
surface_min = 0.01
surface_max = 1

grid_resolution = 1000 #how finely we will slice each property. I have kept this consistent as it makes more sense to do so when sampling
#With a grid_resolution of 1000, the number of possible combinations here is 1000^6 ==> 10^18 so we will randomly sample the array
#I assume you want to get a probability distribution of these combinations.
area_array = np.linspace(area_min,area_max,grid_resolution)
thickness_array = np.linspace(thickness_min,thickness_max,grid_resolution)
porosity_array = np.linspace(porosity_min,porosity_max,grid_resolution)
saturation_array = np.linspace(saturation_min,saturation_max,grid_resolution)
pressure_array = np.linspace(pressure_min,pressure_max,grid_resolution)
surface_array = np.linspace(surface_min,surface_max,grid_resolution)

#it is important to try different sample sizes to be sure your sample is large enough
samples_to_test = [1,2,3,4,5,6] #log10 scale

xmax = 10**8
alpha_array = []
beta_array = []
plt.figure(figsize=(12,10))
for i,s in enumerate(samples_to_test):
    plt.subplot(231+i)
    samples = 10**s
    product_array = []
    for _ in tqdm(range(samples)):
        area = random.choice(area_array)
        thickness = random.choice(thickness_array)
        porosity = random.choice(porosity_array)
        saturation = random.choice(saturation_array)
        pressure = random.choice(pressure_array)
        surface = random.choice(surface_array)
        product_array.append(area*thickness*porosity*saturation*pressure*surface)

    xvals = np.logspace(1,np.log10(xmax),1000)
    [beta,_,alpha] = ss.weibull_min.fit(data=product_array,floc=0)
    alpha_array.append(alpha)
    beta_array.append(beta)
    weibull_yvals = ss.weibull_min.pdf(xvals,beta,scale=alpha)
    plt.plot(xvals,weibull_yvals)
    print('Weibull fit parameters:\nalpha =',alpha,'\nbeta =',beta)
    [mean,variance] = ss.weibull_min.stats(beta, loc=0, scale=alpha, moments='mv')
    median = ss.weibull_min.median(beta, loc=0, scale=alpha)
    print('Mean =',mean)
    print('Median =',median)
    print('Standard deviation =',variance**0.5)

    plt.hist(product_array,bins=1000,density=True)
    plt.yscale('log')
    plt.xscale('log')
    plt.xlabel('Cartesian Product of parameters')
    plt.ylabel('Probability density ($log_{10}$ scale)')
    plt.title(str('Monte Carlo samples = '+str(samples)))
    plt.xlim(10,xmax)
    plt.ylim(10**-8,0.0001)

plt.suptitle('Probability of of a given cartesian product of the specified parameters\nmeasured using different numbers of Monte Carlo samples')
plt.figure(figsize=(12,5))
plt.subplot(121)
plt.semilogx(10**np.array(samples_to_test),alpha_array,label='alpha')
plt.legend()
plt.subplot(122)
plt.semilogx(10**np.array(samples_to_test),beta_array,label='beta')
plt.legend()
plt.suptitle('Test results for alpha and beta')
plt.show()

Output:
100%|██████████| 10/10 [00:00<?, ?it/s]
Weibull fit parameters:
alpha = 86642.0194345818 
beta = 0.4938259951069627
Mean = 177350.7081149186
Median = 41247.66458603765
Standard deviation = 403557.41514732403
100%|██████████| 100/100 [00:00<00:00, 100246.27it/s]
Weibull fit parameters:
alpha = 177861.91287733015 
beta = 0.6310314479279571
Mean = 251385.7124440623
Median = 99503.40459313976
Standard deviation = 415414.97618995525
100%|██████████| 1000/1000 [00:00<00:00, 199131.37it/s]
Weibull fit parameters:
alpha = 171932.22877129668 
beta = 0.5452693527437176
Mean = 296661.14084923535
Median = 87788.61401806296
Standard deviation = 589615.4680695855
100%|██████████| 10000/10000 [00:00<00:00, 179051.70it/s]
Weibull fit parameters:
alpha = 166909.86147776648 
beta = 0.5172460791589029
Mean = 314175.4976503747
Median = 82176.44526800542
Standard deviation = 670314.3944630618
100%|██████████| 100000/100000 [00:00<00:00, 144477.93it/s]
Weibull fit parameters:
alpha = 167711.26073670806 
beta = 0.5194333533253157
Mean = 313393.61873437575
Median = 82817.74728224205
Standard deviation = 664803.5086740599
100%|██████████| 1000000/1000000 [00:07<00:00, 140706.15it/s]
Weibull fit parameters:
alpha = 168089.6178189406 
beta = 0.5186379527889259
Mean = 314930.2501968761
Median = 82914.8108556469
Standard deviation = 669461.6904337168

【讨论】：

感谢您的努力，非常感谢！第一个不明显的东西是'231 + 1' - 那是什么？第二个困难的部分从“for _ in tqdm”开始，直到循环结束。不知道那是什么。第三个与 alpha/beta/weibull 和 matplotlib 有关。我试图通过阅读论文来理解weibull，但没有任何帮助，仍然没有任何线索））最大的问题是我不明白这些计算的一般逻辑。
231+i 用于子图。因此，在要测试的样本循环中，我将是 0、1、2、3、4、5，因此 tubplot 将从 231 开始，然后是 232，然后是 233。查找子图以获取更多信息。 tqdm 是一个进度指示器。这里不需要，但很高兴知道你的 for 循环需要多长时间。使用“for _ in”代替“for x in”，因为我没有使用变量 x，所以我只是使用 _ 告诉 python 不要存储计数器。至于weibull分布，它是一种可以呈现多种形状的概率分布。阅读维基：en.wikipedia.org/wiki/Weibull_distribution
我遵循的步骤是：1）创建输入参数的线性空间数组 2）从每个参数中随机选择一个值 3）将所有这些随机选择的值相乘 4）制作直方图结果 5) 将概率分布拟合到直方图（在本例中为 Weibull 分布） 6) 打印 weibull 分布的平均值、中位数、标准差我还通过查看结果是否检查了 100 万个样本是否足够weibull 分布在 100K 和 100 万之间非常不同，但事实并非如此。所以不需要更多。
正如我之前提到的，我认为将随机样本相乘是无稽之谈，因为您的单位为 km^2.m.pu^3.atm。如果我是你，我希望通过将每个参数的边际概率相乘而不是通过参数本身相乘来获得联合概率分布。为此，您最好从概率分布（例如正态分布）中采样，而不是像我们对线性间隔数组所做的那样从均匀分布中采样。我知道这是很多 python 和统计信息的结合，但你的问题需要它。还有其他问题吗？