Solving hard problems

We help companies solve hard problems in mathematics, statistics, and computing. Let’s explore how we might work together.

PHI and offshore data processing

Posted on 23 October 2018 by John

The US government does not prohibit the transfer of PHI (protected health information) offshore [1], but it does subject offshore data processing to extra reporting [2] and more scrutiny in general. The CMS (Centers for Medicare & Medicaid Services, part of the Department of Health and Human Services) has said

Given the unique risks associated with the use of contractors operating outside the jurisdiction of the United States, CMS encourages sponsors using offshore subcontractors to take extraordinary measures to ensure that offshore arrangements protect beneficiary privacy [3, emphasis added].

Why are the risks unique? For one thing, subcontractors outside the US may not be familiar with US law or held accountable for complying with that law.

What are the necessary “extraordinary measures”? The CMS does not say because the answer depends very much on context. The probability of being able to identify someone from a set of data fields depends on what those fields are. It also can change over time, and so the conclusion needs to be reviewed periodically.

It is legal to give a third party access to PHI if there is a BAA (business associate agreement) in place. For example, if a hospital outsources billing, the company doing the billing must see personal information in order to carry out their business function. But with offshore processing, it seems safest, if practical, to proceed as if there were no BAA in place and deidentify the data before it leaves the US.

If you’d like help with privacy regulation compliance or data de-identification, let’s talk.

[1] Nothing in this blog post is legal advice. Compliance with privacy regulation is both a legal and statistical matter. I address statistical questions and let lawyers address legal questions. I have a lawyer specializing in healthcare law who I recommend if clients don’t have their own legal expert.

[2] CMS memo September 16, 2011. Contract Year (CY) 2012 Medicare Advantage and Part D Readiness Checklist. Available here.

[3] Allen Briskin, Lisa C. Earl, Gerry Hinkley, and Joseph E. Kendall. Offshoring Health Information: Issues and Lingering Concerns. Journal of Health & Life Sciences Law. Vol. 8, No. 1.

Fundamental units

Posted on 23 October 2018 by John

It’s much easier to convert meters to kilometers than to convert yards to miles. You know what’s even easier? Not converting meters to kilometers!

You only need one unit of length, say a meter. You could express all distances in terms of meters and dispense with kilometers, millimeters, etc. And in fact, this is essentially what science does. It hardly makes a difference whether you say the Andromeda galaxy is 2.365 × 10²² meters away or 2.365 × 10¹⁹ kilometers away. Especially if you need scientific notation anyway, you might as well use base units.

Multiple units for the same fundamental quantity are psychologically convenient though not logically necessary. Astronomers, for example, simply use meters for calculation. But they may use other units such as AUs or light years for conceptualization and communication. Scientists do all time calculations in seconds, though it’s easier to understand long time scales in units of years.

It’s commonly said that it’s easier to do science in SI units than in imperial units. This is in some sense true, but it would be more accurate to say it’s easier to do science in terms of fundamental units. The SI prefixes kilo, mega, etc. are a sort of intermediate step to letting go of all but fundamental units.

One could imagine an alternative history in which science standardized on the inch rather than the meter. Distances like yards and miles could have been seen as auxiliary units for conceptual purposes, much like the astronomical unit and light year.

***

This post was motivated by an online discussion of kilobytes and megabytes. I said something about a kilobyte being 1,000 bytes, and someone replied saying technically this should be 1,024. I replied saying that even more technically, it’s 1,000 bytes. This is a good example of coming full circle. The least informed and most informed agree, while those in the middle disagree.

You might naively assume that because 1,000 grams is a kilogram that 1,000 bytes is a kilobyte. But for years, a kilobyte was commonly understood to mean 1,024 bytes. To address the confusion, the IEC introduced new prefixes such as kibi and mebi in 1998. Now a kibibyte is 2¹⁰ = 1,024 bytes, a mebibyte is 2²⁰ = 1,048,576 bytes, etc. And so a kilobyte now refers to 1,000 bytes, just as one might naively have assumed.

