Range trick for larger sample sizes

Posted on by John

The previous post showed that the standard deviation of a sample of size n can be well estimated by multiplying the sample range by a constant dn that depends on n.

The method works well for relatively small n. This should sound strange: typically statistical methods work better for large samples, not small samples. And indeed this method would work better for larger samples. However, we’re not so much interested in the efficiency of the method per se, but its efficiency relative to the standard way of estimating standard deviation. For small samples, both methods are not very accurate, and the two methods appear to work equally well.

If we want to use the range method for larger n, there are a couple questions: how well does the method work, and how do we calculate dn.

Simple extension

Ashley Kanter left a comment on the previous post saying

dn = 3 log n0.75 (where the log is base 10) seems to perform quite well even for larger n.

Ashley doesn’t say where this came from. Maybe it’s an empirical fit.

The constant dn is the expected value of the range of n samples from a standard normal random variable. You could find a (complicated) expression for this and then find a simpler expression using an asymptotic approximation. Maybe I’ll try that later, but for now I need to wrap up this post and move on to client work.

Note that

log10 n0.75 = 0.977 loge(n)

and so we could use

dn = log n

where log is natural log. This seems like something that might fall out of an asymptotic approximation. Maybe Ashley empirically discovered the first term of a series approximation.

Update: See this post for a more detailed exploration of how well log n, square root of n, and another method approximate dn.

Update 2: Ashley Kanter’s approximation was supposed to be 3 (log10 n) 0.75 rather than 3 log10 (n0.75) and is a very good approximation. This is also addressed in the link in the first update.

Simulation

Here’s some Python code to try things out.

    import numpy as np
    from scipy.stats import norm

    np.random.seed(20220309)

    n = 20
    for _ in range(5):
        x = norm.rvs(size=n)
        w = x.max() - x.min()
        print(x.std(ddof=1), w/np.log(n))

And here are the results.

    |   std | w/d_n |
    |-------+-------|
    | 0.930 | 1.340 |
    | 0.919 | 1.104 |
    | 0.999 | 1.270 |
    | 0.735 | 0.963 |
    | 0.956 | 1.175 |

It seems the range method has more variance, though notice in the fourth row that the standard estimate can occasionally wander pretty far from the theoretical value of 1 as well.

We get similar results when we increase n to 50.

    |   std | w/d_n |
    |-------+-------|
    | 0.926 | 1.077 |
    | 0.889 | 1.001 |
    | 0.982 | 1.276 |
    | 1.038 | 1.340 |
    | 1.025 | 1.209 |

Not-so-simple extension

There are ways to improve the range method, if by “improve” you mean make it more accurate. One way is to divide the sample into random partitions, apply the method to each partition, and average the results. If you’re going to do this, partitions of size 8 are optimal [1]. However, the main benefit of the range method [2] is its simplicity.

Related posts

[1] F. E. Grubbs and C. L. Weaver (1947). The best unbiased estimate of a population standard deviation based on group ranges. Journal of the American Statistical Association 42, pp 224–41

[2] The main advantage now is its simplicity, When it was discovered, the method reduced manual calculation, and so it could have been worthwhile to make the method a little more complicated as long as the calculation effort was still less than that of the standard method.

Estimating standard deviation from range

Posted on by John

Suppose you have a small number of samples, say between 2 and 10, and you’d like to estimate the standard deviation σ of the population these samples came from. Of course you could compute the sample standard deviation, but there is a simple and robust alternative.

Let W be the range of our samples, the difference between the largest and smallest value. Think “w” for “width.” Then

W / dn

is an unbiased estimator of σ where the constants dn can be looked up in a table [1].

    |  n | 1/d_n |
    |----+-------|
    |  2 | 0.886 |
    |  3 | 0.591 |
    |  4 | 0.486 |
    |  5 | 0.430 |
    |  6 | 0.395 |
    |  7 | 0.370 |
    |  8 | 0.351 |
    |  9 | 0.337 |
    | 10 | 0.325 |

The values dn in the table were calculated from the expected value of W/σ for normal random variables, but the method may be used on data that do not come from a normal distribution.

Let’s try this out with a little Python code. First we’ll take samples from a standard normal distribution, so the population standard deviation is 1. We’ll draw five samples, and estimate the standard deviation two ways: by the method above and by the sample standard deviation.

    from scipy.stats import norm, gamma

    for _ in range(5):
        x = norm.rvs(size=10)
        w = x.max() - x.min()
        print(x.std(ddof=1), w*0.325)

