Determination of the constant of gravity

The data

Importing the data

We import the data from a csv-List and display them in a grid:

data=Import["/Users/ihn/Documents/Teaching/VP-Leitung/DataAnalysis/2013NewProgram/Lectures given in 2013/7. Lecture/Data_z_t.csv"];
Grid[Prepend[data,{Style["z (cm)",Bold],Style["t (ms)",Bold]}],Background->{None,{Gray,{LightGray,White}}},Dividers->{Black,{Black,Black}},Frame->True]

Converting the data to SI units

In the data list, the lengths (Subscript[z, j]) are given in cm and the times are given in ms. We convert to m and s.

data=Table[{1. (10^(-2))data[[j,1]],(10^(-3))data[[j,2]]},{j,1,Length[data]}];
Grid[Prepend[data,{Style["z (m)",Bold],Style["t (s)",Bold]}],Background->{None,{Gray,{LightGray,White}}},Dividers->{Black,{Black,Black}},Frame->True]

Plotting the data

Since the data are described by a power law, we plot them on double  logarithmic scale and expect a straight line.

g1=ListLogLogPlot[data,PlotRange->{{0.05,1.30},{0.1,0.55}},Ticks->{{0.1,0.2,0.3,0.5,1},{0.1,0.2,0.3,0.5}},FrameLabel->{"z (m)","t (s)"},LabelStyle->24,PlotMarkers->Automatic,ImageSize->600,Frame->True]

Analysis using Model 1 following the lecture

Defining a function for calculating the mean square error

The mean square error is defined as

Q(g) = (1/n)(∑_(j=1)^n)((((Subscript[t, j]^m)-h(g,(Subscript[z, j]))))^2)

We define this function in Mathematica using our specific fitting function
h(g,(Subscript[z, j])) = (Sqrt[((2 (Subscript[z, j]))/g)]).

Q[g_,data_]:=(1/(Length[data]))(∑_(j=1)^(Length[data]))(((data[[j,2]]-(Sqrt[((2 data[[j,1]])/g)])))^2);

Plotting the mean square error function

Plot[Q[g,data],{g,9.759,9.779},PlotRange->{{9.759,9.779},{0.285 (10^(-6)),0.33 (10^(-6))}},LabelStyle->24,PlotStyle->{Thickness[0.005]},ImageSize->600]

Numerical minimization of the mean square error

The Mathematica function FindMinimum does the numerical minimization:

minQ=FindMinimum[Q[g,data],{g,9.78}];

From the result we extract the least mean square error (Subscript[Q, min])...

Qmin=minQ[[1]];

... and the estimated parameter value (OverHat[g]):

gg=g/.minQ[[2]];

The estimate ((OverHat[σ])^2) is then given by

σest2=((Length[data])/(Length[data]-1))Qmin;

Standard error of the parameters

The uncertainty for ((OverHat[σ])^2) is given by

σσ2=((Length[data]Qmin)/(Length[data]-1))(Sqrt[(2/(Length[data]-3))]);

The uncertainty of (OverHat[g]) is

σg=(((((Length[data]-1)/(2 Qmin))(∂_g) (∂_g) Q[g,data]/.{g->gg}))^(-1/2));

Now we can print our final results, properly rounded:

Print["g = (",Round[gg,0.0001]," ± ",Round[σg,0.0001],") m/(s^2)"]
Print["(σ^2) = (", Round[(10^6)σest2,0.001]," ± ",Round[(10^6)σσ2,0.001],") (ms^2)"]

Graphical determination of the standard error of g

Alternatively we can determine (Subscript[σ, g]) using the graphical method. Here is the plot of the normalized mean square error function (N-1)ΔQ/Q:

