gtMath

Preliminaries: How do mathematicians model randomness?, Monte Carlo Simulation - Part 1

In this post, a continuation of the reader-requested Monte Carlo Simulation - Part 1, I will present another application of Monte Carlo methods and a concrete example (the promised $\pi$ approximation), including the C++ source code so that you can see the "guts" of a Monte Carlo simulation.

I will also explain how the $\pi$ example relates to the Law of Large Numbers and Central Limit Theorem, key probabilistic results which apply beyond just Monte Carlo methods.

Monte Carlo approximation of definite integrals (and $\pi$)

In Part 1, I presented a retirement planning asset allocation problem as an example of when one may employ a Monte Carlo method. Given the enormous number of possible scenarios under the many statistical assumptions, this problem was a prime Monte Carlo candidate.

I will now show you another use of Monte Carlo methods: estimating the value of a definite integral $\int_{a}^{b}{f(x) \, dx}$, i.e. the area under a curve $f(x)$ between $x=a$ and $x=b$, when we cannot calculate it analytically (via the antiderivative). This is especially useful for multidimensional integrals, but I will stick to single-dimensional integrals here just to illustrate the method.

For example, suppose we would like to estimate the value of the definite integral $I=\int_{0}^{\pi}{\frac{\sin(x)}{x} \, dx}$. We have a computer choose many random sample points (from the uniform distribution on the rectangle $\left[ 0, \pi \right] \times \left[ 0,1 \right]$) and count how many fall below the curve (the green points in the below image, courtesy of Google):


In the image above, the computer has chosen 314,028 points so far, of which 185,208 were below the curve and 128,820 were above. We thus estimate the value of $I$ as $\frac{{\text 185,208}}{{\text 185,208}+{\text 128,820}} \cdot \pi \approx 1.852854$. This is the percentage of sample points in the rectangle which fell below the curve, multiplied by the total area of the sample rectangle ($1 \cdot \pi$).

To approximate the value of $\pi$, we can use the same technique: we draw a circle of radius 1 centered at $(0,0)$ and choose $N$ random sample points $(x,y)$ from the uniform distribution on the square $\left[ -1, 1 \right] \times \left[ -1,1 \right]$. The total area of the square is 4, and the area of the circle is $\pi \cdot 1^{2} = \pi$. Thus, we would expect the percentage of sample points falling within the circle (the red points in the below animation) to be approximately $\frac{\pi}{4}$.

After running our program, we can count the number $N_1$ of sample points which fell within the circle (a point falls within the circle if its radius is less than 1, i.e. if $\sqrt{x^2 + y^2}<1$). Our $\pi$ approximation is then $\pi \approx 4 \tfrac{N_1}{N}$.

Animation of the Monte Carlo method to estimate $\pi$ (image courtesy of Wikipedia)

The below chart shows three of my simulation runs using the above method, with 100,000 sample points per run ($N= \text{100,000}$). Based on the chart, the approximation seems to generally improve as $N$ increases; while this is not exactly surprising, the Law of Large Numbers, which I will explain in the next section, formalizes this notion. The Law of Large Numbers states that approximations like this one will always yield sample averages close to the true value as $N \rightarrow \infty$. In other words, I didn't just run my simulation 100 times and choose three "lucky" results to show you.

Digging a bit deeper, it appears that $N=\text{10,000}$ would suffice for an approximation within about 0.02-0.03 of the actual value, but if we want precision within 0.01, we need an $N$ somewhere in the neighborhood of 50-100,000 samples. The Central Limit Theorem, which I will also address later in this post, tells us exactly how large we need to make $N$ in order to obtain, with a given level of confidence (e.g. 95%), an approximation with a given level of precision (e.g. 0.01).


I have included the C++ code for my simulation at the end of the post so that interested readers can see the details.

The Law of Large Numbers

The Law of Large Numbers formalizes the notion that sample averages converge to the long-run average. There are actually two versions, known as the Weak Law and Strong Law of Large Numbers, which correspond to different types of random variable convergence. Let's take a look at each and interpret them in the context of the $\pi$ simulation.

For the below definitions, let $X_1, X_2, X_3, \dotsc$ be a sequence of random variables, independent from each other, which all have the same probability distribution. These are called independent and identically distributed (iid for short). Since the distributions are identical, the mean of each $X_i$ is the same value, $\mu$, which we will assume is finite. The symbol $\overline{X}_n$ will denote the sample average, defined as $\frac{\sum_{i=1}^{n}{X_i}}{n}
$.

The Weak Law of Large Numbers states that, as $n \rightarrow \infty$, the sample average converges in probability to $\mu$, which means that, for any threshold $\epsilon > 0$ (no matter how small), we have $$
\lim_{n \rightarrow \infty} {\Bbb P} \left( \left| \overline{X}_n - \mu \right| > \epsilon \right) = 0
$$ In other words, for any $\epsilon$ specified (e.g. one thousandth, one millionth, etc.), the probability that the sample average will be further than $\epsilon$ from $\mu$ can be made arbitrarily small (but not necessarily zero) by choosing a large enough sample size $n$.

The Strong Law of Large Numbers states that the sample means converge almost surely to $\mu$ (almost surely means "with probability 1"), or in symbols, $$
{\Bbb P} \left( \lim_{n \rightarrow \infty}{\overline{X}_n} = \mu \right) = 1
$$ This is equivalent to the statement $$
{\Bbb P} \left( \lim_{n \rightarrow \infty}{\overline{X}_n} \neq \mu \right) = 0
$$ In other words, the probability of a sample outcome in which there is not an $m$ large enough that $\overline{X}_n$ is within $\epsilon$ of $\mu$ whenever $n>m$, is equal to zero. This is a stronger statement than the Weak Law in that it implies the Weak Law.

For the $\pi$ simulation, a single $X_i$ is the outcome of the experiment of choosing a random sample point within the square $\left[ -1, 1 \right] \times \left[ -1, 1 \right]$ and assigning $X_i = 1$ if the point is within the circle of radius 1 or $X_i = 0$ if not. Thus, the expected value of each $X_i$ is $\mu = \pi / 4$, since $$
\begin{align}
{\Bbb E}(X_i) &= 0 \cdot {\Bbb P}(X_i = 0) + 1 \cdot {\Bbb P}(X_i = 1) \\
&= 0 \cdot (1 - \pi / 4) + 1 \cdot (\pi / 4) \\
&= \pi / 4
\end{align}
$$ The Run 1, 2, and 3 lines in the chart above represent the random variable $4 \overline{X}_n$ (with $n$ on the horizontal axis), and $4X_i$ has expected value $\pi$, the dotted red line. Both Laws of Large Numbers essentially tell us that it is unlikely for the simulated values to be too far from $\pi$ when $n$ grows large enough.

More precisely, take $\epsilon = 0.0001$ as an example error threshold and 99.9% as an example probability threshold. The Weak Law tells us that there exists some $m$ large enough that, if we ran a new Run 4 simulation, there would be a probability of at least 99.9% that its line would remain within a band 0.0001 above or below the dotted red line everywhere to the right of $n=m$.

The Strong Law tells us that there is a probability of exactly 100% that, for some large enough $m$, the Run 4 line would remain within the above-mentioned error band to the right of $n=m$. Thus, the Strong Law clearly implies the Weak Law, but while there is a subtle difference in interpretation, they both essentially say that, in some sense, $\overline{X}_n \rightarrow \mu$ as $n \rightarrow \infty$.

