Using EstimatedDistribution with custom circular ProbabilityDistribution

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?

• The choice of the automatic bandwidth selections in SmoothHistogram depends on the sample size. But you have a set of "weighted" data where the sample size is unknown. The sum of the weights is 13777.80435. Maybe explaining the weights would help us better understand how you can end up with a probability distribution.
– JimB
Apr 24, 2018 at 22:22
• Using skd = SmoothKernelDistribution[data, Automatic, {"Bounded", {-180, 180}, "Epanechnikov"}]; Plot[PDF[skd, \[Theta]], {\[Theta], -180, 180}, PlotRange -> {{-180, 180}, {0, 0.006}}] with the Bounded option will more appropriately display the data summary. But I'm still skeptical about the "weights".
– JimB
Apr 24, 2018 at 22:52
• @JimB The plot is just to show that my distribution is okay and otherwise not relevant to the question. I don't get your issue with the weights, but: The angles are measured from a field [theta = f(x,y,z)] defined over a 2-manifold embedded in 3D; the measurements can only be obtained for polygons that tesselate the manifold; the weight is the polygon area and their sum is ~ the total surface area of the manifold. So an observation theta with weight w indicates that a triangle with area w has a measured field value of theta. The distribution should be for an arbitrary point on the manifold.
– nben
Apr 25, 2018 at 13:59
• @JimB More generally, having weights as part of the sample of an underlying distribution is not an unusual concept. See, for example, the wikipedia article on sampling. Mathematica natively extends most of its probability functionality to use weights, and explanations of how to use them are easy to find: en.wikipedia.org/wiki/Sampling_(statistics)#Survey_weights
– nben
Apr 25, 2018 at 14:18
• It is certainly not unusual to have weights when sampling finite populations. However, it is unusual (and at times nonsensical) to have weights when sampling from a continuous probability distribution as you have here. In any event, you should define exactly what these particular weights mean (in the question not just in a comment). For example, using the weights appropriately in a maximum likelihood estimation procedure isn't as simple as adding in an option such as Weights -> w.
– JimB
Apr 25, 2018 at 14:43

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

• Okay, I think I understand the objection with the weights now; thanks!
– nben
Apr 26, 2018 at 0:56