Using EstimatedDistribution with custom circular ProbabilityDistribution

Question

Basic Idea

I'm attempting to estimate the distribution of a set of angle measurements. I have a custom probability distribution that I believe explains the measurements accurately and am attempting to fit the distribution using EstimatedDistribution. However, although I believe I have specified everything correctly, Mathematica seems unable to find a fit at all when the dataset is large (but works fine when the dataset is small).

Sample Data

I've uploaded a sample weighted dataset to tinyupload.com; you can download it here. Assuming you download it to 'dataset.csv', you can import using data = WeightedData@@Transpose@Rest@Import["dataset.csv"]. In this dataset, each observation is an angle in degrees (between -180 and 180).

A simple plot of the sample data looks like this:

SmoothHistogram[data, PlotRange -> {{-180, 180}, {0, Automatic}}]

Probability Distribution

I believe that a good distribution for this data is similar to the sum of two Von Mises distributions; the formula I've settled on is this:

ClearAll[AngleDistribution];
AngleDistribution[t0_, c_] := Block[{t},
  ProbabilityDistribution[
    (Exp[c*Cos[Pi/180 (t - t0)]] + Exp[c*Cos[Pi/180 (t + t0)]])/(720*BesselI[0, c]),
    {t, -180, 180}]];

A plot of this distribution (red) with the data (blue):

With[
  {dist = AngleDistribution[100.0, 2.2]},
  Show[
    {SmoothHistogram[data, PlotRange -> {{-180, 180}, {0, Automatic}}],
     Plot[
       PDF[dist, t], {t, -180, 180},
       PlotRange -> {0, Automatic}, 
       PlotStyle -> {{Red}}]}]]

I've noticed that Mathematica does not seem to be capable of analytically integrating the PDF for this distribution; Integrate[PDF[dist,t], {t,-180,180}] hangs for longer than I'm willing to wait, but NIntegrate[PDF[dist,t], {t,-180,180}] always yields 1. immediately for any reasonable parameterization of dist = AngleDistribution[t0, c].

The Problem

I'm attempting to use EstimatedDistribution to find the parameters for this distribution that best explain the data. In theory, this should be as easy as a call to EstimatedDistribution, perhaps with starting parameters; but this raises a message and returns the distribution without changing the starting paramters:

Block[{t0, c},
 EstimatedDistribution[
   data, 
   AngleDistribution[t0, c],
   {{t0, 105.0}, {c, 2.2}}]]

FindMaximum::fmgz: Encountered a gradient that is effectively zero. The result returned may not be a maximum; it may be a minimum or a saddle point.

ProbabilityDistribution[ 0.000528267 (E^(2.2 Cos[1/180 (-105. + \[FormalX]) \[Pi]]) + E^( 2.2 Cos[1/180 (105. + \[FormalX]) \[Pi]])), {\[FormalX], -180, 180}]

This happens despite the fact that I can observe differences in fit tests for small changes in the parameters:

