creating distributions in mathematica
I have a function which I know to be a multivariate distribution in (x,y), and mathematica is having numerical stability issues when I form the marginal distributions.
For example, marginalizing along y yields the following: 0.e^(154.88-0.5x^2)
Since I know the result must be a distribution, I would like to extract just the e^(-.5x^2) and do a renormalization myself. Alternatively, it would be even nicer if mathematica would let me take a multivariate function and somehow specify it as a probability distribution.
Anyway, does anyone know how to implement either of the above two solutions programatically?
Ok, here's an example of what I mean. Suppose I have the following 2D distribution:
Dist =
3.045975040844157` E^(-(x^2/2) - y^2/
2) (-1 + E^(-1.` (x + 0.1` y) UnitStep[x + 0.1` y]))^2
And I attempt to
Integrate[Dist, {y, -Infinity, Infinity}]
Mathematica does not provide an answer, or at least doesn't do so for quite a while on my computer. Suggestions?
Edit: ok, so it actually does, but takes 5 minutes on my Intel i5 with 4GB ram... I am still hoping theres some way to tap into Mathematica's built in distribution type (though it seems to be single variable only) and make use of their RandomReal[dist]. The best I could hope for is if Mathematica would let me specify this 2D function as a distribution, and be able to call RandomRealVector[dist].
ProbabilityDistribution does take multivariate functions, though your Dist function is a bit too weird for its taste.
Additionaly, it seems that user-defined multivariate distributions currently don't work in combination with RandomVariate (the slightly more versatile V8 version of RandomReal / RandomInteger ). Univariate distributions work. I submitted a bug report to WRI.
Well, dealing with symbolic expressions in Mathematica, it's best to keep things exact, ie avoid approximate numbers:
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 order to attack the problem, it is better to change variables:
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}
Notice that coefficients were chosen so that jacobian is a unity:
In[42]:= jac = Simplify[Det[Outer[D, {x, y} /. cvr, {u, v}]]]
Out[42]= 1
After the change of variables, you see that the density factorizes into a product:
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]
That is, now variables 'u' and 'v' are independent. The 'v' variable is NormalDistribution[0, 1/101] , while the 'u' variable is a little more complicated, but can now be handled by 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])
So you can now define the joint distribution for vector {u,v} :
dist = ProductDistribution[NormalDistribution[0, 1/101],
ProbabilityDistribution[updf, {u, 0, Infinity}]];
Since the relationship between {u,v} and {x,y} is known, generating of {x,y} variates is easy:
XYRandomVariates[len_] :=
RandomVariate[dist, len].{{-1, 10}, {10/101, 1/101}}
You can encapsulate the accumulated knowledge using TransformedDistribution :
origdist =
TransformedDistribution[{(10 u)/101 - v,
u/101 + 10 v}, {Distributed[v, NormalDistribution[0, 1/101]],
Distributed[u, ProbabilityDistribution[updf, {u, 0, Infinity}]]}];
Eg:
In[68]:= Mean[RandomVariate[origdist, 10^4]]
Out[68]= {1.27198, 0.126733}
下一篇: 在mathematica中创建分布