# Fit dataset with customized discontinuous distribution

I have the following parametric CDF:

F[x_, λ_, β_, γ_] = Piecewise[{
{0, x < 0},
{1 - (1 - β) E^(-x λ), 0 <= x < γ },
{1, x >= γ}
}];


that is for example

Plot[F[x, 0.5, 0.2, 5], {x, -10, 10}]


I'm trying to give an estimation of these parameters based on my dataset.

First I build a ProbabilityDistribution object from CDF:

myTruncExp2[λ_, β_, γ_] := ProbabilityDistribution[{"CDF",F[x, λ, β, γ]}, {x, -∞, ∞}, Assumptions -> {λ > 0, β > 0, γ > 0}];


My data list is stored in dt1 list, but these two attempts both fail:

FindDistributionParameters[dt1, myTruncExp2[λ, β, γ]]

DistributionFitTest[dt1, myTruncExp2[λ, β, γ], {"PValue", "FittedDistribution"}]


Perhaps Mathematica is unable to manage distribution object with discontinuous CDF, in fact the plot I got using the object is not what I expect

Plot[Evaluate@CDF[myTruncExp2[0.5, 0.2, 5], x], {x, -20, 20}]


Even when I represent the same function through its PDF, that is

f[x_, λ_, β_, γ_] := Piecewise[{
{0, x < 0},
{(1 - β) λ E^(-x λ), 0 <= x < γ },
{β + (1 - β) E^(-γ λ), x == γ},
{0, γ < x}
}];


putting it into the previous object

myTruncExp2[λ_, β_, γ_] := ProbabilityDistribution[f[x, λ, β, γ], {x, -∞, ∞} ,Assumptions -> {λ > 0, β > 0, γ > 0}];


the integration discards the mass point, and it still comes with the same problem

CDF[myTruncExp2[0.5, 0.2, 5], x]


$$\begin{array}{cc} 0.734332 & x\geq 5. \\ 0.8 e^{-0.5 x} \left(e^{0.5 x}-1.\right) & 0.<x<5. \\ 0. & \text{True} \\ \end{array}$$

There would be a way to manage discontinuous custom distribution objects in Mathematica? That is, there would be a way to let Mathematica know to integrate both function and probability mass points?

May you suggest another strategy to get an estimation of these distribution parameters?

And then another strategy to perform a goodness of fit test of this distribution on my dataset?

• I think it would be helpful if you could give the actual data set dt1 (or create another data set that behaves the same). – demm Sep 23 '15 at 13:37
• From the documentation on ProbabilityDistribution: "For a multivariate ProbabilityDistribution definition, all variables need to be either discrete or continuous; no mixed cases can occur." I would imagine that also means that one can't have a continuous part and a probability mass at a single point for univariate distributions. Also, you can normalize your distribution with the Method->"Normalize" option but if you really have the remaining mass at $\gamma$, then that won't be what you want. – JimB Sep 23 '15 at 18:00
• This is the CDF from your points i.stack.imgur.com/Al1UX.png It doesn't resemble your model too much ... – Dr. belisarius Sep 23 '15 at 19:31
• Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Dr. belisarius Sep 23 '15 at 20:33
• Each row of this file is an instance of dt1 list dropbox.com/s/c9wf46x6w0z90vm/err1.txt?dl=1 Anyway I have to make hypothesis test of that distribution whatever the dataset may contain. – Timothy Sep 23 '15 at 23:48

After looking at this a bit more there are closed-form maximum likelihood estimates for the parameters of this doubly-censored shifted exponential distribution.

Suppose there are $n$ observations with $n_0$ having the minimum value, $n_2$ having the maximum value, and that we label the $n_1=n-n_0-n_2$ values between the minimum and the maximum values as $x_1, x_2, …, x_{n_1}$. The likelihood function is

$$L=\beta^{n_0}\left((1-\beta)e^{-\lambda\gamma}\right)^{n_2}\prod_{i=1}^{n_1} (1-\beta)\lambda e^{-\lambda x_i}$$

with the log likelihood being

$$\log(L)=n_0\log\beta+\left(n_1+n_2\right)\log(1-\beta)+n_1\log\lambda-n_2\lambda\gamma-\lambda\sum_{i=1}^{n_1} x_i$$

The maximum likelihood estimators are

\begin{align*} \hat{\beta}&=\frac{n_0}{n}\\ \hat{\gamma}&=\max\limits_i x_i\\ \hat{\lambda}&=\frac{n_1}{\sum\limits_{i=1}^{n_1} x_i+\hat{\gamma}n_2} \end{align*}