Here’s the output:

    | w/d_n |   std |
    |-------+-------|
    | 1.174 | 1.434 |
    | 1.205 | 1.480 |
    | 1.173 | 0.987 |
    | 1.154 | 1.277 |
    | 0.921 | 1.083 |

Just from this example it seems the range method does about as well as the sample standard deviation.

For a non-normal example, let’s repeat our exercise using a gamma distribution with shape 4, which has standard deviation 2.

    | w/d_n |   std |
    |-------+-------|
    | 2.009 | 1.827 |
    | 1.474 | 1.416 |
    | 1.898 | 2.032 |
    | 2.346 | 2.252 |
    | 2.566 | 2.213 |

Once again, it seems both methods do about equally well. In both examples the uncertainty due to the small sample size is more important than the difference between the two methods.

Update: To calculate dn for other values of n, see this post.

[1] Source: H, A. David. Order Statistics. John Wiley and Sons, 1970.

Varicode

Posted on by John

Varicode is a way of encoding text and control characters into binary using code words of variable length. It was developed as part of the PSK31 protocol for digital communication over amateur radio.

In the spirit of Morse code, it uses short code words for common characters and longer code words for less common characters in the expectation that this will result in shorter encodings.

If you use variable length words, you’ve got to have some way of knowing when one word ends and the next begins. Varicode solves this problem by using only keywords that begin and end with 1 and that do not contain two consecutive zeros. Then 00 is inserted between code words. Since 00 cannot appear inside a code word, these bits unambiguously mark the space between code words.

Synchronization

Varicode is self-synchronizing in the sense that if you jump into a stream of bits produced by Varicode, as soon as you see two zeros, you know you’re at a code word boundary and can start reading from there. You’ve lost any bits that were transmitted before you jumped in, but you can parse everything going forward.

ASCII doesn’t have this problem or this robustness. It doesn’t have the problem of determining code word boundaries because every code word is eight bits long. But then to read a stream of ASCII bits you need to know your position mod 8. If you jump into a stream of bits, you don’t know where the next code word boundary will be, though you may be able to infer it by trying 8 possibilities and seeing which produces the most intelligible results.

Regular expression

Another way to state the rules for forming Varicode code words is to say that 1 is a valid code (the code for a space, ASCII 0x20) and that you can form new codes by prefixing 1 or 10 to a code. In terms of regular expressions, this says a Varicode code word matches

    (1|10)*1

Fibonacci numbers

How many code words are there of length n? Well, there are two ways to make a code word that long: you either put a 1 in front of a code word of length n-1, or you put a 10 in front of a code word of length n-2. So the number of code words of length n equals the number of code words of length n-1 plus the number of code words of length n-2. That is, the number of code words satisfies the same recurrence relation as the Fibonacci numbers.

It’s easy to see that there’s only one code word of length one, and only one code word of length two, so the number of code words of each length satisfies the initial conditions for the Fibonacci sequence as well, and so they are the Fibonacci numbers.

Efficiency

Varicode encodes a lot more than lower case letters—it encodes most ASCII characters— and so it would take some work to discover the relative frequencies of the characters, and this frequency would depend on where Varicode is used. As far as I know Varicode is use primarily (only?) in PSK31, and so the relevant frequency would be the frequencies in messages sent over PSK31, not English more generally.

You can find the code words of each letter here.

To make things easier, let’s suppose messages are limited to lower case letters and spaces, and that the letters follow the same distribution as in English in general.

We can use the letter frequencies here, except these don’t take spaces into account. If we assume words are about 5 letters long, then the probability of a character being a space is 1/6 and the probabilities of the other characters conditional on not being a space are given by the table. This means the letter probabilities need to be multiplied by 5/6.

This gives us a an expected length of 3.89 bits per letter, which is effectively 5.89 bits when you consider the 00 pattern we have to add between letters.

You could represent the 26 letters and a few more characters using 5 bits, but the result would not be self-synchronizing. One way of looking at this to say that the compression provided by variable length encoding nearly pays for the overhead required to make the code self-synchronizing.

Related posts

Catenary kiln

Posted on by John

Will Buckner sent me an email with the following question recently. (I’m sharing this with permission.)

