2
$\begingroup$

I have a collection of sample particles, with (x,y,z) coordinates generated by a simplified Monte Carlo-like code. I expect that these particles will follow an anisotropic diffusion process, which would give me a 3D Gaussian distribution with different variances along different axes (based on the geometry of the problem, I do expect the Gaussian to be axis-aligned; I suppose future development may generate a rotated 3D Gaussian, but I can deal with that generalization when and if it arises).

What I'm looking for is a way to:

  1. Check that my data is consistent with a 3D Gaussian distribution (to confirm that diffusion is a good model);
  2. Get the best-fit mean and variance (along each axis), with error if possible (to measure the drift velocity and the diffusion coefficients).

I thought about binning the data and doing a fit, but I'm not sure that's a great method. I was told to look into the Kolmogorov-Smirnov test, but that appears to be applicable to 1D distributions not 3D. Is there an accepted method for this sort of thing?

I generally work in C++, Fortran, or Python, so I'm hoping to avoid answers along the lines of "Matlab has a function to do that". If Python has a method, I can work with that. Otherwise, I don't mind coding up my own tool, but I don't know what algorithm to use.

Improve this question
$\endgroup$
2
  • $\begingroup$ This sounds like a question for the statistics StackExchange. $\endgroup$
    – user3883
    May 13 '15 at 16:57
  • $\begingroup$ You can consider creating histograms that you feed into a kernel density estimator. $\endgroup$
    – AlexE
    Jul 31 '15 at 14:00
2
$\begingroup$

Your first test should be to compute the mean value and covariance matrix of your point sample. If these converge to the correct values as you increase the number of samples, you are, from a practical perspective on the safe side that your points indeed come from the correct distribution.

Of course, in practice there are many distributions that have the same mean and covariance matrix. If you are still unsure, you could also consider higher order moments of your point set, though I think that in practice mean and covariance are probably sufficient discriminating.

Improve this answer
$\endgroup$
4
  • $\begingroup$ I started by doing that, which tells me the mean and (co)variance, which I can translate to drift velocity and diffusion coefficients. But that doesn't tell me whether or not I have data consistent with a Gaussian, so it gives me no way to judge whether or not the diffusion approximation is correct. I was hoping to find some method of testing my diffusion assumption. $\endgroup$
    – Brendan
    May 11 '15 at 19:36
  • 1
    $\begingroup$ I think you have to ask what it means for a set of samples to be consistent with a give distribution. In essence, this question is equivalent to asking whether the moments of your point set converge to the moments of the continuous distribution. The mean and covariance matrix are simply the first two moments. In the end, your points are consistent with many distributions as long as you only look at single set of points; the real test comes when you observe convergence as you increase the number of points. $\endgroup$
    – Wolfgang Bangerth
    May 11 '15 at 20:31
  • $\begingroup$ I think I see what you're saying: Instead of doing some sort of fitting-like procedure, I should compute a series of moments and compare those to my expected distribution. If enough (for some definition of "enough") of those moments agree to within some tolerance, I can declare that the distributions are consistent. Because while the first and even second moments may be the same for different distributions, an infinite number of moments is equivalent to the distribution, and therefore a long enough list of moments will be discriminating enough. Am I understanding that correctly? $\endgroup$
    – Brendan
    May 11 '15 at 20:39
  • $\begingroup$ Yes, correct. Two distributions with all same moments are identical. $\endgroup$
    – Wolfgang Bangerth
    May 12 '15 at 11:00
3
$\begingroup$

In statistics the widely used test for checking if the distribution is gaussian is the Jarque-Bera test. Koizumi [1] presents an equivalent test for the multivariate case. I don't know if there something ready to use in C++, Fortran, or Python. However, at first sight, the test looks easy to implement.

[1] - Koizumi, Kazuyuki, Naoya Okamoto, and Takashi Seo. "On Jarque-Bera tests for assessing multivariate normality." Journal of Statistics: Advances in Theory and Applications 1.2 (2009): 207-220. (pdf link)

Improve this answer
$\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.