Note: for a proof of the Law of Large Numbers using inequalities about the likely "spread" of a random variable around its mean, see this post.

What the Law of Large Numbers does not tell us is how large this "$m$" value needs to be. This is where the Central Limit Theorem comes in.

The Central Limit Theorem

We know that the expected value of a sample mean $\overline{X}_n$ is the population mean $\mu$, since $$
{\Bbb E}(\overline{X}_n)
= {\Bbb E}\left( \frac{\sum_{i=1}^{n}{X_i}}{n} \right)
= \frac{\sum_{i=1}^{n}{{\Bbb E}(X_i)}}{n}
= \frac{n \mu}{n} = \mu
$$ The Law of Large numbers tells us that, with a very high probability, $\overline{X}_n$ will indeed become close to its expected value for large enough values of $n$. The Central Limit Theorem goes a step further and actually gives us the full probability distribution of $\overline{X}_n$, known as the sampling distribution of the mean.

Central Limit Theorem: If $X_1, X_2, X_3, \dotsc$ are iid with mean $\mu$ and finite variance $\sigma^2$, then, as $n \rightarrow \infty$, the random variable $\overline{X}_n$ converges in distribution to a Normal random variable with mean $\mu$ and variance $\frac{\sigma^2}{n}$.

In this context, convergence in distribution means that the sampling distribution of $\overline{X}_n$ converges pointwise to the applicable Normal distribution.

