ICD-10 chapters and code letters

Posted on by John

I’ve been thinking about ICD-10 codes; they come up a lot in my work.

The ICD-10-CM standard is divided into 21 chapters, which generally correspond to the first letter of a code. However, a chapter may contain blocks beginning with more than one letter, and codes starting with a single letter, namely D, can span two chapters.

Here’s a diagram I made to visualize the relationship between chapters and initial letters of codes.

Notice there’s no letter U on the diagram. That’s because U is reserved for special/provisional codes.

Related posts

Posterior variance

Posted on by John

A few days ago I wrote a post entitled Does additional data always reduce posterior variance?. In a nutshell, the answer is no, not always.

That led the previous post which looked at posterior means for three Bayesian models, showing how the posterior mean is a weighted average of the prior mean and the mean of the new data. The weights are precisions, which means something different for each model.

For the beta-binomial model, variance may increase when seeing unexpected data (details here), but precision always increases.

For the normal-normal model precision is the reciprocal of variance. Every new data point makes precision go up and posterior variance go down.

The Poisson-gamma model may be the most interesting. As stated in the previous post, if data has a Poisson distribution with parameter λ, and λ has a gamma(α0, β0) prior distribution, then the posterior distribution on λ after observing k events over time t has a gamma(α0 + k, β0 + t) posterior distribution. Therefore the posterior variance is

(α0 + k) / (β0 + t)².

Note the posterior variance is an increasing function of k and a decreasing function of t. This means that the posterior variance increases every time an event is observed, and it decreases quadratically between observations.

Here’s an illustration. I simulated data from a Poisson process with λ and used a gamma(1, 1) prior on λ. Here’s a plot of the posterior variance.

Posterior mean

Posted on by John

Common sense says that what you believe after seeing new data should be some sort of compromise between what you believed before and what the new data says. You don’t want to ignore previous information or new information.

How much should new data change your prior beliefs? When prior judgment and new information are in conflict, which one should be given the benefit of the doubt?

Bayesian data models provide a framework for making such decisions quantitative and objective. The choice of a data model is somewhat subjective—whether it’s a Bayesian model or not—but given a Bayesian model, the rules for updating the representation of your beliefs are objective. As some put it, you “turn the Bayesian crank.” A likelihood model and a prior on parameters together specify how new data changes the prior distribution into a posterior distribution.

We will make this more concrete with three examples.

Normal-normal model

Suppose that data X has a normal distribution with unknown mean μ and known variance σ², and we assume that a priori μ has a normal distribution with mean μ0 and variance σ0².

After observing x, the posterior distribution on μ also has a normal distribution, but with a different mean and variance. Its mean is somewhere between the prior mean and x. We will ignore the change in the variance for this post.

The posterior mean of μ is

\mu_{\text{post}} = \frac{\dfrac{\mu_0}{\sigma_0^2} + \dfrac{x}{\sigma^2}}{\dfrac{1}{\sigma_0^2} + \dfrac{1}{\sigma^2}}

This equation becomes more understandable when we introduce precisons τ = 1/σ² and τ0 = 1/σ0².

Then we have

\mu_{\text{post}} = \frac{\mu_0 \cdot \tau_0 + x \cdot \tau}{\tau_0 + \tau}

which you can read as saying the posterior mean is the weighted average of the prior mean and x, with the weights given by the precision. Intuitively, you take the weighted mean of your conclusions from previous data and new data, weighting the mean according to how much confidence you have in each.

Beta-binomial model

Now let’s switch over to a different data model. Now assume X is a binary random variable, with probability of success p and probability of failure 1 − p, and we assume p has a beta(a, b) distribution.

After observing s successes and f failures, the posterior mean of the distribution on p becomes

p_{\text{post}} = \frac{a + s}{a + b + s + f}

We can rewrite this as

p_{\text{post}} = \frac{(a + b) \dfrac{a}{a+b} + (s + f) \dfrac{s}{s+f}}{(a + b) + (s +f)}

This says that the posterior mean is the weighted average of the prior mean a/(a + b) and the mean of the data s/n. The weights are the prior effective sample size a + b and the sample size of the new data n. In this example (effective) sample size is playing the role that precision played in the normal-normal model above.