Show[{Plot[(Length[data]-1)((Q[g+gg,data]-Qmin)/Qmin),{g,9.759-gg,9.779-gg},LabelStyle->24,PlotStyle->{Thickness[0.005]},ImageSize->600,Frame->True,FrameLabel->{"g-(Overscript[g, ⋀]) (m/(s^2))","(N-1)((ΔQ)/Q)"},FrameTicks->{{0,{σg,Round[σg,0.0001]},{-σg,Round[-σg,0.0001]}},{0,1,2,3},{{σg,"(Subscript[σ, g])"},{-σg,"-(Subscript[σ, g])"}},{}}],Graphics[{{Dashed,Line[{{9.759-gg,1},{9.779-gg,1}}],Line[{{-σg,0},{-σg,3.7}}],Line[{{σg,0},{σg,3.7}}]},{LightGray,EdgeForm[Black],Rectangle[{-σg,0},{σg,1}]}}],Plot[{(Length[data]-1)((Q[g+gg,data]-Qmin)/Qmin)},{g,9.759-gg,9.779-gg},Frame->True,FrameTicks->{{},{},{},{}},PlotStyle->{Thickness[0.005]}]}]

We may find the two solutions of (N-1)ΔQ/Q = 1 numerically:

σg1=sg/.FindRoot[(Length[data]-1)((Q[sg+gg,data]-Qmin)/Qmin)==1,{sg,0.005,0,0.1}]

σg2=sg/.FindRoot[(Length[data]-1)((Q[sg+gg,data]-Qmin)/Qmin)==1,{sg,-0.005,-0.1,0}]

These results agree with our previous error estimate.

Plotting the data with the fitted curve

We plot the resulting fit together with the data: this plot shows us graphically, how well the systematic part of the model describes our data.

Show[g1,LogLogPlot[(Sqrt[((2 z)/gg)]),{z,0.05,1.30},PlotStyle->{Red},Frame->True]]

Inspecting the residuals of the fit

We inspect the fitting residuals, in order to see graphically, if the residuals resemble the assumptions that we have made about the noise in our data.

residuals=Table[{data[[j,1]],data[[j,2]]-(Sqrt[((2 data[[j,1]])/gg)])},{j,1,Length[data]}];
Grid[Prepend[residuals,{Style["z (m)",Bold],Style["r (s)",Bold]}],Background->{None,{Gray,{LightGray,White}}},Dividers->{Black,{Black,Black}},Frame->True]

Systematic deviations of the fitted curve from the data can be detected by plotting the residuals versus the (Subscript[z, j]).

ListPlot[residuals,PlotRange->{{0.05,1.30},{-0.001,0.001}},AxesLabel->{"z (m)","Δt (s)"},LabelStyle->24,PlotMarkers->Automatic,ImageSize->600,AxesOrigin->{0.05,0}]

There is a clear trend in the fitting residuals, which indicates a systematic error.

Advanced topic: Further investigation of the residuals with different techniques (beyond the lecture)

For further inspection we generate a correlation plot:

ListPlot[Table[{(10^3)residuals[[j,2]],(10^3)residuals[[j+1,2]]},{j,1,Length[data]-1}],PlotRange->{{-0.5,0.5},{-0.5,0.5}},AxesLabel->{"(Subscript[r, j]) (ms)","(Subscript[r, j + 1]) (ms)"},LabelStyle->24,PlotMarkers->Automatic,ImageSize->600,AxesOrigin->{0,0},AspectRatio->1]

From the correlation plot we do not easily see any suspicious behavior in this case.

Another way to look at the residuals is a histogram plot. We choose histogram bars with a bin width of (Sqrt[((OverHat[σ])^2)])/2.

Histogram[(10^3)residuals[[Range[1,Length[data]],2]],(10^3){-((11 (Sqrt[σest2]))/4),((11 (Sqrt[σest2]))/4),((Sqrt[σest2])/2)},AxesOrigin->{0,0},LabelStyle->24,AxesLabel->{"(Subscript[r, j]) (ms)","rel. frequency"},ImageSize->600]

This distribution of counts should be brought about by a normal distribution, as assumed in our model. We therefore plot the relative frequency (counts normalized to the number of data points) together with the normal distribution which has (Subscript[σ, est]) as its standard deviation.

