michaelnhowes

Maximum likelihood estimation and the method of moments

Maximum likelihood and the method of moments are two ways to estimate parameters from data. In general, the two methods can differ but for one-dimensional exponential families they produce the same estimates.

Suppose that $\{P_\theta\}_{\theta \in \Omega}$ is a one-dimensional exponential family written in canonical form. That is, $\Omega \subseteq \mathbb{R}$ and there exists a reference measure $\mu$ such each distribution $P_\theta$ has a density $p_\theta$ with respect to $\mu$ and,

$p_\theta(x) = h(x)\exp\left(\theta T(x)-A(\theta)\right).$

The random variable $T(X)$ is a sufficient statistic for the model $X \sim P_\theta$ . The function $A(\theta)$ is the log-partition function for the family $\{P_\theta\}_{\theta \in \Omega}$ . The condition, $\int p_\theta(x)\mu(dx)=1$ implies that

$A(\theta) = \log\left(\int h(x)\exp(\theta T(x))\mu(dx) \right).$

It turns out that the function $A(\theta)$ is differentiable and that differentiation and integration are exchangeable. This implies that

$A'(\theta) = \frac{\int h(x)\frac{d}{d\theta}\exp(\theta T(x))\mu(dx)}{\int h(x)\exp(\theta T(x))\mu(dx)} = \frac{\int h(x)\frac{d}{d\theta}\exp(\theta T(x))\mu(dx)}{\exp(A(\theta))}$

Note that $\int h(x)\frac{d}{d\theta}\exp(\theta T(x))\mu(dx) = \int T(x)h(x)\exp(\theta T(x))\mu(dx)$ . Thus,

$A'(\theta) = \int T(x)h(x) \exp(\theta T(x)-A(\theta))\mu(dx) = \int T(x)p_\theta(x)\mu(dx).$

This means that $A'(\theta) = \mathbb{E}_\theta[T(X)]$ , the expectation of $T(X)$ under $X \sim P_\theta$ .

Now suppose that we have an i.i.d. sample $X_1,\ldots, X_n \sim P_\theta$ and we want to use this sample to estimate $\theta$ . One way to estimate $\theta$ is by maximum likelihood. That is, we choose the value of $\theta$ that maximises the likelihood,

$L(\theta|X_1,\ldots,X_n) = \prod_{i=1}^n p_\theta(X_i).$

When using the maximum likelihood estimator, it is often easier to work with the log-likelihood. The log-likelihood is,

$\log L(\theta|X_1,\ldots,X_n) = \sum_{i=1}^n \log\left(p_\theta(X_i)\right) = \sum_{i=1}^n \log(h(X_i))+\theta T(X_i) - A(\theta)$ .

Maximising the likelihood is equivalent to maximising the log-likelihood. For exponential families, the log-likelihood is a concave function of $\theta$ . Thus the maximisers can be found be differentiation and solving the first order equations. Note that,

$\frac{d}{d\theta} \log L(\theta|X_1,\ldots,X_n) =\sum_{i=1}^n T(X_i)-A'(\theta) = -nA'(\theta) + \sum_{i=1}^n T(X_i).$

Thus the maximum likelihood estimate (MLE) $\widehat{\theta}$ solves the equation,

$-nA'(\widehat{\theta}) + \sum_{i=1}^n T(X_i) = 0.$

But recall that $A'(\widehat{\theta}) = \mathbb{E}_{\widehat{\theta}}[T(X)]$ . Thus the MLE is the solution to the equation,

$\mathbb{E}_{\widehat{\theta}}[T(X)] = \frac{1}{n}\sum_{i=1}^n T(X_i)$ .

Thus the MLE is the value of $\theta$ for which the expectation of $T(X)$ matches the empirical average from our sample. That is, the maximum likelihood estimator for an exponential family is a method of moments estimator. Specifically, the maximum likelihood estimator matches the moments of the sufficient statistic $T(X)$ .

A counter example

It is a special property of maximum likelihood estimators that the MLE is a method of moments estimator for the sufficient statistic. When we leave the nice world of exponential families, the estimators may differ.

Suppose that we have data $X_1,\ldots,X_n \sim P_\theta$ where $P_\theta$ is the uniform distribution on $[0,\theta]$ . A minimal sufficient statistic for this model is $X_{(n)}$ – the maximum of $X_1,\ldots, X_n$ . Given what we saw before, we might imague that the MLE for this model would be a method of moments estimator for $X_{(n)}$ but this isn’t the case.

The likelihood for $X_1,\ldots,X_n$ is,

$L(\theta|X_1,\ldots,X_n) = \begin{cases} \frac{1}{\theta^n} & \text{if } X_{(n)} \le \theta,\\ 0 & \text{if } X_{(n)} > \theta. \end{cases}$

Thus the MLE is $\widehat{\theta} = X_{(n)}$ . However, under $P_\theta$ , $\frac{1}{\theta}X_{(n)}$ has a $\text{Beta}(n,1)$ distribution. Thus, $\mathbb{E}_\theta[X_{(n)}] = \theta \frac{n}{n+1}$ so the method of moments estimator would be $\widehat{\theta}' = \frac{n+1}{n}X_{(n)} \neq X_{(n)} = \widehat{\theta}$ .

A not so simple conditional expectation

It is winter 2022 and my PhD cohort has moved on the second quarter of our first year statistics courses. This means we’ll be learning about generalised linear models in our applied course, asymptotic statistics in our theory course and conditional expectations and martingales in our probability course.

In the first week of our probability course we’ve been busy defining and proving the existence of the conditional expectation. Our approach has been similar to how we constructed the Lebesgue integral in the previous course. Last quarter, we first defined the Lebesgue integral for simple functions, then we used a limiting argument to define the Lebesgue integral for non-negative functions and then finally we defined the Lebesgue integral for arbitrary functions by considering their positive and negative parts.

Our approach to the conditional expectation has been similar but the journey has been different. We again started with simple random variables, then progressed to non-negative random variables and then proved the existence of the conditional expectation of any arbitrary integrable random variable. Unlike the Lebesgue integral, the hardest step was proving the existence of the conditional expectation of a simple random variable. Progressing from simple random variables to arbitrary random variables was a straight forward application of the monotone convergence theorem and linearity of expectation. But to prove the existence of the conditional expectation of a simple random variable we needed to work with projections in the Hilbert space $L^2(\Omega, \mathbb{P},\mathcal{F})$ .

