我有一个函数，我知道它是 (x,y) 中的多元分布，当我形成边际分布时，mathematica 存在数值稳定性问题。
例如，沿 y 边缘化会产生以下结果： 0.e^(154.88-0.5x^2)
因为我知道结果必须是分布，所以我想只提取 e^(-.5x^2) 并自己进行重整化。或者，如果mathematica 允许我采用多元函数并以某种方式将其指定为概率分布，那就更好了。
无论如何，有人知道如何以编程方式实现上述两种解决方案吗？
【问题讨论】：
标签： distribution wolfram-mathematica probability
好的，这是我的意思的一个例子。假设我有以下二维分布：
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 允许我将这个 2D 函数指定为分布，并且能够调用 RandomRealVector[dist]。
【讨论】：
ProbabilityDistribution 确实采用了多元函数，尽管您的 Dist 函数对于它的口味来说有点太奇怪了。
此外，用户定义的多元分布目前似乎无法与RandomVariate（RandomReal/RandomInteger 的功能稍多的 V8 版本）结合使用。单变量分布有效。我向世界资源研究所提交了一份错误报告。
【讨论】：
好吧，在 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])
所以你现在可以定义向量{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}
【讨论】：