Python:为变量创建高斯分布并使用高斯值循环运行程序

Question

我仍在尝试使用 python，并且正在努力解决这个问题:

我正在使用名为 Fatiando a Terra 的工具包在 python 中运行脚本；我在这里定义了一个物体的参数，然后用重力/梯度法对其进行建模，如下所示:

import numpy as np
from fatiando.vis import mpl
from fatiando.mesher import Prism
from fatiando.gravmag import prism
from fatiando.constants import G
from fatiando import utils
import fatiando
import matplotlib.pyplot as plt

model = Prism(-1000, 1000, -1000, 1000, 1000, 1100, {'density': 300})
n = 500
x = np.zeros(n)
y = np.zeros(n)
z = np.linspace(0, 2000, n)
data = np.array([prism.gx(x, y, z, [model]),
                 prism.gy(x, y, z, [model]),
                 prism.gz(x, y, z, [model]),
                 prism.gxx(x, y, z, [model]),
                 prism.gxy(x, y, z, [model]),
                 prism.gxz(x, y, z, [model]),
                 prism.gyy(x, y, z, [model]),
                 prism.gyz(x, y, z, [model]),
                 prism.gzz(x, y, z, [model])])

在此之后，我正在使用 python 绘制它，但是这段代码不相关。

我手动执行了一个微型高斯阵列，使用:

from numpy.random import multivariate_normal
multivariate_normal([300]. [[300]], 10)

它为我提供的 10 个值是我手动输入到脚本中并生成输出的。

我想做的是对密度值为 300 执行高斯分布 +-100。我想生成 1000 个值并将它们循环回脚本，运行程序 1000 次，一次用于每个密度变体。

对于输出，我目前使用:

titles = ['gx', 'gy', 'gz', 'gxx', 'gxy', 'gxz', 'gyy', 'gyz', 'gzz']
np.savetxt(title, np.vstack((z.ravel(), d.ravel())).T)

理想情况下，对于高斯数组，我想以相同的方式输出文本，但使用 1000 个样本的平均值而不是任何单个值。标准偏差也会很大。

---

如果我问得太多，我深表歉意，非常感谢您提供的任何帮助。干杯!

Answer 1

这更像是一个扩展的评论而不是一个答案，尽管它也可以作为后者。

将您的代码包装在一个函数中，然后在循环中调用该函数，最后对结果进行平均，可能看起来像这样(部分符号):

##import statements

NDENSITIES = 1000

def dowork(density):
    model = Prism(-1000, 1000, -1000, 1000, 1000, 1100, {'density': 300})
    ##other code
    return data   ## or whatever is important

densities = multivariate_normal([300], [[300]], NDENSITIES)
results = []
for density in range(densities):
    results.append(dowork(density))
results = np.vstack(results)   ## or hstack or dstack, depending on the dimensions of `data`
mean = results.mean()   # possibly: results.mean(axis=0) or with axis=1 etc. Again dependent on the dimensions
std = results.std()