Unlike the Lebesgue integral, defining the conditional expectation of a simple random variable is not straight forward. One reason for this is that the conditional expectation of a random variable need not be a simple random variable. This comment was made off hand by our Professor and sparked my curiosity. The following example is what I came up with. Below I first go over some definitions and then we dive into the example.

A simple random variable with a conditional expectation that is not simple

Let $(\Omega, \mathbb{P}, \mathcal{F})$ be a probability space and let $\mathcal{G} \subseteq \mathcal{F}$ be a sub- $\sigma$ -algebra. The conditional expectation of an integrable random variable $X$ is a random variable $\mathbb{E}(X|\mathcal{G})$ that satisfies the following two conditions:

The random variable $\mathbb{E}(X|\mathcal{G})$ is $\mathcal{G}$ -measurable.
For all $B \in \mathcal{G}$ , $\mathbb{E}[X1_B] = \mathbb{E}[\mathbb{E}(X|\mathcal{G})1_B]$ , where $1_B$ is the indicator function of $B$ .

The conditional expectation of an integrable random variable is unique and always exists. One can think of $\mathbb{E}(X|\mathcal{G})$ as the expected value of $X$ given the information in $\mathcal{G}$ .

A simple random variable is a random variable $X$ that take only finitely many values. Simple random variables are always integrable and so $\mathbb{E}(X|\mathcal{G})$ always exists but we will see that $\mathbb{E}(X|\mathcal{G})$ need not be simple.

Consider a random vector $(U,V)$ uniformly distributed on the square $[-1,1]^2 \subseteq \mathbb{R}^2$ . Let $D$ be the unit disc $D = \{(u,v) \in \mathbb{R}^2:u^2+v^2 \le 1\}$ . The random variable $X = 1_D(U,V)$ is a simple random variable since $X$ equals $1$ if $(U,V) \in D$ and $X$ equals $0$ otherwise. Let $\mathcal{G} = \sigma(U)$ the $\sigma$ -algebra generated by $U$ . It turns out that

$\mathbb{E}(X|\mathcal{G}) = \sqrt{1-U^2}$ .

Thus $\mathbb{E}(X|\mathcal{G})$ is not a simple random variable. Let $Y = \sqrt{1-U^2}$ . Since $Y$ is a continuous function of $U$ , the random variable is $\mathcal{G}$ -measurable. Thus $Y$ satisfies condition 1. Furthermore if $B \in \mathcal{G}$ , then $B = \{U \in A\}$ for some measurable set $A\subseteq [-1,1]$ . Thus $X1_B$ equals $1$ if and only if $U \in A$ and $V \in [-\sqrt{1-U^2}, \sqrt{1-U^2}]$ . Since $(U,V)$ is uniformly distributed we thus have

$\mathbb{E}[X1_B] = \int_A \int_{-\sqrt{1-u^2}}^{\sqrt{1-u^2}} \frac{1}{4}dvdu = \int_A \frac{1}{2}\sqrt{1-u^2}du$ .

The random variable $U$ is uniformly distributed on $[-1,1]$ and thus has density $\frac{1}{2}1_{[-1,1]}$ . Therefore,

$\mathbb{E}[Y1_B] = \mathbb{E}[\sqrt{1-U^2}1_{\{U \in A\}}] = \int_A \frac{1}{2}\sqrt{1-u^2}du$ .

Thus $\mathbb{E}[X1_B] = \mathbb{E}[Y1_B]$ and therefore $Y = \sqrt{1-U^2}$ equals $\mathbb{E}(X|\mathcal{G})$ . Intuitively we can see this because given $U=u$ , we know that $X$ is $1$ when $V \in [-\sqrt{1-u^2},\sqrt{1+u^2}]$ and that $X$ is $0$ otherwise. Since $V$ is uniformly distributed on $[-1,1]$ the probability that $V$ is in $[-\sqrt{1-u^2},\sqrt{1+u^2}]$ is $\sqrt{1-u^2}$ . Thus given $U=u$ , the expected value of $X$ is $\sqrt{1-u^2}$ .

An extension

The previous example suggests an extension that shows just how “complicated” the conditional expectation of a simple random variable can be. I’ll state the extension as an exercise:

Let $f:[-1,1]\to \mathbb{R}$ be any continuous function with $f(x) \in [0,1]$ . With $(U,V)$ and $\mathcal{G}$ as above show that there exists a measurable set $A \subseteq [-1,1]^2$ such that $\mathbb{E}(1_A|\mathcal{G}) = f(U)$ .

Sports climbing at the 2020 Olympics

Last year sports climbing made its Olympic debut. My partner Claire is an avid climber and we watched some of the competition together. The method of ranking the competitors had some surprising consequences which are the subject of this blog.

In the finals, the climbers had to compete in three disciplines – speed climbing, bouldering and lead climbing. Each climber was given a rank for each discipline. For speed climbing, the rank was given based on how quickly the climbers scaled a set route. For bouldering and lead climbing, the rank was based on how far they climbed on routes none of them had seen before.

To get a final score, the ranks from the three disciplines are multiplied together. The best climber is the one with the lowest final score. For example, if a climber was 5th in speed climbing, 2nd in bouldering and 3rd in lead climbing, then they would have a final score of $5\times 2 \times 3 = 30$ . A climber who came 8th in speed, 1st in bouldering and 2nd in lead climbing would have a lower (and thus better) final score of $8 \times 1 \times 2 = 16$ .

One thing that makes this scoring system interesting is that your final score is very dependent on how well your competitors perform. This was most evident during Jakob Schubert’s lead climb in the men’s final.

Lead climbing is the final discipline and Schubert was the last climber. This meant that ranks for speed climbing and bouldering were already decided. However, the final rankings of the climbers fluctuated hugely as Schubert scaled the 15 metre lead wall. This was because if a climber had done really well in the boulder and speed rounds, being overtaken by Schubert didn’t increase their final score by that much. However, if a climber did poorly in the previous rounds, then being overtaken by Schubert meant their score increased by a lot and they plummeted down the rankings. This can be visualised in the following plot:

Along the x-axis we have how far up the lead wall Schubert climbed (measured by the “hold” he had reached). On the y-axis we have the finalists’ ranking at different stages of Schubert’s climb. We can see that as Schubert climbed and overtook the other climbers the whole ranking fluctuated. Here is a similar plot which also shows when Schubert overtook each competitor on the lead wall:

