Maxwell fit in 2 dimensions

Question

The histogram below shows the 2D-velocities v of a micrometer sized particle in a plasma measured from image to image (60 Hz) over 33,33 sec.

I would like to fit this data by the 2D Maxwell-Boltzmann formula, to determine the temperature T:

p(v) = m/(k*T)*v*Exp(-m*v^2/(2*k*T)); 

with:

m = 5.99736*10^-13; (*particle mass in kg*)
k = 1.3806488*10^-23; (*Boltzmann constant in J/K*)

To solve that in Mathematica I used:

binSize = 0.000025; (*m/sec*)
hList = HistogramList[v, {binSize}];

Here are the 2d velocity data: http://pastebin.com/rtteUpk4

The result of this is:

vHist=hList[[1]]={0., 0.000025, 0.00005, 0.000075, 0.0001, 0.000125, 0.00015, 0.000175, 0.0002, 0.000225, 0.00025, 0.000275, 0.0003, 0.000325, 0.00035, 0.000375, 0.0004, 0.000425, 0.00045, 0.000475, 0.0005, 0.000525, 0.00055, 0.000575, 0.0006, 0.000625, 0.00065, 0.000675, 0.0007, 0.000725, 0.00075, 0.000775, 0.0008, 0.000825, 0.00085}

and

numHist=hList[[2]]={5, 16, 18, 36, 54, 59, 81, 107, 120, 110, 134, 135, 111, 97, 111, 115, 99, 95, 91, 81, 58, 58, 48, 41, 44, 24, 16, 18, 6, 6, 2, 1, 1, 1}

What do I have to correctly use as data (since vHist has 1 more element than numHist)? If they would have the same dimension I would try:

data=Transpose[{vHist,numHist}]

Now I thought to do the following as a test to plot the fit:

lenhthvHist = Length[vHist];
vHistHelp=vHist[[1 ;; lengthvHist - 1]];
data = Transpose[{vHistHelp, numHist}];
model = m/(k*T)*vHistHelp*Exp(-m*vHistHelp^2/(2*k*T));
fit = FindFit[data, model, {T}, vHistHelp];
modelf = Function[{t}, Evaluate[model /. fit]];
Plot[modelf[t], {t, 0, 0.0009}, Epilog -> Map[Point, data]]

Unfortunately this does not work.

How do I obtain the fitting curve?

Answer 1

numHist = numHist/Tr@numHist/vHist[[2]] // N;
data = Transpose[{MovingAverage[vHist, 2], numHist}];
m = 5.99736*10^-13;(*particle mass in kg*)
k = 1.3806488*10^-23;(*Boltzmann constant in J/K*)
model = m/(k*T)*v*Exp[-m*v^2/(2*k*T)];
fit = FindFit[data, model, {{T, 10}}, v]
modelf = Function[{v}, Evaluate[model /. fit]];
Plot[modelf[t], {t, 0, 0.0009}, Epilog -> Map[Point, data]]

(* {T -> 3315.85} *)

Answer 2

