0
$\begingroup$

I would be tremendously grateful for help with this one:

I wrote a function that, for a given frequency distribution (observed data) that i call fbar, and a probability distribution function of the Laplace family with parameters mu and sigma, minimizes the Kullback-Leibler divergence.

Then I use Minimize[] to find the KBD minimizing parameters mu and sigma. When I did this for a very small (fake) dataset and a simple exponential function, the minimization gives a numerical value.

However, for the more complicated Laplace distribution and the larger dataset (200 entries), it only gives me {mu, sigma} as output, which is rather useless.

I attach the code, I am pretty sure, I did something wrong.

Laplace PDF:

laplace[x_, mu_, sigma_] := 1/(2*sigma) E^(-Abs[x - mu]/sigma)

KB Divergence:

KLlaplace[x_, mu_, sigma_] :=  Sum[If[fbar[[i]] == 0, 0,
     laplace[x[[i]], mu, sigma] * Log[laplace[x[[i]], mu, sigma] fbar[[i]]]], {i, 1, 200}]

Minimization:

Minimize[KLlaplace[x, mu, sigma], {mu, sigma}]

Edit: x is just an interval between -.99 and 1 with 200 bins, so Laplace[ ] calculates the PDF for the corresponding bins. I am afraid, I would not know how to attach the fbar dataset. It is a list of 200 values, corresponding frequencies to 200 bins of the same domain.

$\endgroup$
  • $\begingroup$ Are you aware that LaplaceDistribution[] is built-in? $\endgroup$ – J. M. will be back soon Mar 10 '18 at 8:58
  • $\begingroup$ What is x? Does it have a fixed value or is it meant as further optimization variable? $\endgroup$ – Henrik Schumacher Mar 10 '18 at 10:00
  • $\begingroup$ I tried to answer the questions in an edit. I did not know about the built-in Laplace function, will try to implement it. $\endgroup$ – Patrick M. Mar 11 '18 at 17:27
1
$\begingroup$

Following @J.M. 's suggestion about the LaplaceDistribution[] and modifying your definition of the Kullback-Leibler divergence...

(* Generate data from a Laplace distribution *)
SeedRandom[12345];
n = 1000;
y = RandomVariate[LaplaceDistribution[0, 0.16], n];

(* Put into bins *)
x = Table[i/100, {i, -99, 100}];
(* Proportion of observations in each bin *)
fbar = BinCounts[y, {-99/100, 1, 1/100}]/n;

(* Function that gives Kullback-Leibler divergence of Laplace distribution
 with binned version of observed distribution *)
kl[mu_, sigma_] := 
 Sum[If[fbar[[i]] == 0, 0, 
   fbar[[i]] (Log[fbar[[i]]] - 
      Log[CDF[LaplaceDistribution[mu, sigma], x[[i + 1]]] - 
        CDF[LaplaceDistribution[mu, sigma], x[[i]]]])], 
   {i, Length[fbar]}]

(* Find mu and sigma that minimize Kullback-Leibler divergence statistic *)
FindMinimum[kl[mu, sigma], {{mu, 0}, {sigma, 0.15}}]
(* {0.10119387738464815,{mu -> 0.0018434507409612973, sigma -> 0.1505939050029673}} *)

(* Maximum likelihood estimation on original data *)
FindDistributionParameters[y, LaplaceDistribution[mu, sigma]]
(* {mu -> 0.0013829360457713437, sigma -> 0.15054298262072455} *)

Update: Minimizing the Kullback-Leibler divergence between the observed data and a specific distributional form is the same as performing maximum likelihood (Choi, date unknown). So if all you have is the binned data (rather than the raw data), the maximum likelihood estimates can be obtained a little bit more directly with the following:

FindMaximum[
 Sum[fbar[[i]] Log[CDF[LaplaceDistribution[mu, sigma], x[[i + 1]]] -
     CDF[LaplaceDistribution[mu, sigma], x[[i]]]],
  {i, Length[fbar]}], {{mu, 0}, {sigma, 0.15}}]
(* {-4.405331024529697, {mu -> 0.0018434507409612552, sigma -> 0.15059390500296513}} *)
$\endgroup$
  • $\begingroup$ Thank you, this runs. I will take some time to get into the dynamics of it, since the duality of maximum likelihood and minimum KBD does not reproduce in the results yet, but this is excellent. Thanks again. $\endgroup$ – Patrick M. Mar 12 '18 at 17:08
  • $\begingroup$ ??? But it does reproduce the results exactly (both in theory and in practice). $\endgroup$ – JimB Mar 12 '18 at 17:26

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.