What's amazing about the Central Limit Theorem is that the limiting distribution of $\overline{X}_n$ is Normal regardless of the distribution of the $X_i$'s themselves. In the $\pi$ example, the distribution of $X_i$ is obviously not Normal: it's the discrete distribution $$
p(x) = \begin{cases} 1 - \pi / 4 & \text{if $x=0$} \\ \pi / 4 & \text{if $x=1$}
\end{cases}
$$ Nonetheless, the CLT tells us that the sampling distribution of $\overline{X}_n$ will still be approximately Normal for large enough values of $n$. In practice, the rule of thumb is that the Normal approximation is valid when $n$ is at least 30; this "For Dummies" article provides a bit more color on this topic. The below image from Wikipedia shows the sampling distribution of the mean for differing sizes of $n$, using the example of coin flips (the $X_i$'s here are the same as the $\pi$ approximation except that they have "success" probability $1/2$ instead of $\pi/4$):

Sampling distribution of the mean $\frac{1}{n}\sum_{i=1}^{n}{X_i}$ for various values of $n$, where each $X_i$ is 1 for heads and 0 for tails. Notice that the CLT's Normal approximation seems to be most accurate starting around $n=30$. (image courtesy of Wikipedia)

We can use the CLT-provided Normal approximation to determine how many simulations we need in order to have at least, say, 95% probability of obtaining a $\pi$ estimate within, say, 0.01 of the true value of 3.1415[...]. In order to do this, we need to find the variance of the random variable $4 X_i$. This is given by the computation $$
\begin{align}
\sigma^2_{4 X_i} &= 4^2 \cdot \sigma^2_{X_i} = 16 \cdot {\Bbb E}\left[ (X_i - \mu_{X_i})^2 \right] \\[2mm]

&= 16 \bigg[
\underbrace{\left( 1 - \frac{\pi}{4} \right)}_{p(0)} \left(0 - \frac{\pi}{4} \right)^2
+ \underbrace{\frac{\pi}{4}}_{p(1)} \left( 1 - \frac{\pi}{4} \right)^2
\bigg] \\[2mm]

&= 16 \left[
\left(\frac{\pi}{4} \right)^2 - \left( \frac{\pi}{4} \right)^3 + \frac{\pi}{4} - 2 \left( \frac{\pi}{4} \right)^2 + \left( \frac{\pi}{4} \right)^3
\right] \\[2mm]

&= 16 \left[
\frac{\pi}{4} - \left( \frac{\pi}{4} \right)^2
\right] \\[2mm]

&= 4 \pi \left[ 1 - \frac{\pi}{4} \right] \\[2mm]

&= 4 \pi - \pi^2
\end{align}
$$ The below Excel screenshot shows the math needed to obtain a suitable value of $N$, which turns out to be just over the value of 100,000 which I used for my simulations.

$\pi$ approximation source code (C++)

Below is the full source code for my program, as well as the result and the first few lines of the output file "Program Output.txt". Underneath, I have added some commentary about what the different sections do.

C++ program output

"Program Output.txt" output

My notes on the program:

• Lines 1-7: C++ standard libraries included for various purposes
• Lines 9-20: definition of the "point" class: in the "main" program, the class is instantiated to create a new "point" object $p$, which has members called $x$ and $y$. Since these members are declared using the "private" keyword in the class definition, the "main" program cannot access them directly. Instead, it will populate their values using the constructor (e.g. when they are initiated with 0's in the variable declaration in line 43) or the "assign_vals" function (in line 51). The "len" function returns a number representing the "radius", $\sqrt{x^2+y^2}$, of the point.
• Lines 22-28: "seed" the random number generator using the current clock time, which ensures we will get different random numbers each time we run the program. I also define the "rand_num()" function to generate random numbers from the uniform distribution on $[-1,1] \times [-1,1]$.
• "main" program, lines 37-43: variable declarations
• "main" program, lines 48-66: looping from $j=1$ to $N$, populate the point $p$ with random values of $x$ and $y$, and if the radius $r<1$, populate the $j$-th entry of the "results" vector with a 1. If $r \geq 1$, then "results[j]" is populated with a 0 instead.
• "main" program, line 69: the $\pi$ approximation is calculated as 4 times the sum of the "results" vector (i.e. the number of "successes") divided by $N$. Since both $N$ and "sum" are integer variables, "static_cast<double>" is included to ensure that the quotient is calculated as a decimal number rather than being rounded off to the nearest integer.

Preliminary: How do mathematicians model randomness?

In this two-part post, I will introduce Monte Carlo methods, an umbrella term for techniques involving many random simulations of a phenomenon of interest to obtain a numerical result (or range of results). This post is a response (albeit a late one- my apologies for that) to a reader request from Anonymous dated August 2016.

In Part 1, I will explain at a high level what a Monte Carlo simulation is, what kinds of typical inputs and outputs we may expect, and the benefits and limitations of Monte Carlo methods. I use the example of retirement/investment planning software, which Anonymous also mentioned in his request.

In Part 2, I will walk through a considerably simpler example in detail: a Monte Carlo method to approximate the value of $\pi$. I will also explain two well known statistical results, the Law of Large Numbers and the Central Limit Theorem, which justify the use of Monte Carlo methods and allow us to quantify their accuracy. The $\pi$ simulation will directly illustrate the LLN and CLT.

Monte Carlo methods - high-level overview

In the most general terms, a Monte Carlo simulation is a method whereby we simulate one or many random inputs many times (usually thousands or millions), making assumptions as to their probability distributions, and combine the results in a deterministic (i.e. non-random) way to obtain a numerical output. In this section, I will explain in layman's terms why and how we may do this as well as some limitations of these methods.

When to employ a Monte Carlo method

In the reader request, Anonymous mentioned a typical situation in which an analyst may employ a Monte Carlo method: a wealth manager wishes to build a portfolio of various assets for a client with the goal of providing sufficient income for the client after his retirement. More precisely, for a prospective portfolio to be considered acceptable, it must have a sufficiently high probability of providing enough cash to cover the client's estimated living costs at all times after his retirement date.

Suppose we would like to assess the viability of a proposed portfolio consisting of a number of equities (stocks) and fixed-income assets (bonds). We will need to make static assumptions about the client's age, retirement age, initial portfolio size, and post-retirement income and living costs. We will also need to make certain assumptions regarding the probability distributions of interest rates, equity returns, debt yields, dividend rates, inflation, mortality, etc. Note that we may re-categorize some of the "static" assumptions as variable and vice-versa, depending on the goals of the analysis.

In any case, with so many variable inputs, it is no simple task to determine the probability that the proposed portfolio will be acceptable at any given time, let alone at all times. The number of random inputs could be enormous, and we do not necessarily have a tractable way of combining them all to arrive at a probability distribution for the portfolio value at some time $t$, so we can't solve this out analytically. However, we do have clear assumptions for all the inputs, and we have computers.

Instead of an analytical solution, we can use a Monte Carlo simulation to arrive at a numerical solution.

Typical inputs, outputs, and interpretation

I mentioned before that a Monte Carlo simulation consists of many "runs" of a deterministic procedure based on randomized inputs.

To perform a single run of a Monte Carlo simulation, we program a computer to simulate each input by drawing random numbers from an assumed probability distribution. In the retirement example, we would randomly generate a sample path over time for each relevant interest rate, stock price, bond price, etc. Given these paths, we can (deterministically!) compute the value of the portfolio at each time, and then it is a matter of arithmetic to determine whether the portfolio covered the living costs at each time or failed. This run represents one of the many possible scenarios which could occur.

We repeat this process for, say, 100,000 runs (in Part 2, I will elaborate further on how many runs we need to use), calculating the simulated probability of the portfolio's success as the number of successes divided by 100,000. If that probability exceeds a predetermined threshold (e.g. 90%), then we deem the portfolio "acceptable".

A typical software would also provide more sophisticated chart outputs such as the below (which I conveniently found on Google). The first image shows median portfolio values at each time based on certain fixed input parameters, with the median taken over all the simulation runs.

The second image shows a "heat map" of success probabilities based on different values of the input parameters. Each pixel on the heat map summarizes a full probability distribution (presumably of portfolio shortfall) for a fixed set of values of the input parameters: a green pixel indicates that the input parameter values (the $x$- and $y$-axis values) lead to a portfolio which is likely to cover costs in each year after retirement, while a red pixel indicates that the parameter values lead to a portfolio which is more likely to fail, i.e. to not cover costs (even with previous years' excesses) in at least one year after retirement.



The software must make assumptions about the statistical properties of the portfolio in order to generate chart outputs like the above. With respect to the probabilistic inputs (these do not include the deterministic inputs such as retirement age, tax rates, initial portfolio size, etc.), the images show only a mean and standard deviation of the investment returns. Therefore, this software likely assumes Normally distributed returns.

It's also worth noting that, for example, Life Expectancy seems to be an input here with a value of 95. This suggests that the simulations may not be randomizing that input. In interpreting this output, it would be important to read the software documentation and understand the major assumptions employed.


Benefits and limitations of Monte Carlo methods

The most obvious benefit of a Monte Carlo simulation is that it allows us to run millions of complex scenarios in just seconds or minutes, a task which would be impossible without a computer. Another key benefit is a bit more subtle.

If we didn't have computers available, we may attempt to answer our retirement question using "what-if scenarios": we would choose a few sets of very simple assumptions about equity and debt returns which allow us to use "back-of-the-envelope" arithmetic to calculate the portfolio value. For example, we may assume fixed, non-random annual returns with three cases as follows:


We may even go a step further and assume probability distributions for these returns (e.g. Normal distributions) which allow us to arrive at an analytical solution for the portfolio value's probability distribution over time. This is certainly a reasonable approach; however, these what-if scenarios give us no indication of how likely the three different cases are. A Monte Carlo simulation, on the other hand, gives us a distribution of possible outcomes and the probabilities of those outcomes. Unlikely scenarios will not be over-represented among the simulation runs.

Of course, no method is without limitations. On the practical side, without readymade software, we typically need to write our own code to perform a Monte Carlo simulation; this takes time and effort which may not be necessary. For very complex simulations, we may need to worry about the integrity of the random number generator(s) and, in some cases, the speed of the code, though these are beyond the scope of this post.

From a more qualitative standpoint, we have the "garbage in, garbage out" principle: a Monte Carlo simulation, like any analysis, is based on certain input assumptions. We must always keep in mind which assumptions went into a simulation when interpreting its results. These assumptions include both input probability distributions and the input parameter values used to calibrate those distributions.

For example, the chart outputs above probably assumed "normal" market conditions and Normally distributed returns. The above simulation may therefore underestimate the probability of portfolio failure due to a financial crisis (a non-normal and certainly non-Normal event). Similarly, if life expectancy was indeed a static input, then the analysis may underestimate the probability of failure for people who are very healthy and thus likely to live longer than average.

Finally, we must determine "reasonable" values of the input parameters. More often than not, we base these on historical data as our best guess. This implies that the simulation will weight possible scenarios based on their historical relative frequencies, while future frequencies may differ. On the other hand, assuming different parameter values based on our own informed estimates clearly introduces a different bias into the model.

In conclusion, Monte Carlo simulation is a powerful tool, but it is still just that: a tool. We must remain critical of our analysis and the assumptions that underlie it.

Stay tuned for Part 2, in which we will walk through a simple Monte Carlo method to estimate $\pi$ and introduce the Law of Large Numbers and Central Limit Theorem.

In 1827, botanist Robert Brown observed through his microscope that pollen particles in resting water exhibited motion which he was unable to explain. This motion came to be known as Brownian motion, which now also has a precise mathematical definition as a stochastic (i.e. random) process.

Fast forward to 1905: the scientific community had not yet accepted, as we do today, that matter consists of atoms. Though analysis based on atoms and molecules had proven useful in some applications, many physicists regarded atoms as a hypothetical tool rather than a physical reality. Alternative theories postulated that the structure of the world consisted, for example, of different forms of energy and energy transformations. Albert Einstein's 1905 paper "On the movement of particles suspended in stationary liquids required by the molecular-kinetic theory of heat" (kind of a mouthful, but hey, the guy was German, so what did you expect?) explained Brownian motion as arising from the random bombardments of the pollen particles by moving water molecules, like in the animation below:

Einstein was able to derive formulas for the mean square displacement of a suspended particle as well as Avogadro's number, and when Jean Perrin experimentally verified these results in 1908, the remaining skeptics accepted this as evidence of the physical reality of atoms.

In this post, I will walk you through Einstein's pivotal analysis of Brownian motion.

Note: nothing in this post is new- it's all copied from various sources, including Einstein's paper, Robert E. Kennedy's book on the same (as well as other Einstein papers), and other sources from around the internet. I found it very difficult to find any single source that provided a satisfactory explanation of every step in Einstein's paper, so I hope this post, which has taken me quite a while, achieves that and is clear enough to follow easily without the need to scour the internet for additional information. If not, please post any questions in the comments section!

The setup

We have a container with volume $V$ filled with liquid (the solvent- let's just call it water going forward) and a solute (e.g. sugar) dissolved therein. The container is divided into two sections by a wall through which the water can flow freely, but which is impermeable to the sugar, and the latter is confined to the walled-off section of the container, which has volume $V^*$. The volume of the remainder of the container is then $V-V^*$.

Since the sugar molecules cannot pass through the wall, they bounce back when they hit it and thus exert a pressure on the wall (the water passes right through the wall and thus does not exert a pressure on it). Van 't Hoff's law states that if the number of moles of sugar in the volume $V^*$ is $z$, then this pressure $\Pi$, known as the osmotic pressure, is given by $\Pi V^* = zRT$, where $R$ is the gas constant and $T$ is the temperature in Kelvin. This is the same as the ideal gas law, which is not totally unexpected since the pressure is due to the sugar molecules bouncing around in the volume $V^*$ the same way a gas in a closed container would.

If, instead of sugar, some larger particles are suspended (but not dissolved, the difference lying solely in the particle size) in the water, according to Einstein, the osmotic pressure on the wall should be given by the same formula as that for solutions above; in other words, the osmotic pressure only depends on the number (or more precisely, the concentration, or number per unit volume) of dissolved/suspended particles. The next section shows his derivation based on the "molecular-kinetic theory of heat."

The osmotic pressure due to suspended particles

In the last post, I showed you Einstein's derivation of a formula for the entropy of $n$ particles in the water bouncing around due to random molecular motion and collisions: $$
S = \frac{\bar{E}}{T} + k \ln \left[ \int e^{\frac{-E}{kT}}dp_1 \, dp_2 \, ... \, dp_{3n} \, dq_1 \, dq_2 \, ... \, dq_{3n} \right] \tag{1}
$$ Here, $\bar E$ is the total energy of the $n$ particles, $T$ is the absolute temperature in Kelvin, $k$ is the Boltzmann constant (equal to the gas constant $R$ divided by Avogadro's number $N$), and the $p_i$'s and $q_i$'s are the state variables of the system of $n$ particles, representing the 3 components of momentum and position of each particle, hence the indices' running up to $3n$. The integral is taken over the different configurations of the system (i.e. different possible values of position and momentum of the particles). The $E$ inside the integral is the energy of the system as a function of the state variables.

The free energy (Helmholtz free energy) of the system is defined as $F = E - TS \ \ \text{(2)}$. If you read the Wikipedia article on Helmholtz free energy, you can see the derivation of the formula $dF = -S \, dT - P \, dV$ via basic calculus (in particular, the product rule for differentiation). From this formula, it's clear that $-\dfrac{\partial F}{\partial V} = P$ when the volume increases by an infinitesimal amount at constant pressure. We'll use this relation later on.

Plugging in equation (1) for $S$ in equation (2), we obtain $$
\begin{align}
F &= -kT \ln \left[ \int e^{\frac{-E}{kT}}dp_1 \, dp_2 \, ... \, dp_{3n} \, dq_1 \, dq_2 \, ... \, dq_{3n} \right] \\[3mm]
&= -kT \ln B \tag{3}
\end{align}
$$ The key insight that makes the nearly-impossible calculation of the integral $B$ unnecessary is that $B$ can be shown to have the form $B = J(V^{*})^n \ \ \text{(4)}$, where $J$ is some function that does not depend on $V^*$. I'll go through the details of that derivation at the end of this post so as not to deviate too far from the more interesting line of reasoning.

Combining equations (3) and (4), we recover van 't Hoff's law for the osmotic pressure $\Pi$: $$
\begin{align}
\Pi &= -\dfrac{\partial F}{\partial V^*} \\[3mm]
&= \dfrac{\partial}{\partial V^*} kT \ln [ J(V^{*})^n ] \\[3mm]
&= \dfrac{\partial}{\partial V^*} kT \, [\ln J + n \ln V^{*} ] \\[3mm]
&= kT\dfrac{n}{V^*} \\[3mm]
&\Downarrow \\[3mm]
\Pi V^* &= (Nk)T\dfrac{n}{N} \\[3mm]
\Pi V^* &= zRT
\end{align}
$$ In the last line, we used the fact that $k = \dfrac{R}{N}$ and denote the number of moles of particles by $z$. We have shown that van 't Hoff's law is a consequence of the molecular-kinetic theory. This law can equivalently be expressed as $\Pi = kT \nu \tag{$\spadesuit$}$ where $\nu = \dfrac{n}{V^*}$ is the concentration of particles in the volume $V^*$.

The diffusion coefficient

In deriving a formula for the diffusion coefficient (i.e. the constant in the differential equation governing the diffusion of the particles over time- if you aren't familiar with the diffusion equation, you'll see presently...),  Einstein introduced a fictitious force and used the fact that the free energy is minimized when the system is in dynamic equilibrium. I couldn't fully understand this argument (if you do, please click here and answer my questions!), so instead, I will show a slightly more direct derivation by letting the force due to osmotic pressure be balanced by friction in the fluid.

Suppose again we have $n$ particles in the volume $V^*$, which has a cross-sectional area $A$ perpendicular to the $x$-axis, and that the concentration of the particles $\nu (x)$ may vary with $x$. Consider the volume $\Delta V^*$ between $x$ and $x+ \Delta x$, which contains $N_1$ particles. We have: $$
\frac{\Pi (x) - \Pi (x+ \Delta x)}{\Delta x} = \frac{F(x) - F(x+ \Delta x)}{A \Delta x} = \frac{N_{1} F_{\Pi}}{A \Delta x}
$$ where $\Pi$ is the osmotic pressure. $F(x)$ is the force on the cross-sectional area at $x$ (positive $F$ would be to the right), and $F_{\Pi}$ is the average force on a particle in the region $\Delta V^*$. Taking the limit as $\Delta x \rightarrow 0$ gives: $$
\frac{\partial \Pi}{\partial x} = - \nu F_{\Pi} \tag{5}
$$ Note that both sides are functions of $x$ and time, but I am leaving the $x$'s and $t$'s out to simplify the notation a bit. The $\frac{N_1}{A \Delta x} = \frac{N_1}{\Delta V^*}$ became $\nu$ in the limit since that was the number of particles per unit volume.

The force $F_{\Pi}$ acts on the particles and "competes" with the force of friction in the fluid. The drag force due to friction is proportional to the particle velocity since a faster-moving particle hits more molecules that slow it down; if the particles are spheres, the drag force is given by Stokes' Law: $$
F_d = -6 \pi \mu R_{p} v \tag{6}
$$ where $\mu$ is the viscosity coefficient of the fluid (measurable by simple experiments), $R_p$ is the particle radius, and $v$ is the particle velocity. This lower-case $\pi$ is the usual 3.14. The negative sign reflects the fact that the drag force points in the opposite direction of the velocity. I'm not going to present a derivation of Stokes' Law since it would be very long, so let's take that as a given and plug on.

Before combining all this, we need one more tidbit. Taking  the partial derivative with respect to $x$ on both sides of equation $(\spadesuit)$ from above, we obtain: $$
\frac{\partial \Pi}{\partial x} = kT \frac{\partial \nu}{\partial x}
$$ In dynamic equilibrium, $F_{\Pi} + F_{d} = 0$, so combining (5) and (6) with the above gives: $$
\begin{align}
- \frac{1}{\nu} \frac{\partial \Pi}{\partial x} &= 6 \pi \mu R_{p} v \\[3mm]
\implies -kT \frac{\partial \nu}{\partial x} &= 6 \pi \mu R_{p} \nu v \tag{7}
\end{align}
$$ Now $\nu v$, the density times the velocity, is the flux of particles (at position $x$ and time $t$), i.e. the number of particles passing through a unit area perpendicular to the $x$-axis, per second. Notice that the units are $\dfrac{\text{particles}}{\text{m}^3} \times \dfrac{\text{m}}{\text{s}} = \dfrac{\text{particles}}{\text{m}^2 \ \text{s}}$, as one would expect of such a quantity.

Fick's First Law states that this flux equals $-D \frac{\partial \nu}{\partial x}$, where $D$ is the diffusion coefficient referred to above. The negative sign means that if the particle density is higher on the left, particles tend to diffuse to the right. $D$ has units of $\frac{\text{m}^2}{\text{s}}$ and measures the mean squared displacement of particle diffusion per unit time. This Stack Exchange thread has a detailed explanation of the physical interpretation of the diffusion coefficient, and this Wikipedia article presents a quick and simple derivation of Fick's First Law.

Combining Fick's First Law with equation (7), we see that $$
D = \frac{kT}{6 \pi \mu R_{p}} = \frac{RT}{N} \frac{1}{6 \pi \mu R_{p}} \tag{8}
$$

Root mean squared displacement

Assume the $n$ particles' movements are independent and introduce a time interval $\tau$. This time interval is much shorter than the time intervals of observation but long enough that we can consider the movements of a single particle in consecutive time intervals of length $\tau$ to be independent.

Since the movements of the $n$ particles are independent of each other, we can think of them as $n$ different stochastic processes. In other words, we can consider them $n$ separate observations of results of the same random experiment. In a time interval $\tau$, the $x$-coordinate of the position of any given particle will change by some amount $\Delta$. If a particle moves to the right, $\Delta > 0$, and if it moves to the left, $\Delta < 0$. The values of $\Delta$ follow some probability distribution $\phi (\Delta)$ so that after a time period $\tau$ elapses, the proportion of particles $\frac{dn}{n}$ which have experienced a displacement in the $x$ direction between $\Delta$ and $\Delta + d \Delta$ satisfies the equation $$
\frac{dn}{n} = \phi(\Delta) \, d \Delta
$$ There are a few conditions that the probability distribution $\phi$ should satisfy. First, all probabilities should add up to 1, so $$
\int \limits_{-\infty}^{\infty}{\phi(\Delta) \, d \Delta} = 1
$$ Also, since $\tau$ is small, $\phi( \Delta)$ should be zero except for very small absolute values of $\Delta$. Finally, a particle should be equally likely to have been displaced to the left or right, so $\phi$ should be an even function, i.e. $$
\phi(\Delta) = \phi( - \Delta)
$$ If we denote the particle density by $\nu = f(x,t)$, then by the definition of $\phi$, we have: $$
f(x,t + \tau) = \int \limits_{-\infty}^{\infty}{f(x+ \Delta , t) \phi(\Delta) \, d \Delta} \tag{9}
$$ This equation states that the number of particles at position $x$ at time $t + \tau$ is the sum (integral) of the numbers of particles that were at positions $x + \Delta$ at time $t$ and were displaced by $-\Delta$ in the time $\tau$ (recall that $\phi(-\Delta) = \phi(\Delta)$), summed over the possible values of $\Delta$. Since $\phi(\Delta)$ is assumed to be non-zero only for very small values of $\Delta$, extending the integral out from $-\infty$ to $\infty$ doesn't change the value, but it will be useful below since it will allow us to use well known results for Gaussian integrals.

Now we can expand $f$ into a Taylor series on both sides of equation (9). On the left, we have $$
\begin{align}
f(x, t + \tau) &= f(x,t) + \tau \frac{\partial f}{\partial t} + \frac{1}{2} \tau^{2} \frac{\partial^{2} f}{\partial t^2} + ... \\[3mm]
&\approx f(x,t) + \tau \frac{\partial f}{\partial t}
\end{align}
$$ On the right side, we need to expand out to second order (you'll see why in a second)- note that it is justified to do the Taylor expansion inside the integral since we are expanding around $\Delta = 0$, and only small values of $\Delta$ (i.e. those for which the Taylor expansion is accurate) contribute anything to the integral: $$
\begin{align}
f(x,t) + \tau \frac{\partial f}{\partial t} &\approx
\int \limits_{-\infty}^{\infty}{\left[ f(x,t) + \Delta \frac{\partial f}{\partial x} + \frac{1}{2} \Delta^{2} \frac{\partial^{2} f}{\partial x^2} \right] \phi(\Delta) \, d \Delta} \\[3mm]

&= f(x,t) \int \limits_{-\infty}^{\infty}{\phi(\Delta) \, d \Delta}
+ \frac{\partial f}{\partial x} \int \limits_{-\infty}^{\infty}{\Delta \phi(\Delta) \, d \Delta}
+ \frac{1}{2} \frac{\partial^{2} f}{\partial x^2} \int \limits_{-\infty}^{\infty}{\Delta^{2} \phi(\Delta) \, d \Delta}
\end{align}
$$ The first integral on the right-hand side is equal to 1, and the second is zero since $\phi(\Delta)$ is an even function (so $\Delta \phi(\Delta)$ is an odd function). So the above simplifies to $$
\frac{\partial f}{\partial t} = \frac{1}{2 \tau} \frac{\partial^{2} f}{\partial x^2} \int \limits_{-\infty}^{\infty}{\Delta^{2} \phi(\Delta) \, d \Delta}
$$ For any choice of the distribution $\phi$, the quantity $\frac{1}{2 \tau} \int_{-\infty}^{\infty}{\Delta^{2} \phi(\Delta) \, d \Delta}$ will be a constant, which we call $D$. It follows that the particle density $f$ satisfies the well known diffusion equation (or Fick's Second Law): $$
\frac{\partial f}{\partial t} = D \frac{\partial^{2} f}{\partial x^2} \tag{10}
$$ This $D$ is the same one referred to in Fick's First Law from above (in fact, the same Wikipedia article shows a derivation of the diffusion equation from Fick's First Law). This means that the formula derived above for the diffusion coefficient is valid here, and equation (10) completely describes the evolution of the particle density.

If we consider a separate coordinate system for each particle whose origin is at the particle's position at $t=0$, then equation (10) still holds, except now $f(x,t)$, instead of describing number of particles at position $x$ at time $t$, would describe the number of particles experiencing a displacement of $x$ from their respective initial positions in an elapsed time $t$. Given the conditions $f(x,0)=0$ for $x \neq 0$ and $\int \limits_{-\infty}^{\infty}{f(x,t) \, dx} = n$, the differential equation (10) is that of diffusion from a point, and its solution is: $$
f(x,t) = \frac{n}{\sqrt{4 \pi D t}} \exp \left(\frac{-x^2}{4Dt} \right)
$$ Since $f(x,t)$ represents a number of particles, $\frac{1}{n} f(x,t)$ is the probability that a given particle experiences a displacement $x$ in time $t$. Indeed, using the Gaussian integral formula $\int_{-\infty}^{\infty}{e^{-ax^2}dx} = \sqrt{\frac{\pi}{a}}$ with $a=4Dt$, we see that $\int_{-\infty}^{\infty}{\frac{1}{n}f(x,t) \, dx}=1$. It follows that the mean value of the squared displacement is: $$
\begin{align}
\left< x^2 \right> &= \int \limits_{-\infty}^{\infty}{x^{2} \frac{1}{n} f(x,t) \, dx} \\
&= \frac{1}{\sqrt{4 \pi Dt}}\int \limits_{-\infty}^{\infty}{x^{2} \exp \left( \frac{-x^2}{4Dt} \right) \, dx} \\[3mm]
&= 2Dt
\end{align}
$$ where the last step utilized the Gaussian integral formula $\int_{-\infty}^{\infty}{x^{2} e^{-ax^2} \, dx} = \frac{1}{2} \sqrt{\frac{\pi}{a^3}}$. Finally, it we see that the root mean sqaured displacement, a measure of the average distance a particle is expected to travel in time $t$, is $\sqrt{\left<  x^2 \right>} = \sqrt{2Dt}$, which is proportional to the square root of $t$.

This result was verified by experiment not long after Einstein published his paper, and this is often considered the defining result that convinced the remaining skeptics at the time that atoms existed.

I hope you enjoyed reading and please post any questions in the comments section.

For those interested, I promised an appendix for one of the details that was glossed over above:

Appendix: derivation of $B = J (V^*)^n$

Consider the $n$ particles at positions $(x_1 , y_1, z_1 ), ... , (x_n , y_n, z_n )$ and each surrounded by an infinitesimal parallelpiped region $dx_i \, dy_i \, dz_i$ which is contained in $V^*$. Consider the same $n$ particles, but now at different positions $(x_i ', y_i ', z_i ')$ and surrounded by parallelpiped regions of the same size as before (thus no label change necessary on those). Writing $$
dB  = J \, dx_1 \, dy_1 \, ... \, dz_n
$$ and $$
dB'  = J' \, dx_1 \, dy_1 \, ... \, dz_n
$$ we see that $\dfrac{dB}{dB'} = \dfrac{J}{J'}$.

The probability that the system is in the first configuration is $dB$ divided by the integral over all configurations, $B$, and likewise, the probability the system is in the second configuration is $\frac{dB'}{B}$. Since the particles move independently and the parallelpiped regions are the same size, the two probabilities must be equal, and so $$
\frac{dB}{B} = \frac{dB'}{B} \implies dB = dB' \implies J = J'
$$ This shows that $J$ is independent of the particle positions and of the volume $V^*$, and so $$
\begin{align}
B = \int dB &= \int \limits_{V^*}{J \, dx_1 \, dy_1 \, ... \, dz_n} \\[3mm]
&= J \int \limits_{V^*}{dx_1 \, dy_1 \, ... \, dz_n} \\[3mm]
&= J (V^*)^n \ \tag*{$\square$}
\end{align}
$$

Entropy is often considered a measure of "disorder" which, as you may recall from chemistry, is supposed to increase over time. A physical system tends to evolve in such a way that its useful energy dissipates. The value of entropy measures how far along the system is in that process. There are a few different formulations which each capture this idea; each is a formula that must be calculated as opposed to a physically observable property of the system (like temperature or pressure would be).

In this post, I'm going to show you Einstein's derivation of a formula for entropy, which will also shed some light on what exactly this quantity represents. This formula will be an ingredient in the forthcoming post on Einstein's Brownian motion paper, so pay attention!

Imagine a system, consisting of $n$ atoms of an ideal gas in a closed container ($n$ will be a big number, like on the order of $10^{23}$ big). Actually, it doesn't need to be a gas, but that just seems to be the easiest to picture. An ideal gas means that the molecules are monatomic, and we can ignore rotational energy of the atoms, as well as interactions between them. In order words, we are concerned only with their translational motion and assume all collisions are elastic. Finally, we assume the container of gas is surrounded by an ambient reservoir, with which it can exchange heat, so that the system's temperature $T$ remains pretty much constant.



Each atom has a position and a momentum, each a 3-dimensional vector, i.e. a vector with 3 components. At any given time, the $2 \times 3n = 6n$ components of position and momentum, denoted $q_1, q_2, ..., q_{3n}$ and $p_1, p_2, ..., p_{3n}$ respectively, for all the atoms determine the configuration of the system, and the set of all possible configurations is called configuration space, a subset of $\Bbb R ^ {6n}$. The $q_i$'s and $p_i$'s are called the state variables of the system.

The First Law of Thermodynamics states that energy can be converted into different forms but not created or destroyed, and thus as the system evolves, its change in energy is the work done on the system plus the heat supplied to the system: $$
dE = dW + dQ
$$ If the system is described by a set of parameters $\lambda_1, \lambda_2, ... , \lambda_m$ and the state variables mentioned above, then if it undergoes a small change over a time interval $dt$, the resulting change in energy is given by: $$
dE = \sum_{i=1}^{m}{\frac{\partial E}{\partial \lambda_i}\frac{d\lambda_i}{dt}dt}
+ \sum_{j=1}^{3n}{\left( \frac{\partial E}{\partial q_j}\frac{dq_j}{dt}dt + \frac{\partial E}{\partial p_j}\frac{dp_j}{dt}dt \right)} \tag{$\spadesuit$}
$$ To be a bit more concrete, the "parameters" above would be things like volume of the container, so the first sum is identified with the work done on the system (usually just $dW = P \ dV$), and the second sum is identified with the heat supplied to the system; more heat in the system (all else equal) means higher temperature, i.e. the atoms bounce around faster, and so the second sum is our $dQ$. In the example we're talking about, the energy is only due to translational kinetic energy of the atoms, and thus the energy would be $E=\sum_{j=1}^{3n}{\frac{p_i^2}{2m}}$. There would also be a term involving the $q_i$'s if we took into account gravitational potential energy, which depends on the particles' positions.

The above formula doesn't depend on which path the system takes through state space, so in particular, it holds for an adiabatic change, i.e. one in which the system does not gain/lose any heat so $dQ=0$. Before such a change occurs, the probability of finding the system in a state with energy $E$ is given by the volume of a little box in state space times its probability density: $$
d{\Bbb P} = Ce^{\frac{-E}{kT}}dp_1 \, dp_2 \, ... \, dp_{3n} \, dq_1 \, dq_2 \, ... \, dq_{3n}
$$ The probability density $Ce^{\frac{-E}{kT}}$ deserves a bit of attention. $T$ is the temperature of the system, considered to be constant as mentioned above, $k$ is the Boltzmann constant, $1.38 \times 10^{-23}$ Joules/Kelvin (units of energy per temperature), which makes the exponent dimensionless, and $C$ is a constant which makes all the probabilities add up to 1. We'll solve for $C$ in a moment, but why is the probability density given by an exponential?

If we have two pockets of gas in the container with energies $\epsilon_1$ and $\epsilon_2$, since the model is that they are independent, the probability of finding pocket 1 at energy $\epsilon_1$ and finding pocket 2 at energy $\epsilon_2$ should be the product of the individual probabilities of finding the respective pockets at those energy levels. Furthermore, the energies add, so that the total energy of the two pockets combined is $\epsilon_1 + \epsilon_2$. The exponential function has the desired property: $$
\exp \left( \frac{-(\epsilon_1 + \epsilon_2)}{kT} \right) = \exp \left( \frac{-\epsilon_1}{kT} \right) \exp \left( \frac{-\epsilon_2}{kT} \right)
$$ The fact that the "drop-off factor" in the denominator of the exponent is proportional to $T$ is a consequence of the fact that the average kinetic energy an atom in the gas is $\frac{3}{2}kT$, which in turn follows from the ideal gas law $PV = NkT$. For a simple and insightful derivation of the $e^{\frac{-E}{kT}}$ based on Maxwell's analysis, click here.

Back to the main line: in order to solve for $C$, we note that the probabilities of all possible configurations must equal 1. Probabilities are always non-negative, so we can assume that $C$ is of the form $e^c$, and thus: $$
\begin{align}
1 &= \int{d{\Bbb P}} \\[3mm]
&= \int e^c e^{\frac{-E}{kT}}dp_1 \, dp_2 \, ... \, dp_{3n} \, dq_1 \, dq_2 \, ... \, dq_{3n} \\[3mm]
&= e^c \int e^{\frac{-E}{kT}}dp_1 \, dp_2 \, ... \, dp_{3n} \, dq_1 \, dq_2 \, ... \, dq_{3n} \\[3mm]
& \Downarrow \\[3mm]
c &= - \ln \left[ \int e^{\frac{-E}{kT}}dp_1 \, dp_2 \, ... \, dp_{3n} \, dq_1 \, dq_2 \, ... \, dq_{3n} \right]
\end{align}
$$ The integrals above are taken over the entire range of possible values of the $q_i$'s and $p_i$'s, i.e. over all state space.

After our adiabatic change in the system, there will be a similar expression for $d{\Bbb P}$, except that $c$ now may have shifted a bit from $c$ to $c + dc$, as may have $\beta$ to $\beta + d\beta$ (where we are now defining $\beta := \frac{1}{2kT}$ for notational convenience). The energy $E$ will also shift to $E + dE = E + \sum_{i=1}^{m}{\frac{\partial E}{\partial \lambda_i}\frac{d\lambda_i}{dt}dt}$. Here, we've used equation $(\spadesuit)$ and the fact that $dQ$ = 0, so only the $dW$ term comes into play. To save on the symbols, I'll also start referring to the $\frac{d\lambda_i}{dt}dt$'s simply as $d \lambda$.

Now similar to the above, we have: $$
\begin{align}
1 = &\int{d{\Bbb P}} \\[3mm]

=  &\int \exp \left( (c+dc)-2(\beta+d\beta)\left(E + \sum{\frac{\partial E}{\partial \lambda} d\lambda}\right) \right) dp_1 \, dp_2 \, ... \, dp_{3n} \, dq_1 \, dq_2 \, ... \, dq_{3n} \\[3mm]

 = &\int \exp
\left(
dc - 2
\left(
E \, d\beta + \beta \sum{\frac{\partial E}{\partial \lambda} d\lambda} +d\beta \sum{\frac{\partial E}{\partial \lambda} d\lambda}
\right)
\right) \\[1mm]
&\times \exp \left( c-\frac{E}{kT} \right)
dp_1 \, dp_2 \, ... \, dp_{3n} \, dq_1 \, dq_2 \, ... \, dq_{3n}
\end{align}

$$ We can expand the first exponential into a Taylor series and then neglect the terms past first order: $$
\begin{align}
1 = \int
&\left[
1 + dc - 2
\left(
E \, d\beta + \beta \sum{\frac{\partial E}{\partial \lambda} d\lambda} +d\beta \sum{\frac{\partial E}{\partial \lambda} d\lambda}
\right)
+ \frac{1}{2} (dc - 2 (...))^2 + ... \right] \\[1mm]

& \times \exp \left( c-\frac{E}{kT} \right)
dp_1 \, dp_2 \, ... \, dp_{3n} \, dq_1 \, dq_2 \, ... \, dq_{3n}\\[3mm]

\approx \int
&\left[
1 + dc - 2
\left(
E \, d\beta + \beta \sum{\frac{\partial E}{\partial \lambda} d\lambda} +d\beta \sum{\frac{\partial E}{\partial \lambda} d\lambda}
\right) \right] \\[1mm]

& \times \exp \left( c-\frac{E}{kT} \right)
dp_1 \, dp_2 \, ... \, dp_{3n} \, dq_1 \, dq_2 \, ... \, dq_{3n}\\[3mm]

= \int &\exp \left( c-\frac{E}{kT} \right)
dp_1 \, dp_2 \, ... \, dp_{3n} \, dq_1 \, dq_2 \, ... \, dq_{3n}\\[1mm]

+ \int & \left[
dc - 2
\left(
E \, d\beta + \beta \sum{\frac{\partial E}{\partial \lambda} d\lambda} +d\beta \sum{\frac{\partial E}{\partial \lambda} d\lambda}
\right) \right] \\[1mm]

& \times \exp \left( c-\frac{E}{kT} \right)
dp_1 \, dp_2 \, ... \, dp_{3n} \, dq_1 \, dq_2 \, ... \, dq_{3n}\\[3mm]


= \ \ \ & 1 \\[1mm]

+ \int & \left[
dc - 2
\left(
E \, d\beta + \beta \sum{\frac{\partial E}{\partial \lambda} d\lambda} +d\beta \sum{\frac{\partial E}{\partial \lambda} d\lambda}
\right) \right] \\[1mm]

& \times \exp \left( c-\frac{E}{kT} \right)
dp_1 \, dp_2 \, ... \, dp_{3n} \, dq_1 \, dq_2 \, ... \, dq_{3n}\\[3mm]


\Downarrow \\[3mm]

0  \ \approx \int & \left[
dc - 2
\left(
E \, d\beta + \beta \sum{\frac{\partial E}{\partial \lambda} d\lambda}
\right) \right]  \exp \left( c-\frac{E}{kT} \right)
dp_1 \, dp_2 \, ... \, dp_{3n} \, dq_1 \, dq_2 \, ... \, dq_{3n}\\

\end{align}

$$ Note: in the line before the $\Downarrow$, the first integrand is the probability density, so its integral over all state space must equal 1. Also, after the $\Downarrow$, we dropped the last term in the square brackets because we ignored terms above first order, i.e. terms containing products of two or more differentials.

Since $\exp \left( c-\frac{E}{kT} \right)$ is never negative, the only way the integral in the last line above can equal 0 is if the expression in square brackets equals 0. Thus:$$
dc -2E \, d\beta - 2 \beta \sum{\frac{\partial E}{\partial \lambda} d\lambda} = 0 \tag{1}
$$ On the other hand, multiplying the equation $dE = \sum{\frac{\partial E}{\partial \lambda}d\lambda} + dQ$ by $2 \beta$ and rearranging gives: $$
-2 \beta \, dE + 2 \beta \sum{\frac{\partial E}{\partial \lambda} d \lambda} + 2 \beta \, dQ = 0 \tag{2}
$$ Adding $(1)$ and $(2)$ eliminates the $\lambda$'s and yields: $$
\begin{align}
0 &= dc - 2E \, d\beta - 2\beta \, dE + 2\beta \, dQ \\[2mm]
&= dc -2 (E \, d\beta + \beta \, dE) + 2 \beta \, dQ \\[2mm]
&= dc -2 \, d(\beta E) + 2 \beta \, dQ \\[2mm]
\implies 2 \beta \, dQ &= d(2 \beta E - c) \\[3mm]
\implies \frac{dQ}{T} &= d \left( \frac{E}{T}-kc \right) := dS
\end{align}
$$ In the last step, we plugged in $\beta = \frac{1}{2kT}$ and then multiplied through by the constant $k$.

We have shown that $\frac{dQ}{T}$ is the total differential of some quantity related to energy and temperature, which we call entropy and denote by $S$. Evidently, $S$ is given by: $$
S = \frac{E}{T} + k \ln \left[ \int e^{\frac{-E}{kT}}dp_1 \, dp_2 \, ... \, dp_{3n} \, dq_1 \, dq_2 \, ... \, dq_{3n} \right]
$$ where we have used the formula for $c$ which was derived above. This entropy equation will be used in the next post on Brownian motion- stay tuned.

Thanks for reading, and please post any questions in the comments section.

Prerequisites: Vectors

In this two-part post, we'll go through some of the basic geometry of pool/billiards.

In Part 2, we'll derive simple methods to make one-rail bank and kick shots (to be defined below). We'll also go briefly into an example of a two-rail kick, with which it's extremely difficult to actually pocket the target ball, but at least you can work your way out of some sticky situations and avoid giving your opponent ball-in-hand.

This analysis also works for mini-golf on flat surfaces, by the way.

How to pocket a ball

This section involves a bit of physics, which I'll explain for those who don't know it already, but you may need to quickly read up on vectors here before proceeding.

Suppose we have the following set-up:
The white ball is the cue ball, and we want to hit it into the 2 ball (the blue one, also known as the object ball) to pocket the latter in the side pocket on the top of the image.

In order to accomplish this, we need the cue ball, upon contact with the 2, to impart a force upon the latter which makes it move in the direction of the pocket. The way to make the force be in that direction is to hit such that the point of contact of the two balls lies along the line between the center of the pocket (where we want the object ball to go), and the center of the object ball as in the next diagram:
It's a subtle difference, but please note that you are aiming for the point of contact, and not the center of the cue ball, to lie along the yellow line upon contact with the 2. The translucent cue ball in the above diagram shows the desired cue ball position upon contact.

If you get this right, and hit the cue ball at center with a reasonable speed, then the force imparted on the 2 will be along the yellow line, and thus it will roll along the yellow line and into the pocket. This is Newton's second law at work, which states that if the mass (i.e. how many kilograms) of the 2 ball is $m$, and the force the cue ball imparts on the 2 is ${\bf F}$ (note that this is a vector quantity, which is why it has a direction, while the mass is a scalar), then ${\bf F} = m {\bf a}$, i.e. the 2 will gain an acceleration ${\bf a}$ due to the force ${\bf F}$. The units of the acceleration are meters per second squared, and thus the units of the force are killograms*meters per second squared, also called Newtons after the same Isaac Newton we were just talking about.

Acceleration is the change, both of magnitude and direction, in velocity per unit time (thus change over one second, in how many meters per second the ball is traveling at that moment). Velocity is just the vector quantity whose magnitude is the speed (in units of meters per second) and whose direction is the direction of motion of the ball. These quantities can all vary over time, as can the force. The ball's mass $m$ is a scalar quantity that is constant over time and is a measure of how much matter is contained in the ball. The heavier the ball, the more force it takes to accelerate the ball by an equivalent amount. That's what the magnitude part of the vector equation ${\bf F} = m {\bf a}$ tells us. The direction part tells us that the acceleration is in the same direction as the force.

Make sense? Ok good- if there are questions on that, they can go in the comments section or maybe I can do a separate post, but my point was that since the cue ball is round and thus contacts the (also round) 2 at exactly one point, the cue ball must impart a force on the 2 i the direction of the yellow line in the diagram above, and thus the 2 will accelerate along that line after the contact. Since no forces act on the ball that would cause it to deviate off of that line after the contact, it will continue along that line and thus into the side pocket.

Where does the cue ball go after the contact?

That's an important question, and good players need to take this into account when planning a series of shots.

As it turns out, the cue ball bounces off perpendicular to the yellow line as in the next diagram:

To see why this is the case, we need to use conservation of momentum. What does this mean? Well, momentum is the vector quantity $m{\bf v}$ where $m$ and ${\bf v}$ are mass and velocity as above. To say that momentum is conserved means that the momentum vector of a system of objects (for a system of multiple objects, this would be the vector sum of the individual momenta) remains constant in the absence of a net external force. In our system of two balls, the force between the balls would not qualify as external. Gravity would, but it is counteracted by the force of the table pushing back up on the balls, which causes them to not fall to the ground. Thus, ignoring friction, energy loss due to the sound of the balls' hitting, etc., the momentum of the system of the two balls is the same right before and right after the collision.

In the diagram above, we have labeled the velocities, and let's assume the cue and 2 have the same mass $m$. Momentum is conserved before and after the collision, which means: $$m{\bf v}_0 = m{\bf v}_1 + m{\bf v}_2
$$
Note that the 2 ball had no velocity initially, so the left-hand side of the equation has only the cue ball's momentum. The $m$'s cancel out to give the vector equation $${\bf v}_0 = {\bf v}_1 + {\bf v}_2 $$ which actually comprises 2 algebraic equations, one in the $x$-component and one in the $y$-component:
$$
\begin{align}
v_{0x} &= v_{1x} + v_{2x} \tag{1}\\[2mm]
v_{0y} &= v_{1y} + v_{2y} \tag{2}
\end{align}
 $$Now, we can use $(1)$ and $(2)$ to obtain: $$
\begin{align}
\| {\bf v}_0 \|^2 &= v_{0x}^2 + v_{0y}^2 \\[2mm]
&= (v_{1x} + v_{2x})^2 + (v_{1y} + v_{2y})^2 \\[2mm]
&= \| {\bf v}_1 \|^2 + \| {\bf v}_2 \|^2 + 2v_{1x}v_{2x} + 2v_{1y}v_{2y} \\[2mm]
&= \| {\bf v}_1 \|^2 + \| {\bf v}_2 \|^2 + 2({\bf v}_1 \cdot {\bf v}_2) \tag{3}
\end{align}
$$ We also know that the energy of the system is conserved. The energy of an object of mass $m$ and speed $v$ is $\frac{1}{2}mv^2$. Technically, this is only the kinetic energy (energy due to motion of a massive particle), but there is no potential energy in this system (e.g. an object high up about to fall and gain speed, and thus kinetic energy, would have potential energy).

Conservation of energy tells us that $$
\begin{align}
&\frac{1}{2} m \| {\bf v}_0 \|^2 = \frac{1}{2} m \| {\bf v}_1 \|^2 + \frac{1}{2} m \| {\bf v}_2 \|^2 \\[2mm]
\Longrightarrow \ &\| {\bf v}_0 \|^2 = \| {\bf v}_1 \|^2 + \| {\bf v}_2 \|^2 \tag{4}
\end{align}
$$ Subtracting equation $(4)$ from equation $(3)$ shows that ${\bf v}_1 \cdot {\bf v}_2 = 0$, i.e. ${\bf v}_1$ and ${\bf v}_2$ are perpendicular. This means that the cue ball indeed bounces off at a right angle to the direction of the 2 after contact.

In part 2 of this post, we'll explore bank and kick shots...