3
$\begingroup$

What I want to do:
1) Define a joint distribution function using CopulaDistribution
2) then transform it with TransformedDistribution
3) to get the new Culmulative Distribution Function (CDF).

The problem:
If I choose correlation 0 or the Product Copula, everything is working perfectly, I obtain a CDF. If I choose a copula different from independence, I don't get a result anymore. In the code below, change the parameter for the copula from 0. to e.g. 0.2, then there's no result for the CDF anymore. (I'm using Mathematica Version 8.0)

The code:

dist = CopulaDistribution[{"Binormal", 0.}, {UniformDistribution[{0, 10}], 
    UniformDistribution[{0, 10}]}]; 

f[x_, y_] = 
  Piecewise[{{x - 3, x > 3 && y > 3}, {x - 6 + y, 
     x > 3 && y <= 3 && x + y - 6 >= 0}}]; 

transdist = 
  TransformedDistribution[f[x, y], {x, y} \[Distributed] dist];

CDF[transdist, 2.]
$\endgroup$

2 Answers 2

2
$\begingroup$

Copulas (with dependent structures) can result in complicated functional forms. Your example highlights why working with functions as BLACK BOXES can rapidly get one into trouble. It is always a good idea to check both the input and output from such functions at every stage, so that we have some idea of what we are actually working with, and what we could reasonably hope to attain.

For your Binormal copula, for general parameter $a$, the pdf has a complicated form ... check it using PDF[dist, {x, y}] ... which yields the joint pdf of $(X,Y)$ as say $g(x,y)$:


(source: tri.org.au)

Your transformation equation is:

  f[x_, y_] = Piecewise[{{x-3, x>3 && y>3}, {x-6+y, x>3 && y<=3 && x+y-6 >= 0}}]


(source: tri.org.au)

You seek the cdf of $Z = f(x,y)$, i.e. $P(f(x,y) < z)$:


(source: tri.org.au)

where I am using the Prob function from the mathStatica package for Mathematica. Although Mathematica cannot obtain a closed-form, we now have the functional form for the cdf, and it is easy to obtain solutions for given parameter values. For example, here is the solution expressed as a numerical cdf (just swapping Integrate for NIntegrate):


(source: tri.org.au)

[To zoom in on the above, just open the picture in a new tab or window.]

For your problem, parameter $a=1/5$ and you seek $P(Z<2)$:

 numCDF[2, 1/5]

0.540261

Here is a quick plot of the CDF (when $a=\frac15$):

Plot[numCDF[z, 1/5], {z, -2, 8}, PlotRange -> All, AxesLabel -> {"z", "cdf"}]


(source: tri.org.au)

It looks almost linear! :) Note the discrete mass at $Z = 0$ that arises due to the structure of your transformation equation f[x,y].

$\endgroup$
1
$\begingroup$

I post this just for version 10. It appears to handle this. I would guided by @wolfies as to correctness. Using the definitions from OP:

Plot[Evaluate[CDF[transdist, z]], {z, 0, 10}]

enter image description here

The symbolic definitions:

enter image description here

Note all the following yield 0.545:

CDF[transdist, 2]
Integrate[PDF[transdist, x], {x, 0, 2}] /.HeavisideTheta[x_] -> UnitStep[x]
Probability[f[x, y] < 2, {x, y} \[Distributed] dist]
Probability[z < 2, z \[Distributed] transdist]

and the following yielded: 0.5449

N[Length[Cases[#, _?(# < 2 &)]]/Length[#]] &@RandomVariate[transdist, 100000]

and finally as this just treating as uniform distribution:

reg = ImplicitRegion[
  Reduce[f[x, y] < 2 && 0 < x < 10 && 0 < y < 10, {x, y}], {x, y}];
rp = RegionPlot[reg]
NIntegrate[1, {x, y} \[Element] reg]/100.

The region:

enter image description here

and the integral: 0.545

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.