The volatility of the rankings can really be seen by following Adam Ondra’s ranking (in purple). When Schubert started his climb, Ondra was ranked second. After Schubert passed Albert Gines Lopez, Ondra was ranked first. But then Schubert passed Ondra, Ondra finished the event in 6th place and Gines Lopez came first. If you go to the 4 hour and 27 minutes mark here you can watch Schubert’s climb and hear the commentator explain how both Gines Lopez and Ondra are in the running to win gold.

Similar things happened in the women’s finals. Janja Garnbret was the last lead climber in the women’s final. Here is the same plot which shows the climbers’ final rankings and lead ranking.

Garnbret was a favourite at the Olympics and had come fifth in the speed climbing and first in bouldering. This meant that as long as she didn’t come last in the lead climbing she would at least take home silver and otherwise she’ll get the gold. Garnbret ended up coming first in the lead climbing and we can see that as she overtook the last few climbers, their ranking fluctuated wildly.

Here is one more plot which shows each competitor’s final score at different points of Garnbret’s climb.

In the plot you can really see that, depending on how they performed in the previous two events, each climber’s score changed by a different amount once they were overtaken. It also shows how Garnbret is just so ahead of the competition – especially when you compare the men’s and women’s finals. Here is the same plot for the men. You can see that the men’s final scores were a lot closer together.

Before I end this post, I would like to make one comment about the culture of sport climbing. In this post I wanted to highlight how tumultuous and complicated the sport climbing rankings were but if you watched the athletes you’d have no idea the stakes were so high. The climbers celebrate their personal bests as if no one was watching, they trade ideas on how to tackle the lead wall and the day after the final, they returned to the bouldering wall to try the routes together. Sports climbing is such a friendly and welcoming sport and I would hate for my analysis to give anyone the wrong idea.

Visualising Strava data with R

I’ve recently had some fun downloading and displaying my running data from Strava. I’ve been tracking my runs on Strava for the last five years and I thought it would be interesting to make a map showing where I run. Here is one of the plots I made. Each circle is a place in south east Australia where I ran in the given year. The size of the circle corresponds to how many hours I ran there.

I think the plot is a nice visual diary. Looking at the plot, you can see that most of my running took place in my hometown Canberra and that the time I spend running has been increasing. The plot also shows that most years I’ve spent some time running on the south coast and on my grandfather’s farm near Kempsey. You can also see my trips in 2019 and 2020 to the AMSI Summer School. In 2019, the summer school was hosted by UNSW in Sydney and in 2020 it was hosted by La Trobe University in Melbourne. You can also see a circle from this year at Mt Kosciuszko where I ran the Australian Alpine Ascent with my good friend Sarah.

I also made a plot of all my runs on a world map which shows my recent move to California. In this plot all the circles are the same size and I grouped all the runs across the five different years.

I learnt a few things creating these plots and so I thought I would document how I made them.

Creating the plots

Strava lets you download all your data by doing a bulk export. The export includes a zipped folder with all your activities in their original file format.
My activities where saved as .gpx files and I used this handy python library to convert them to .csv files which I could read into R. For the R code I used the packages “tidyverse”, “maps” and “ozmaps”.
Now I had a .csv files for each run. In these files each row corresponded to my location at a given second during the run. What I wanted was a single data frame where each row corresponded to a different run. I found the following way to read in and edit each .csv file:

files <- list.files(path = ".", pattern = "*.csv")
listcsv <- lapply(files, function(x) read_csv(paste0(x)) %>% 
                    select(lat, lon, time) %>% 
                    mutate(hours = n()/3600) %>% 
                    filter(row_number() == 1)
                  )

The first line creates a list of all the .csv files in the working directory. The second line then goes through the list of file names and converts each .csv file into a tibble. I then selected the rows with the time and my location and added a new column with the duration of the run in hours. Finally I removed all the rows except the first row which contains the information about where my run started.

Next I combined these separate tibbles into a single tibble using rbind(). I then added some new columns for grouping the runs. I added a column with the year and columns with the longitude and latitude rounded to the nearest whole number.

runs <- do.call(rbind, listcsv) %>% 
  mutate(year = format(time,"%Y"),
         approx_lat = round(lat),
         approx_lon = round(lon))

To create the plot where you can see where I ran each year, I grouped the runs by the approximate location and by year. I then calculated the total time spent running at each location each year and calculated the average longitude and latitude. I also removed the runs in the USA by only keeping the runs with a negative latitude.

run_counts_year <- runs %>% 
  group_by(approx_lat, approx_lon, year) %>% 
  summarise(hours = sum(hours),
            lat = mean(lat),
            lon = mean(lon),
            .groups = "drop") %>% 
  select(!approx_lat & !approx_lon)

oz_counts_year <- run_counts_year %>% 
  filter(lat < 0)

I then used the package “ozmaps” to plot my running locations on a map of the ACT, New South Wales and Victoria.

oz_states <- ozmaps::ozmap_states %>% 
  filter(NAME == "New South Wales" |
         NAME == "Victoria" |
         NAME == "Australian Capital Territory")

ggplot() + 
  geom_sf(data = oz_states) +
  coord_sf() +
  geom_point(data = oz_counts_year,
             mapping = aes(x = lon, 
                           y = lat,
                           size = hours),
            color = "blue",
            shape = 21) +
  facet_wrap(~ year) +
  theme_bw()

Creating the world map was similar except I didn’t group by year and I kept the runs with positive latitude.

run_counts <- runs %>% 
  group_by(approx_lat, approx_lon) %>% 
  summarise(hours = sum(hours),
            lat = mean(lat),
            lon = mean(lon),
            .groups = "drop") %>% 
  select(!approx_lat & !approx_lon)

world <- map_data("world")
ggplot() +
  geom_polygon(data = world,
               mapping = aes(x = long, y = lat, group = group),
               fill = "lightgrey",
               color = "black",
               size = 0.1) +
  geom_point(data = run_counts,
             mapping = aes(x = lon, 
                           y = lat),
             size = 2,
             shape = 1,
             color = "blue") +
  theme_bw()

Extremal couplings

This post is inspired by an assignment question I had to answer for STATS 310A – a probability course at Stanford for first year students in the statistics PhD program. In the question we had to derive a few results about couplings. I found myself thinking and talking about the question long after submitting the assignment and decided to put my thoughts on paper. I would like to thank our lecturer Prof. Diaconis for answering my questions and pointing me in the right direction.

What are couplings?

Given two distribution functions $F$ and $G$ on $\mathbb{R}$ , a coupling of $F$ and $G$ is a distribution function $H$ on $\mathbb{R}^2$ such that the marginals of $H$ are $F$ and $G$ . Couplings can be used to give probabilistic proofs of analytic statements about $F$ and $G$ (see here). Couplings are also are studied in their own right in the theory optimal transport.

