créer des distributions dans mathematica

Question

J'ai une fonction que je sais être une distribution multivariée dans (x, y), et mathematica a des problèmes de stabilité numérique quand je forme les distributions marginales.

Par exemple, marginalisant le long y donne le résultat suivant: 0.e^(^ 154.88-0.5x 2)

Depuis que je sais que le résultat doit être une distribution, je voudrais extraire simplement l'e^(-.5x^2) et faire une renormalisation moi-même. Alternativement, ce serait encore mieux si mathematica me permettait de prendre une fonction multivariée et de le spécifier en quelque sorte comme une distribution de probabilité.

De toute façon, quelqu'un sait comment implémenter l'une des deux solutions ci-dessus par programmation?

Answer 1

Ok, voici un exemple de ce que je veux dire. Supposons que j'ai la distribution 2D suivante:

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

Et je tente de

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

Mathematica ne fournit pas de réponse, ou tout au moins ne le fait pas pour un certain temps sur mon ordinateur. Suggestions?

Edit: ok, donc ça le fait, mais ça prend 5 minutes sur mon Intel i5 avec 4 Go de RAM ... J'espère toujours avoir un moyen de puiser dans le type de distribution intégré de Mathematica (bien qu'il semble être une seule variable seulement) et faire usage de leur RandomReal [dist]. Le mieux que je puisse espérer est si Mathematica me laisse spécifier cette fonction 2D comme une distribution, et être capable d'appeler RandomRealVector [dist].

Answer 2

ProbabilityDistribution prend des fonctions multivariées, bien que votre fonction Dist soit un peu trop bizarre à son goût.

En outre, il semble que les distributions multivariées définies par l'utilisateur ne fonctionnent actuellement pas en combinaison avec RandomVariate (la version V8 légèrement plus polyvalente de RandomReal/RandomInteger). Les distributions univariées fonctionnent. J'ai soumis un rapport de bogue à WRI.

Answer 3

Eh bien, traitant des expressions symboliques dans Mathematica, il est préférable de garder les choses exactement, à savoir éviter le nombre approximatif:

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] 

Pour attaquer le problème, il est préférable de modifier les 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} 

Notez que les coefficients ont été choisis de telle sorte que jacobian est une unité:

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

Out[42]= 1 

Après le changement de variables, vous voyez que la densité factorise en un produit:

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] 

C'est, les variables maintenant « u » et « v » sont indépendants. La variable 'v' est NormalDistribution[0, 1/101], tandis que la variable 'u' est un peu plus compliquée, mais peut maintenant être traitée par 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]) 

Vous pouvez maintenant définir la distribution conjointe pour le vecteur {u,v}:

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

Depuis la relation entre {u,v} et {x,y} est connu, générant des {x,y} Taxipost est facile:

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

Vous pouvez encapsuler les connaissances accumulées en utilisant TransformedDistribution:

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

.: par exemple

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

Out[68]= {1.27198, 0.126733}