Creating Distributions in Math

Question

I have a function that, as I know, is a multidimensional distribution in (x, y), and mathematics has problems with numerical stability when I form marginal distributions.

For example, marginalization along y gives the following: 0.e ^ (154.88-0.5x ^ 2)

Since I know that the result should be a distribution, I would like to extract only e ^ (-. 5x ^ 2) and do the renormalization myself. Alternatively, it would be even better if mathematics allowed me to take a multidimensional function and somehow indicate it as a probability distribution.

Anyway, does anyone know how to programmatically implement either of these two solutions?

Answer 1

Well, here is an example of what I mean. Suppose I have the following two-dimensional 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'm trying

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

Mathematica does not give an answer or at least does not do this for a while on my computer. Suggestions?

Editing: ok, so it really is, but it takes 5 minutes on my Intel i5 processor with 4 GB of RAM ... I still hope they have some way to connect to Mathematica, built by distribution type (although it seems that this is only one variable) and use their RandomReal [dist]. The best I can hope for is that Mathematica will allow me to specify this 2D function as a distribution and be able to call RandomRealVector [dist].

Answer 2

ProbabilityDistribution , Dist .

, , RandomVariate ( V8 RandomReal/RandomInteger). . WRI.

Answer 3

, 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])

So now you can define the joint distribution for the vector {u,v}:

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

Since the connection between {u,v}and is {x,y}known, the generation of variations {x,y}is simple:

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

You can encapsulate 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}