I am building a kiln using a catenary arch. The rear wall and front wall/door will be vertical and fill in the space under the arch, which has the dimensions of 41″W x 39.5″H. I need the area within this arch in order to calculate how many bricks I need to construct the two walls. Can you help me calculate that area?

Here’s my solution. I fit a parabola and a catenary in part to check my work and in part to see how different they are.

Parabola

We’ll start with a parabola as our warmup because it’s easier.

Let w = 41 and h = 39.5. We want a quadratic function y(x) such that y(0) = h and y(±w/2) = 0.

Assume y has the form

y(x) = b − ax²

y(0) = h leads to the conclusion that b = h, and y(±w/2) = 0. leads to

a(w/2) = h

And this leads to the final form.

y = h(1 – (2x/w)²)

Catenary

Our catenary will have the form

y(x) = b – a cosh(x/a)

and again y(0) = h and y(±w/2) = 0. These two requirements lead to the equations

b − a = h

and

b − a cosh(w/2a) = 0.

Then

b = h + a

and

h + a (1 − cosh(w/2a) = 0.

The latter equation cannot be solved in closed form, but it’s easy enough to solve numerically. When I asked Mathematica to compute the root with

    FindRoot[39.5 + a - a Cosh[20.5/a] == 0, a]

it failed. Then I plotted the function, saw that the root was around 8.5, then gave Mathematica a hint.

    FindRoot[39.5 + a - a Cosh[20.5/a] == 0, {a, 8.5}]

Then Mathematica came back with

a = 8.475359432102444.

So the equation of our catenary is

y = h + a(1 – cosh(0.5 x/a))

where h and a are given above.

Plots

Here’s a plot of both curves. The blue curve inside is the parabola, and the gold curve on the outside is the catenary.

Here’s the code that made the plot.

    Plot[{h (1 - (2 x/w)^2), h + a - a Cosh[x/a] }, {x, -w/2, w/2}]

Area

And now for the area, the thing I was asked to find. For the parabola

    Integrate[h (1 - (2 x/w)^2), {x, -w/2, w/2}]

returns 1079.67. For the catenary,

    Integrate[h + a (1 -  Cosh[x/a]), {x, -w/2, w/2}]

returns 1166.56.

The number of bricks this will take is roughly the area of the arch divided by the area of a brick, but being more precise is a complicated question that depends on how much mortar you use, and how you use plan to use rectangular bricks to make a curved arch.

More catenary posts

Connecting powers of two and decibels

Posted on by John

Colin Wright pointed out a pattern in my previous post that I hadn’t seen before. He wrote about it here a couple years ago.

Start with the powers of 2 from the top of the post:

21 = 2
22 = 4
23 = 8
24 = 16
25 = 32
26 = 64
27 = 128
28 = 256
29 = 512

and sort the numbers on the right hand side in lexical order, i.e. sort them as you’d sort words. This is almost never what you want to do, but here it is [1].

128
16
2
256
32
4
512
64
8

Next, add a decimal point after the first digit in each number.

1.28
1.6
2
2.56
3.2
4
5.12
6.4
8

Now compare the decibel values listed at the bottom of the previous post.

1.25
1.6
2
2.5
3.2
4
5
6.3
8

Five of the numbers are the same and the remaining 4 are close. And where the numbers differ, the exact decibel value is between the the two approximate values.

Here’s a table to compare the exact value rounded to 3 decimals, the approximation given before, and the value obtained by sorting powers of 2 and moving the decimal.

    |---+-----------+--------+-------|
    | n | 10^(n/10) | approx | power |
    |---+-----------+--------+-------|
    | 1 |     1.259 |   1.25 |  1.28 |
    | 2 |     1.585 |   1.60 |  1.60 |
    | 3 |     1.995 |   2.00 |  2.00 |
    | 4 |     2.512 |   2.50 |  2.56 |
    | 5 |     3.162 |   3.20 |  3.20 |
    | 6 |     3.981 |   4.00 |  4.00 |
    | 7 |     5.012 |   5.00 |  5.12 |
    | 8 |     6.310 |   6.30 |  6.40 |
    | 9 |     7.943 |   8.00 |  8.00 |
    |---+-----------+--------+-------|

So the numbers at the top and bottom of my list are practically the same, but in a different order.

Related post: 2s, 5s, and decibels.

[1] One reason I use ISO dates (YYYY-MM-DD) in my personal work is that that lexical order equals chronological order. Otherwise you can get weird things like December coming before March because D comes before M or because 1 (as in 12) comes before 3. Using year-month-day and padding days and months with zeros as needed eliminates this problem.

100 digits worth memorizing

Posted on by John

I was thinking today about how people memorize many digits of π, and how it would be much more practical to memorize a moderate amount of numbers to low precision.

So suppose instead of memorizing 100 digits of π, you memorized 100 digits of other numbers. What might those numbers be? I decided to take a shot at it. I exclude things that are common knowledge, like multiplication tables up to 12 and familiar constants like the number of hours in a day.

There’s some ambiguity over what constitutes a digit. For example, if you say the speed of light is 300,000 km/s, is that one digit? Five digits? My first thought was to count it as one digit, but then what about Avagadro’s number 6×1023? I decided to write numbers in scientific notation, so the speed of light is 2 digits (3e5 km/s) and Avagadro’s number is 3 digits (6e23).

Here are 40 numbers worth remembering, with a total of 100 digits.

Powers of 2

23 = 8
24 = 16
25 = 32
26 = 64
27 = 128
28 = 256
29 = 512
210 = 1024

Squares

13² = 169
14² = 196
15² = 225

Probability

P(|Z| < 1) ≈ 0.68
P(|Z| < 2) ≈ 0.95

Here Z is a standard normal random variable.

Music

Middle C ≈ 262 Hz
27/12 ≈ 1.5

The second fact says that seven half steps equals one (Pythagorean) fifth.

Mathematical constants

π ≈ 3.14
√2 ≈ 1.414
1/√2 ≈ 0.707
φ ≈ 1.618
loge 10 ≈ 2.3
log10 e ≈ 0.4343
γ ≈ 0.577
e ≈ 2.718

Here φ is the golden ratio and γ is the Euler-Mascheroni constant.

I included √2 and 1/√2 because both come up so often.

Similarly, loge 10 and log10 e are reciprocals, but it’s convenient to know both.

The number of significant figures above is intentionally inconsistent. It’s at least as easy, and maybe easier, to remember √2 as 1.414 than as 1.41. Similarly, if you’re going to memorize that log10 e  is 0.43, you might as well memorize 0.4343. Buy two get two free.

Each constant is truncated before a digit less than 5, so all the figures are correct and correctly rounded. For φ and loge 10 the next digit is a zero, so you get an implicit extra digit of precision.

The requirement that truncation = rounding means that you have to truncate e at either 2.7 or 2.718. If you’re going to memorize the latter, you could memorize six more digits with little extra effort since these digits are repetitive:

e = 2.7 1828 1828

Measurements and unit conversion

c = 3e8 m/s
g = 9.8 m/s²
NA = 6e23
Earth circumference = 4e7m
1 AU = 1.5e8 km
1 inch = 2.54 cm

Floating point

Maximum double = 1.8e308
Epsilon = 2e-16

These numbers could change from system to system, but they rarely do. See Anatomy of a floating point number.

Decibels

1 dB = 100.1 = 1.25
2 dB = 100.2 = 1.6
3 dB = 100.3 = 2
4 dB = 100.4 = 2.5
5 dB = 100.5 = 3.2
6 dB = 100.6 = 4
7 dB = 100.7 = 5
8 dB = 100.8 = 6.3
9 db = 100.9 = 8

These numbers are handy, even if you don’t work with decibels per se.

Update: The next post points out a remarkable pattern between the first and last sets of numbers in this post.

Related posts

Morse code in musical notation

Posted on by John

Maybe this has been done before, but I haven’t seen it: Morse code in musical notation.

Here’s the Morse code alphabet, one letter per measure; in practice there would be less space between letters [1]. A dash is supposed to be three times as long as a dot, so a dot is a sixteenth note and a dash is a dotted eighth note.

Morse code is often at a frequency between 600 and 800 Hz. I picked the E above middle C (660 Hz) because it’s in that range.

Rhythm

Officially a dash is three times as long as a dot. But there’s also a space equal to the length of a dot between parts of a letter. So the sheet music above would be more accurate if you imagined all the sixteenth notes are staccato and the dotted eighth notes are really eighth notes followed by a sixteenth rest.

This doesn’t make much difference because individual operators have varying “fists,” styles of sending Morse code, and won’t exactly follow the official length and spacing rules.

You could rewrite the music above as follows, but it’s all an approximation.

Tempo

According to Wikipedia, “the dit length at 20 words per minute is 50 milliseconds.” So if a sixteenth note has a duration of 50 milliseconds, this would mean five quarter notes per second, or 300 beats per minute. But according to this video, the shortest duration people can distinguish is about 50 milliseconds.

That would imply that copying Morse code at 20 wpm is pushing the limits of human hearing. But copying at 20 wpm is common. Some people can copy Morse code at more than 50 words per minute or more, but at that speed they’re not hearing individual dits and dahs. An H, for example, four dits in a row, sounds like a single rough sound. In fact, they’re not really hearing letters at all but recognizing the shape of words.

How the image was made

I made the image above with LaTeX and Lilypond.

Adding the letters above each measure was kind of a hack. I used rehearsal markings to label the measures, but there was one problem: the software skips from letter H to letter J. That meant that the labels I and all subsequent letters were one ahead of what they should be, and the final letter Z was labeled AA. I tried several tricks, and Lilypond steadfastly refused to label a measure with ‘I’ even though I’ve seen such a label in the documentation.

My way around this was to make it label two consecutive measures with H, then in image editing software I turned the second H into an I. No doubt there’s a better way, but this worked.

I may play around with this and try to improve it a bit. If you have any suggestions, particularly related to Lilypond, please let me know.

Related posts

[1] You could think of the musical score above as a sort of transcription of the Farnsworth method of teaching Morse code. Students learn the letters at full speed, but with extra space between the letters at first. The faster speed discourages consciously counting the dits and dahs, forcing the student to listen to the overall rhythm of the letters.

Find log normal parameters for given mean and variance

Posted on by John

Earlier today I needed to solve for log normal parameters that yield a given mean and variance. I’m going to save the calculation here in case I needed in the future or in case a reader needs it. The derivation is simple, but in the heat of the moment I’d rather look it up and keep going with my train of thought.

NB: The parameters μ and σ² of a log normal distribution are not the mean and variance of the distribution; they are the mean and variance of its log.

If m is the mean and v is the variance then

\begin{align*} m &= \exp(\mu + \sigma^2/2) \\ v &= (\exp(\sigma^2) - 1) \exp(2\mu + \sigma^2) \end{align*}

Notice that the square of the m term matches the second part of the v term.

Then

\frac{v}{m^2} = \exp(\sigma^2) -1

and so

\sigma^2 = \log\left(\frac{v}{m^2} + 1 \right)

and once you have σ² you can find μ by

\mu = \log m - \sigma^2/2

Here’s Python code to implement the above.

    from numpy immport log
    def solve_for_log_normal_parameters(mean, variance):
        sigma2 = log(variance/mean**2 + 1)
        mu = log(mean) - sigma2/2
        return (mu, sigma2)

And here’s a little test code for the code above.

    mean = 3.4
    variance = 5.6

    mu, sigma2 = solve_for_log_normal_parameters(mean, variance)

    X = lognorm(scale=exp(mu), s=sigma2**0.5)
    assert(abs(mean - X.mean()) < 1e-10)
    assert(abs(variance - X.var()) < 1e-10)

Related posts

Q codes in Seveneves

Posted on by John

The first time I heard of Q codes was when reading the novel Seveneves by Neal Stephenson. These are three-letter abbreviations using in Morse code that all begin with Q.

Since Q is always followed by U in native English words, Q can be used to begin a sort of escape sequence [1].

There are dozens of Q codes used in amateur radio [2], and more used in other contexts, but there are only 10 Q codes used in Seveneves [3]. All begin with Q, followed by R, S, or T.

Tree[Q, {Tree[R, {A, K, N, S, T}], Tree[S, {B, L, O}], Tree[T, {H, X}]}]

Each Q code can be used both as a question and as an answer or statement. For example, QRS can mean “Would you like me to slow down” or “Please slow down.” I’ll just give the interrogative forms below.

Here are the 10 codes that appear in Stephenson’s novel.

QRA
What is your call sign?
QRK
Is my signal intelligible?
QRN
Is static a problem?
QRS
Should I slow down?
QRT
Should I stop sending?
QSB
Is my signal fading?
QSL
Are you still there?
QSO
Could you communicate with …?
QTH
Where are you?
QTX
Will you keep your station open for talking with me?

Related posts

[1] Some Q codes have a U as the second letter. I don’t know why—there are plenty of unused TLAs that begin with Q—but it is what it is.

[2] You can find a list here.

[3] There is one non-standard code in the novel: QET for “not on planet Earth.”