7
$\begingroup$

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*)

enter image description here

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?

$\endgroup$
  • 2
    $\begingroup$ I must leave in 5 min, but a few comments: 1. Wouldn't it be easier to determine the temperature from the mean velocity? 2. {bins, counts} = HistogramList[...] is much easier to write. 3. Element -1 of a list is the same as element Length[list] - 1. 4. I didn't have time to understand your model fully, but it looks weird. It should be a simple formula, why do you mix in the data? 5. Get the bin centres with MovingAverage[bins,2]. Then you have the same number of values as in counts. $\endgroup$ – Szabolcs Nov 5 '15 at 16:04
  • 2
    $\begingroup$ 6. It seems you are trying to fit a probability density to unnormalized particle counts. That doesn't work. Finally: since I wrote this in a hurry, there's probably something I misunderstood, for which I apologize in advance. Hope some of the comments helped. $\endgroup$ – Szabolcs Nov 5 '15 at 16:04
7
$\begingroup$
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} *)

Mathematica graphics

$\endgroup$
  • $\begingroup$ Probably you want your velocities in m/s instead of in mm/s ... so to get your T result in reasonable units ... $\endgroup$ – Dr. belisarius Nov 5 '15 at 16:47
  • $\begingroup$ Yes, the unit in the plot is m/sec. I put the wrong label and will correct it ... thanks ... and thanks a lot for the solution ... Is it possible to use a second fit parameter (amplitude) to get the correct count values? What happens exactly in the first line: numHist = numHist/Tr@numHist/vHist[[2]] // N; And what is the meaning of {{T,10}}? $\endgroup$ – mrz Nov 5 '15 at 18:11
12
$\begingroup$

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]]

Histogram and fit

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:

enter image description here

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]

Periodicity

The periodicity is even more evident in the complete series:

Periodicity in 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.)

$\endgroup$
  • 1
    $\begingroup$ I have some suggestions on your code. Log is listable, so use Total@Log[v] instead of an explicit Sum. Similarly, v.v or Total[v^2] for the sum of squares. This enables Mathematica to invoke the vector engine which will vastly speed up your code. $\endgroup$ – rcollyer Nov 5 '15 at 17:19
  • 1
    $\begingroup$ @rcollyer. Thanks! I'll update the answer with those suggestions. (Such suggestions is why I really like this forum.) $\endgroup$ – JimB Nov 5 '15 at 17:24
  • 1
    $\begingroup$ BTW: RayleighDistribution[] is built-in. This can now be used with FindDistributionParameters[] to perform maximum likelihood fitting. $\endgroup$ – J. M. is away Nov 5 '15 at 21:46
  • 2
    $\begingroup$ A few additional notes: 1. the Boltzmann constant is known to Mathematica: QuantityMagnitude[UnitConvert[Quantity["BoltzmannConstant"]]]. 2. If a goodness-of-fit is also needed, then one uses DistributionFitTest[] for that. $\endgroup$ – J. M. is away Nov 6 '15 at 5:39
  • 1
    $\begingroup$ Here are the 2d velocity data: pastebin.com/rtteUpk4 $\endgroup$ – mrz Nov 9 '15 at 13:51

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.