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.)
{bins, counts} = HistogramList[...]
is much easier to write. 3. Element-1
of a list is the same as elementLength[list] - 1
. 4. I didn't have time to understand yourmodel
fully, but it looks weird. It should be a simple formula, why do you mix in the data? 5. Get the bin centres withMovingAverage[bins,2]
. Then you have the same number of values as incounts
. $\endgroup$ – Szabolcs Nov 5 '15 at 16:04