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 ParameterEstimator
s 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?
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. $\endgroup$skd = SmoothKernelDistribution[data, Automatic, {"Bounded", {-180, 180}, "Epanechnikov"}]; Plot[PDF[skd, \[Theta]], {\[Theta], -180, 180}, PlotRange -> {{-180, 180}, {0, 0.006}}]
with theBounded
option will more appropriately display the data summary. But I'm still skeptical about the "weights". $\endgroup$Weights -> w
. $\endgroup$