In practice the difference hardly matters. If you want to indicate exactly 2²⁰ bytes, you can say “one mebibyte.” But in my opinion it’s even better to simply say “2²⁰ bytes” because then there is no chance of being misunderstood.

To my mind, terms like “megabyte” are something like light years, conceptual units. Whether you think a megabyte is exactly a million bytes or approximately a million bytes doesn’t matter in informal communication. If it’s necessary to be precise, simply state the number of bytes as an integer, using powers of 10 or 2 if that helps.

The purpose of communication is to be understood, not to be correct. If you use megabyte to mean exactly 1,000,000 bytes, you may be correct according to the IEC, but you run a significant risk of being misunderstood. Better, in my opinion, to just use the fundamental unit bytes.

I made a similar point in my post on names for extremely large numbers. When you use these names, there’s a good chance your listener will either not know what you’re talking about, or worse, misunderstand what you mean.

Milk and wine

Posted on 22 October 2018 by John

gallon of milk and bottle of wine

The US uses a mix of imperial and metric units of measure. Some people, almost all outside the US, are quite agitated by this. In practice, the two systems coexist peacefully.

Americans buy milk by the gallon and wine by the milliliter. Milk typically comes in gallon jugs, and wine typically comes in 750 milliliter bottles. The inconsistency is curious but harmless.

For practical purposes, milk is sold in units of jugs, not gallons. And wine is sold in units of bottles, not milliliters. When I go to the store to buy a jug of milk or a bottle of wine, I’m not that aware of the units. If overnight the size of a milk jug changed from one gallon to four liters, or if wine went from 750 milliliter bottles to 25 fluid ounce bottles, no one would know.

If you went to a store determined to buy exactly equal volumes of milk and wine, you couldn’t do it. You’d need to buy more containers of each than any store could possibly sell [1]. But nobody ever does this.

And if for some bizarre reason you really did want to buy equal volumes of milk and wine, the biggest problem would not be gallons versus milliliters but rather the uncertainty in both. A jug of milk isn’t exactly one gallon, nor is a bottle of wine exactly 750 milliliters.

[1] A US gallon is 231 cubic inches, and an inch is 2.54 centimeters. From this can work out that 31,250,000 milk jugs contain the same volume as 157,725,491 wine bottles.

Integer odds and prime numbers

Posted on 17 October 2018 by John

For every integer m > 1, it’s possible to choose N so that the proportion of primes in the sequence 1, 2, 3, … N is 1/m. To put it another way, you can make the odds against one of the first N natural numbers being prime any integer value you’d like [1].

For example, suppose you wanted to find N so that 1/7 of the first N positive integers are prime. Then the following Python code shows you could pick N = 3059.

    from sympy import primepi

    m = 7

    N = 2*m
    while N / primepi(N) != m:
        N += m
    print(N)

[1] Solomon Golomb. On the Ratio of N to π(N). The American Mathematical Monthly, Vol. 69, No. 1 (Jan., 1962), pp. 36-37.

The proof is short, and doesn’t particularly depend on the distribution of primes. Golomb proves a more general theorem for any class of integers whose density goes to zero.

My most popular posts on Reddit

Posted on 16 October 2018 by John

There are only three posts on this top 10 list that are also on the top 10 list for Hacker News.

Comparing trig functions and Jacobi functions

Posted on 13 October 2018 by John

My previous post looked at Jacobi functions from a reference perspective: given a Jacobi function defined one way, how do I relate that to the same function defined another way, and how would you compute it?

This post explores the analogy between trigonometric functions and Jacobi elliptic functions.

Related basic Jacobi functions to trig functions

In the previous post we mentioned a connection between the argument u of a Jacobi function and the amplitude φ:

$u = \int_0^{\varphi} \frac{d\theta}{\sqrt{1-m\sin^2\theta}}$

We can use this to define the functions sn and cn. Leaving the dependence on m implicit, we have

$\begin{align*} \mathrm{sn}(u) &= \sin(\varphi) \\ \mathrm{cn}(u) &= \cos(\varphi) \end{align*}$

If m = 0, then u = φ and so sn and cn are exactly sine and cosine.

