Best Choice for the Expected Error

  In the last blog we found two values for the estimate of σ, the expected value, μ₁, and the rms error, √μ₂ = √n, and we have to ask which is the best choice. The expected value is the mean that one would get for an estimate of the error but it has some uncertainty, sd, associated with it. If we want to combine these two uncertainties we have to add their squares and the result is the second moment, μ₂ = n.

If the sum involves a short string of numbers then we are most likely to random values for z near the peak of the distribution where
z_peak = √n-1.

The rms error, s, for the sum is then equal to σ√n-1. We can turn things around to get an estimate of σ for the probability distribution from that of s and we conclude that σ = s/√n-1.

Expected Error For the Sum of a Number of Trials

  A scientist like any other witness can be impeached for using poor practices.  So it is important for him to maintain his reputation for accuracy. When one publishes the results of an experiment it is customary to include the error bounds for numerical estimates. Averaging a number of trials gives a better estimate of an observation and involves the sum of a number of individual observations, x_k. What error error would one expect for this sum? The formula that one usually uses is that the error for the sum is equal σ√n, where σ is the standard deviation for x_k. The following proof is based on the expected values involved. The x_k are assumed to be equal to the mean value, μ, plus a random error which can be expressed in terms of z-values.

The square of the difference between the sum and its mean value, nμ, involves a sum of the products of the individual z-values which can be split up into parts where the indices j and k are equal and unequal. The term that is crossed out makes a small contribution to the sum and on average is zero since the two z-values are uncorrelated.

  A more sophisticated analysis makes better use of probability theory.  We start by determining the probability distribution for the sum of two and three trials assuming normal distributions for each of the random variables and we need only consider the sum of the z-values, x + y and x + y + z.  The square of the difference between the sum and its mean involves a sum of squares which is also found in the exponential of the joint probability distribution and so we can just consider the distribution for the radial length of a vector whose components are the z-values involved in the sum.

Notice that the area and volume contribute powers of r to the joint probability distribution and for n terms in the sum the contribution of the spatial unit space's differential element will be some constant times r^n-1 which can be determined from the condition that the integral of the probability density is 1. Having found the probability function we can determine formulas the 1st and 2nd moments and the standard deviation for r.

To test these formulas we can we can generate two sets of random numbers, create a set of sums and observe the resulting distributions. The formulas give a very good fit for the generated data. The mean value for r (or z below) is slightly less than √n.

For a sum of more numbers we can combine a new set of random numbers with the previous sum. For the sum of three numbers we get:

Continuing the process we get for the sum of 10 numbers:

We got n for the sum of the squares vindicating the previous result. Note that the joint probability functions are similar in shape to the Poisson distribution and may be considered continuous analogs. The standard deviation of the mean values is small and levels off after the sum of about 20 numbers and becomes approximately equal to √2/2.

Supplemental (Dec 2): The definite integral above is evaluated in Reif, Fundamentals of Statistical and Thermal Physics (p. 608). It is also found in Boltzmann, Lectures on Gas Theory (p. 64).

Edit (Dec 3): The differential elements are surface lamina of nD hyperspheres.

Using Grid Searches to Optimize a Fit

  We were able to estimate the reaction rate constants from functions of three unknown parameters and to improve the fit with a simple search for k₀. For a better fit we need to search for all three unknown parameters. We can choose X₀, X_eq and k₀ as the unknown parameters with X = A - A₀ being the amount of product at some particular time. X simplifies the formulas which is the reason why it was chosen. We have estimates for X_eq and k₀ determined by the data points and can search the neighborhood of the values for the best fit. So a three dimensional grid search will allow us to improve on our initial estimates. Here is a sample of how one might start out. The first three lines determine the bounds of the search and the next three compute the grid points. X₀ was initially assumed to be a small positive number.

I wrote a Mathcad program to do the search for the integers λ, μ and ν which give the lowest value for the variance, V, of A and A'.

Three such searches zoomed in on values very close to the original values.

They give a relative rms error for the fit of about three tenth of one percent. The fit of the "lower left corner" is also much better.

An error propagation formula will help us determine how errors in the grid search parameters will affect the final results. We can estimate how two independent variables, x and y, affect a third dependent variable, z, starting with the chain rule. If dx and dy represent small errors then the chain rule gives us a good estimate of dz. The formula can be modified to give the rms error of dz but we need to use an inequality to express it in terms of the rms errors of dx and dy. One can choose the right hand side of the last equation to represent a bound for the rms error of dz. One usually sees this propagation of error formula written as an equality.

The formula for the propagation of error can be used to give the expected relative errors of the computed parameters for the fit based on the relative errors of the fit variables. The matrix Λ was found by applying the propagation of error formula to each of the computed variables and evaluating it using estimated values. The differences between the fit values and original values can be used to get approximations for the relative errors of the fit variables. The bottom left shows the expected bounds and can be thought of as standard deviations. The bottom right shows the z-scores for the the actual relative errors.