Show[Histogram[(10^3)residuals[[Range[1,Length[data]],2]],(10^3){-((11 (Sqrt[σest2]))/4),((11 (Sqrt[σest2]))/4),((Sqrt[σest2])/2)},"Probability",AxesOrigin->{0,0},LabelStyle->24,AxesLabel->{"(Subscript[r, j]) (ms)","rel. frequency"},ImageSize->600],Plot[(((10^3)(Sqrt[σest2]))/2)PDF[NormalDistribution[0,(10^3)(Sqrt[σest2])],x],{x,-(11/4)(10^3)(Sqrt[σest2]),(11/4)(10^3)(Sqrt[σest2])},PlotStyle->{Red,Thickness[0.005]}]]

However, what is now missing, is a way to compare the relative frequency with the red distribution. To this end we may use the number of counts in each bin for estimating the probabilities that a particular measured value of the noise arises within this bin. For this purpose we generate a list of counts in the individual bins of the histogram:

hl=HistogramList[(10^3)residuals[[Range[1,Length[data]],2]],(10^3){-((11 (Sqrt[σest2]))/4),((11 (Sqrt[σest2]))/4),((Sqrt[σest2])/2)}][[2]]

For estimating a probability from the number of counts in each bin, we use the mean of the Beta-Distribution BetaDistribution[k+1,n-k+1]. Below we generate a list of these values, which will later be plotted.

ml=Table[{(j-6)(((10^3)(Sqrt[σest2]))/2),Mean[BetaDistribution[hl[[j]]+1,Length[data]-hl[[j]]+1]]},{j,1,11}];
Grid[Prepend[ml,{Style["bin center (ms)",Bold],Style["⟨p⟩",Bold]}],Background->{None,{Gray,{LightGray,White}}},Dividers->{Black,{Black,Black}},Frame->True]

Below we generate a list of lines that will be used as the error bars for the probability estimate. The error bars span the range from the 5% to the 95% quantile.

eb=Table[Line[{{(j-6)(((10^3)(Sqrt[σest2]))/2),N[Quantile[BetaDistribution[hl[[j]]+1,Length[data]-hl[[j]]+1],0.95]]},{(j-6)(((10^3)(Sqrt[σest2]))/2),N[Quantile[BetaDistribution[hl[[j]]+1,Length[data]-hl[[j]]+1],0.05]]}}],{j,1,11}]

Below we plot the all these results in one plot. The plot comprises the histogram of the data, this time plotted as a relative frequency (counts normalized to the number of data points), the normal distribution that results from our estimated value (Subscript[σ, est]) (red curve with blue dots), and the estimated probabilities with their error bars (blue dots with black error bars).

Show[Histogram[(10^3)residuals[[Range[1,Length[data]],2]],(10^3){-((11 (Sqrt[σest2]))/4),((11 (Sqrt[σest2]))/4),((Sqrt[σest2])/2)},"Probability",AxesOrigin->{0,0},LabelStyle->24,AxesLabel->{"(Subscript[r, j]) (ms)","rel. frequency"},ImageSize->600],ListPlot[ml,PlotMarkers->Automatic],Plot[(((10^3)(Sqrt[σest2]))/2)PDF[NormalDistribution[0,(10^3)(Sqrt[σest2])],x],{x,-(11/4)(10^3)(Sqrt[σest2]),(11/4)(10^3)(Sqrt[σest2])},PlotStyle->{Red,Thickness[0.005]}],ListPlot[Table[{k(((10^3)(Sqrt[σest2]))/2),(CDF[NormalDistribution[0,(Sqrt[σest2])],((Sqrt[σest2])/4)+k((Sqrt[σest2])/2)]-CDF[NormalDistribution[0,(Sqrt[σest2])],-((Sqrt[σest2])/4)+k((Sqrt[σest2])/2)])},{k,-5,5,1}],LabelStyle->24,PlotMarkers->Automatic],Graphics[{Thickness[0.005],eb}],AspectRatio->1]

There is one histogram bar suggesting a probability for the bin that can hardly be produced by the normal distribution. (We have to be careful here: since we are dealing with probabilities, such a situation can occur; there is no reason that a normal distribution would not bring about such an exceptional bin. However, we would have to call this a rare event.)

Analysis using Model I1 following the lecture

Defining a new mean square error function

The new fitting function is

(Subscript[t, j])=(Sqrt[((2 ((Subscript[z, j])-(Subscript[z, 0])))/g)])=:h(g,(Subscript[z, 0]),(Subscript[z, j]))

