Fit dataset with customized discontinuous distribution

Question

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?

Answer 1

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 *)

Answer 2

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 ...

Answer 3

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.