There’s a third Jacobi function we didn’t discuss much last time, and that’s dn. It would be understandable to expect dn might be analogous to tangent, but it’s not. The function dn is the derivative of φ with respect to u, or equivalently

$\mathrm{dn}(u) = \sqrt{1 - m \sin^2\varphi}$

The rest of the Jacobi functions

Just as there are more trig functions than just sine and cosine, there are more Jacobi functions than sn, cn, and dn. There are two ways to define the rest of the Jacobi functions: in terms of ratios, and in terms of zeros and poles.

Ratios

I wrote a blog post one time asking how many trig functions there are. The answer depends on your perspective, and I gave arguments for 1, 3, 6, and 12. For this post, lets say there are six: sin, cos, tan, sec, csc, and ctn.

One way to look at this is to say there are as many trig functions as there are ways to select distinct numerators and denominators from the set {sin, cos, 1}. So we have tan = sin/cos, csc = 1/sin, etc.

There are 12 Jacobi elliptic functions, one for each choice of distinct numerator and denominator from {sn, cn, dn, 1}. The name of a Jacobi function is two letters, denoting the numerator and denominator, where we have {s, c, d, n} abbreviating {sn, cn, dn, 1}.

For example, cs(u) = sn(u) / cn(u) and ns(u) = 1 / sn(u).

Note that to take the reciprocal of a Jacobi function, you just reverse the two letters in its name.

Zeros and poles

The twelve Jacobi functions can be classified [1] in terms of their zeros and poles over a rectangle whose sides have length equal to quarter periods. Let’s look at an analogous situation for trig functions before exploring Jacobi functions further.

Trig functions are periodic in one direction, while elliptic functions are periodic in two directions in the complex plane. A quarter period for the basic trig functions is π/2. The six trig functions take one value out of {0, 1, ∞} at 0 and different value at π/2. So we have one trig function for each of the six ways to chose an permutation of length 2 from a set of 3 elements.

In the previous post we defined the two quarter periods K and K‘. Imagine a rectangle whose corners are labeled

s = (0, 0)
c = (K, 0)
d = (K, K‘)
n = (0, K‘)

Each Jacobi function has a zero at one corner and a pole at another. The 12 Jacobi function correspond to the 12 ways to chose a permutation of two items from a set of four.

The name of a Jacobi function is two letters, the first letter corresponding to where the zero is, and the second letter corresponding to the pole. So, for example, sn has a zero at s and a pole at n. These give same names as the ratio convention above.

Identities

The Jacobi functions satisfy many identities analogous to trigonometric identities. For example, sn and cn satisfy a Pythagorean identity just like sine and cosine.

$\mathrm{sn}^2 u + \mathrm{cn}^2 u = 1$

Also, the Jacobi functions have addition theorems, though they’re more complicated than their trigonometric counterparts.

$\begin{align*} \mathrm{sn}(u + v) &= \frac{\mathrm{sn}\,u\, \mathrm{cn}\, v\, \mathrm{dn}\,v\, + \mathrm{sn}\,v\, \mathrm{cn}\, u\, \mathrm{dn}\,u\,}{1 - m\, \mathrm{sn}^2 u\, \, \mathrm{sn^2} v\,} \\ \\ \mathrm{\mathrm{cn}\,}(u + v) &= \frac{\mathrm{cn}\, u\, \mathrm{cn}\, v\, - \mathrm{sn}\,u\, \mathrm{dn}\,u\, \mathrm{sn}\,v\, \mathrm{dn}\,v\,}{1 - m\, \mathrm{sn}^2 u\, \, \mathrm{sn^2} v\,} \\ \\ \mathrm{dn}(u + v) &= \frac{\mathrm{dn}\,u\, \mathrm{dn}\,v\, - m\, \mathrm{sn}\,u\, \mathrm{cn}\, u\, \mathrm{sn}\,v\, \mathrm{cn}\, v\,}{1 - m\, \mathrm{sn}^2 u\, \, \mathrm{sn^2} v\,} \end{align*}$

Derivatives