The parameters of the model are now g, (Subscript[z, 0]), and σ. The mean square error function of model II is now

Q[g_,z0_,data_]:=(1/(Length[data]))(∑_(j=1)^(Length[data]))(((data[[j,2]]-(Sqrt[((2 (data[[j,1]]-z0))/g)])))^2);

Plotting the mean square error function

We use a contour plot to look at the mean square error function around its minimum.

ContourPlot[Q[g,z0,data],{g,9.78,9.82},{z0,-0.003,-0.001}]

The contour lines are ellipses with their main axes tilted with respect to the coordinate axes. This indicates that the parameters g and (Subscript[z, 0]) are correlated.

Numerical minimization of the mean square error

minQ=FindMinimum[Q[g,z0,data],{g,9.78},{z0,0}];
Qmin=minQ[[1]];
gg=g/.minQ[[2]];
zz0=z0/.minQ[[2]];
σest2=((Length[data])/(Length[data]-2))Qmin;

Standard error of the parameters

σσ2=((Length[data])/(Length[data]-2))Qmin(Sqrt[(2/(Length[data]-4))]);

In order to find the uncertainties in the estimated parameters g and (Subscript[z, 0]) we calculate the Hesse Matrix:

hessematrix=(({{(∂_g) (∂_g) Q[g,z0,data]/.{g->gg,z0->zz0}, (∂_z0) (∂_g) Q[g,z0,data]/.{g->gg,z0->zz0}}, {(∂_g) (∂_z0) Q[g,z0,data]/.{g->gg,z0->zz0}, (∂_z0) (∂_z0) Q[g,z0,data]/.{g->gg,z0->zz0}}}));

Inversion of the Hesse matrix gives the matrix of the variances and covariances of the parameters:

correlationmatrix=Inverse[((Length[data]-2)/(2 Qmin))hessematrix];
σg=(Sqrt[correlationmatrix[[1,1]]]);
σz0=(Sqrt[correlationmatrix[[2,2]]]);
ρ=((correlationmatrix[[1,2]])/(σg σz0));

Here we print our final results:

Print["g = (",Round[gg,0.0001],"±" ,Round[σg,0.0001],") m/(s^2)"]
Print["(Subscript[z, 0]) = (",Round[(10^3)zz0,0.01],"±" ,Round[(10^3)σz0,0.01],") mm"]
Print["(σ^2) = (",Round[(10^6)σest2,0.001],"±" ,Round[(10^6)σσ2,0.001],") (ms^2)"]
Print["(Subscript[ρ, g, (Subscript[z, 0])]) = ",Round[ρ,0.01]]

Graphical determination of the standard error of g and (Subscript[z, 0])

We replot the mean square error function in order to illustrate the estimated values (white spot) and their uncertainties. The innermost (red) contour line is at (Subscript[Q, min])(1+(1/(n-2))). Its extent in g-direction indicated by the vertical black dashed lines gives the ±(Subscript[σ, g]) range around the estimate. Its extent in (Subscript[z, 0])-direction indicated by the horizontal black dashed lines gives the ±(Subscript[σ, z0]) range around the estimate of (Subscript[z, 0]). The green dashed line indicates the intersection of the red contour line with the white horizontal line. It has a distance (Subscript[σ, g])(Sqrt[1-(ρ^2)]) from the white vertical dashed line and can therefore serve for the graphical determination of ρ.

