Re: equiprobable intervals with triangular pdf
- To: mathgroup at smc.vnet.net
- Subject: [mg41474] Re: [mg41443] equiprobable intervals with triangular pdf
- From: Bobby Treat <drmajorbob+MathGroup3528 at mailblocks.com>
- Date: Wed, 21 May 2003 08:03:39 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
Here's a faster, simpler way, though perhaps not as accurate.
<< Statistics`ContinuousDistributions`
Unprotect[PDF, CDF, Quantile, InterpolatedQuantile];
Clear[PDF, CDF, Quantile, InterpolatedQuantile]
PDF[TriangularDistribution[a_, b_], x_] := ((4*UnitStep[x -
a]*UnitStep[b - \
x])/(b - a)^2)*((x - a)*UnitStep[(a + b)/2 - x] + (b - x)*UnitStep[
x - (a + b)/2])
CDF[TriangularDistribution[(a_)?NumericQ, (
b_)?NumericQ], (x_)?NumericQ] :=
Block[{t}, NIntegrate[PDF[TriangularDistribution[a, b], t],
{t, a,
x}]]
Quantile[dist :
TriangularDistribution[(a_)?NumericQ, (b_)?NumericQ], (q_)?
NumericQ] := Block[{x, $DisplayFunction = Identity},
If[Head[InterpolatedQuantile[dist]] ===
InterpolatedQuantile, \
Unprotect[InterpolatedQuantile];
InterpolatedQuantile[dist] = Interpolation[Reverse /@ \
First[Cases[Plot[CDF[dist, x], {x, a, b}],
Line[pts_] :> pts, Infinity]]]; Protect[InterpolatedQuantile]];
InterpolatedQuantile[dist][q]]
Protect[PDF, CDF, Quantile, InterpolatedQuantile];
dist = TriangularDistribution[0, 1];
Plot[{Quantile[dist, x]}, {x, 0, 1}];
That can be made more accurate by adding PlotPoints->100 (for instance)
to the Plot used in defining Quantile. Error in the following can be
attributed to inaccuracy in either CDF or Quantile, but it's mostly in
Quantile:
dist = TriangularDistribution[0, 1];
Plot[{x - CDF[dist, Quantile[dist, x]]}, {x, 0, 1}, PlotRange -> All];
The more complicate version, with InterpolatedRoot, is about three
orders of magnitude more accurate than the above with PlotPoints->100.
Bobby
-----Original Message-----
From: Bobby Treat <drmajorbob+MathGroup3528 at mailblocks.com>
To: mathgroup at smc.vnet.net
Subject: [mg41474] Re: [mg41443] equiprobable intervals with triangular pdf
Here's one way to do it:
<< "Statistics`ContinuousDistributions`"
Unprotect[PDF, CDF, Quantile];
Clear[PDF, CDF, Quantile]
PDF[TriangularDistribution[a_, b_], x_] :=
((4*UnitStep[x - a]*UnitStep[b - x])/(b - a)^2)*
((x - a)*UnitStep[(a + b)/2 - x] +
(b - x)*UnitStep[x - (a + b)/2])
CDF[TriangularDistribution[(a_)?NumericQ,
(b_)?NumericQ], (x_)?NumericQ] :=
Block[{t}, NIntegrate[
PDF[TriangularDistribution[a, b], t],
{t, a, x}]]
Quantile[dist:TriangularDistribution[(a_)?NumericQ,
(b_)?NumericQ], (q_)?NumericQ] :=
Which[q >= 1., b, q <= 0., a, True,
Block[{x}, x /. FindRoot[CDF[dist, x] - q,
{x, (2*a + b)/3, (a + 2*b)/3},
MaxIterations -> 30]]]
Protect[PDF, CDF, Quantile];
dist = TriangularDistribution[0, 1];
Plot[{PDF[dist, xx], Quantile[dist, xx]}, {xx, 0, 1}];
Then, to get equiprobably intervals, use Quantile as before.
Bobby
-----Original Message-----
From: S White <susanlcw at aol.com>
To: mathgroup at smc.vnet.net
Subject: [mg41474] [mg41443] equiprobable intervals with triangular pdf
Hello all,
I posted a couple of weeks ago about dividing a normal distribution
into n equiprobable intervals and received some great responses. I am
now working with a triangular pdf and need to do the same thing.
I am defining the triangular pdf on the interval [a,b] with mean at
(a+b)/2 as follows:
triPdf[x_,a_,b_]:=(2/(b-a))^2*(x-a)/;a<=x<(b+a)/2;
triPdf[x_,a_,b_]:=(2/(b-a))^2*(b-x)/;(b+a)/2<=x<=b;
triPdf[x_,a_,b_]:=0/;a>x||x>b;
triCdf[x_,a_,b_]:=N[Integrate[triPdf[y,a,b],{y,a,x}]]
When working with the normal command, the following function gave me
output in the form I need:
equiprob[dist_,n_]:=Partition[Table[Quantile[dist,k/n],{k,0,n}],2,1]
However, the Quantile function doesn't work on this triangular pdf.
I have worked on defining some function that would do the same thing
as the Quantile function does but I keep running into numerous error
messages and it only works for certain a and b.
Does anyone have a suggestion of a function that will give me
equiprobable intervals in the output form
{{a,x1},{x1,x2},...,{xn,b}},
where a,x1,...,xn,b are the endpoints of the equiprobable intervals?
I really appreciate any help...
Susan