在python中动态构建随机变量的函数

Question

我有一些使用scipy.stats的随机变量如下：

import scipy.stats as st
x1 = st.uniform()
x2 = st.uniform()

现在我想根据以前的随机变量制作另一个随机变量，并对新的随机变量进行一些计算，例如var。假设我希望新的随机变量类似于max(2, x1) + x2。如何动态定义它？

Answer 1

我的旧答案如下：

（当编辑引用 SO 文档的答案以删除这些引用时，系统提示我再次查看此问题。我认为无论如何这是一个更好的答案。）

首先，据我所知，对于两个或多个变量的非线性函数的方差，没有一般的方法可以得到一个很好的封闭形式的表达式。可能大多数人都会采用某种蒙特卡洛策略来近似这样的数量。

这里有一些代码可以生成针对这种特定情况执行此操作的绘图。它适用于许多其他人。

从单位均匀随机变量中生成两个伪随机样本，然后根据这些样本的元素计算伪随机变量Y。

>>> import scipy.stats as stats
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> X1 = stats.uniform.rvs(0,1, 5000)
>>> X2 = stats.uniform.rvs(0,1, 5000)
>>> Y = [max(2,x1)+x2 for (x1,x2) in zip(X1,X2)]

现在，希望确定这个函数的密度函数，绘制它的直方图。

>>> plt.hist(Y)
(array([ 501.,  526.,  490.,  481.,  513.,  488.,  525.,  490.,  521.,  465.]), array([ 2.00012599,  2.10007992,  2.20003386,  2.2999878 ,  2.39994173,
        2.49989567,  2.59984961,  2.69980354,  2.79975748,  2.89971141,
        2.99966535]), <a list of 10 Patch objects>)
>>> plt.show()

我们很幸运，因为它很容易识别。在这里。

它是一个以闭区间 [2,3] 为支撑的制服。我们可以再次使用 scipy，这次是为了获得它的方差。其他时刻可用；请参阅文档。

>>> stats.uniform.stats(2,1, moments='v')
array(0.08333333333333333)

这些都不是真的必要，不是吗？

作为一个 U(0,1) 随机变量 X1 永远不会超过 1。因此，max(X1,2) 必须是 2。那么 2+X2 必须是 U(2,3)。这个随机变量的尺度与 X2 相同；只是它的位置发生了变化。所以它的方差一定是一样的，一个U(0,1)的方差是0.0833333。

编辑“下一天”：

刚刚了解到（来自https://stackoverflow.com/a/46383333/131187）sympy 现在支持随机变量，我很想尝试解决这个问题。

>>> from sympy.stats import Uniform, Variance
>>> from sympy import symbols, Integral
>>> X1 = Uniform('X1', 0, 1)
>>> X2 = Uniform('X2', 0, 1)

唉，作为其他答案的作者，它似乎无法处理涉及max 的表达式。

>>> Variance(max(2, X1) + X2)
Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
  File "C:\Python34\lib\site-packages\sympy-1.0.1.dev0-py3.4.egg\sympy\core\relational.py", line 195, in __nonzero__
    raise TypeError("cannot determine truth value of Relational")
TypeError: cannot determine truth value of Relational

但在这个问题的情况下，这不是必需的。它很容易被淘汰。我们有，这会产生方差积分的确切值。

>>> Variance(2 + X2)
Variance(X2 + 2)
>>> Variance(2 + X2).evaluate_integral()
1/12

“旧答案”从这里开始：

我认为不是直接的。但是，这种方法可能对您有用。

首先假设您知道感兴趣的随机变量函数的 pdf 或 cdf。然后你可以使用 scipy.stats 中的 rv_continuous 来计算该函数的方差和其他矩。

显然，“乐趣”从这里开始。通常您会尝试定义 cdf。对于随机变量的任何给定值，这是您给出的表达式不超过给定值的概率。因此，确定 cdf 简化为解决两个变量中的（无限）不等式集合。当然，通常有一种强大的模式可以大大降低执行这项任务的复杂性和难度。

Answer 2

在 OpenTURNS 中，使用 Symbolic functions 允许您进行更多种类的操作。

在您的情况下，x1 和 x2 将代表独立分布

import openturns as ot
x1 = ot.Uniform() 
x2 = ot.Uniform()

因此，边际为 x1 和 x2 的复合分布将是：

dist = ot.ComposedDistribution([x1, x2], ot.IndependentCopula(2))
dist.setDescription(["x1", "x2"])  # labels 

# note the use of "IndependentCopula of dimension 2" as second argument 

如果你想要一个 size = 5 的样本

sample = dist.getSample(5)
print(sample)
>>>     [ x1        x2        ]
0 : [ -0.752141 -0.897212 ]
1 : [  0.850966  0.857914 ]
2 : [ -0.340213 -0.344882 ]
3 : [ -0.166526  0.458643 ]
4 : [  0.378453 -0.908958 ]

如前所述，您可以将基于 (x1, x2) 的模型定义为符号函数。在您的示例中： y = max(2, x1) + x2

model = ot.SymbolicFunction(["x1", "x2"], ["max(2, x1) + x2"])

你可以申请model(sample)

    [ y0      ]
0 : [ 1.10279 ]
1 : [ 2.85791 ]
2 : [ 1.65512 ]
3 : [ 2.45864 ]
4 : [ 1.09104 ]

但是您的模型可以是多维的。例如：

model = ot.SymbolicFunction(["x1", "x2"], ["x1^2+x2", "x2^2+x1"])

应用于样本将给出二维样本

>>>    [ y0         y1      ]
0   -0.331496   0.05284813
1   1.582057    1.586982
2   -0.2291374  -0.2212693
3   0.4863741   0.04382738
4   -0.7657314  1.204657

这在创建更高级的模型时非常有趣。在最后一种情况下，绘制大小为 10,000 out = model(dist.getSample(10000)) 的输出给出

import matplotlib.pyplot as plt
plt.scatter(out.getMarginal(0),out.getMarginal(1), s=0.5)