We can think of $F$ and $G$ as being the cumulative distribution functions of some random variables $X$ and $Y$ . A coupling $H$ of $F$ and $G$ thus corresponds to a random vector $(\widetilde{X},\widetilde{Y})$ where $\widetilde{X}$ has the same distribution as $X$ , $\widetilde{Y}$ has the same distribution as $Y$ and $(\widetilde{X},\widetilde{Y}) \sim H$ .

The independent coupling

For two given distributions function $F$ and $G$ there exist many possible couplings. For example we could take $H = H_I$ where $H_I(x,y) = F(x)G(y)$ . This coupling corresponds to a random vector $(\widetilde{X}_I,\widetilde{Y}_I)$ where $\widetilde{X}_I$ and $\widetilde{Y}_I$ are independent and (as is required for all couplings) $\widetilde{X}_I \stackrel{\text{dist}}{=} X$ , $\widetilde{Y}_I \stackrel{\text{dist}}{=} Y$ .

In some sense the coupling $H_I$ is in the “middle” of all couplings. This is because $\widetilde{X}$ and $\widetilde{Y}$ are independent and so $\widetilde{X}$ doesn’t carry any information about $\widetilde{Y}$ . As the title of the post suggests, there are couplings were this isn’t the case and $\widetilde{X}$ carries “as much information as possible” about $\widetilde{Y}$ .

The two extremal couplings

Define two function $H_L, H_U :\mathbb{R}^2 \to [0,1]$ by

$H_U(x,y) = \min\{F(x), G(y)\}$ and $H_L(x,y) = \max\{F(x)+G(y) - 1, 0\}$ .

With some work, one can show that $H_L$ and $H_U$ are distributions functions on $\mathbb{R}^2$ and that they have the correct marginals. In this post I would like to talk about how to construct random vectors $(\widetilde{X}_U, \widetilde{Y}_U) \sim H_U$ and $(\widetilde{X}_L, \widetilde{Y}_L) \sim H_L$ .

Let $F^{-1}$ and $G^{-1}$ be the quantile functions of $F$ and $G$ . That is,

$F^{-1}(c) = \inf\{ x \in \mathbb{R} : F(x) \ge c\}$ and $G^{-1}(c) = \inf\{ x \in \mathbb{R} : G(x) \ge c\}$ .

Now let $V$ be a random variable that is uniformly distributed on $[0,1]$ and define

$\widetilde{X}_U = F^{-1}(V)$ and $\widetilde{Y}_U = G^{-1}(V)$ .

Since $F^{-1}(V) \le x$ if and only if $V \le F(x)$ , we have $\widetilde{X}_U \stackrel{\text{dist}}{=} X$ and likewise $\widetilde{Y}_U \stackrel{\text{dist}}{=} Y$ . Furthermore $\widetilde{X}_U \le x, \widetilde{Y}_U \le y$ occurs if and only if $V \le F(x), V \le G(y)$ which is equivalent to $V \le \min\{F(x),G(y)\}$ . Thus

$\mathbb{P}(\widetilde{X}_U \le x, \widetilde{Y}_U \le y) = \mathbb{P}(V \le \min\{F(x),G(y)\})= \min\{F(x),G(y)\}.$

Thus $(\widetilde{X}_U,\widetilde{Y}_U)$ is distributed according to $H_U$ . We see that under the coupling $H_U$ , $\widetilde{X}_U$ and $\widetilde{Y}_U$ are closely related as they are both increasing functions of a common random variable $V$ .

We can follow a similar construction for $H_L$ . Define

$\widetilde{X}_L = F^{-1}(V)$ and $\widetilde{Y}_L = G^{-1}(1-V)$ .

Thus $\widetilde{X}_L$ and $\widetilde{Y}_L$ are again functions of a common random variable $V$ but $\widetilde{X}_L$ is an increasing function of $V$ and $\widetilde{Y}_L$ is a decreasing function of $V$ . Note that $1-V$ is also uniformly distributed on $[0,1]$ . Thus $\widetilde{X}_L \stackrel{\text{dist}}{=} X$ and $\widetilde{Y}_L \stackrel{\text{dist}}{=} Y$ .

Now $\widetilde{X}_L \le x, \widetilde{Y}_L \le y$ occurs if and only if $V \le F(x)$ and $1-V \le G(y)$ which occurs if and only if $1-G(y) \le V \le F(x)$ . If $1-G(y) \le F(x)$ , then $F(x)+G(y)-1 \ge 0$ and $\mathbb{P}(1-G(y) \le V \le F(x)) =F(x)+G(y)-1$ . On the other hand, if $1 - G(y) > F(x)$ , then $F(x)+G(y)-1< 0$ and $\mathbb{P}(1-G(y) \le V \le F(x))=0$ . Thus

$\mathbb{P}(\widetilde{X}_L \le x, \widetilde{Y}_L \le y) = \mathbb{P}(1-G(y) \le V \le F(x)) = \max\{F(x)+G(y)-1,0\}$ ,

and so $(\widetilde{X}_L,\widetilde{Y}_L)$ is distributed according to $H_L$ .

What makes $H_U$ and $H_L$ extreme?

Now that we know that $H_U$ and $H_L$ are indeed couplings, it is natural to ask what makes them “extreme”. What we would like to say is that $\widetilde{Y}_U$ is an increasing function of $\widetilde{X}_U$ and $\widetilde{Y}_L$ is a decreasing function of $\widetilde{X}_L$ . Unfortunately this isn’t always the case as can be seen by taking $X$ to be constant and $Y$ to be continuous.

However the intuition that $\widetilde{Y}_U$ is increasing in $\widetilde{X}_U$ and $\widetilde{Y}_L$ is decreasing in $\widetilde{X}_L$ is close to correct. Given a coupling $(\widetilde{X},\widetilde{Y}) \sim H$ , we can look at the quantity

$C(x,y) = \mathbb{P}(\widetilde{Y} \le y | \widetilde{X} \le x) -\mathbb{P}(\widetilde{Y} \le y) = \frac{H(x,y)}{F(x)}-G(y)$

This quantity tells us something about how $\widetilde{Y}$ changes with $\widetilde{X}$ . For instance if $\widetilde{X}$ and $\widetilde{Y}$ were positively correlated, then $C(x,y)$ would be positive and if $\widetilde{X}$ and $\widetilde{Y}$ were negatively correlated, then $C(x,y)$ would be negative.