With[ (* Note: this takes some time to run *)
  {dist1 = AngleDistribution[105.0, 2.2],
   dist2 = AngleDistribution[106.0, 2.2],
   sample = RandomVariate[EmpiricalDistribution[data], 2000]},
  DistributionFitTest[sample, #] & /@ {dist1, dist2}]

{0.0000852185, 0.0000134059}

Another odd observation: in the above code-block, I used a sample of 2000 data-points from the empirical distribution of data; if I use a smaller sample, the fits are consistently much better; for a sample of 100 instead of 2000 points, an example result is {0.327189, 0.298757}; though I assume that this is a feature related to precision and not a bug.

Currently, I'm fairly uncertain as to whether this is an issue related to my specification of the distribution, how I'm calling EstimatedDistribution, or something else entirely.

What I've Tried

My first instinct was to change the ParameterEstimator. This, however, only seems to cause EstimatedDistribution to hang; it's possible that the estimator is working and just taking a very long time, but since I don't see a documented way to monitor steps (like with FindFit and StepMonitor) I can't tell. This seems to be the case with all ParameterEstimators except "MaximumLikelihood" (which has the behavior detailed above).

The next thing I tried was to use a sample of the data instead of the full weighted dataset:

estimateWithSample[data_, n_] := Block[{t0,ct},
  With[
    {sample = RandomVariate[EmpiricalDistribution[data], n]},
    EstimatedDistribution[
      sample,
      AngleDistribution[t0, ct],
      {{t0, 100.0}, {ct, 2.2}}]]]

Confusingly, this seems to yield an identical result if the sample is large enough (> 1500 or so) and seems to work as I would expect if the sample is small:

estimateWithSample[data, 2000]

FindMaximum::fmgz: Encountered a gradient that is effectively zero. The result returned may not be a maximum; it may be a minimum or a saddle point.

ProbabilityDistribution[ 0.000528267 (E^(2.2 Cos[1/180 (-105. + \[FormalX]) \[Pi]]) + E^( 2.2 Cos[1/180 (105. + \[FormalX]) \[Pi]])), {\[FormalX], -180, 180}]

estimateWithSample[data, 100]

ProbabilityDistribution[ 0.000784007 (E^(-1.61967 Cos[1/180 (-91.1175 + \[FormalX]) \[Pi]]) + E^(-1.61967 Cos[ 1/180 (91.1175 + \[FormalX]) \[Pi]])), {\[FormalX], -180, 180}]

Unfortunately, with such small samples, the distribution parameters are not great fits for the entire dataset. I suppose one solution is to bootstrap over many fits to estimate the parameters, but this seems like a clunky work around.

Does anyone know why EstimatedDistribution is having such a hard time with my dataset and distribution? Or how to get a good fit for a dataset like the one I've uploaded?

Answer 1

If the distribution is a mixture of two von Mises distributions, then the maximum likelihood estimates of the associated parameters can be found.

(* Define a mixture of two von Mises distributions *)
mix = ProbabilityDistribution[w PDF[VonMisesDistribution[μ1 π/180, k1], x][[1, 1, 1]] + 
  (1 - w) PDF[VonMisesDistribution[μ2 π/180, k2], x][[1, 1, 1]], {x, -∞, ∞}];

(* Construct the log likelihood *)
logL = LogLikelihood[mix, data[[2, 1]] π/180];

(* Find the maximum likelihood estimates of the 5 parameters *)
sol = FindMaximum[{logL, k1 > 0 && k2 > 0}, {{w, 0.54}, {μ1, -98}, {μ2, 98},
  {k1, 1.56}, {k2, 2.1}}]
(* {-39119.944639131994, {w -> 0.5458155533472334, μ1 -> -1.709486155119687,
   μ2 -> 1.728812322769151, k1 -> 1.5671613366157455,k2 -> 2.0834567826989288}} *)

(* Nonparametric density estimate *)
skd = SmoothKernelDistribution[data, Automatic, {"Bounded", {-180, 180}, "Epanechnikov"}];

(* Show results *)
Show[Histogram[data[[2, 1]], "Scott", "PDF"],
 Plot[{PDF[skd, x], PDF[mix /. sol[[2]], x π/180] π/180}, {x, -180, 180}, 
  PlotRange -> {Automatic, {0, Automatic}},
  PlotLegends -> {"Nonparametric density estimate", 
    "Mixture of two von Mises distributions"}]]

The mixture does not look very good compared to the histogram and the nonparametric density estimate. Given the large amount of data that you have, you shouldn't restrict yourself to a poor-fitting parametric mixture.

I have ignored the "weights" as I don't see how they are relevant unless they are associated with the precision of an individual measurement. In that case you have the observed random variable $Y$ but it is the unobserved and desired distribution of the random variable $X$ contaminated with a varying error (which might be an estimated error or associated standard error - it's not clear) that is of interest.

If the weights somehow modify what distribution that you're sampling from, then simply using the option Weights -> w that some of the Mathematica procedures have won't account properly for such deviations from the hypothesized probability distribution. (Or at least you need to define what the weights really mean and specifically how they affect the random selection from the hypothesized probability distribution.)