Re: Re: equiprobable intervals with triangular pdf [LONG]
- To: mathgroup at smc.vnet.net
- Subject: [mg41549] Re: Re: equiprobable intervals with triangular pdf [LONG]
- From: bobhanlon at aol.com (Bob Hanlon)
- Date: Mon, 26 May 2003 05:46:27 -0400 (EDT)
- References: <baiar6$dks$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Rather than define separate names for the limiting cases you can handle special
parameters separately, e.g.,
PDF[TriangularDistribution[a_, b_, b_], x_] := ...
Needs["Statistics`ContinuousDistributions`"];
(* called to avoid shadowing function names *)
Needs["Graphics`Colors`"];
Clear[TriangularDistribution];
TriangularDistribution::usage =
"TriangularDistribution[a, b, c] represents the triangular distribution \
on interval [a, b] with mode c. The default for the mode is the midpoint \
(a+b)/2.";
TriangularDistribution[a_, b_] :=
TriangularDistribution[a, b, (a+b)/2];
(* sets default for mode to the midpoint (a+b)/2 *)
TriangularDistribution /: PDF[
TriangularDistribution[a_, b_, a_],
x_] :=
(2*(b-x)*(UnitStep[x-a]-UnitStep[x-b]))/(a-b)^2;
(* Limit[PDF[TriangularDistribution[a, b, c], x], c->a] *)
TriangularDistribution /: PDF[
TriangularDistribution[a_, b_, b_],
x_] :=
(2*(x-a)*(UnitStep[x-a]-UnitStep[x-b]))/(b-a)^2;
(* Limit[PDF[TriangularDistribution[a, b, c], x], c->b] *)
TriangularDistribution/:PDF[
TriangularDistribution[a_,b_,c_],
x_]:=-((2*((b-c)*(a-x)*UnitStep[x-a]-(a-c)*(b-x)*UnitStep[x-b]+
(a-b)*(c-x)*UnitStep[x-c]))/((a-b)*(a-c)*(b-c)));
TriangularDistribution /: CDF[
TriangularDistribution[a_, b_, a_],
x_] :=
((-UnitStep[a-b])*(a-b)^2+(a-x)*(a-2*b+x)*UnitStep[x-a]+
(b-x)^2*UnitStep[x-b])/(a-b)^2;
TriangularDistribution /: CDF[
TriangularDistribution[a_, b_, b_],
x_] :=
((-UnitStep[a-b])*(a-b)^2+(a-x)^2*UnitStep[x-a]+
(2*a-b-x)*(x-b)*UnitStep[x-b])/(a-b)^2;
(* Limit[CDF[TriangularDistribution[a, b, c], x], c->b] *)
TriangularDistribution /: CDF[
TriangularDistribution[a_, b_, c_],
x_] :=
(1/((a-b)*(a-c)*(b-c)))*(b*UnitStep[x-a]*a^2-
c*UnitStep[x-a]*a^2-
2*b*x*UnitStep[x-a]*a+2*c*x*UnitStep[x-a]*a-
b^2*UnitStep[x-b]*a-x^2*UnitStep[x-b]*a+
2*b*x*UnitStep[x-b]*a+(a-b)^2*(a-c)*
UnitStep[a-b]-(a-b)*(a-c)^2*UnitStep[a-c]+
b*x^2*UnitStep[x-a]-c*x^2*UnitStep[x-a]+
c*x^2*UnitStep[x-b]+b^2*c*UnitStep[x-b]-
2*b*c*x*UnitStep[x-b]+(a-b)*(c-x)^2*UnitStep[x-c]);
TriangularDistribution /: MGF[
TriangularDistribution[a_, b_, a_],
t_] :=
(2*(E^(b*t)+E^(a*t)*(-1+a*t-b*t)))/((a-b)^2*t^2);
(* Limit[MGF[TriangularDistribution[a, b, c]], c->a] *)
TriangularDistribution /: MGF[
TriangularDistribution[a_, b_, b_],
t_] :=
(2*(E^(a*t)+E^(b*t)*(-1-a*t+b*t)))/((a-b)^2*t^2);
(* Limit[MGF[TriangularDistribution[a, b, c]], c->b] *)
TriangularDistribution /: MGF[
TriangularDistribution[a_, b_, c_],
t_] :=
(2*(c*(-E^(a*t)+E^(b*t))+b*(E^(a*t)-E^(c*t))+
a*(-E^(b*t)+E^(c*t))))/((a-b)*(a-c)*(b-c)*t^2);
(* <Exp[X*t]> *)
TriangularDistribution /: Moment[
TriangularDistribution[a_, b_, a_],
n_] :=
(2*((a*(n+1)-b*(n+2))*a^(n+1)+b^(n+2)))/
((a-b)^2*(n+1)*(n+2));
(* Limit[Moment[TriangularDistribution[a, b, c], n], c->a] *)
TriangularDistribution /: Moment[
TriangularDistribution[a_, b_, b_],
n_] :=
(2*((n*b+b-a*(n+2))*b^(n+1)+a^(n+2)))/
((a-b)^2*(n+1)*(n+2));
(* Limit[Moment[TriangularDistribution[a, b, c], n], c->b ] *)
TriangularDistribution /: Moment[
TriangularDistribution[a_, b_, c_],
n_] :=
(2*(b*a^(n+2)-c*a^(n+2)-b^(n+2)*a+
(a-b)*c^(n+2)+b^(n+2)*c))/
((a-b)*(a-c)*(b-c)*(n+1)*(n+2));
TriangularDistribution /: HarmonicMean[
TriangularDistribution[a_, b_,
a_]] :=
(a-b)^2/(2*(a-b-b*Log[a]+b*Log[b]));
(* Limit[HarmonicMean[TriangularDistribution[a, b, c]], c->a] *)
TriangularDistribution /: HarmonicMean[
TriangularDistribution[a_, b_,
b_]] :=
(a-b)^2/(2*(-a+b+a*Log[a]-a*Log[b]));
(* Limit[HarmonicMean[TriangularDistribution[a, b, c]], c->b] *)
TriangularDistribution /: HarmonicMean[
TriangularDistribution[a_, b_, c_]] :=
((a-b)*(a-c)*(b-c))/
(2*(a*(b-c)*Log[a]+b*(c-a)*Log[b]+(a-b)*c*Log[c]));
TriangularDistribution /: GeometricMean[
TriangularDistribution[a_, b_, a_]] :=
a^((a*(a-2*b))/(a-b)^2)*b^(b^2/(a-b)^2)*
E^((a-3*b)/(2*b-2*a));
(* Limit[GeometricMean[TriangularDistribution[a, b, c]], c->a] *)
TriangularDistribution /: GeometricMean[
TriangularDistribution[a_, b_, b_]] :=
a^(a^2/(a-b)^2)*b^((b*(b-2*a))/(a-b)^2)*
E^((b-3*a)/(2*(a-b)));
(* Limit[GeometricMean[TriangularDistribution[a, b, c]], c->b] *)
TriangularDistribution /: GeometricMean[
TriangularDistribution[a_, b_,
c_]] :=
(a^(a^2/((a-b)*(a-c)))*c^(c^2/((a-c)*(b-c))))/
(b^(b^2/((a-b)*(b-c)))*E^(3/2));
(* Exp[<Log[X]>] *)
TriangularDistribution /: Mean[
TriangularDistribution[a_, b_, c_]] :=
(a+b+c)/3;
TriangularDistribution /: RMS[
TriangularDistribution[a_, b_, c_]] :=
Sqrt[(a^2+(b+c)*a+b^2+c^2+b*c)/6];
TriangularDistribution /: Variance[
TriangularDistribution[a_, b_,
c_]] :=
(1/18)*(a^2-(b+c)*a+b^2+c^2-b*c);
TriangularDistribution /: StandardDeviation[
TriangularDistribution[a_, b_, c_]] :=
Sqrt[(1/18)*(a^2-(b+c)*a+b^2+c^2-b*c)];
TriangularDistribution /: Mode[
TriangularDistribution[a_, b_, c_]] := c;
TriangularDistribution /: Median[
TriangularDistribution[a_, b_,
c_]] :=
(1/
2)*(UnitStep[c-(a+b)/2]*(2*a+Sqrt[2]*Sqrt[(a-b)(a-c)])+
(1-UnitStep[c-(a+b)/2])*(2*b-Sqrt[2]*Sqrt[(b-a)(b-c)]));
TriangularDistribution/:Quantile[
TriangularDistribution[a_,b_,c_],q_]:=
(a+Sqrt[(-a+b)*(-a+c)*q])*
(1-UnitStep[-((a-c)/(a-b))+q])+
(b-Sqrt[(-a+b)*(b-c)*(1-q)])*
UnitStep[-((a-c)/(a-b))+q]/;0<=q<=1;
TriangularDistribution/:equiProb[
TriangularDistribution[a_,b_,c_],n_Integer?Positive]:=
Partition[Table[Quantile[
TriangularDistribution[a,b,c],k/n],{k,0,n}],2,1];
With[{a=2,b=5, n=4},
dist = TriangularDistribution[a, b, c];
d = (b-a)/n;
Plot[Evaluate[Transpose[
Table[{PDF[dist,x],CDF[dist,x]},
{c,a,b,d}]]],
{x,a-.1,b+.1},
PlotStyle->Table[Hue[(5+3k/n)/10],{k,0,n}],
Frame->True, Axes->False,
FrameLabel -> {"\nx", None,
" PDF and CDF with Median for\nTriangularDistribution[" <>
ToString[a] <> ", " <> ToString[b] <> ", c]", None},
Prolog->{AbsolutePointSize[4],Red,
Point /@ Table[{Median[dist],1/2},
{c,a,b,d}]}, ImageSize->400]];
With[{a=2,b=5,n=4},
dist = TriangularDistribution[a, b, c];
d = (b-a)/n;
Plot[Evaluate[
Table[
Quantile[dist,q],
{c,a,b,d}]],
{q,0,1},
PlotStyle->Table[Hue[(5+3k/n)/10],{k,0,n}],
AspectRatio->GoldenRatio, ImageSize->247,
Frame->True, Axes->False,
FrameLabel -> {"\nq", "x\n",
" Quantile and Median for\nTriangularDistribution[" <>
ToString[a] <> ", " <> ToString[b] <> ", c]", None},
Prolog->{Text["c = a = " <> ToString[a], {.6,
Quantile[TriangularDistribution[a, b, a],.6]}, {-1,1}],
Text["c = b = " <> ToString[b], {.4,
Quantile[TriangularDistribution[a, b, b],.4]}, {1,-1}],
AbsolutePointSize[4],Red,
Point /@ Table[{1/2,Median[dist]},
{c,a,b,d}]}]];
With[{a = 2, b = 5, n = 4},
dist = TriangularDistribution[a, b, c];
d = (b - a)/n;
Plot[Evaluate[Table[Moment[dist, t]^(1/t),
{c, a, b, d}]], {t, -3/2, 5/2},
Frame -> True,
Axes -> False,
PlotStyle -> Table[Hue[(5 + 3*(k/n))/10], {k, 0, n}],
FrameLabel -> {"t", "<X^t>^(1/t)\n", StringJoin[
"Generalized Mean & Special Cases\nTriangularDistribution[",
ToString[a], ", ", ToString[b], ", c]"], None},
Prolog -> {
Text["\!\(\[Mu]\_H\)",
{-1, HarmonicMean[dist2 = TriangularDistribution[a, b, a + d/2]]},
TextStyle -> {FontSize -> 14}],
Text["\!\(\[Mu]\_G\)", {0, GeometricMean[dist2]},
TextStyle -> {FontSize -> 14}],
Text["\[Mu]", {1, Mean[dist2]}, TextStyle -> {FontSize -> 14}],
Text["RMS", {2, RMS[dist2]}],
Text[StringJoin["c = a = ", ToString[a]],
{1.5, GeometricMean[TriangularDistribution[a, b, a]]}],
Text[StringJoin["c = b = ", ToString[b]],
{-0.5, Mean[TriangularDistribution[a, b, b]]}],
AbsolutePointSize[4], Red, Table[Point /@
{{-1, HarmonicMean[dist]},
{0, GeometricMean[dist]},
{1, Mean[dist]},
{2, RMS[dist]}},
{c, a, b, d}]},
ImageSize -> 400]];
With[{a=2,b=5,c=3, n=5},
dist= TriangularDistribution[a, b, c];
Plot[{CDF[dist,x],PDF[dist,x]},{x,a,b},
PlotStyle->{{Dashing[{}],Blue},Green},
PlotLabel->" Equi-Probable Intervals for\nTriangularDistribution[" <>
ToString[a] <> ", " <> ToString[b] <> ", " <> ToString[c] <> "]",
Prolog->{Red,AbsoluteDashing[{5,5}],
Line[{{a,y=CDF[dist,#]},{#,y},{#,0}}]& /@
(equiProb[dist,n]//Flatten//Union)},
ImageSize->400]];
Bob Hanlon
In article <baiar6$dks$1 at smc.vnet.net>,
<drmajorbob+MathGroup3528 at mailblocks.com> wrote:
<<
Subject: Re: Re: equiprobable intervals with triangular pdf
From: Bobby Treat <drmajorbob+MathGroup3528 at mailblocks.com>
To: mathgroup at smc.vnet.net
Date: Thu, 22 May 2003 11:04:06 +0000 (UTC)
Bob,
I corrected a couple of errors (dist undefined in equiProb, x arguments
not needed in Mean, Variance, Mode, etc.) and I extended the code to
asymmetric triangular distributions. I learned a lot in the process,
but not enough, I'm thinking!
I even started handling the two special cases, where the mode is at one
of the endpoints. That got really complicated if I defined those to be
TriangularDistribution[a_,a_,c_] or TriangularDistribution[a_,c_,c_],
and it got really tedious (hence unfinished) when I gave those
distributions their own names. I wonder if there's a best way to
handle that? (The tedious way, I'm thinking.)
OVER and OVER I ran into problems getting things to Simplify properly.
For instance CDF applied to the results of equiProb doesn't simplify,
if we leave things symbolic.
The behavior of UnitStep precisely at 0 was often a problem, too. That
made the following line necessary:
TriangularDistribution /: Quantile[TriangularDistribution[a_, b_, c_],
1] := c
Very annoying. I was free to use Which or If in some cases but held to
UnitStep instead -- and that may have caused some of my problems. But,
for the PDF at least, I needed UnitStep.
ClearAll[TriangularDistribution, DescendingTriangularDistribution,
AscendingTriangularDistribution, a, b, c, x, dist];
TriangularDistribution::usage = "TriangularDistribution[a, b, c]
represents
the triangular distribution defined on the interval [a, c] with
peak at \
b.\n TriangularDistribution[a, c] is TriangularDistribution[a, (a+c)/2,
c].";
DescendingTriangularDistribution::usage = \
"DescendingTriangularDistribution[a, b] is equivalent to \
TriangularDistribution[a, a, b]";
AscendingTriangularDistribution::usage =
"AscendingTriangularDistribution[a, \
b] is equivalent to TriangularDistribution[a, b, b]";
TriangularDistribution[args__] /; ! OrderedQ[{args}] := \
TriangularDistribution @@ Sort[{args}]
TriangularDistribution[a_, b_, b_] :=
AscendingTriangularDistribution[a, b]
TriangularDistribution[a_, a_, b_] :=
DescendingTriangularDistribution[a, b]
TriangularDistribution[a_, c_] := TriangularDistribution[a, (a + c)/2,
c]
AscendingTriangularDistribution /:
PDF[AscendingTriangularDistribution[a_,
c_], x_] :=
Evaluate[Simplify[((2*(-a + x))/((a - c)*(a - c)))*(UnitStep[x - a]
-
UnitStep[x - c])]]
DescendingTriangularDistribution /:
PDF[DescendingTriangularDistribution[a_, \
c_], x_] :=
Evaluate[Simplify[(-((2*(c - x))/((a - c)*(-a + \
c))))*(UnitStep[x - a] - UnitStep[x - c])]]
TriangularDistribution /: PDF[TriangularDistribution[a_, b_, c_], x_]
:=
Evaluate[Block[{one, two}, one = (2*(-a + x))/((a - b)*(a - c));
two = -((2*(c - x))/((b - c)*(-a + c)));
Simplify[UnitStep[x - a]*one + UnitStep[x - b]*(two -
one) - UnitStep[x - c]*two]]]
TriangularDistribution /: CDF[TriangularDistribution[a_, b_,
c_], x_] :=
Evaluate[FullSimplify[Integrate[PDF[TriangularDistribution[
a, b, c], x], x], a < b < c]]
TriangularDistribution /: Mean[TriangularDistribution[a_,
b_, c_]] := Simplify[(a + b + c)/3];
TriangularDistribution /: Mode[TriangularDistribution[a_, b_, c_]] := b;
TriangularDistribution /: Variance[TriangularDistribution[a_, b_, c_]]
:= \
Simplify[((a + b + c)^2 - 3*(a*b + b*c + a*c))/18];
TriangularDistribution /: StandardDeviation[TriangularDistribution[a_,
b_, \
c_]] := Sqrt[Simplify[((a + b + c)^2 - 3*(a*b + b*c + a*c))/18]];
TriangularDistribution /: Moment[
TriangularDistribution[a_, b_,
c_], (n_Integer)?NonNegative] :=
Simplify[(2*(a^(2 + n)*b - a*
b^(2 + n) - a^(2 + n)*c + b^(2 + n)*c + (a -
b)*c^(2 + \
n)))/((a - b)*(a - c)*(b - c)*(1 + n)*(2 + n))]
TriangularDistribution /: Quantile[TriangularDistribution[a_,
b_, c_], 1] := c
TriangularDistribution /: Quantile[TriangularDistribution[
a_, b_, c_], q_] := Evaluate[Block[{cutoff = (a - b)/(a - c),
one = a - Sqrt[a - b]*Sqrt[a - c]*Sqrt[q], two =
c - Sqrt[(-(a - c))*(b - c)*(-1 + q)]},
UnitStep[q]*one + UnitStep[q - cutoff]*(two - one) - UnitStep[q -
1]*two]]
TriangularDistribution /: equiProb[TriangularDistribution[a_, b_,
c_], (n_Integer)?
Positive] := Partition[Table[Quantile[TriangularDistribution[a, b,
c], k/
n], {k, 0, n}], 2, 1];
dist = TriangularDistribution[0, 1, 1/7];
Plot[Evaluate[#[dist,
x] & /@ {PDF, CDF, Quantile, x &}], {x, 0, 1}, PlotStyle ->
{Yellow, \
Red, Blue, Black}, AspectRatio -> Automatic, ImageSize -> 150];
Bobby
-----Original Message-----
From: Bob Hanlon <bobhanlon at aol.com>
To: mathgroup at smc.vnet.net
Subject: [mg41549] Re: equiprobable intervals with triangular pdf
It is generally better to define piecewise continuous functions with
UnitStep.
It is also recommended that you follow the conventions used in
Statistics`ContinuousDistributions` for defining distributions so that
use of
your distribution is consistent with use of the standard distributions.