The Mathematica code to generate samples from this distribution and estimate the parameters follows:

(* Set some values for the parameters to be used in the example *)
parms = {λ -> 0.5, β -> 0.2, γ -> 5};

(* Take a random sample from a shifted exponential *)
x = RandomVariate[ExponentialDistribution[λ], 500] + Log[1 - β]/λ /. parms;
(* Censor values to be between 0 and γ *)
y = Map[Min[γ, Max[0, #]] &, x] /. parms;

(* Construct some summary statistics to create the maximum likelihood estimates *)
n = Length[y]; (* Total number of observations *)
nMinimum = Total[Map[If[Min[y] == #, 1, 0] &, y]]; (* Number of observations equaling the minimum *)
nMaximum = Total[Map[If[Max[y] == #, 1, 0] &, y]]; (* Number of observations equaling the maximum *)
z = Map[If[Min[y] < # < Max[y], 1, 0] &, y];
sumY = z.y; (* Sum of the obervations between 0 and γ *)
nMiddle = Total[z];  (* Total number of observations between 0 and γ *)

(* Maximum likelihood estimates *)
β0 = N[nMinimum/n] (* Proportion of observations equaling the minimum *)
(* 0.206 *)
γ0 = Max[y] (* Minimum possible value of γ is maximum likelihood estimate *)
(* 5 *)
λ0 = nMiddle/(sumY + γ0 nMaximum) (* Similar to reciprocal of mean of non-zero observations *)
(* 0.5276208413626468 *)
aic = (-2 Sum[Log[f[y[[i]], λ, β, γ]], {i, Length[y]}] /. {λ -> λ0, β -> β0, γ -> γ0}) + 2*3
(* 1727.7477639507404 *)

(* Use FindMaximum as a partial check *)
(* Define probability function for a doubly-censored shifted exponential *)
f[x_, λ_, β_, γ_] := Piecewise[{
{β, x == 0},
{λ (1 - β) Exp[-λ x], 0 < x < γ},
{(1 - β) Exp[-λ γ], x == γ},
{0, x > γ}}]

mle = FindMaximum[{Sum[Log[f[y[[i]], λ, β, γ]], {i, Length[y]}]}, {{λ, λ0}, {β, β0}, {γ, γ0}}]
(* {-860.8738819753705,{λ->0.5276208413626468,β->0.206,γ->5.}} *)
aicFindMaximum = -2 mle[[1]] + 2*3
(* 1727.747763950741 *)


This is an extended comment.

I gathered all your data in l and git this CDF

l1 = Join @@ l;
hl = HistogramList[l1, "Knuth", "CumulativeCount"];
data = Transpose[{MovingAverage[hl[[1]], 2], Rescale@hl[[2]]}];
ListPlot@data


This plot strongly suggest that your model isn't a good fit because:

1) There are negative values {-1 ...1} and your model starts at zero
2) Your model can't account for the positive second derivative near zero
3) There is no a clear "massive" point (discontinuity)


I believe you should rethink your model ...

• Thanks for your reply. I understand your points and I agree with them. However that model has been proposed in the '70s by an emeritus professor to approximate this error distribution, not by me. It may be the case that the dataset hasn't been accurately generated by me, I know. Anyway I have to prove that a given dataset doesn't fit with this model, then I have to "formally" perform a goodness of fit test of this distribution. If possible using Mathematica... – Timothy Sep 23 '15 at 20:31
• @Timothy You can do that with NonLinearModelFit[ ]. But not with this data set ... – Dr. belisarius Sep 23 '15 at 20:58
• @Timothy. I'm going to suggest that you don't want NonLinearModelFit if the proposed model is really a probability distribution and what you have are random samples from that distribution. (If the observations are measurements in which their values match the shape of a probability distribution, then NonLinearModelFit is what you want.) What you'd need to construct is the likelihood function (and not the cumulative distribution function) and use maximum likelihood. That would also get you proper standard errors for the coefficients which least squares with NonLinearModel Fit won't. – JimB Sep 23 '15 at 22:01

2nd Update:

Here's brute force method to obtain maximum likelihood estimates and the AIC statistic which you could use in comparing the fit to other distributions.

What you have is a doubly censored shifted exponential distribution. The mixed density/probability function is given by

(* Define probability function *)
f[x_, λ_, β_, γ_] := Piecewise[{
{1 - Exp[λ Log[1 - β]/λ], x == 0},
{λ Exp[-λ (x - Log[1 - β]/λ)], 0 < x < γ},
{Exp[-λ (γ - Log[1 - β]/λ)], x == γ},
{0, x > γ}}]


One can create a random sample from this distribution:

(* Set some values for the parameters to be used in the example *)
parms = {λ -> 0.8, β -> 0.5, γ -> 5};

(* Take a random sample from a shifted exponential *)
x = RandomVariate[ExponentialDistribution[0.5], 100] +
Log[1 - β]/λ /. parms;

(* Censor values to be between 0 and λ *)
y = Map[Min[γ, Max[0, #]] &, x] /. parms;


And now find the maximum likelihood estimates of the parameters:

(* Find maximum likelihood estimate given the starting values *)
λ0 = 1/Mean[y];  (* Estimate for a non-censored exponential *)
β0 = Total[Map[If[Min[y] == #, 1, 0] &, y]]/Length[y];  (* Proportion of observations equaling the minimum *)
γ0 = Max[y];  (* Start with minimum possible value of γ *)
mle = FindMaximum[{Sum[Log[f[y[[i]], λ, β, γ]], {i, Length[y]}], γ >= γ0}, {{λ, λ0}, {β, β0}, {γ, γ0}}]
aic = -2 mle[[1]] + 2 3


(* {-167.736, {λ -> 0.530281, β -> 0.38, γ -> 5.00107}} ) ( 341.472 *)

End of 2nd update

Update: While I still think the issue is that the ProbabilityDistribution function does not use the information from the Piecewise function as you want, I now see that the distribution intended is a doubly censored shifted exponential. As such it would seem hopeful to use the built-in probability functions to obtain estimates of the parameters from data. (This does not address the goodness-of-fit issue, however.)

So first define a doubly censored exponential and plot it:

(* Define an doubly censored exponential distribution *)
d[λ_, β_, γ_] := CensoredDistribution[{-Log[1 - β]/λ, γ - Log[1 - β]/λ},
ExponentialDistribution[λ]]

(* Plot an example *)
Plot[CDF[d[0.5, 0.2, 5], x - Log[1 - 0.2]/0.5], {x, -10, 10}, PlotRange -> {All, {0, 1}}]


Generate a random sample from that distribution and then shift.

data = RandomVariate[d[0.5, 0.2, 5], 100] - Log[1 - 0.2]/0.5;


Then estimate the parameters:

FindDistributionParameters[data, d[λ, β, γ], {{λ, 0.5}, {β, 0.2}, {γ, 5.}},
ParameterEstimator -> {"MethodOfMoments", MaxIterations -> 500, AccuracyGoal -> 8}]

(* {λ -> 0.401674, β -> 0.329015, γ -> 4.79587} *)


I haven't been successful with the MaximumLikelihood option and I need to check a bit more on the "shifting" adjustment. It does seem that because the distribution function at zero is above zero that there is also a probability mass intended at zero as well as at the upper value at γ hence the double censoring.

End of update

I think the issue is that you don't have a probability density that integrates to 1.

As to why the last plot looks wrong here is what Plot sees:

F[x_, λ_, β_, γ_] :=  Piecewise[{{0, x < 0},
{1 - (1 - β) E^(-x λ),  0 <= x <= γ},
{1.0, x > γ}}];
F[x, 0.5, 0.2, 5]

myTruncExp2[λ_, β_, γ_] := ProbabilityDistribution[{"CDF", F[x, λ, β, γ]}, {x, -∞, ∞},
Assumptions -> {λ > 0, β > 0, γ > 0}]
myTruncExp2[0.5, 0.2, 5]

Simplify[CDF[myTruncExp2[0.5, 0.2, 5], x]]


But this is all because you don't have a probability density function that integrates to 1.

• If you divide 1 - (1 - β) E^(-x λ) by $1 + E^(-\gamma \lambda) (-1 + \beta) - \beta$, then you'll end up with a proper density. – JimB Sep 23 '15 at 14:49
• Thanks for reply. By dividing for this expression, of course I get a CDF that it sums up to 1, but it differently behaves from the one I defined at the beginning of the question: it doesn't star from β and it continues to raise after γ. And anyway it doesn't work with DistributionFitTest. – Timothy Sep 23 '15 at 18:06
• I believe that what you want to define isn't what you get with the Mathematica statements defining F. I assume that you want a probability mass of $(\exp(-\gamma \lambda) (-1 + \beta + \exp(\gamma \lambda) (1 - \beta + \lambda - \gamma \lambda)))/\lambda$ at $\gamma$ and the continuous density between 0 and $\gamma$. – JimB Sep 23 '15 at 18:14