Show[ContourPlot[Q[g,z0,data],{g,gg-(3/2)σg,gg+(3/2)σg},{z0,zz0-(3/2)σz0,zz0+(3/2)σz0},Frame->True,Axes->True,AxesOrigin->{gg-(3/2)σg,zz0-(3/2)σz0},Contours->Table[Qmin(1+(k/(Length[data]-2))),{k,1,10}],ContourStyle->{{Red,Thickness[0.007]},Blue,Blue,Blue,Blue,Blue,Blue,Blue,Blue,Blue},FrameTicks->{{{gg,Round[gg,0.0001]},{gg+σg,Round[gg+σg,0.0001]},{gg-σg,Round[gg-σg,0.0001]}},{{zz0,Round[(10^3)zz0,0.01]},{zz0-σz0,Round[(10^3)(zz0-σz0),0.01]},{zz0+σz0,Round[(10^3)(zz0+σz0),0.01]}},{{gg,"(Overscript[g, ⋀])"},{gg+σg,"(Overscript[g, ⋀])+(Subscript[σ, g])"},{gg-σg,"(Overscript[g, ⋀])-(Subscript[σ, g])"}},{{zz0,"(Subscript[(Overscript[z, ⋀]), 0])"},{zz0-σz0,"(Subscript[(Overscript[z, ⋀]), 0])-(Subscript[σ, (Subscript[z, 0])])"},{zz0+σz0,"(Subscript[(Overscript[z, ⋀]), 0])+(Subscript[σ, (Subscript[z, 0])])"}}},FrameLabel->{"g (m/(s^2))","(Subscript[z, 0]) (mm)"}],Graphics[{{White,Dashed,Line[{{gg,zz0-(3/2)σz0},{gg,zz0+(3/2)σz0}}],Line[{{gg-(3/2)σg,zz0},{gg+(3/2)σg,zz0}}]},{White,Disk[{gg,zz0},{σg/30,σz0/30}]},{Black,Dashed,Line[{{gg+σg,zz0-(3/2)σz0},{gg+σg,zz0+(3/2)σz0}}],Line[{{gg-σg,zz0-(3/2)σz0},{gg-σg,zz0+(3/2)σz0}}],Line[{{gg-(3/2)σg,zz0-σz0},{gg+(3/2)σg,zz0-σz0}}],Line[{{gg-(3/2)σg,zz0+σz0},{gg+(3/2)σg,zz0+σz0}}]},{Green,Dashed,Line[{{gg+σg(Sqrt[1-(ρ^2)]),zz0-(3/2)σz0},{gg+σg(Sqrt[1-(ρ^2)]),zz0+(3/2)σz0}}]},{Green,Text[Style[Round[gg+σg(Sqrt[1-(ρ^2)]),0.0001],24],{gg+σg(Sqrt[1-(ρ^2)]),zz0-(3/2)σz0},{0,-1},{1,0},Background->White]},{Green,Text[Style["(Overscript[g, ⋀])+(Subscript[σ, g])(Sqrt[1 - (ρ^2)])",24],{gg+σg(Sqrt[1-(ρ^2)]),zz0+(8/7)σz0},{0,-1},{1,0},Background->White]}}],AxesLabel->{"g (m/(s^2))","(Subscript[z, 0]) (m)"},LabelStyle->24,ImageSize->600,Ticks->{{Round[gg-σg,0.0001],Round[gg,0.0001],Round[gg+σg,0.0001]},{Round[zz0-σz0,0.00001],Round[zz0,0.00001],Round[zz0+σz0,0.00001]}}]

Plotting the data with the fitted curve

Show[g1,LogLogPlot[(Sqrt[((2 (z-zz0))/gg)]),{z,0.05,1.30},PlotStyle->{Red}]]

Inspecting the residuals of the fit

Again we inspect the fitting residuals:

residuals=Table[{data[[j,1]],data[[j,2]]-(Sqrt[((2 (data[[j,1]]-zz0))/gg)])},{j,1,Length[data]}];
Grid[Prepend[residuals,{Style["z (m)",Bold],Style["r (s)",Bold]}],Background->{None,{Gray,{LightGray,White}}},Dividers->{Black,{Black,Black}},Frame->True]

ListPlot[residuals,PlotRange->{{0.05,1.30},{-0.001,0.001}},AxesLabel->{"L (m)","Δt (s)"},LabelStyle->24,PlotMarkers->Automatic,ImageSize->600,AxesOrigin->{0.05,0}]

We see that the improved model has removed the trend in the residuals.

Advanced topic : Further investigation of the residuals with different techniques (beyond the lecture)