Gamma-Poisson model

Suppose data have a Poisson distribution with parameter λ, and λ has a gamma(α0, β0) prior distribution [1]. And suppose you observe k events over time t. Then the posterior distribution of λ given the data has a gamma(α0 + k, β0 + t) prior distribution and the mean of the posterior distribution is given by

\lambda_{\text{post}} = \frac{\alpha_0 + k}{\beta_0 + t} = \frac{\beta_0 (\alpha_0 / \beta_0) + t (k / t)}{\beta_0 + t}

As before, the posterior mean is a weighted average of the prior mean and new data, and the weights are interpretable as some sort of measure of confidence, namely time. The variable t is directly time and the parameter β0 is sort of an effective time, just as a + b is an effective sample size for the beta distribution.

Common thread

In each example the posterior mean is the weighted average of the prior mean and the mean of the data, with the weights given by a precision. However, precision means something different in each example. In the normal-normal model, precision is the reciprocal of variance, but in the beta-binomial model precision is sample size and in the Poisson-gamma model precision is time.

What all three examples have in common is that they are conjugate models using distributions from the “exponential family” of probability distributions. In technical terms, precision is the multiplicative factor on the sufficient statistic in the exponent of the posterior kernel.

Related posts

[1] There are multiple conventions for parameterizing the gamma distribution. Here we’re using the shape-rate parameterization, where the mean is α/β.

Progress on Gilbreath’s conjecture

Posted on by John

Years ago I wrote about Gilbreath’s conjecture. It’s a simple conjecture; you could explain it to anyone who understands what prime numbers are. See the linked post for a description of the problem.

Gilbreath’s conjecture is simple, but it’s also kinda weird. As I wrote before,

Paul Erdős speculated that Gilbreath’s conjecture is true but it would be 200 years before anyone could prove it. I find Erdős’s conjecture more interesting than Gilbreath’s conjecture.

The conjecture is hard in a way that, say, solving a nasty-looking differential equation is not. Over the last three centuries, mathematics has developed quite a toolbox for solving differential equations. But Gilbreath’s conjecture is just odd enough that it’s not at all clear what kind of tool might be useful in approaching it.

Terence Tao has a new blog post announcing a paper that he and two coauthors wrote on a random model intended to mimic Gilbreath’s calculation on primes. This random model is more sophisticated than the little game Gilbreath was playing, but it’s also much more amenable to analysis by established techniques. Tao’s post gives a heuristic explanation for why Gilbreath’s conjecture is plausible, but then adds

However, it seems well beyond current technology to try to make these heuristics rigorous; even the first step … is far out of reach.

Related posts

Reproducing a geometry theorem diagram

Posted on by John

I ran across a geometry theorem with the following diagram.

The theorem corresponding to the diagram is interesting, but I found reproducing the diagram more interesting.

The segment AB is a diameter and the line CD is perpendicular to the diameter.

Assume the outer circle is a unit circle. I guessed C = (cos(1), sin(1)) and made the following diagram.

I guessed the value of C by eyeballing it, but in retrospect this would have been a convenient value for the creator of the original diagram to have chosen.

Drawing the blue circle inscribed in the triangle was easy using the equations for the center and radius from this post. Drawing the other two circles, the green and orange circles, was harder. They are also inscribed circles, but not inscribed in a triangle. They’re inscribed in a three-sided figure with two perpendicular sides and a circular arc.

The radius r of the green circle is the distance from the center of the circle to each of its tangent lines. Also, the distance from the origin to the center of the circle must be 1 − r. This is enough information to set up a quadratic equation for r. The same reasoning applies to the orange circle.

The original diagram comes from [1] and the theorem it illustrates says the diameter of the blue circle equals the sum of the radii of the green and orange circles.

Python code

In case you’re interested, here’s the code that created the diagram.

#!/usr/bin/env -S uv run --script

# /// script
# dependencies = ["numpy", "matplotlib"]
# ///

import numpy as np
import matplotlib.pyplot as plt

def connect(A, B, color='gray'):
    plt.plot([A[0], B[0]], [A[1], B[1]], color=color, linewidth=2)

