4
$\begingroup$

I would like to estimate the parameters of a distribution that I obtain via TransformedDistribution:

bDist[σ_, λ_] := TransformedDistribution[real + noise, 
                  {real \[Distributed] ExponentialDistribution[λ], 
                  noise \[Distributed] NormalDistribution[0, σ]}]

(This new object is an exponentially-modified Gaussian for reference.) So, I generate some data using the newly-created distribution:

data = RandomVariate[bDist[2, 1/3], {10000}];

Then try to estimate its parameters, using FindDistributionParameters:

FindDistributionParameters[data, bDist[σ, λ]]

Which produces an error:

"The support of the distribution TransformedDistribution.... could not be determined. 
The validity of the data for TransformedDistribution.... could not be determined"

Where "..." is the full detail of the TransformedDistribution I form above.

It seems to me that these errors suggest that I need to "manually" tell Mathematica what the support of this transformed distribution is. This seems odd - I would have thought that Mathematica could determine that the distribution has support over the entire real line.

I have tried placing assumptions on the parameters (constraining them to be positive), but this doesn't seem to work. I have also tried using EstimatedDistribution, which produces the same errors as above (makes sense since I think they are the same thing under the hood.)

I have also tried using the Laplace distribution which naturally has support over all the real line (opposed to the exponential which is only over the positive reals), but this produces the same error.

There aren't any examples where Mathematica uses derived distributions formed from TransformedDistribution on the help page for the FindDistributionParameters function. Neither is there any information on how one could define support on a function.

Finally, I tried using Maximize on a LogLikelihood object, but this takes far longer (as this question taught me), and only works for small data samples.

Does anyone have an idea as to how I can solve this problem?

$\endgroup$
  • 2
    $\begingroup$ The problem seems to be related to the fact that DistributionDomain[bDist[σ, λ]] (used under the hood IIRC) remains unevaluated, instead of returning the expected Interval[{-∞, ∞}]. I'm not sure if it's a bug or a limitation of DistributionDomain[]. $\endgroup$ – J. M. will be back soon Feb 21 '16 at 3:59
  • $\begingroup$ Thanks for your comment. I wonder if I can get round this issue by censoring a Laplace distribution over the negatives? I will try tomorrow to see if this helps. Thanks again! Best, Ben $\endgroup$ – ben18785 Feb 21 '16 at 4:09
5
$\begingroup$

The problem

Let $Z = X + Y$ where: $X \sim \text{Exponential}(\frac1\lambda)$ and $Y \sim N(0,\sigma^2)$ are independent. We observe $(Z_1, \dots, Z_n)$, and wish to estimate the parameters $(\lambda, \sigma)$.

Comment (too long for the comment box ...)

A neater and simpler starting point is to observe that:

$$E[Z] = E[X] + E[Y] = \frac1\lambda$$

Since $\lambda = \frac{1}{E[Z]}$, a natural choice to estimate $\lambda$ is: $$\hat{\lambda}= \frac{n}{\sum_{i=1}^n Z_i} = \frac{1}{\bar{Z}}$$

... as Bob Hanlon has intimated in his answer (since he starts the estimation search for $\lambda$ here). Note that $\hat{\lambda}$ is the MLE for $\lambda$ (biased but consistent). If you are happy to use the MLE for $\lambda$, then we don't need to start the search here ... we do not need to search for it at all.

Your problem then reduces to estimation of a single parameter, so you can set up the observed loglikelihood purely as function of $\sigma$, which makes life (and the optimisation problem) much simpler.

$\endgroup$
  • $\begingroup$ Many many thanks here. I can't tell you how long I've spent trying to find the ML estimator here. I hadn't seen the woods for the trees it seems! Best, Ben $\endgroup$ – ben18785 Feb 21 '16 at 17:37
3
$\begingroup$
$Version

(*  "10.3.1 for Mac OS X x86 (64-bit) (December 9, 2015)"  *)

bDist[σ_, λ_] := 
  TransformedDistribution[
   real + noise, {real \[Distributed] ExponentialDistribution[λ], 
    noise \[Distributed] NormalDistribution[0, σ]}];

data = RandomVariate[bDist[2, 1/3], {10000}];

Although the FindDistributionParameters::nvsprt warning message is still generated, you can get decent estimates by providing initial estimates for the parameters.

FindDistributionParameters[data, bDist[σ, λ], 
 {{σ, Mean[data]}, {λ, 1/Mean[data]}}]

(*  FindDistributionParameters::nvsprt: The support ... could not be determined. >>

{σ -> 2.00001, λ -> 0.335188}  *)

Alternatively, instead of using the TransformedDistribution you can use ProbabilityDistribution with the PDF given in your cited reference for the exponentially modified Gaussian distribution.

dist[σ_, λ_] = ProbabilityDistribution[
   λ/2 Exp[λ/2 (λ σ^2 - 2 x)]*
    Erfc[(λ σ^2 - x)/(Sqrt[2] σ)],
   {x, -∞, ∞}, 
   Assumptions -> {λ > 0, σ > 0}];

Again, by providing initial estimates for the parameters

FindDistributionParameters[data, dist[σ, λ], 
 {{σ, Mean[data]}, {λ, 1/Mean[data]}}]

(*  {σ -> 2.00001, λ -> 0.335188}  *)

This provides the same estimates but without the warning message.

$\endgroup$
1
$\begingroup$

Just to confirm if a bug was introduced in v10, this is the result (reasonable in spite of detecting the same error) in v9:

SeedRandom[42];
bd[l_, s_] := TransformedDistribution[real + noise, 
               {real \[Distributed] ExponentialDistribution[l], 
               noise \[Distributed] NormalDistribution[0, s]}]
data = RandomVariate[bd[2, 1], {100}];
FindDistributionParameters[data, bd[l, s]]

FindDistributionParameters::nvsprt: The support of the distribution TransformedDistribution[[FormalX]1+[FormalX]2,{[FormalX]1[Distributed]ExponentialDistribution[l],[FormalX]2[Distributed]NormalDistribution[0,s]}] could not be determined. The validity of the data for TransformedDistribution[[FormalX]1+[FormalX]2,{[FormalX]1[Distributed]ExponentialDistribution[l],[FormalX]2[Distributed]NormalDistribution[0,s]}] could not be determined. >>

(* {l -> 2.36135, s -> 1.01611} *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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