The derivatives of the basic Jacobi functions are given below.

$\begin{align*} \mathrm{sn}'(u) &= \mathrm{cn}(u)\, \mathrm{dn}(u) \\ \\ \mathrm{cn}'(u) &= -\mathrm{sn}(u) \,\mathrm{dn}(u) \\ \\ \mathrm{dn}'(u) &= -m\,\mathrm{sn}(u)\, \mathrm{cn}(u) \\ \end{align*}$

Note that the implicit parameter m makes an appearance in the derivative of dn. We will also need the complementary parameter m‘ = 1 – m.

The derivatives of all Jacobi functions are summarized in the table below.

$\begin{table} \centering \begin{tabular}{l|rrrr} \multicolumn{1}{l}{} & s & n & d & c \\ \cline{2-5} s & & dn cn & nd cd & nc dc \\ n & $-$ ds cs & & $m$ sd cd & sc dc \\ d & $-$ ns cs & $-m$ sn cn & & $m'$ sc nc \\ c & $-$ ns ds & $-$ sn dn & $-m'$ sd nd & \end{tabular} \end{table}$

The derivatives of the basic Jacobi functions resemble those of trig functions. They may look more complicated at first, but they’re actually more regular. You could remember them all by observing the patterns pointed out below.

Let w, x, y, z be any permutation of {s, n, d, c}. Then the derivative of wx is proportional to yx zx. That is, the derivative of every Jacobi function f is a constant times two other Jacobi functions. The names of these two functions both end in the same letter as the name of f, and the initial letters are the two letters not in the name of f.

The proportionality constants also follow a pattern. The sign is positive if and only if the letters in the name of f appear in order in the sequence s, n, d, c. Here’s a table of just the constants.

$\begin{table} \centering \begin{tabular}{l|rrrr} \multicolumn{1}{l}{} & s & n & d & c \\ \cline{2-5} s & & 1 & 1 & 1 \\ n & $-1$ & & $m$ & 1 \\ d & $-1$ & $-m$ & & $m'$ \\ c & $-1$ & $-1$ & $-m'$ & \end{tabular} \end{table}$

Note that the table is skew symmetric, i.e. its transpose is its negative.

[1] An elliptic function is determined, up to a constant multiple, by its periods, zeros, and poles. So not only do the Jacobi functions have the pattern of zeros and poles described here, these patterns uniquely determine the Jacob functions, modulo a constant. For (singly) periodic functions, the period, zeros, and poles do not uniquely determine the function. So the discussion of zeros and poles of trig functions is included for comparison, but it does not define the trig functions.

Clearing up the confusion around Jacobi functions

Posted on 12 October 2018 by John

The Jacobi elliptic functions sn and cn are analogous to the trigonometric functions sine and cosine. The come up in applications such as nonlinear oscillations and conformal mapping. Unfortunately there are multiple conventions for defining these functions. The purpose of this post is to clear up the confusion around these different conventions.

Plot of Jacobi sn

The image above is a plot of the function sn [1].

Modulus, parameter, and modular angle

Jacobi functions take two inputs. We typically think of a Jacobi function as being a function of its first input, the second input being fixed. This second input is a “dial” you can turn that changes their behavior.

There are several ways to specify this dial. I started with the word “dial” rather than “parameter” because in this context parameter takes on a technical meaning, one way of describing the dial. In addition to the “parameter,” you could describe a Jacobi function in terms of its modulus or modular angle. This post will be a Rosetta stone of sorts, showing how each of these ways of describing a Jacobi elliptic function are related.

The parameter m is a real number in [0, 1]. The complementary parameter is m‘ = 1 – m. Abramowitz and Stegun, for example, write the Jacobi functions sn and cn as sn(u | m) and cn(u | m). They also use m₁ = rather than m‘ to denote the complementary parameter.

The modulus k is the square root of m. It would be easier to remember if m stood for modulus, but that’s not conventional. Instead, m is for parameter and k is for modulus. The complementary modulus k‘ is the square root of the complementary parameter.

The modular angle α is defined by m = sin² α.