def circle(c, r, color='gray'):
    t = np.linspace(0, 2*np.pi)
    plt.plot(c[0] + r*np.cos(t), c[1] + r*np.sin(t), color=color, linewidth=2)

def quadratic(a, b, c):
    det = b**2 - 4*a*c
    return ((-b - det**0.5)/(2*a), (-b + det**0.5)/(2*a))

A = np.array([-1, 0])
B = np.array([ 1, 0])
C = np.array([np.cos(1), np.sin(1)])
a = np.linalg.norm(B - C)
b = np.linalg.norm(A - C)
c = np.linalg.norm(B - A)
s = (a + b + c)/2

circle([0,0], 1)
connect(A, B,)
connect(A, C)
connect(C, B)
connect(C, C*np.array([1, -1]))

center = (a*A + b*B + c*C)/(2*s)
radius = 0.5*a*b/s
circle(center, radius, 'C0')

Ex = C[0]
roots = quadratic(1, 2 + 2*Ex, Ex**2 - 1)
r = roots[1] # Smaller root is negaive
print(roots)
center = (r + Ex, -r)
circle(center, r, 'C1')

roots = quadratic(1, 2 - 2*Ex, Ex**2 - 1)
r = roots[1] # Smaller root is negaive
center = (Ex - r, -r)
circle(center, r, 'C2')

plt.gca().set_aspect("equal")
plt.axis("off")
plt.show()

[1] Leon Bankoff. A Geometrical Coincidence. Mathematics Magazine, Vol. 37, No. 5 (Nov., 1964), p. 324.

e approximation

Posted on by John

I ran across the approximation

e ≈ 2721/1001

recently. What makes this remarkable is its accuracy relative to the size of the denominator.

You can create a trivial approximation just by truncating a decimal expansion

e ≈ 2718/1000

but this is only good to four significant figures, but 2721/1001 is good to seven, almost eight, significant figures.

e         = 2.71828182… 
2721/1001 = 2.71828171…

The comparison is more impressive in binary.

$ bc -l
>>> obase=2
>>> 2721/1001
10.10110111111000010100…
>>> e(1)
10.10110111111000010101…

The denominator is a 10-bit number but the approximation is accurate to 21 bits.

Does additional data always reduce posterior variance?

Posted on by John

A discussion over lunch today brought up the fact that additional data does not always decrease the size of a confidence interval. This post will look at this from a Bayesian perspective.

In general, new information reduces your uncertainty regarding whatever you’re estimating. The posterior distribution becomes more concentrated as more data are collected.

That’s what happens “in general” but does it necessarily happen every time you get new data? Conceivably if you get surprising data, data that is very unlikely given your current prior, posterior uncertainty might increase.

Beta-binomial model

To show that this is the case, suppose the probability of success in some binary trial has parameter θ and that θ has a beta prior. You could imagine this prior to be the posterior after having made some number of previous observations. Can a new observation increase the posterior variance in θ? If so, under what conditions?

The variance of a beta(a, b) random variable is

ab / (a + b)²(a + b + 1).

After observing a successful trial, the posterior distribution on θ is beta(a + 1, b). We can calculate the ratio of the posterior variance to the prior variance and ask under what circumstances, if any, the ratio is greater than 1.

If 2a ≥ b the posterior variance will be strictly less than the prior variance. This says if the prior mean odds against a success are no more than 2 : 1, observing a success will reduce the variance. (So will observing a failure.) But for any value of b, you can find a small enough value of a that observing a success will increase the variance.

Normal-normal model

Whether an observation can increase the posterior variance depends on the data model. If your data have a normal likelihood function with known variance and a normal prior on the mean θ, the posterior variance is always less than the prior observation, and it reduces by the same amount, independent of the observation x. If x is very unlikely a priori then it will pull the posterior mean toward itself more than an observation that is more concordant with the prior would have, but the change in the posterior variance is the same.

Proof of beta theorem

Here is a proof in Lean 4 of the statement above that if 2a ≥ b the posterior variance will be strictly less than the prior variance.

import Mathlib

set_option linter.style.header false

noncomputable def f (a b : ℝ) : ℝ := a * b / ((a + b) ^ 2 * (a + b + 1))