ListPlot[Table[{(10^3)residuals[[j,2]],(10^3)residuals[[j+1,2]]},{j,1,Length[data]-1}],PlotRange->{{-0.5,0.5},{-0.5,0.5}},AxesLabel->{"(Subscript[r, j]) (ms)","(Subscript[r, j + 1]) (ms)"},LabelStyle->24,PlotMarkers->Automatic,ImageSize->600,AxesOrigin->{0,0},AspectRatio->1]

The same histogram analysis done above leads to

σest=(Sqrt[σest2]);
hl=HistogramList[(10^3)residuals[[Range[1,Length[data]],2]],(10^3){-((11 σest)/4),((11 σest)/4),((σest)/2)}][[2]]

ml=Table[{(j-6)(((10^3)σest)/2),Mean[BetaDistribution[hl[[j]]+1,Length[data]-hl[[j]]+1]]},{j,1,11}];
Grid[Prepend[ml,{Style["bin center (ms)",Bold],Style["⟨p⟩",Bold]}],Background->{None,{Gray,{LightGray,White}}},Dividers->{Black,{Black,Black}},Frame->True]

eb=Table[Line[{{(j-6)(((10^3)σest)/2),N[Quantile[BetaDistribution[hl[[j]]+1,Length[data]-hl[[j]]+1],0.95]]},{(j-6)(((10^3)σest)/2),N[Quantile[BetaDistribution[hl[[j]]+1,Length[data]-hl[[j]]+1],0.05]]}}],{j,1,11}]

Show[Histogram[(10^3)residuals[[Range[1,Length[data]],2]],(10^3){-((11 σest)/4),((11 σest)/4),((σest)/2)},"Probability",AxesOrigin->{0,0},LabelStyle->24,AxesLabel->{"(Subscript[r, j]) (ms)","rel. frequency"},ImageSize->600],ListPlot[ml,PlotMarkers->Automatic],Plot[(((10^3)σest)/2)PDF[NormalDistribution[0,(10^3)σest],x],{x,-(11/4)(10^3)σest,(11/4)(10^3)σest},PlotStyle->{Red,Thickness[0.005]}],ListPlot[Table[{k(((10^3)σest)/2),(CDF[NormalDistribution[0,σest],((σest)/4)+k((σest)/2)]-CDF[NormalDistribution[0,σest],-((σest)/4)+k((σest)/2)])},{k,-5,5,1}],LabelStyle->24,PlotMarkers->Automatic],Graphics[{Thickness[0.005],eb}],AspectRatio->1]

The red curve runs reasonably within the error bars of the probability estimate, confirming qualitatively that the assumption of a normal distribution for the scatter of the time measurement is reasonable. The exceptional bin encountered in the analysis with the previous model has disappeared. However, there is nothing in this comparison that would tell us under which circumstances we better reject the hypothesis of a normal distribution. The situation could be made to look much worse, if we were to plot error bars spanning only the range between the 30% and the 70% quantile.

Model I, a shortcut offered by Mathematica

The command NonlinearModelFit performs most of the tasks that we did above:

nlm=NonlinearModelFit[data,(Sqrt[((2 z)/g)]),{g},z];

The functional form of the fitted model is:

Normal[nlm]

We can plot the data together with the result of the fit:

Show[ListPlot[data],Plot[nlm[z],{z,0.1,1.25}],Frame->True]

The fitting parameters and their standard errors are also available:

nlm["ParameterTable"]

The list of fit residuals can be extracted...

nlm["FitResiduals"]

... and plotted.

ListPlot[nlm["FitResiduals"],Filling->Axis]

Model II, using the shortcut offered by Mathematica

We use a constraint for the offset making sure that the values under the square root are certainly positive. Furthermore, we use the constraint that g is positive.

nlm=NonlinearModelFit[data,{(Sqrt[((2 (z-offs))/g)]),{-0.1<offs<0.1,g>0}},{g,offs},z];

The functional form of the fitted model is:

Normal[nlm]

We can plot the data together with the result of the fit:

Show[ListPlot[data],Plot[nlm[z],{z,0.1,1.25}],Frame->True]

The fitting parameters and their standard errors are also available:

nlm["ParameterTable"]

The list of fit residuals can be extracted...

nlm["FitResiduals"]

... and plotted.

ListPlot[nlm["FitResiduals"],Filling->Axis]