For the independent coupling $(\widetilde{X}_I,\widetilde{Y}_I) \sim H_I$ , the quantity $C(x,y)$ is constantly $0$ . It turns out that the above probability is maximised by the coupling $(\widetilde{X}_U, \widetilde{Y}_U) \sim H_U$ and minimised by $(\widetilde{X}_L,\widetilde{Y}_L) \sim H_L$ and it is in this sense that they are extremal. This final claim is the two dimensional version of the Fréchet-Hoeffding Theorem and checking it is a good exercise.

You could have proved the Neyman-Pearson lemma

The Neyman-Pearson lemma is foundational and important result in the theory of hypothesis testing. When presented in class the proof seemed magical and I had no idea where the ideas came from. I even drew a face like this 😲 next to the usual $\square$ in my book when the proof was finished. Later in class we learnt the method of undetermined multipliers and suddenly I saw where the Neyman-Pearson lemma came from.

In this blog post, I’ll first give some background and set up notation for the Neyman-Pearson lemma. Then I’ll talk about the method of undetermined multipliers and show how it can be used to derive and prove the Neyman-Pearson lemma. Finally, I’ll write about why I still think the Neyman-Pearson lemma is magical despite the demystified proof.

Background

In the set up of the Neyman-Pearson lemma we have data $x$ which is a realisation of some random variable $X$ . We wish to conclude something about the distribution of $X$ from our data $x$ by doing a hypothesis test.

In the Neyman-Pearson lemma we have simple hypotheses. That is our data either comes from the distribution $\mathbb{P}_0$ or from the distribution $\mathbb{P}_1$ . Thus, our null hypothesis is $H_0 : X \sim \mathbb{P}_0$ and our alternative hypothesis is $H_1 : X \sim \mathbb{P}_1$ .

A test of $H_0$ against $H_1$ is a function $\phi$ that takes in data $x$ and returns a number $\phi(x) \in [0,1]$ . The value $\phi(x)$ is the probability under the test $\phi$ of rejecting $H_0$ given the observed data $x$ . That is, if $\phi(x) = 1$ , we always reject $H_0$ and if $\phi(x)=0$ we never reject $H_0$ . For in-between values $\phi(x) = \gamma \in (0,1)$ , we reject $H_0$ with probability $\gamma$ .

An ideal test would have two desirable properties. We would like a test that rejects $H_0$ with a low probability when $H_0$ is true but we would also like the test to reject $H_0$ with a high probability when $H_1$ is true. To state this more formally, let $\mathbb{E}_0[\phi(X)]$ and $\mathbb{E}_1[\phi(X)]$ be the expectation of $\phi(X)$ under $H_0$ and $H_1$ respectively. The quantity $\mathbb{E}_0[\phi(X)]$ is the probability we reject $H_0$ when $H_0$ is true. Likewise, the quantity $\mathbb{E}_1[\phi(X)]$ is the probability we reject $H_0$ when $H_1$ is true. An optimal test would be one that minimises $\mathbb{E}_0[\phi(X)]$ and maximises $\mathbb{E}_1[\phi(X)]$ .

Unfortunately the goals of minimising $\mathbb{E}_0[\phi(X)]$ and maximising $\mathbb{E}_1[\phi(X)]$ are at odds with one another. This is because we want $\phi$ to be small in order to minimise $\mathbb{E}_0[\phi(X)]$ and we want $\phi$ to be large to maximise $\mathbb{E}_1[\phi(X)]$ . In nearly all cases we have to trade off between these two goals and there is no single test that simultaneously achieves both.

To work around this, a standard approach is to focus on maximising $\mathbb{E}_1[\phi(X)]$ while requiring that $\mathbb{E}_0[\phi(X)]$ remains below some threshold. The quantity $\mathbb{E}_1[\phi(X)]$ is called the power of the test $\phi$ . If $\alpha$ is a number between $0$ and $1$ , we will say that $\phi$ has level– $\alpha$ if $\mathbb{E}_1[\phi(X)] \le \alpha$ . A test $\phi$ is said to be most powerful at level- $\alpha$ , if

The test $\phi$ is level- $\alpha$ .
For all level- $\alpha$ tests $\phi'$ , the test $\phi$ is more powerful than $\phi'$ . That is,

