Avogadro Constant

The Avogadro constant is a fundamental physical constant that plays a crucial role in the field of chemistry. In our first lecture, we discussed one way to determine the Avogadro constant from Brownian motion. This technique was very relevant at the beginning of the 20th century. A list of historic measurements can be found in table 1 of Peter Becker, Rep. Prog. Phys. 64 (2001) 1945–2008.

In this exercise we will test if applying the Kalman filter to the measurement enhances the accuracy of the estimated constant.

• Select at least 15 values with corresponding uncertainties from the table of historic measurements (do not consider the 10 latest measurements) and plot them over time in a plt.errorbar()-plot in Python.
• Plot the histogram of the measured values, calculate the mean and the variance, and display both in the plot
• Implement a 1-dimensional Kalman filter in Python and apply the filter to the series of measurements. Plot the estimates for the mean and the variance at each time step display the final mean value as plt.text()
• How does the error converge?
• Is this realistic?
• Please explain possible reasons for this observation
• Rescale the error to get a realistic error convergence and plot the estimate for the mean and the variance at each timestep

Metropolis-Algoritmus

In thermodynamics, the probability density of paricles in a Potential V(x) is given by the Boltzmann distribution. It reads (if constants are set to 1) \[ \rho(x)=e^{-V}. \] We want to study the case of a double-well potential \[ V(x) = x^4-x^2. \]

The Metropolis algorithm is a special case of the Metropolis-Hastings algorithm. It is a Markov Chain Monte Carlo method that uses one artificially generated series (a Markov chain) to approximate the integral \[ \langle f(x) = \int f(x)\rho(x)\mathrm{d}x. \]

The algorithm consists of the following steps:

• Create an array x of length 100000 to save all the states
• Choose an initial condition, e.g. x[0]=3

• At each step a new state \[x_t+q_t\] is proposed. The acceptance probability of the new state is given by \[ \min\left[ 1, \frac{\rho(x_t+q_t)}{\rho(x_t)} \right]. \] If the proposal is rejected the new state is identical to the old state xt+1=xt. Repeat the following steps for t=0 to t=99999
• create a new proposal using qt=np.random.normal(0,1)
• create a random varible at=np.random.uniform(0,1)
• calculate the acceptance probability pt of the proposed state xt+qt; if a>pt, reject the state and define the new state as the old state; otherwise, accept the new state and define the new state as \[x_{t+1}=x_t+q_t\]

Plot the generated series x against t. What is the problem with the first values of the series? Remove the initial part of the series and only keep the part, where states are within the region of interest (x=x[1000:])

Plot the histogram of the series (plt.hist(x,bins=40,density=True)); also plot the function rho in a suitable range and compare both plots; was the approximation of the density successful?

Calculate the first, second, and fourth moment of the distribution, using the series x