Note that as the parameter m goes to zero, so does the modulus k and the modular angle α. As any one of these three goes to zero, the Jacobi functions sn and cn converge to their counterparts sine and cosine. So whether your dial is the parameter, modulus, or modular angle, sn converges to sine and cn converges to cosine as you turn the dial toward zero.

As the parameter m goes to 1, so does the modulus k, but the modular angle α goes to π/2. So if your dial is the parameter or the modulus, it goes to 1. But if you think of your dial as modular angle, it goes to π/2. In either case, as you turn the dial to the right as far as it will go, sn converges to hyperbolic secant, and cn converges to the constant function 1.

Quarter periods

In addition to parameter, modulus, and modular angle, you’ll also see Jacobi function described in terms of K and K‘. These are called the quarter periods for good reason. The functions sn and cn have period 4K as you move along the real axis, or move horizontally anywhere in the complex plane. They also have period 4iK‘. That is, the functions repeat when you move a distance 4K‘ vertically [2].

The quarter periods are a function of the modulus. The quarter period K along the real axis is

$K(m) = \int_0^{\pi/2} \frac{d\theta}{\sqrt{1-m\sin^2\theta}}$

and the quarter period K‘ along the imaginary axis is given by K‘(m) = K(m‘) = K(1 – m).

The function K(m) is known as “the complete elliptic integral of the first kind.”

Amplitude

So far we’ve focused on the second input to the Jacobi functions, and three conventions for specifying it.

There are two conventions for specifying the first argument, written either as φ or as u. These are related by

$u = \int_0^{\varphi} \frac{d\theta}{\sqrt{1-m\sin^2\theta}}$

The angle φ is called the amplitude. (Yes, it’s an angle, but it’s called an amplitude.)

When we said above that the Jacobi functions had period 4K, this was in terms of the variable u. Note that when φ = π/2, u = K.

Jacobi elliptic functions in Mathematica

Mathematica uses the u convention for the first argument and the parameter convention for the second argument.

The Mathematica function JacobiSN[u, m] computes the function sn with argument u and parameter m. In the notation of A&S, sn(u | m).

Similarly, JacobiCN[u, m] computes the function cn with argument u and parameter m. In the notation of A&S, cn(u | m).

We haven’t talked about the Jacobi function dn up to this point, but it is implemented in Mathematica as JacobiDN[u, m].

The function that solves for the amplitude φ as a function of u is JacobiAmplitude[um m].

The function that computes the quarter period K from the parameter m is EllipticK[m].

Jacobi elliptic functions in Python

The SciPy library has one Python function that computes four mathematical functions at once. The function scipy.special.ellipj takes two arguments, u and m, just like Mathematica, and returns sn(u | m), cn(u | m), dn(u | m), and the amplitude φ(u, m).

The function K(m) is implemented in Python as scipy.special.ellipk.

[1] The plot was made using JacobiSN[0.5, z] and the function ComplexPlot described here.

[2] Strictly speaking, 4iK‘ is a period. It’s the smallest vertical period for cn, but 2iK‘ is the smallest vertical period for sn.

Hadamard product

Posted on 10 October 2018 by John

The first time you see matrices, if someone asked you how you multiply two matrices together, your first idea might be to multiply every element of the first matrix by the element in the same position of the corresponding matrix, analogous to the way you add matrices.

But that’s not usually how we multiply matrices. That notion of multiplication hardly involves the matrix structure; it treats the matrix as an ordered container of numbers, but not as a way of representing a linear transformation. Once you have a little experience with linear algebra, the customary way of multiplying matrices seems natural, and the way that may have seemed natural at first glance seems kinda strange.

The componentwise product of matrices is called the Hadamard product or sometimes the Schur product. Given two m by n matrices A and B, the Hadamard product of A and B, written A ∘ B, is the m by n matrix C with elements given by

c_ij = a_ij b_ij.

Because the Hadamard product hardly uses the linear structure of a matrix, you wouldn’t expect it to interact nicely with operations that depend critically on the linear structure. And yet we can give a couple theorems that do show a nice interaction, at least when A and B are positive semi-definite matrices.