Because you are implying that you have a random sample of values of $v$ from a specified probability distribution, you should consider using maximum likelihood rather than least squares. (In fact, I must comment that in this forum least squares seems to be often used for the estimation of parameters given a random sample which I would argue is almost never the way to do it. At best it is inefficient and maybe the only use is to obtain starting values for an iterative maximum likelihood procedure but that's about it. And there are other methods superior to least-squares.)

In this case there is an explicit maximum likelihood estimator for $T$. If the observed velocities are store in a list named v, then the log of the likelihood function is

logL = n Log[m/(k T)] + Total@Log@v - (m v.v)/(2 k T)

The maximum likelihood estimator of T is That ($\hat{T}$)

That = m v.v/(2 k n)

found with

That = T/.Solve[D[logL,T]==0,T][[1]]

The estimate of the asymptotic standard error is

sigma=FullSimplify[1/((-D[logL,{T,2}]/.T->That))^(1/2)]

or

sigma=1/(2 Sqrt[(k^2 n^3)/(m^2 (v.v)^2)])

Update: Given the good suggestions by @rcollyer and @J.M. here is an expanded answer.

This distribution is named a Rayleigh distribution which is typically parameterized as

$$p(v)={{v}\over{\sigma^2}} e^{-v^2/(2 \sigma^2)}$$

where for this question $T=\sigma^2 m/k$. We define some parameters so that a random sample of size $n$ can be taken:

m = 5.99736*10^-13;(*particle mass in kg*)
k = 1.3806488*10^-23;(*Boltzmann constant in J/K*)
T = 3316;
n = 1999;

Determine the parameter required by Mathematica's RayleighDistribution and take a random sample:

σ = (k T/m)^(1/2)
(* 0.000276292 *)
v = RandomVariate[RayleighDistribution[σ], n];

The maximum likelihood estimate of $T$ is found with

sol = FindDistributionParameters[v, RayleighDistribution[s]]
{s -> 0.000274438}
σhat = s /. sol
(* 0.000274438 *)
That = σhat^2  m/k
(* 3271.63 *)

The direct formula for $\hat T$ is

That = m v.v/(2 k n)
(* 3271.63 *)

One shortcoming of FindDistributionParameters is that no estimates of standard errors seem to be available. For this maximum likelihood estimator there is a closed-form for the estimate of the asymptotic standard error:

m v.v/(2 k n^(3/2))
(* 73.1741 *)

One can bin the values to create a histogram of the data but the estimation should be performed with the raw data. And with this much data using SmoothKernelDistribution should be used instead of Histogram. Using the bins and counts from above we can approximate the maximum likelihood estimate:

counts = {5, 16, 18, 36, 54, 59, 81, 107, 120, 110, 134, 135, 111, 97,
    111, 115, 99, 95, 91, 81, 58, 58, 48, 41, 44, 24, 16, 18, 6, 6, 2,
    1, 1, 1};

vHist = {0., 0.000025, 0.00005, 0.000075, 0.0001, 0.000125, 0.00015, 
   0.000175, 0.0002, 0.000225, 0.00025, 0.000275, 0.0003, 0.000325, 
   0.00035, 0.000375, 0.0004, 0.000425, 0.00045, 0.000475, 0.0005, 
   0.000525, 0.00055, 0.000575, 0.0006, 0.000625, 0.00065, 0.000675, 
   0.0007, 0.000725, 0.00075, 0.000775, 0.0008, 0.000825, 0.00085};

midPoints = Table[(vHist[[i]] + vHist[[i + 1]])/2, {i, Length[counts]}];

That = m counts.(midPoints^2) /(2 k Total[counts])
(* 3043.82 *)
sigma = m counts.(midPoints^2)/(2 k Total[counts]^(3/2))
(* 68.0789 *)

A plot of the fit is shown below:

  Show[ListPlot[Flatten[Table[{{vHist[[i]], 0}, {vHist[[i]], 
      counts[[i]]/0.049975}, {vHist[[i + 1]], 
      counts[[i]]/0.049975}, {vHist[[i + 1]], 0}}, {i, 
     Length[counts]}], 1], Joined -> True, Frame -> True, 
  PlotRange -> {Automatic, {0, 3000}}, 
  FrameLabel -> {{Style["Probability density", Large, Bold, Black], 
     ""}, {"", ""}}],
 Plot[(m/(k That)) vv Exp[-m vv^2/(2 k That)], {vv, 0, 0.001}, 
  PlotStyle -> Red]]

If only binned data is available, a better approach would be to used the censored likelihood approach. And after all of this goodness-of-fit should be examined but I've run out of time for now.

Second update

The raw data has been made available and below is the fit using that data. The 1,999 data points have been assigned to the variable named v.

(* Set constants and Rayleigh distribution parameter *)
m = 5.99736*10^-13;(*particle mass in kg*)
k = QuantityMagnitude[
  UnitConvert[
   Quantity["BoltzmannConstant"]]];(*Boltzmann constant in J/K*)
n = Length[v];
sigma = (k T/m)^(1/2);

(* Find maximum likelihood estimate of T *)
sol = FindDistributionParameters[v, RayleighDistribution[s]];
sigmahat = s /. sol
(* 0.00026481579806363645 *)
That = sigmahat^2 m/k
(* 3046.2439475575343 *)
(* As a check using the closed-form maximum likelihood solution *)
That = m v.v/(2 k n)
(* 3046.2439475575347 *)

(* Standard error of That *)
seT = m v.v/(2 k n^(3/2))
(* 68.13312083776144` *)

(* Approximate 95% confidence interval for T *)
{That - 1.96 seT, That + 1.96 seT}
(* {2912.7030307155223`,3179.784864399547`} *)

A plot of the results:

d = SmoothKernelDistribution[v];
Show[Histogram[v, Automatic, "PDF"],
 Plot[PDF[d, v], {v, 0, 0.00085}],
 Plot[PDF[RayleighDistribution[sigmahat], v], {v, 0, 0.00085}, PlotStyle -> Red], ImageSize -> Large]
DistributionFitTest[v, RayleighDistribution[sigmahat]]

With the SmoothKernelDistribution in gray/blue and the maximum likelihood fit in red:

It certainly does not look like a perfect fit. Using DistributionFitTest results in a very, very small P-value suggesting a lack of fit:

DistributionFitTest[v, RayleighDistribution[sigmahat]]
(* 8.790774019828973`*^-8 *)

Part of the problem using P-values is that small and practically inconsequential differences can become "statistically significant" with large sample sizes. So an examination of the data is order. Plotting the raw data in sequence results in a very non-random arrangement of values in time order. This lack of expected randomness is critical because the rationale for being able to estimate $T$ relies on random and independent samples from the Rayleigh distribution.

Here we see high autocorrelation

tsm = TimeSeriesModelFit[v]["BestFit"]
(* ARProcess[0.00003833580283247393`,0.5734763681300956`,0.4570403165209849`,
0.08967306679823858`,-0.018997998762076446`,-0.12439400350041019`,
-0.08876573889266684`},3.3484025317002574`*^-9] *)

and maybe some periodicity as evidenced by showing the first 200 data points:

ListLinePlot[Table[{i, v[[i]]}, {i, 1, 200}], Frame -> True]

The periodicity is even more evident in the complete series:

It might be best to address this issue at https://stats.stackexchange.com/. One possibility is to only use every 25th sample point as distasteful as that sounds. (But the assumptions should still be examined even after doing that.)