$\mathbb{E}_1[\phi'(X)] \le \mathbb{E}_1[\phi(X)]$ .

Thus we can see that finding a most powerful level- $\alpha$ test is a constrained optimisation problem. We wish to maximise the quantity

$\mathbb{E}_1[\phi(X)]$

subject to the constraint

$\mathbb{E}_0[\phi(X)] \le \alpha$

With this in mind, we turn to the method of undetermined multipliers.

The method of undetermined multipliers

The method of undetermined multipliers (also called the method of Lagrange multipliers) is a very general optimisation tool. Suppose that we have a set $U$ and two function $f,g : U \to \mathbb{R}$ and we wish to maximise $f(u)$ subject to the constraint $g(u) \le 0$ .

In the context of hypothesis testing, the set $U$ is the set of all tests $\phi$ . The objective function $f$ is given by $f(\phi) = \mathbb{E}_1[\phi(X)]$ . That is, $f(\phi)$ is the power of the test $\phi$ . The constraint function $g$ is given by $g(\phi)=\mathbb{E}_1[\phi(X)]-\alpha$ so that $g(\phi) \le 0$ if and only if $\phi$ has level- $\alpha$ .

The method of undetermined multipliers allows us to reduce this constrained optimisation problem to a (hopefully easier) unconstrained optimisation problem. More specifically we have the following result:

Proposition: Suppose that $u^* \in U$ is such that:

$g(u^*) = 0$ ,
There exists $k \ge 0$ , such that $u^*$ maximises $f(u)-kg(u)$ over all $u \in U$ .

Then $u$ maximises $f(u)$ under the constraint $g(u) \le 0$ .

Proof: Suppose that $u^*$ satisfies the above two dot points. We need to show that for all $u \in U$ , if $g(u) \le 0$ , then $f(u) \le f(u^*)$ . By assumption we know that $f(u^*)=0$ and $u^*$ maximises $f(u)-kg(u)$ . Thus, for all $u \in U$ ,

$f(u^*)=f(u^*)-kg(u^*) \ge f(u)-kg(u)$ .

Now suppose $g(u) \le 0$ . Then, $kg(u) \le 0$ and so $f(u)-kg(u)\ge f(u)$ and so $f(u^*) \ge f(u)$ as needed. $\square$

The constant $k$ is the undetermined multiplier. The way to use the method of undetermined is to find values $u_k^*$ that maximise $h_k(u) = f(u)-kg(u)$ for each $k \ge 0$ . The multiplier $k$ is then varied so that the constraint $g(u^*_k) = 0$ is satisfied.

Proving the Neyman-Pearson lemma

Now let’s use the method of undetermined multipliers to find most powerful tests. Recall the set $U$ which we are optimising over is the set of all tests $\phi$ . Recall also that we wish to optimise $f(\phi) = \mathbb{E}_1[\phi(X)]$ subject to the constraint $g(\phi) = \mathbb{E}_0[\phi(X)] - \alpha \le 0$ . The method of undetermined multipliers says that we should consider maximising the function

$h_k(\phi) = \mathbb{E}_1[\phi(X)] - k\mathbb{E}_0[\phi(X)]$ ,

where $k \ge 0$ . Suppose that both $\mathbb{P}_0$ and $\mathbb{P}_1$ have densities₁ $p_0$ and $p_1$ with respect to some measure $\mu$ . We can we can write the above expectations as integrals. That is,

$\mathbb{E}_0[\phi(X)] = \int \phi(x)p_0(x)\mu(dx)$ and $\mathbb{E}_1[\phi(X)] = \int \phi(x)p_1(x)\mu(dx)$ .

Thus the function $h_k$ is equal to

$h_k(\phi) = \int \phi(x)p_1(x)\mu(dx)- k\int \phi(x)p_0(x)\mu(dx)=\int \phi(x)(p_1(x)-kp_0(x))\mu(dx)$ .

We now wish to maximise the function $h_k(\phi)$ . Recall that $\phi$ is a function that take values in $[0,1]$ . Thus, the integral

$\int \phi(x)(p_1(x)-kp_0(x))\mu(dx)$ ,

is maximised if and only if $\phi(x)=1$ when $p_1(x)-kp_0(x) > 0$ and $\phi(x)=0$ when $p_1(x)-kp_0(x) < 0$ . Note that $p_1(x)-kp_0(x) > 0$ if and only if $\frac{p_1(x)}{p_0(x)} > k$ . Thus for each $k \ge 0$ , a test $\phi_k$ maximises $h_k(\phi)$ if and only if

$\phi_k(x) = \begin{cases} 1 & \text{if } \frac{p_1(x)}{p_0(x)} > k,\\ 0 &\text{if } \frac{p_1(x)}{p_0(x)} < k. \end{cases}$

The method of undetermined multipliers says that if we can find $k$ so that the above is satisfied and $g(\phi_k) = 0$ , then $\phi_k$ is a most powerful test. Recall that $g(\phi_k)=0$ is equivalent to $\mathbb{E}_1[\phi(X)]=\alpha$ . By summarising the above argument, we arrive at the Neyman-Pearson lemma,

Neyman-Pearson Lemma₂: Suppose that $\phi$ is a test such that

$\mathbb{E}_0[\phi(X)] = \alpha$ , and
For some $k \ge 0$ , $\phi(x) = \begin{cases} 1 & \text{if } \frac{p_1(x)}{p_0(x)} > k,\\ 0 & \text{if } \frac{p_1(x)}{p_0(x)}< k.\end{cases}$

then $\phi$ is most powerful at level- $\alpha$ for testing $H_0 : X \sim \mathbb{P}_0$ against $H_1 : X \sim \mathbb{P}_1$ .

The magic of Neyman-Pearson

By learning about undetermined multipliers I’ve been able to better understand the proof of the Neyman-Pearson lemma. I now view it is as clever solution to a constrained optimisation problem rather than something that comes out of nowhere.

There is, however, a different aspect of Neyman-Pearson that continues to surprise me. This aspect is the usefulness of the lemma. At first glance the Neyman-Pearson lemma seems to be a very specialised result because it is about simple hypothesis testing. In reality most interesting hypothesis tests have composite nulls or composite alternatives or both. It turns out that Neyman-Pearson is still incredibly useful for composite testing. Through ideas like monotone likelihood ratios, least favourable distributions and unbiasedness, the Neyman-Pearson lemma or similar ideas can be used to find optimal tests in a variety of settings.

Thus I must admit that the title of this blog post is a little inaccurate and deceptive. I do believe that, given the tools of undetermined multipliers and the set up of simple hypothesis testing, one is naturally led to the Neyman-Pearson lemma. However, I don’t believe that many could have realised how useful and interesting simple hypothesis testing would be.

Footnotes

The assumption that $\mathbb{P}_0$ and $\mathbb{P}_1$ have densities with respect to a common measure $\mu$ is not a restrictive assumption since one can always take $\mu = \mathbb{P}_0+\mathbb{P}_1$ and the apply Radon-Nikodym. However there is often a more natural choice of $\mu$ such as Lebesgue measure on $\mathbb{R}^d$ or the counting measure on $\mathbb{N}$ .
What I call the Neyman-Pearson lemma is really only a third of the Neyman-Pearson lemma. There are two other parts. One that guarantees the existence of a most powerful test and one that is a partial converse to the statement in this post.

An art and maths collaboration

Over the course of the past year I have had the pleasure to work with the artist Sanne Carroll on her honours project at the Australian National University. I was one of two mathematics students that collaborated with Sanne. Over the course of the project Sanne drew patterns and would ask Ciaran and I to recreate them using some mathematical or algorithmic ideas. You can see the final version of project here: https://www.sannecarroll.com/ (best viewed on a computer).

I always loved the patterns Sanne drew and the final project is so well put together. Sanne does a great job of incorporating her drawings, the mathematical descriptions and the communication between her, Ciaran and me. Her website building skills also far surpass anything I’ve done on this blog!

It was also a lot of fun to work with Sanne. Hearing about her patterns and talking about maths with her was always fun. I also learnt a few things about GeoGebra which made the animations in my previous post a lot quicker to make. Sanne has told me that she’ll be starting a PhD soon and I’m looking forward to any future collaborations that might arise.

Leverages Scores

I am very excited to be writing a blog post again – it has been nearly a year! This post marks a new era for the blog. In September I started a statistics PhD at Stanford University. I am really enjoying my classes and I am learning a lot. I might have to change the name of the blog soon but for now let’s stick with “Maths to Share” although you will undoubtedly see more and more statistics here.

Today I would like to talk about leverages scores. Leverages scores are a way to quantify how sensitive a model is and they can be used to explain the different behaviour in these two animations

Linear Models

I recently learnt about leverage scores in the applied statistics course STATS 305A. This course is all about the linear model. In the linear model we assume with have $n$ data points $(x_i,y_i)$ where $x_i$ is a vector in $\mathbb{R}^d$ and $y_i$ is a number in $\mathbb{R}$ . We model $y_i$ as a linear function of $x_i$ plus noise. That is we assume $y_i = x_i^T\beta + \varepsilon_i$ , where $\beta \in \mathbb{R}^d$ is a unknown vector of coefficients and $\varepsilon_i$ is a random variable with mean $0$ and variance $\sigma^2$ . We also require that for $i \neq j$ , the random variable $\varepsilon_i$ is uncorrelated with $\varepsilon_j$ .

We can also write this as a matrix equation. Define $y$ to be the vector with entries $y_i$ and define $X$ to be the matrix with rows $x_i^T$ , that is

$y = \begin{bmatrix} y_1\\ y_2 \\ \vdots \\ y_n \end{bmatrix} \in \mathbb{R}^n$ and $X = \begin{bmatrix} -x_1^T-\\-x_2^T-\\ \vdots \\ -x_n^T-\end{bmatrix} \in \mathbb{R}^{n \times d}.$

Then our model can be rewritten as

$y = X\beta + \varepsilon,$

where $\varepsilon \in \mathbb{R}^n$ is a random vector with mean $0$ and covariance matrix $\sigma^2 I_n$ . To simplify calculations we will assume that $X$ contains an intercept term. This means that the first column of $X$ consists of all 1’s.

In the two animations at the start of this post we have two nearly identical data sets. The data sets are an example of simple regression when each vector $x_i$ is of the form $(1,z_i)$ where $z_i$ is a number. The values $z_i$ are on the horizontal axis and the values $y_i$ are on the vertical axis.

Estimating the coefficients

In the linear model we wish to estimate the parameter $\beta$ which contains the coefficients of our model. That is, given a sample $(y_i,x_i)_{i=1}^n$ , we wish to construct a vector $\widehat{\beta}$ which approximates the true parameter $\beta$ . In ordinary least square regression we choose $\widehat{\beta}$ to be the vector $b \in \mathbb{R}^d$ that minimizes the quantity

$\sum_{i=1}^n (x_i^T b - y_i)^2=\left \Vert Xb - y \right \Vert_2^2$ .

Differentiating with respect to $b$ and setting the derivative equal to $0$ shows that $\widehat{\beta}$ is a solution to the normal equations:

$X^TXb = X^T y.$

We will assume that the matrix $X^TX$ is invertible. In this case then the normal equations have a unique solution $\widehat{\beta} = (X^TX)^{-1}X^T y$ .

Now that we have our estimate $\widehat{\beta}$ , we can do prediction. If we are given a new value $x' \in \mathbb{R}^d$ we would use $x'^T\widehat{\beta}$ to predict the corresponding value of $y'$ . This was how the straight lines in the two animations were calculated.

We can also calculate the model’s predicted values for the data $x_i$ that we used to fit the model. These are denoted by $\widehat{y}$ . Note that

$\widehat{y} = X\widehat{\beta} = X(X^TX)^{-1}X^Ty = Hy,$

where $H = X(X^TX)^{-1}X^T$ is called the hat matrix for the model (since it puts the hat $\widehat{ }$ on $y$ .

Leverage scores

We are now ready to talk about leverage scores and the two animations. For reference, here they are again:

In both animations the stationary line corresponds to an estimator $\widehat{\beta}$ that was calculated using only the black data points. The red points are new data points with different $x$ values and varying $y$ values. The moving line corresponds to an estimator $\widehat{\beta}$ calculated using the red data point as well as all the black points. We can see immediately that if the red point is far away from the “bulk” of the other $x$ points, then the moving line is a lot more sensitive to the $y$ value of the red point.

The leverage score of a data point $(x_i,y_i)$ is defined to be $\frac{\partial \widehat{y}_i}{\partial y_i}.$ That is, the leverage score tells us how much does the prediction $\widehat{y}_i$ change if we change $y_i$ .

Since $\widehat{y} = Hy$ , the leverage score of $(x_i,y_i)$ is $H_{ii}$ , the $i^{th}$ diagonal element of the hat matrix $H$ . The idea is that if a data point $(x_i,y_i)$ has a large leverage score, then the model is more sensitive to changes in that value of $y_i$ .

This can be seen in a leave one out calculation. This calculation tells us what we should expect if we make a leave-one-out model – a model that uses all the data points apart from one. In our animations, this corresponds to the stationary line.

The leave one out calculation says that the predicted value using all the data is always between the true value and the predicted value from the leave-one-out model. In our animations this can be seen by noting that the moving line (the full model) is always between the red point (the true value) and the stationary line (the leave-one-out model).

Furthermore the leverage score tells us exactly how close the predicted value is to the true value. We can see that the moving line is much closer to the red dot in the high leverage example on the right than the low leverage example on the left.

Mahalanobis distance

We now know that the two animations are showing the sensitivity of a model to two different data points. We know that a model is more sensitive to point with high leverage than to points with low leverage. We still haven’t spoken about why some point have higher leverage and why the point on the right has higher leverage.

It turns out that leverage score are measuring how far away a data point is from the “bulk” of the other $x_i$ ‘s. More specifically in a one dimensional example like what we have in the animations

$H_{ii} = \frac{1}{n}\left(\frac{1}{S^2}(x_i-\bar{x})^2+1\right),$

where $n$ is the number of data points, $\bar{x} = \frac{1}{n}\sum_{j=1}^n x_j$ is the sample mean and $S^2 = \frac{1}{n}\sum_{j=1}^n (x_j-\bar{x})^2$ is the sample variance. Thus high leverage scores correspond to points that are far away from the centre of our data $x_i$ . In higher dimensions a similar result holds if we measure distance using Mahalanobis distance.

The mean of the black data points is approximately 2 and so we can now see why the second point has higher leverage. The two animations were made in Geogebra. You can play around with them here and here.

Commuting in Groups

In July 2020 I submitted my honours thesis which was titled “Commuting in Groups” and supervised by Professor Anthony Licata. You can read the abstract below and the full thesis here.

Abstract

In this thesis, we study four different classes of groups coming from geometric group theory. Each of these classes are defined in terms of fellow traveller conditions. First we study automatic groups and show that such groups have a solvable word problem. We then study hyperbolic groups and show that a group is hyperbolic
if and only if it is strongly geodesically automatic. We also show that a group is hyperbolic if and only if it has a divergence function. We next study combable groups and examine some examples of groups that are combable but not automatic. Finally we introduce biautomatic groups. We show that biautomatic groups have a solvable conjugacy problem and that many nilpotent groups cannot be subgroups of biautomatic groups. Finally we introduce Artin groups of finite type and show that all such groups are biautomatic.

Sums of Squares – Part 2: Motzkin

In the previous blog post we observed that if a polynomial $g$ could be written as a sum $\sum_{i=1}^n h_i^2$ where each $h_i$ is a polynomial, then $g$ must be non-negative. We then asked if the converse was true. That is, if $g$ is a non-negative polynomial, can $g$ be written a sum of squares of polynomials? As we saw, a positive answer to this question would allow us to optimize a polynomial function $f$ by looking at the values of $\lambda$ such that $g = f - \lambda$ can be written as a sum of squares.

Unfortunately, as we shall see, not all non-negative polynomials can be written as sums of squares.

The Motzkin Polynomial

The Motzkin polynomial is a two variable polynomial defined by

$g(x,y) = x^4y^2 +x^2y^4-3x^2y^2+1$ .

We wish to prove two things about this polynomial. The first claim is that $g(x,y) \ge 0$ for all $x,y \in \mathbb{R}$ and the second claim is that $g$ cannot be written as a sum of squares of polynomials. To prove the first claim we will use the arithmetic mean geometric mean inequality. This inequality states that for all $n \in \mathbb{N}$ and all $a_1,a_2,\ldots, a_n \ge 0$ , we have that

$\left(a_1a_2\ldots a_n\right)^{1/n} \le \frac{a_1+a_2+\ldots+a_n}{n}.$

We will apply this inequality with $n = 3$ , $a_1 = x^4 y^2$ , $a_2 = x^2 y^4$ and $a_3 = 1$ . This gives that

$\left(x^4 y^2 x^2 y^4 1\right)^\frac{1}{3} \le \frac{x^4 y^2 + y^2 x^4 +1 }{3}$ .

Simplifying the left hand side and multiplying both sides by $3$ gives that

$3x^2 y^2 \le x^4 y^2 + y^2 x^4 + 1$ .

Since our Motzkin polynomial $g$ is given by $g(x,y) = x^4 y^2 + y^2 x^2 - 3x^2 y^2 +1$ , the above inequality is equivalent to $g(x,y) \ge 0$ .

Newton Polytopes

Thus we have shown that the Motzkin polynomial is non-negative. We now wish to show that it is not a sum of squares. To do this we will first have to define the Newton polytope of a polynomial. If our polynomial $f$ has $n$ variables, then the Newton polytope of $f$ , denoted $N(f)$ is a convex subset of $\mathbb{R}^n$ . To define $N(f)$ we associate a point in $\mathbb{R}^n$ to each term of the polynomial $f$ based on the degree of each variable. For instance the term $3x^2y^4$ is assigned the point $(2,4) \in \mathbb{R}^2$ and the term $4xyz^3$ is assigned the point $(1,1,3) \in \mathbb{R}^2$ . Note that the coefficients of the term don’t affect this assignment.

We then define $N(f)$ to be the convex hull of all points assigned to terms of the polynomial $f$ . For instance if $f(x,y) = x^4y^3+5x^3y - 7y^3+8$ , then the Newton polytope of $f$ is the following subset of $\mathbb{R}^2$ .

Note again that changing the non-zero coefficients of the polynomial $f$ does not change $N(f)$ .

Now suppose that our polynomial is a sum of squares, that is $f = \sum_{i=1}^n h_i^2$ . It turns out the the Newton polytope of $f$ can be defined in terms of the Newton polytopes of $h_i$ . In particular we have that $N(f) =2X$ where $X$ is the convex hull of the union of $N(h_i)$ for $i = 1,\ldots, n$ .

To see why this is true, note that if $a$ and $b$ are the points corresponding to two terms $p$ and $q$ , then $a+b$ corresponds to $pq$ . Thus we can see that every term of $f$ corresponds to a point that can be written as $a+b$ for some $a,b \in N(h_i)$ . It follows that $N(f) \subseteq 2X$ .

Conversely note that if we have terms $p, q, r$ corresponding to points $a,b, c$ and $pq = r^2$ , then $a+b = 2c$ . It follows that if $c$ is a vertex in $X$ corresponding to the term $r$ in some polynomial $h_i$ and $p,q$ are two terms in a polynomial $h_j$ such that $pq = r^2$ , then $p=q=r$ . This is because if $r$ was not equal to either $p$ or $q$ , then the point $c$ would not be a vertex and would instead lie on the line between the points assigned to $p$ and $q$ .

It follows that if $c$ is a vertex of $X$ with corresponding term $r$ , then $r^2$ appears with a positive coefficient in $f = \sum_{i=1}^n h_i$ . It follows that $2X \subseteq N(f)$ and so $2X = N(f)$ .

Not a sum of squares

Let’s now take another look at the Motzkin polynomial which is defined to be

$g(x,y) = x^4y^2 + x^2 y^4 - 3x^2y^2 + 1$ .

Thus the Newton polytope of $g$ has corners $(4,2), (2,4), (0,0)$ and looks like this

Thus if the Motzkin polynomial was a sum of squares $g = \sum_{i=1}^n h_i^2$ , then the Newton polytope of each $h_i$ must be contained in $\frac{1}{2}N(g)$ . Now $\frac{1}{2}N(g)$ looks like this

The only integer points contained in $\frac{1}{2}N(g)$ are $(0,0), (1,1), (1,2)$ and $(2,1)$ . Thus each $h_i$ must be equal to $h_i(x,y) = a_i x^2 y+b_i xy^2 + c_i xy + d_i$ . If we square this polynomial we see that the coefficient of $x^2 y^2$ is $c_i^2$ . Thus the coefficient of $x^2y^2$ in $\sum_{i=1}^n h_i^2$ must be $c_1^2+c_2^2+\ldots +c_n^2 \ge 0$ . But in the Motzkin polynomial, $x^2y^2$ has coefficient $-3$ .

Thus the Motzkin polynomial is a concrete example of a non-negative polynomial that is not a sum of squares.

References

As stated in the first part of this series, these two blog posts are based off conversations I had with Abhishek Bhardwaj last year. I also used these slides to remind myself why the Motzkin polynomial is not a sum of squares.