theorem f_ratio_lt_one' (a b : ℝ) (ha : 0 < a) (hb : 0 < b) (hab : b ≤ 2 * a) :
    f (a + 1) b / f a b < 1 := by
  have hs : 0 < a + b := by linarith
  have h2ab : 0 ≤ 2 * a - b := by linarith
  have hprod : 0 ≤ (a + b) * (2 * a - b) := mul_nonneg hs.le h2ab
  -- key polynomial inequality (∗)
  have key : (a + 1) * (a + b) ^ 2 < a * ((a + b + 1) * (a + b + 2)) := by
    nlinarith [hprod, ha]
  -- nonzero facts needed to clear denominators
  have ha' : a ≠ 0 := ne_of_gt ha
  have hb' : b ≠ 0 := ne_of_gt hb
  have hs' : a + b ≠ 0 := ne_of_gt hs
  have hs1' : a + b + 1 ≠ 0 := by positivity
  have hs2' : a + b + 2 ≠ 0 := by positivity
  have ha1' : a + 1 ≠ 0 := by positivity
  -- express the ratio as a single closed-form fraction
  have hratio : f (a + 1) b / f a b
      = ((a + 1) * (a + b) ^ 2) / (a * ((a + b + 1) * (a + b + 2))) := by
    unfold f
    have e : a + 1 + b = a + b + 1 := by ring
    rw [e]
    field_simp
    ring
  rw [hratio, div_lt_one (by positivity)]
  exact key

DNA Sequence Alignment and Kings

Posted on by John

This morning I wrote a post that included the central Delannoy numbers. The nth central Delannoy number Dn counts the number of ways a king can move from one corner of a chessboard to the diagonally opposite corner without backtracking.

The more general Delannoy numbers Dm,n are the analogy for an m × n rectangular board, not necessarily square.

Dm,n is also the number of possible sequence alignments for a strand of DNA with m base pairs and a strand with n base pairs [1]. At each step in the alignment process, you can introduce a gap in the first strand, the second strand or neither, which is analogous to the king who can move N, E, or NE at each step.

The Delannoy numbers can be computed recursively:

def D(m, n):
    if m == 0 or n == 0:
        return 1
    return D(m - 1, n) + D(m, n - 1) + D(m - 1, n - 1)

The code above can be sped up tremendously by adding the decorator

@lru_cache(maxsize=None)

above the function definition to turn on memoization. I did an experiment computing D12,15 with and without memoization and the times were 77.1805 seconds and 0.000062 seconds respectively, i.e. memoization made the code over a million times faster.

Incidentally, D12,15 = 2653649025 and so there are a lot of ways to align even short sequences unless you place some restriction on the permissible alignments.

Update: Here’s a heatmap plotting log10(Dm,n). Obviously the function increases with m and n: bigger chessboards have more possible paths. Moreover, it’s larger along the diagonal (i.e. the central Delannoy numbers). If you look along northeast to southwest diagonals, the function is largest in the middle where m = n.

[1] Torres, A., Cabada, A., & Nieto, J. J. (2003). An exact formula for the number of alignments between two DNA sequences. DNA Sequence, 14(6), 427–430. https://doi.org/10.1080/10425170310001617894

Distinguishing variables from parameters

Posted on by John

Imagine the following dialog.

Professor: f is a function of a real variable x that takes a real parameter k.

Student: What’s a parameter?

Professor: It’s a constant that can vary.

Student: Then if it can vary, isn’t it a variable?

Professor: Sorta, but no not really.

This conversation plays out over and over, and unfortunately it often ends as it does above, with the student confused. Here’s how I believe the conversation should continue.

Professor: You’re absolutely right that f is a function of two variables, x and k. But usually k is fixed in the context of a specific application and x is not. A different application might have a different, but also fixed, value of k. So it is helpful to think of f(x; k), a function of x with a parameter k, rather than f(x, k), a function of two variables. The former carries more information, giving a hint as to how the numbers are used.

Is there really a difference between a parameter and a variable? In a reductionistic sense, no. But in a practical sense, yes, absolutely.

