# How to fit off axis 2D Lorentzian Distribution

I have a 2D dataset. I can easily fit an on axis Lorentzian by assuming the independence of the x and y coordinates. In the case where no such assumption can be made, what method should I use to fit 2D Lorentzian distributed data?

Note, I am restricted to using mathematica in this project. I am wondering what sorts of tools are available in mathematica to perform this task.

• I am asking how to do it with mathematica. Perhaps I should have stated that. Apr 8, 2017 at 21:00
• What bivariate probability density function are you considering? (And not being a physicist I don't know what "off axis" and "on axis" means.) ${1 \over {2\pi}} [{\gamma \over {((x-x_0)^2+(y-y_0)^2+\gamma^2)^{1.5}}}]$ ?
– JimB
Apr 8, 2017 at 21:03
• @Kuba I think that link deals with regression-type issues rather than fitting probability distributions (assuming that the question above is actually dealing with estimating parameters from a bivariate probability distribution). But I agree: supply some data would go a long way to help define the need.
– JimB
Apr 8, 2017 at 21:08
• @JimBaldwin Just in the simpliest way--not in line with the xy-axis. Yes that distribution function is appropriate. I've tried using it, but using it with nonlinear fit appears to end up in a local minimum. Apr 8, 2017 at 21:12
• NonlinearModelFit is for regression problems not for fitting distributions. (Unless the desired relationship in the regression problem has the same shape as a probability density function - but even then the available error structures currently in Mathematica are not usually appropriate.) I'll enter an answer for your consideration.
– JimB
Apr 8, 2017 at 21:30

Suppose your two variables have marginal Cauchy (if you're a statistician) or a Lorentzian (if you're a physicist) distributions but are "correlated" (which is a loose term for "not independent"). Consider a CopulaDistribution which joins two marginal distributions and imposes some form of non-independence:

d = CopulaDistribution[{"FGM", α},
{CauchyDistribution[a1, b1], CauchyDistribution[a2, b2]}]


Take a random sample from that distribution:

parms = {a1 -> 1, a2 -> 3, b1 -> 1, b2 -> 4, α -> 1/2};
sample = RandomVariate[d /. parms, 100];


Find maximum likelihood estimators of the 5 parameters:

FindMaximum[{LogLikelihood[d, sample],
b1 > 0 && b2 > 0 && 0 <= α <= 1}, {{a1, 1}, {b1, 1}, {a2, 3}, {b2, 4}, {α, 1/2}}]
(* {-642.961757205368, {a1 -> 1.167266125882441, b1 -> 1.1267968189006063,
a2 -> 3.4697597908386792, b2 -> 3.4687845216307216, α -> .4876710819301035}} *)


Update

It is important to have good starting values but unlike the above example, one doesn't know the true values and must obtain starting values from the data. For this maximum likelihood estimation one can use the sample medians for a1 and a2 and the semi-interquartile range for b1 and b2. I really have no good idea in finding a good starting value for $\alpha$. It could start at 0 (implying idependence) or below I've used a completely arbitrary and little tested starting value for $\alpha$: the square root of the Spearman's rho statistic. (Your mileage may vary.)

αInit = SpearmanRho[sample[[All, 1]], sample[[All, 2]]]
αInit = Sign[αInit] Abs[αInit]^0.5
sol = FindMaximum[{LogLikelihood[d, sample],
b1 > 0 && b2 > 0 && -1 <= α <= 1},
{{a1, Median[sample[[All, 1]]]}, {b1, InterquartileRange[sample[[All, 1]]]/2},
{a2, Median[sample[[All, 2]]]}, {b2, InterquartileRange[sample[[All, 2]]]/2},
{α, αInit}}]


Another important consideration is to obtain estimates of precision. Here are the estimates of the asymptotic variances and covariances of the maximum likelihood estimators:

cov = -Inverse[(D[LogLikelihood[d, sample], {{a1, b1, a2, b2, α}, 2}]) /. sol[]]


The correlation matrix can be found in the following manner:

cor = Table[cov[[i, j]]/(cov[[i, i]]^0.5 cov[[j, j]]^0.5), {i, 5}, {j, 5}]