The first is the Schur product theorem. It says that if A and B are positive semi-definite n by n matrices, then

det(A ∘ B) ≥ det(A) det(B)

where det stands for determinant.

Also, there is the following theorem of Pólya and Szegö. Assume A and B are symmetric positive semi-definite n by n matrices. If the eigenvalues of A and B, listed in increasing order, are α_i and β_i respectively, then for every eigenvalue λ of A ∘ B, we have

α₁ β₁ ≤ λ ≤ α_n β_n.

Python implementation

If you multiply two (multidimensional) arrays in NumPy, you’ll get the componentwise product. So if you multiply two matrices as arrays you’ll get the Hadamard product, but if you multiply them as matrices you’ll get the usual matrix product. We’ll illustrate that below. Note that the function eigvalsh returns the eigenvalues of a matrix. The name may look a little strange, but the “h” on the end stands for “Hermitian.” We’re telling NumPy that the matrix is Hermitian so it can run software specialized for that case [1].

    
    from numpy import array, matrix, array_equal, all
    from numpy.linalg import det, eigvalsh
    
    A = array([
        [3, 1],
        [1, 3]
    ])
    
    B = array([
        [5, -1],
        [-1, 5]
    ])
    
    H = array([
        [15, -1],
        [-1, 15]
    ])
    
    AB = array([
        [14,  2],
        [ 2, 14]
    ])
    
    # Hadamard product
    assert(array_equal(A*B, H))
    
    # Ordinary matrix product
    assert(array_equal(A@B, AB))
    
    # Schur product theorem
    assert(det(H) >= det(A)*det(B))
    
    # Eigenvalues
    eigA = eigvalsh(A)
    eigB = eigvalsh(B)
    eigH = eigvalsh(A*B)
    
    lower = eigA[0]*eigB[0]
    upper = eigA[1]*eigB[1]
    assert(all(eigH >= lower))
    assert(all(eigH <= upper))

The code above shows that the eigenvalues of A are [2, 4], the eigenvalues of B are [4, 6], and the eigenvalues of A ∘ B are [14, 16].

[1] For complex matrices, Hermitian means conjugate symmetric, which in the real case reduces to simply symmetric. The theorem of Pólya and Szegö is actually valid for Hermitian matrices, but I simplified the statement for the case of real-valued matrices.

Prime interruption

Posted on 8 October 2018 by John

Suppose you have a number that you believe to be prime. You start reading your number aloud, and someone interrupts “Stop right there! No prime starts with the digits you’ve read so far.”

It turns out the person interrupting you shouldn’t be so sure. There are no restrictions on the digits a prime number can begin with. (Ending digits are another matter. No prime ends in 0, for example.) Said another way, given any sequence of digits, it’s possible to add more digits to the end and make a prime. [1]

For example, take today’s date in ISO format: 20181008. Obviously not a prime. Can we find digits to add to make it into a prime? Yes, we can add 03 on to the end because 2018100803 is prime.

What about my work phone number: 83242286846? Yes, just add a 9 on the end because 832422868469 is prime.

This works in any base you’d like. For example, the hexadecimal number CAFEBABE is not prime, but CAFEBABE1 is. Or if you interpret SQUEAMISH as a base 36 number, you can form a base 36 prime by sticking a T on the end. [2]

In each of these example, we haven’t had to add much on the end to form a prime. You can show that this is to be expected from the distribution of prime numbers.

[1] Source: R. S. Bird. Integers with Given Initial Digits. The American Mathematical Monthly, Vol. 79, No. 4 (Apr., 1972), pp. 367-370

[2] CAFEBABE is a magic number used at the beginning of Java bytecode files. The word “squeamish” here is an homage to “The Magic Words are Squeamish Ossifrage,” an example cleartext used by the inventors of RSA encryption.

Probability of winning the World Series

Posted on 6 October 2018 by John

Astros win 2018 World Series

Suppose you have two baseball teams, A and B, playing in the World Series. If you like, say A stands for Houston Astros and B for Milwaukee Brewers. Suppose that in each game the probability that A wins is p, and the probability of A losing is q = 1 – p. What is the probability that A will win the series?

