在mathematica中创建分布
我有一个函数，我知道它是（x，y）中的多元分布，当我形成边际分布时，数学是有数值稳定性问题的。
例如，沿y的边际化产生以下结果：0.e ^（154.88-0.5x ^ 2）
既然我知道结果必须是一个分布，我想只提取e ^（ - 。5x ^ 2）并且自己做一个重整化。 另外，如果数学可以让我采取多元函数并以某种方式将其指定为概率分布，那将会更好。
无论如何，有没有人知道如何以编程方式实现上述两种解决方案之一？
好的，这里是我的意思的一个例子。 假设我有以下2D分布：
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不提供答案，或者至少在我的计算机上不会这么做。 建议？
编辑：好的，所以它实际上是这样，但需要5分钟，我的英特尔I5与4GB内存...我仍然希望能够采用Mathematica的内置分配类型（尽管它似乎是单变量），并使使用他们的RandomReal [dist]。 我希望的最好的方式是，如果Mathematica允许我将这个2D函数指定为分布，并且能够调用RandomRealVector [dist]。
ProbabilityDistribution确实需要多元函数，尽管你的Dist函数对它的味道有点太古怪了。
另外，似乎用户定义的多变量分布目前不能与RandomVariate （稍微更通用的V8版本的RandomReal / RandomInteger ）结合使用。 单变量分布工作。 我向WRI提交了一个错误报告。
那么，在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}
上一篇: creating distributions in mathematica
下一篇: Calculating expectation for a custom distribution in Mathematica