The grid search we used isn't optimal but it allows us to improve on the original estimate which a plot showed to be slightly off. The final result is what can be expected from a grid search.

Supplemental (Nov 12): The program p_min(V) puts the rows of a table containing the pointer p and corresponding V values in the order of the V values. The resulting location for the minimum value of V is the first p value on the new table.

Fit For When Two Reactants are Consumed

  The fit for problem worked in the blog Increasing the Rate of Dissociation by Adding Acid is similar to that in the last blog. It's what one would have to do for the case where the two reactants are consumed in the reaction. There was a minor complication with the complex κ.

The estimate value for κ turned out to be close to the complex conjugate of the original κ. For z = x + i y the complex conjugate is z* = x - i y. It turns out that tanh(κ*) = tanh(κ) when the result is real. There is another formula for the tanh of a complex number involving ordinary and hyperbolic sines and cosines that helps explain the situation. Either the cosine or the sine of the imaginary part in the result has to be equal to zero.

Some data sets are more difficult to fit than others. One has to play with the numbers of terms in the empirical fit for best results. The empirical fit for the initial rate of change might be good but that doesn't guarantee that the fit for the formula will be good. The rounded corner at the lower left of the curve below seems to be hardest to get right. The procedure is sensitive to random errors in the first few data points so some method of adjusting errors might prove useful.

Supplemental (Nov 8): A weakness of the fit procedure is that an error in k'₀ propagates through as one calculates all the other fit constants so if one searches the neighborhood of k'₀ for the value which gives a minimum for the variance of the difference between A' and A then one will get consistently good fits to the data points. The resulting fit is a least squares fit of the data.

Fitting a Simple Catalyzed Reaction

  In a catalyzed reaction the catalyst is not used up in the reaction. With the dissociation of sucrose in an acid solution water takes part in the reaction but the quantity of acid and its ions remain constant. The quantity of hydrogen ions in solution is maintained by the dissociation of water. We can take this into account and it simplifies the equations for the constants in the formula for the amount of sucrose as a function of time. We proceed as we did before determining the coefficients for the quadratic equation for the reaction rate but let B remain unchanged. Integrating the reaction rate equation again produces a formula with a hyperbolic tangent function but the constants in it are slightly different due to changes in relative magnitude. We can assume some initial conditions for the sake of argument.

We then generate some random data to work with and allow sufficient time for equilibrium to be reached.

If one checks the initial reaction rate one gets dA/dt = −k₀B₀A which agrees with the formula that Wilhelmy used so the initial solution is approximately an exponential decay function. To make the change easier to fit we can use the natural log of A instead of A and fit a polynomial equation. F(t) below is a vector function containing the powers of t and is needed by linfit to find the coefficients of the polynomial. We can estimate the initial value for A by averaging the values of A divided by powers of e with x values corresponding to the data points. A check of the standard deviation shows that the rms error is about one unit.

Once we have the fit we can draw a smooth curve through the data points using the empirical power law.

By taking the derivative of the fit polynomial at t = 0 we get the initial rate constant and this gives us an estimate of the first rate constant, k₀, on division by B₀. To find the second rate constant, k₁, we need an estimate of the equilibrium value for A and the equilibrium condition k₀B₀A_eq = k₁(A₀ - A_eq)² which states that the amount of product that is produced is canceled by the reverse reaction.

With the estimates for the rate constants we can now estimate the constants needed to fit the formula to the data. Using an auxiliary constant γ makes it easier to determine an estimate for κ. The fit agrees fairly well with the initial assumptions.

Increasing the Rate of Dissociation by Adding Acid

  Sucrose is a disaccharide produced by the hydrolysis of glucose and fructose molecules in which an hydroxyl group (-OH) in each sugar molecule join producing an oxygen bridge between them by releasing an H₂O water molecule. One can see the oxygen atom linking the hexagonal ring (glucose) to the pentagonal ring (fructose) in the .gif below.

animation credit: Wikipedia article Sucrose

Dissociating a disaccharide molecule requires a hydration reaction to replace the missing atoms of the simple sugars. There are always hydrogen and hydroxyl radicals normally present in a sugar solution which can react with the sucrose molecules but adding an acid as Wilhelmy did acts as a catalyst to increase the reaction rate.

  To study the effect of adding the acid to the solution we need to modify the equation for the reaction and deduce the corresponding reaction rate equation. In the reaction below we can let A correspond to the sucrose molecule, B to the hydrogen molecules present in the solution resulting from the addition of acid and C and D to the glucose and fructose products. We can also let X represent the progress of the reaction in order to determine the amount of sucrose and other molecules present at any given time. The acid molecules easily dissociate so the amount of hydrogen present is essentially equal to the concentration of acid. Plugging the modified amounts into the reaction rate formula again reduces it to a first order differential equation involving the sucrose concentration [A].

The reaction again follows a binomial rate law and we can simplify it by "completing the square" and rewriting the equation using the auxiliary parameters u and v.