The World Series is a best-of-seven series, so the first team to win 4 games wins the series. Once one team wins four games there’s no point in playing the rest of the games because the series winner has been determined.

At least four games will be played, so if you win the series, you win on the 4th, 5th, 6th, or 7th game.

The probability of A winning the series after the fourth game is simply p⁴.

The probability of A winning after the fifth game is 4 p⁴ q because A must have lost one game, and it could be any one of the first four games.

The probability of A winning after the sixth game is 10 p⁴ q² because A must have lost two of the first five games, and there are 10 ways to choose two items from a set of five.

Finally, the probability of A winning after the seventh game is 20 p⁴ q³ because A must have lost three of the first six games, and there are 20 ways to choose three items from a set of six.

The probability of winning the World Series is the sum of the probabilities of winning after 4, 5, 6, and 7 games which is

p⁴(1 + 4q + 10q² + 20q³)

Here’s a plot:

Probability of winning the world series as a function of the probability of winning one game

Obviously, the more likely you are to win each game, the more likely you are to win the series. But it’s not a straight line because the better team is more likely to win the series than to win any given game.

Now if you only wanted to compute the probability of winning the series, not the probability of winning after different numbers of games, you could pretend that all the games are played, even though some may be unnecessary to determine the winner. Then we compute the probability that a Binomial(7, p) random variable takes on a value greater than or equal to 4, which is

35p⁴q³ + 21p⁵q² + 7p⁶q + p⁷

While looks very different than the expression we worked out above, they’re actually the same. If you stick in (1 – p) for q and work everything out, you’ll see they’re the same.

Physical constants in Python

Posted on 4 October 2018 by John

You can find a large collection of physical constants in scipy.constants. The most frequently used constants are available directly, and hundreds more are in a dictionary physical_constants.

The fine structure constant α is defined as a function of other physical constants:

$\alpha = \frac{e^2}{4 \pi \varepsilon_0 \hbar c}$

The following code shows that the fine structure constant and the other constants that go into it are available in scipy.constants.

    import scipy.constants as sc

    a = sc.elementary_charge**2
    b = 4 * sc.pi * sc.epsilon_0 * sc.hbar * sc.c
    assert( abs(a/b - sc.fine_structure) < 1e-12 )

Eddington’s constant

In the 1930’s Arthur Eddington believed that the number of protons in the observable universe was exactly the Eddington number

$N_{\mathrm{Edd}} = \frac{2^{256}}{\alpha}$

Since at the time the fine structure constant was thought to be 1/136, this made the number of protons a nice even 136 × 2²⁵⁶. Later he revised his number when it looked like the fine structure constant was 1/137. According to the Python code above, the current estimate is more like 1/137.036.

Eddington was a very accomplished scientist, though he had some ideas that seem odd today. His number is a not a bad estimate, though nobody believes it could be exact.

The constants in scipy.constants have come up in a couple previous blog posts.

The post on Koide’s coincidence shows how to use the physical_constants dictionary, which includes not just the physical constant values but also their units and uncertainty.

The post on Benford’s law shows that the leading digits of the constants in scipy.constants follow the logarithmic distribution observed by Frank Benford (and earlier by Simon Newcomb).

Technological context

Posted on 3 October 2018 by John

As you read this blog post, you probably have a hierarchy of contexts in the back of your mind. It comes so naturally to you that you don’t consciously think of it.

If you’re reading this in a web browser, you probably know what browser you’re using. If not, you’re at least aware that you are using a browser, even if you forget momentarily which one you have open. And you probably know what operating system is hosting your browser. You understand that you are reading a page on my site, that this page is not your browser, but content hosted by your browser. If you’ve subscribed via RSS or email, you know what email or RSS client you’re using and understand how this post is organized with respect to your other content.

Some people do not have this kind of context. Anything on their computer is simply “The Computer.” They don’t really understand what an operating system, web browser, or email client are. And they don’t need to know, most of the time. They can get their work done, but then occasionally they have inexplicable problems.

