MathGroup Archive 2009

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: piecewice pdf, problems with cdf

  • To: mathgroup at smc.vnet.net
  • Subject: [mg105373] Re: [mg105365] piecewice pdf, problems with cdf
  • From: DrMajorBob <btreat1 at austin.rr.com>
  • Date: Mon, 30 Nov 2009 06:10:05 -0500 (EST)
  • References: <200911291012.FAA16385@smc.vnet.net>
  • Reply-to: drmajorbob at yahoo.com

I didn't normalize that correctly and it's a poor example when plotted, so  
here's a slightly different example:

Clear[f, tmp, x, cdf]
tmp[x_] =
   Sqrt[2] + UnitStep[x + 2] Cos@x + UnitStep[x - 3] Sin@x +
    UnitStep[x - 5] Exp[x - 7] - UnitStep[x - 10] Exp[x - 7];
f[x_] = tmp@x/Integrate[tmp@v, {v, -2, 10}];
cdf[x_] = Integrate[f@v, {v, -3, x}];
Plot[{f@x, cdf@x}, {x, -3, 11}, PlotRange -> All,
  AxesOrigin -> {0, 0}]

This won't do for every problem, however. For instance:

Clear[f, tmp, x, cdf]
tmp[x_] =
   Sqrt[2] + UnitStep[x + 2] Cos@x + UnitStep[x - 3] Sin@x +
    UnitStep[x - 5] Log@x - UnitStep[x - 10] Log@x;
f[x_] = tmp@x/Integrate[tmp@v, {v, -2, 10}]
Plot[f@v, {v, -3, 10}]

(Sqrt[2] - Log[x] UnitStep[-10 + x] + Log[x] UnitStep[-5 + x] +
    Sin[x] UnitStep[-3 + x] + Cos[x] UnitStep[2 + x])/(-5 +
    12 Sqrt[2] + Cos[3] - Cos[10] + 5 Log[20] + Sin[2] + Sin[10])

So far so good, but Integrate can't do the next step!

You could use NIntegrate in theory, but it's too slow in a Plot, as in  
your original code.

NDSolve does the trick instead:

Clear[cdf]
cdf[x_] /; x <= -2 = 0;
cdf[x_] /; x >= 10 = 1;
cdf[x_] =
   g[x] /. First@NDSolve[{g'[x] == f@x, g[-2] == 0}, g[x], {x, -2, 10}];
Plot[{f@x, cdf@x}, {x, -2, 10}, PlotRange -> All,
  AxesOrigin -> {0, 0}]

If the first Integrate (used to normalize "tmp") happens to fail, you can  
try NIntegrate, since it's only calculated once.

NIntegrate[tmp@v, {v, -2, 10}]

27.1636

Bobby

On Sun, 29 Nov 2009 21:32:04 -0600, DrMajorBob <btreat1 at austin.rr.com>  
wrote:

> I really meant to post this solution:
>
> Clear[x, f1, cdf]
> f1[x_] = Piecewise[{{x^2/9, 0 <= x <= 3}}];
> cdf[x_?Negative] = 0;
> cdf[x_] /; x >= 3 = 1;
> cdf[x_] = Assuming[0 < x < 3, Integrate[f1@v, {v, -Infinity, x}]]
> Plot[cdf@x, {x, -5, 5}, PlotRange -> All]
>
> If there are many pieces, UnitStep is your friend. For instance, here's  
> a function that changes at -2, +3, +5, and +10, where I start with  
> something arbitrary and normalize by dividing by the integral:
>
> Clear[f, tmp, x, cdf]
> tmp[x_] =
>    UnitStep[x + 2] Sin@x + UnitStep[x - 3] Cos@x +
>     UnitStep[x - 5] Exp@x - UnitStep[x - 10] Exp@x;
> f[x_] = tmp@x/Integrate[tmp@v, {v, 0, 10}]
> cdf[x_] = Integrate[f@v, {v, 0, x}]
>
> Integrate is a SYMBOLIC solver, so it's inherently limited. It doesn't  
> understand functions built with If.
>
> Bobby
>
> Plot[cdf@x, {x, -3, 11}, PlotRange -> All]
>
> On Sun, 29 Nov 2009 20:07:42 -0600, michael partensky  
> <partensky at gmail.com> wrote:
>
>> Thanks, Bobby.
>>
>> Two questions though.
>> (1) Why some other definitions  would not work with Integrate.
>> I know the derivative has an issue with Delta function, but what's the
>> problem with the  Integrate given that all values of the integrand  
>> (pdf) are
>> defined.
>> At least, it should work if each piece is handled individually.
>> (2) What if we  have many pieces? Should we enter all them via  
>> "Assuming"
>> (which could be inconvenient). After all, the assumptions reproduce  
>> x-ranges
>> that are already in  the function definition. Can we directly   
>> Integrate pdf
>> to produce  a piecewise CDF function that will properly return for any  
>> value
>> of x (automatically locating it in a proper piece) ?
>>
>> Thanks
>> Michael
>>
>>
>>
>> On Sun, Nov 29, 2009 at 7:49 PM, DrMajorBob <btreat1 at austin.rr.com>  
>> wrote:
>>
>>> If it's a piecewise function... use Piecewise!!
>>>
>>> Clear[x, f1, cdf]
>>> f1[x_] = Piecewise[{{x^2/9, 0 <= x <= 3}}];
>>> cdf[x_] =
>>>  Piecewise[{{Assuming[0 < x <= 3,
>>>      Integrate[f1@v, {v, -Infinity, x}]], 0 < x <= 3}}];
>>>
>>> Plot[cdf@x, {x, -5, 5}]
>>>
>>> Bobby
>>>
>>>
>>> On Sun, 29 Nov 2009 04:12:35 -0600, michael partensky  
>>> <partensky at gmail.com>
>>> wrote:
>>>
>>>  Hi! Teaching the continuous distributions, I needed to introduce the
>>>> piecewise functions.
>>>> Here is the example that did not work well:
>>>>
>>>> In[56]:= f1[x_] /; 0 < x <= 3 := 1/9  x ^2;
>>>> f1[x_] := 0;
>>>>
>>>> Plot[f1[x],{x,-1,4}] works fine. However, the results for cdf are
>>>> ambiguous
>>>> In[57]:= cdf[x_] := Integrate[f1[v], {v, -\[Infinity], x}]
>>>>
>>>> In[59]:= cdf[1]
>>>> Out[59]= 0
>>>>
>>>> I thought that may be the second definition (for some reason)  
>>>> overwrote
>>>> the
>>>> first, but apparently this was not the case.
>>>>
>>>> Then I tried using Which,
>>>>
>>>> f1[x_] := Which[0 < x <= 3, x^2/9, x <= 0 || x > 3, 0];
>>>>
>>>> Plot[f2[x], {x, -1, 4}] worked fine.
>>>>
>>>> However, Plotting CDF is very slow.
>>>>
>>>> What is the reason for the first error and how to accelerate  
>>>> (compile?)
>>>>  the
>>>> second?
>>>>
>>>> Thanks
>>>> Michael
>>>>
>>>> PS: I was aware about the issues with the derivatives of Piecewise
>>>> functions, but expected  integration to be safe. What did i do wrong?
>>>>
>>>> On Tue, Nov 24, 2009 at 5:50 AM, Bill Rowe <readnews at sbcglobal.net>
>>>> wrote:
>>>>
>>>>  On 11/23/09 at 6:53 AM, lshifr at gmail.com (Leonid Shifrin) wrote:
>>>>>
>>>>> >Hi, Michael.
>>>>>
>>>>> >Your solution is indeed very memory -  hungry due to the mapping of
>>>>> >Permutations, as you mentioned. The total number of permutations can
>>>>> >be easily deduced from the list of multiplicities of elements in a
>>>>> >given partition: n!/(n1!n2!...nk!), where n1, ..., nk are
>>>>> >multiplicities of elements, and n is the length of the partition:
>>>>> >n=n1+...+nk. The multiplicities can be obtained by Tally. The
>>>>> >following modification can be realistically used in a much wider
>>>>> >region of the problem's parameter space than your original one, and
>>>>> >may possibly serve your needs.
>>>>>
>>>>> >In[1]:=
>>>>> >Clear[outsNew];
>>>>> >outsNew[sum_, thr_] :=
>>>>> >Total[Factorial[Length[#]]/
>>>>> >Times @@ Factorial[Tally[#][[All, 2]]] & /@
>>>>> >Cases[IntegerPartitions[sum, thr, {1, 2, 3, 4, 5, 6}],
>>>>> >Table[_, {thr}]]];
>>>>>
>>>>> The above can be made more efficient by using the form of
>>>>> IntegerPartitions that only generates partitions of a given
>>>>> length. That is
>>>>>
>>>>> outs[sum_, thr_] :=
>>>>>   Total[Factorial[Length[#]]/
>>>>>      Times @@ Factorial[Tally[#][[All, 2]]] & /@
>>>>>     IntegerPartitions[sum, {thr}, {1, 2, 3, 4, 5, 6}]]
>>>>>
>>>>> does the same thing a bit more efficiently
>>>>>
>>>>>
>>>>>
>>>>>
>>>>
>>>
>>> --
>>> DrMajorBob at yahoo.com
>>>
>
>


-- 
DrMajorBob at yahoo.com


  • Prev by Date: Re: piecewice pdf, problems with cdf
  • Next by Date: Re: piecewice pdf, problems with cdf
  • Previous by thread: Re: piecewice pdf, problems with cdf
  • Next by thread: Re: piecewice pdf, problems with cdf