3

我有一个函数，我知道它是 (x,y) 中的多元分布，当我形成边际分布时，mathematica 存在数值稳定性问题。

例如，沿 y 边缘化会产生以下结果：0.e^(154.88-0.5x^2)

因为我知道结果必须是分布，所以我想只提取 e^(-.5x^2) 并自己进行重新归一化。或者，如果mathematica 允许我采用多元函数并以某种方式将其指定为概率分布，那就更好了。

无论如何，有谁知道如何以编程方式实现上述两种解决方案中的任何一种？

4

3 回答 3

1

好吧，在 Mathematica 中处理符号表达式，最好保持精确，即避免使用近似数字：

In[36]:= pdf = PiecewiseExpand[Rationalize[E^(-(x^2/2) - y^2/2)*
         (-1 + E^(-1.*(x + 0.1*y)*UnitStep[x + 0.1*y]))^2], 
  Element[{x, y}, Reals]]

Out[36]= Piecewise[{{E^(-2*x - x^2/2 - y/5 - y^2/2)*(-1 + 
       E^(x + y/10))^2, 10*x + y >= 0}}, 0]

为了解决问题，最好改变变量：

In[56]:= cvr = 
 First[Solve[{10 x + y == u, (10 y - x)/101 == v}, {x, y}]]

Out[56]= {x -> (10 u)/101 - v, y -> u/101 + 10 v}

请注意，选择系数以使雅可比是一个单位：

In[42]:= jac = Simplify[Det[Outer[D, {x, y} /. cvr, {u, v}]]]

Out[42]= 1

更改变量后，您会看到密度分解为乘积：

In[45]:= npdf = FullSimplify[jac*pdf /. cvr]

Out[45]= Piecewise[{{E^(-(u/5) - u^2/202 - (101*v^2)/2)*(-1 + 
       E^(u/10))^2, u >= 0}}, 0]

也就是说，现在变量“u”和“v”是独立的。'v' 变量是NormalDistribution[0, 1/101]，而 'u' 变量稍微复杂一些，但现在可以由ProbabilityDistribution.

In[53]:= updf = 
 Refine[npdf/nc, u >= 0]/PDF[NormalDistribution[0, 1/Sqrt[101]], v]

Out[53]= (E^(-(u/5) - u^2/202)*(-1 + E^(u/10))^2*Sqrt[2/(101*Pi)])/
   (1 - 2*E^(101/200)*Erfc[Sqrt[101/2]/10] + 
   E^(101/50)*Erfc[Sqrt[101/2]/5])

因此，您现在可以定义 vector 的联合分布{u,v}：

dist = ProductDistribution[NormalDistribution[0, 1/101], 
   ProbabilityDistribution[updf, {u, 0, Infinity}]];

{u,v}由于和之间的关系{x,y}是已知的，因此很容易生成{x,y}变量：

XYRandomVariates[len_] := 
 RandomVariate[dist, len].{{-1, 10}, {10/101, 1/101}}

您可以使用以下方法封装积累的知识TransformedDistribution：

origdist = 
  TransformedDistribution[{(10 u)/101 - v, 
    u/101 + 10 v}, {Distributed[v, NormalDistribution[0, 1/101]], 
    Distributed[u, ProbabilityDistribution[updf, {u, 0, Infinity}]]}];

例如：

In[68]:= Mean[RandomVariate[origdist, 10^4]]

Out[68]= {1.27198, 0.126733}
于 2011-02-19T19:52:27.870 回答
1

好的，这是我的意思的一个例子。假设我有以下二维分布：

Dist = 
3.045975040844157` E^(-(x^2/2) - y^2/
2) (-1 + E^(-1.` (x + 0.1` y) UnitStep[x + 0.1` y]))^2

我试图

Integrate[Dist, {y, -Infinity, Infinity}]

Mathematica 没有提供答案，或者至少在我的计算机上很长一段时间内没有提供答案。建议？

编辑：好的，确实如此，但在我的 Intel i5 上需要 5 分钟，4GB 内存...我仍然希望有一些方法可以利用 Mathematica 的内置分发类型（尽管它似乎只是单个变量）并制作使用他们的 RandomReal[dist]。我能希望的最好结果是 Mathematica 是否允许我将此二维函数指定为分布，并能够调用 RandomRealVector[dist]。

于 2009-10-19T00:36:01.353 回答
1

ProbabilityDistribution确实采用多元函数，尽管您的 Dist 函数对于它的口味来说有点太奇怪了。

此外，用户定义的多元分布目前似乎不能与RandomVariate（稍微通用的 V8 版本RandomReal/ RandomInteger）结合使用。单变量分布有效。我向世界资源研究所提交了一份错误报告。

于 2011-02-14T23:05:45.907 回答