Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to use Mathematica for more of my reliability engineering work and am having trouble using SplicedDistribution.

I have a Weibull distribution (dist1) associated with a low stress condition, β1 = 2.59 and η1 = 25889.

I have a different Weibull distribution (dist2) associated with a high stress condition with the characteristic life being 1/3 that for dist1, β2 = 2.59 and η2 = 8630.

The need is to create a spliced distribution such that for the first 5000 hours only dist1 applies and thereafter only dist2 applies. Think of a device that operates in a relatively clean, stress-free environment for the first 5000 hours of its life, then it is abruptly placed in a stressful, corrosive atmosphere thereafter.

My attempt in using SplicedDistribution does not seem to work as shown in the comparison plot below. SplicedDistribution appears to be “mixing” the two distributions for all time, t as opposed to letting only dist1 govern for the first 5000 hours then allowing only dist2 to govern thereafter. I would expect that dist3 would be equivalent to dist1 for the first 5000 hours then its survival (reliability) curve would depart from dist1.

Can someone show me how to accomplish what I need using Mathematica ?

enter image description here

Here is my version 9 code:

{β1 = 2.59, η1 = 25889.};
dist1 = WeibullDistribution[β1, η1]
plot1 = Plot[SurvivalFunction[dist1][t], {t, 0, 18000}, 
  PlotStyle -> Black, Frame -> True, 
  FrameLabel -> {"Time (hours)", "Reliabilty,  R[t]"}, 
  PlotRange -> {0, 1.05}, 
  PlotLabel -> "Distribution 1 Reliability Plot"]

{β2 = 2.59, η2 = 8630.};
dist2 = WeibullDistribution[β2, η2]
plot2 = Plot[SurvivalFunction[dist2][t], {t, 0, 18000}, 
  PlotStyle -> Blue, Frame -> True, 
  FrameLabel -> {"Time (hours)", "Reliabilty,  R[t]"}, 
  PlotRange -> {0, 1.05}, 
  PlotLabel -> "Distribution 2 Reliability Plot"]

plot3 = Show[plot1, plot2, PlotRange -> All, 
  PlotLabel -> "Distribution 1 & 2 Comparison Plot"]

dist3 = SplicedDistribution[{1, 1}, {0, 5000, ∞}, {dist1, dist2}]
plot4 = Plot[SurvivalFunction[dist3][t], {t, 0, 18000}, 
  PlotStyle -> Red, Frame -> True, 
  FrameLabel -> {"Time (hours)", "Reliabilty,  R[t]"}, 
  PlotLabel -> "Distribution 3 (Spliced) Reliability Plot"]

plot5 = Show[plot3, plot4, PlotRange -> All, AxesOrigin -> {0, 0}, 
  PlotLabel -> "Comparison Plot", 
  Epilog -> {Text["Distribution 1 (low stress)", {15000, .9}],
    Text["Distribution 2 (high stress)", {10500, .6}],
    Text["Distribution 3 (spliced)", {6000, .2}]}]
share|improve this question

1 Answer 1

(I don't have v9 so I'll be using Truncated/Mixture Distributions)

One way to do it would be to truncate dist1 at 5000 and shift dist2 by 5000, then mixing them together giving weight CDF[dist1, 5000] to truncated dist1.

{β1 = 2.59, η1 = 25889.};
dist1 = WeibullDistribution[β1, η1];
{β2 = 2.59, η2 = 8630.};
dist2 = WeibullDistribution[β2, η2]

trunc1 = TruncatedDistribution[{0, 5000}, dist1];
shift2 = TransformedDistribution[x + 5000, x \[Distributed] dist2];

dist3shift = MixtureDistribution[
   {CDF[dist1, 5000], 1 - CDF[dist1, 5000]},
   {trunc1, shift2}];

plotshift = Plot[{
  SurvivalFunction[dist1, t],
  SurvivalFunction[dist2, t],
  SurvivalFunction[dist3shift, t]
  }, {t, 0, 20000}]

plot

This assumes that the surviving part that is placed in the bad environment hasn't been damadged, this might not be the case.

Show[plotshift,
 Plot[
  SurvivalFunction[dist1, 5000] SurvivalFunction[dist2, t - 5000],
  {t, 0, 20000},
  PlotStyle -> Directive[Black, Dashed]]]

plot

Closeup around 5000:

closeup

So it might be prudent to start dist2 at some corresponding damage level. Here this is done by finding where the survival rate of dist2 corresponds to the survival rate dist1 has at 5000. A new distribution is created that acts like dist1 up to 5000, where it jumps to the corresponding point in dist2:

betterShift = t /. FindRoot[
  SurvivalFunction[dist2, t] == SurvivalFunction[dist1, 5000],
  {t, 5000}];

trunc2 = TruncatedDistribution[{betterShift, Infinity}, dist2];
shift2Better =
  TransformedDistribution[
    x + 5000 - betterShift,
    x \[Distributed] trunc2];

dist3truncShift = 
  MixtureDistribution[
    {CDF[dist1, 5000], 1 - CDF[dist1, 5000]},
    {trunc1, shift2Better}];

Plot[{
  SurvivalFunction[dist1, t],
  SurvivalFunction[dist2, t],
  SurvivalFunction[dist3truncShift, t]}
 , {t, 0, 20000},
 GridLines -> {{5000}, None}]

better shift

Another closeup:

closeup

share|improve this answer
    
Very nice indeed and many thanks! It will take me some time to fully explore your solution but I think this is what I need. You are right in assuming that the part will be degraded when dist2 kicks in. It will be degraded by an amount caused by 5000 hours of being exposed to dist1. One thing I will try to confirm is whether the shifted/transformed distribution truly is a valid statistical distribution, that is, it's PDF integrates to unity from t zero to infinity. As a side note I do hope Wolfram Research continues adding Reliability & Survival Analysis capability to future versions. –  Steve Dec 11 '13 at 17:39

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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