I’m not saying this to criticize or make fun of the people who don’t have the kind of technological context I’m talking about. It’s a remarkable achievement that software has gotten so easy to use that people can get along without knowing much technological context. But if this sort of thing is second nature to you, you might have a hard time understanding how a large number of people work.

You probably take it for granted that you can access the same web site from different computers. Some people do not. Their desktop at work is one computer, and their iPhone is a different computer. They don’t really understand what a web site is.

I know what a web browser is because I have been using computers since before there was a web. Old timers know what various technological components are because they’ve seen them develop. And “digital natives” know to get things done because they’ve grown up with computers, though their gaps in context show occasionally. Seems like the people in the middle would have the hardest time, not having grown up with the technology but not having watched it develop either.

I’m writing this because I’m becoming increasingly aware of how difficult life can be for people who don’t have adequate mental models for technology. I imagine most of my readers are tech savvy, and may have a hard time seeing some of the same things that I’ve had a hard time seeing, that a lot of people don’t understand things we take for granted.

Source http://dilbert.com/strip/1995-06-24

It used to be that anybody who used a computer had to know some basic things. If you were a Unix user a generation ago, you might not know anything about the internals of Unix, but you at least knew that you were a Unix user. There were wizards and mere mortals, but the two groups shared more context than the most tech savvy and least tech savvy share today.

It’s good that people don’t need to know as much context, but occasionally it produces bewildering situations, both for the user and the person trying to help them.

Interlaced zeros of ODEs

Posted on 2 October 2018 by John

Sturm’s separation theorem says that the zeros of independent solutions to an equation of the form

alternate. That is, between any two consecutive zeros of one solution, there is exactly one zero of the other solution. This is an important theorem because a lot of differential equations of this form come up in applications.

If we let p(x) = 0 and q(x) = 1, then sin(x) and cos(x) are independent solutions and we know that their zeros interlace. The zeros of sin(x) are of the form nπ, and the zeros of cos(x) are multiples of (n + 1/2)π.

What’s less obvious is that if we take two different linear combinations of sine and cosine, as long as they’re not proportional, then their zeros interlace as well. For example, we could take f(x) = 3 sin(x) + 5 cos(x) and g(x) = 7 sin(x) – 11 cos(x). These are also linearly independent solutions to the same differential equation, and so the Sturm separation theorem says their roots have to interlace.

If we take p(x) = 1/x and q(x) = 1 – (ν/x)² then our differential equation becomes Bessel’s equation, and the Bessel functions J_ν and Y_ν are independent solutions. Here’s a little Python code to show how the zeros alternate when ν = 3.

    import matplotlib.pyplot as plt
    from scipy import linspace
    from scipy.special import jn, yn

    x = linspace(4, 30, 100)
    plt.plot(x, jn(3, x), "-")
    plt.plot(x, yn(3, x), "-.")
    plt.legend(["$J_3$", "$Y_3$"])
    plt.axhline(y=0,linewidth=1, color="k")
    plt.show()

Plotting Bessel functions J_3 and Y_3

Thirty days hath September …

Posted on 1 October 2018 by John

Here’s today’s exponential sum, 10/1/2018:

And here’s what it would look like if September had 31 days, 9/31/2018:

Solving hard problems

PHI and offshore data processing

Related

Fundamental units

Milk and wine

Integer odds and prime numbers

Related posts

My most popular posts on Reddit

Comparing trig functions and Jacobi functions

Related basic Jacobi functions to trig functions

The rest of the Jacobi functions

Ratios

Zeros and poles

Identities

Derivatives

Clearing up the confusion around Jacobi functions

Modulus, parameter, and modular angle

Quarter periods

Amplitude

Jacobi elliptic functions in Mathematica

Jacobi elliptic functions in Python

Related posts

Hadamard product

Python implementation

Related posts

Prime interruption

Related posts

Probability of winning the World Series

Physical constants in Python

Eddington’s constant

Related posts

Technological context

Interlaced zeros of ODEs

Related posts

Thirty days hath September …