It might sound pedantic to distinguish a variable from a parameter, and it is, in the best sense of the word. Pedant literally means teacher. Usually pedantic carries a negative connotation, such as making a distinction without a difference. But here the pedant would be making a helpful distinction.

For example, we might write a probability density function as f(x; μ, σ). The function gives the probability density at a point x. The density depends on parameters μ and σ, and these parameters change between applications, but for a given application they have fixed values.

You find the probability of a random variable taking on values in an interval [a, b] by integrating f over that interval. When I say that, you know that I mean you’d integrate with respect to x, because f is a function of x. It is also, in an abstract sense, a function of μ and σ, but it’s typically not useful to think of it that way.

Hypergeometric functions have two sets of parameters, and so you may see two semicolons, such as f(x; a, b; c). This denotes a function of the variable x, with upper parameters a and b, and a lower parameter c. In some abstract sense this is a function of four variables, but it acts very differently with respect to x than with respect to a, b, and c. There’s also a difference between a and b on the one hand and c on the other, one worth paying attention to, though it is less of a difference than between x and the parameters collectively.

Sometimes you’ll see a vertical bar rather than a semicolon to separate variables from parameters. This works out even better for probability densities because then f(x | μ, σ) suggests the probability density of x given μ and σ since the vertical bar is also used for conditional probability. You might also see f(x | a, b; c) for hypergeometric functions, with the vertical bar separating variables from parameters and the semicolon separating two kinds of parameters.

When I first saw a semicolon separating variables from parameters, no explanation was given, and I figured I could mentally replace the semicolon with a comma. Then later I realized that the semicolon was an act of kindness by the author giving the reader additional information.

Silver Rectangles and the Ways of Kings

Posted on by John

Golden rectangles

The defining property of golden rectangle is that if you stick a square on its longer side, you get another golden rectangle.

The smaller vertical rectangle is similar to the larger horizontal rectangle. This means

φ / 1 = (1 + φ) / φ

which tells us φ² = 1 + φ and so the golden ratio φ equals (1 + √5)/2.

Silver rectangles

A silver rectangle is one that if you stick two squares on its longer side you get another rectangle with the same aspect ratio.

This tells us

σ / 1 = (1 + 2σ) / σ

and so σ² = 1 + 2σ and the silver ratio is σ = 1 + √2.

Just as you can define a golden ratio and a silver ratio, there’s an analogous way to define a sequence of metallic ratios.

Kings and Delannoy numbers

The silver ratio has several connections to the ways of ways kings. By that I mean the number of ways a king can go from one corner of a chessboard to the diagonally opposite corner without backtracking.

A king can move one space in any direction. If we start with a king in the bottom left corner of the board, the no-backtracking requirement means the king can move up, right, or up and right.

The number of paths a king can take from one corner to the opposite corner of an n × n chessboard is the nth central Delannoy number Dn. more generally Delannoy numbers are defined for an m × n chessboard, but I’ll stick to the case m = n called the central Delannoy number, or just Delannoy numbers for short.

The first Delannoy number is 1 because there’s only one way for a king to get from one corner to the other: do nothing, because the opposite corner is the same corner. The second Delannoy number is 3 because the king can move up then right, or right then up, or move diagonally up and right.

For a 3 × 3 grid things are significantly more complicated, and D3 = 13. For an 8 × 8 grid the number of paths is 48,639.

Generating function

How would you estimate the number of paths on an n × n board for large values of n without calculating it exactly? You might start by finding a generating function for the Delannoy numbers, which works out to be

(x² − 6x + 1)−1/2

The radius of convergence r for the generating function series is the distance from 0 to the closest singularity of the generating function, which is the smaller root of

x² − 6x + 1

which is

3 − √8 = (3 + √8)−1 = (1 + √2)−2 = 1/σ²

i.e. the radius of convergence is the reciprocal of the silver ratio squared.

Asymptotic estimate

The radius of convergence gives us a first approximation to the asymptotic size of the series coefficients. Since we’re working with the generating function of the Delannoy numbers, these coefficients are the Delannoy numbers. That is,

Dn ~ r−n = (σ2)n = σ2n.

That’s as good as you can do just knowing the radius of convergence. A more careful analysis would refine this estimate by dividing by a factor proportional to √n.

Related posts