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.
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[[2]]]
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}]
NonlinearModelFitis 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. $\endgroup$