The result is a similar formula involving the hyperbolic tangent function but here λ is replaced by βc so the net effect is that the reaction rate constants are multiplied by the factor β increasing the reaction rate considerably. We can check our math by working an example. From the positive polynomial coefficients a, b and c in the rate equation we can compute the constant parameters α and β and κ. Here κ is complex but that isn't a problem since we can use tan(iθ) = i tanh(θ) with a complex θ. The complex κ is easily handled by Mathcad. The complex value is the results from ([A]₀ + α)/β being greater than 1.

The results from an iterative evaluation of the reaction rate formula and those from the derived formula produce identical curves for the reaction, justifying the use of the derived formula. Setting the reaction rate polynomial to zero allows us to compute the equilibrium sucrose concentration A_eq which agrees with the result of the numerical computation.

Supplemental (Nov 4): Try k₀=1×10^-4, k₁=1.234567901234567×10<sup>-6

Suppemental (Nov 5): For the hyperbolic tangent of a complex number one can use the following identity which is similar to the identity for tanh(u+v).</sup>

Even though κ is complex, tanh(κ) is real and consequentially, A is real.

Note also that since tanh(κ + βct) > 0, the new formula for A subtracts 1 from the tanh instead of subtracting the tanh from 1.

Wilhelmy's 1850 Paper on the Dissociation of Cane Sugar

  Wilhelmy published a paper, Ueber das Gesetz, nach welchem die Einwirkung der Säuren auf den Rohrzucker stattfindet (On the law, which takes place after the action of acids on cane sugar) in Annalen der Physik und Chemie in 1850 in which he gives the amount of dissociation of cane sugar (sucrose) into its simpler components (glucose and fructose) as a function of time. One can observe the change in the optical rotation of the sugar solution and determine the rate at which the amount of sugar decreases by means of an inversion coefficient, μ. Wilhelmy uses the optical rotation angle corresponding to the amount of sugar present to also specify its quantity, Z, present in the solution. When a portion, X, of the amount of sugar in solution is inverted it no longer rotates the light passing through the polarimeter X degrees to the right but instead rotates it μX degrees to the left. Wilhelmy explains as follows,

 "The herein stated task will now be to determine, depending on whether and in what way M depends on the different physical conditions of the process, i.e., whether and in what manner M is a function of time, the amount of sugar, the amount of acid, the amount of Resolution agent, the quality of the acid, the temperature and air pressure. Answering these questions should be attempted in the following order.
  The above however is yet an important point to discuss for the calculation of the experiments. The action of the acid on the sugar not only removes a portion of the dextrorotatory sugar but transforms it into opposite rotation. What therefore was originally the rotation =Z°, is, after a rotating X° amount of sugar, converted to =(Z - X - μX)°. To find out that given by observation quantity X itself, one needs to know μ. Hr. Biot has introduced for μ, i.e., for the amount of reverse rotation to the left, produced by a quantity of sugar which had a clockwise rotation of 1°, the name of the inversion coefficient. X then also gives Z = Z₀ - X. The determination of the inversion coefficients has caused some initial difficulties for me. In addition Hr. Biot has been wavering in his statements about the same. In his work on the sugar content of corn ¹) he gives for μ
     for hydrochloric acid  = 0.38
     for sulfuric acid         = 0.3867,
It can be said however that he himself noticed that he received very different results, which he thinks he can ascribe to different purities of the sugar. In a later essay ²) he then gives:
     for sulfuric acid         μ = 0.417
     for nitric acid             μ = 0.394
     for hydrochloric acid  μ = 0.38.
  Thus it is initially to be noted that Hr. Biot, so far as I know, did not state for which temperature these coefficients apply; but since the rotatory power of galactose is dependent on the temperature, as well as the value of μ must be also differ according to the temperature of the reading. If one therefore wishes to know the value of μ for each temperature, one must first determine the law according to which the rotatory power of the galactose is dependent on the temperature."

¹) Comptes rend. 15, 529.
²) Comptes rend. 17, 757.

  The results of Wilhelmy's observations are presented in Table III of his paper. The column headings are time, T, optical rotation, Drehung, the log of the relative change in the amount of sugar, Z, and the temperature. There is a footnote at the bottom of the page which gives a formula for computing the amount of sucrose, Z, present given the rotation angle, D. One can derive the formula as follows,

If we plot Wilhelmy's data we get a reasonably good fit to a decaying exponential function.

Supplemental (Nov 1): The value of the inversion coefficient used for the experimental data is μ = 0.3966. Wilhelmy refers to Biot's corn sugar in German as Schleimzucker which is literally translated as "slime sugar" but in the German Wikipedia it is another name given for galactose. Wilhelmy's definition of the inversion coefficient appears to be based on the relation Z₀ - D = (1 + μ)(Z₀ - Z) where the quantity or concentration of Z is given in units of 1° rotation. For problems in measuring the optical rotation associated with dissociation of complex molecules see the Wikipedia article on mutarotation.