Probability that notes are in the same scale

Question

What is the probability that a random sequence of notes (on the 12-note chromatic scale) of length n is in the same major scale?

Quite some time ago, I came across a song written by converting pi to base-12. (The song didn’t go on forever, obviously, but it was really cool.) The guy who wrote the song claimed that pi is a “musical number”.

I wanted to make it a project to figure out whether pi is truly a musical number. That is, if pi in base-12 begins with a string of digits that, when converted to musical notes, stay in the same scale longer than statistically expected.

This is a big topic depending on what types of scale(s) one could choose and what key transitions would be allowable, but I thought major would be a good place to start.

Haven’t gotten very far mathematically, mainly because I can’t figure out exactly what to do or whether this is an easy or difficult problem.

Edit: Since I saw this in the comments—which scale is not specified in advance. So you have the notes, and then you see if they fit ANY scale.

Answer 1

Having run a simulation, I am finding an average run of $4.69086$ notes all able to be said to stay in the same major key having run over 10 million trials... with a standard deviation of around $2.327$. That said, the initial string of digits of $\pi$ in base 12 is not significantly more than the average.

Here is my javascript code. Note, I use a few custom methods for looping over objects and cloning and removing. It should be clear what they do

//Generate the scales
majorscale = [0,2,4,5,7,9,11];
majorscales = [];
for (i=0; i<12;i++){cur = []; 
  majorscale.$each(function(i){cur.push((i+1)%12)}); 
  majorscale = cur; majorscales.push(cur)}

//Define a function to generate a random song, length 25 is probably fine
generateSong = function(){song = []; for(i=0;i<25; i++){song.push(Math.floor(Math.random()*12))}; 
  return song}

//Define a function to find the length of a run in the same key
testLength = function(song,dbug){if (dbug){console.log(song)}; 
  //Scales not yet disqualified stored in remainingScales
  remainingScales = majorscales.$clone();  
  len = 0; 
  //Loop until all scales disqualified
  while(remainingScales.length){
    clone = remainingScales.$clone(); 
    if(dbug){console.log(song[len]); 
      console.table(clone)} 
    remainingScales.$each(function(scale){if(!scale.contains(song[len])){
      //If current note in song does not appear in a remaining scale, disqualify it
      //modify the clone so as to not mess up where we are in the loop and skip a scale
      clone.$remove(scale)}}); 
    remainingScales = clone; 
    len=len+1} 
  return len-1} //Off by one

//Run simulation
tot=0; for(i=0; i<10000000; i++){tot = tot + testLength(generateSong())};
console.log("Average: " + (tot/10000000))
tot2=0; for(i=0; i<10000000; i++){tot2 = tot2 + (testLength(generateSong())-tot/10000000)**2};
console.log("Std dev: " + Math.sqrt(tot2/10000000))

Answer 2

Number the 12 notes, starting from Middle C = $0$, $(0,1,2,3,4,5,6,7,8,9,10,11)$, and consider mod 12. Then the notes in C Major are $(0,2,4,5,7,9,11)$, D Major are $(2,4,6,7,9,11,1)$, etc for all 12 Major Scales.

Represent the notes in each scale by a vector of length 12 with $1$s and $0$s. Where the $i$th value is $1$ if note $i-1$ is contained in that scale and $0$ otherwise. Eg, C Major has the representation $(1,0,1,0,1,1,0,1,0,1,0,1)$.

Consider the bijection $(0,1,2,3,4,5,6,7,8,9,10,11) \to (5,0,7,2,9,4,11,6,1,8,3,10)$. This sends all 12 Major scales to a representation where the $1s$ and $0s$ are consecutive (allowing wrapping at the ends). Eg, C Major has the transformed representation $(1,1,1,1,1,1,1,0,0,0,0,0)$, C Sharp Major has the transformed representation $(1,1,0,0,0,0,0,1,1,1,1,1)$.

Therefore the probability that a set of $n$ notes are all in the same major scale is then equal to the probability that $n$ numbers selected with replacement from the numbers 1 to 12 have a consecutive gap of at least 5 numbers, allowing wrapping at the ends.

---

To find this probability we use Markov Chains. WLOG, let the first note played be Middle C, and we enter the state 1, with representation $(1,0,0,0,0,0,0,0,0,0,0,0)$ indicating that note 0 has been hit and no other notes have been hit. At the next note we either stay in this state with probability 1/12, or move to another state where a different note is hit, each with probability 1/12. With each note we continue to either stand in the same state, with probability of the sum of that state's representation / 12, or move to another state by hitting a new note. If ever we reach a state where there is not a sequence of 5 consecutive $0$'s we instead go to the absorbing state of $(1,1,1,1,1,1,1,1,1,1,1,1)$. There are 256 valid states. Denote the $256\times256$ matrix of transition probabilities as $T$. The probability of hitting the absorbing state on the nth note is the vector $(1,0,0,\ldots) \times T^{n-1}$. (Recall that we are assuming the first note hit is Middle C)

The probabilities of hitting the absorbing state on the nth note and by at least the nth note are (rounded to 3 decimal places, code to calculate is below).

n	P(X=n)	P(X <= n)
1	0.000	0.000
2	0.000	0.000
3	0.139	0.139
4	0.226	0.365
5	0.201	0.565
6	0.150	0.715
7	0.103	0.818
8	0.068	0.886
9	0.044	0.930
10	0.027	0.957
11	0.017	0.974
12	0.010	0.984
13	0.006	0.991
14	0.004	0.994
15	0.002	0.997

Note that this is the time to hit a note that is not in the major key, so we must subtract 1 to find the number of notes that can be played within the major key.

The mean number of notes to hit the absorbing state is $5.690909$. Subtracting 1 from this gives $4.690909$ notes before hitting a note which is not in the major key. This is in excellent agreement with @JMoravitz' simulation

---

R Code:

library(tidyverse)
# Construct list of scales and show bijection works
CMajor = c(0,2,4,5,7,9,11)
majors = matrix(0,nrow=12,ncol=12)
for(i in 0:11) {
  for(j in 0:11) {
    if(j %in% ((CMajor+i)%%12)) {
      majors[i+1,j+1] <- 1
    }
  }
}
perm_mat <- matrix(data = c(
  0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
  0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
  0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
  1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
  0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
  0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
  0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0),
  nrow = 12, ncol = 12, byrow = TRUE)
permuted_maj <- majors %*% perm_mat
print(permuted_maj)
print(matrix(1:12, nrow = 1,ncol = 12) %*% perm_mat)

# construct a list of valid states ------------------------
# 12 x T will be the absorbing state.
# Any time there is no longer a string of 5+
# Falses, it will go to the absorbing state,
# WLOG, let the first hit strike position 1

# Create the list of states and simplify it
states <- expand.grid(
  x1 = c(T, F), x2 = c(T, F), x3 = c(T, F), x4 = c(T, F),
  x5 = c(T, F), x6 = c(T, F), x7 = c(T, F), x8 = c(T, F),
  x9 = c(T, F), x10 = c(T, F), x11 = c(T, F), x12 = c(T, F))
states <- states %>%
  subset(x1 == T)
states$rowsum <- rowSums(states)
states <- states %>%
  filter((rowsum < 8)|(rowsum == 12))
states <- as.matrix(states)
for(i in nrow(states):2) {
  drop <- TRUE
  for(j in 1:12) {
    if(sum(states[i,(((j:(j+4)%%12)+1))]) == 0) {
      drop <- FALSE
    }
  }
  if(drop == TRUE) {
    states <- states[-c(i),]
  }
}
states <- states[order(
  states[, 13],
  states[, 1], states[, 2], states[, 3], states[, 4],
  states[, 5], states[, 6], states[, 7], states[, 8],
  states[, 9], states[, 10], states[, 11], states[, 12]), ]
print(head(states))
print(tail(states))

# construction transition matrix --------------------------------
transition <- matrix(data = 0, nrow = 256, ncol = 256)
for(i in 1:nrow(transition)) {
  rowcount <- states[i,13]
  transition[i,i] <- states[i,13]/12
  for(j in 1:nrow(transition)) {
    if(i == j) next
    if((states[i,13]+1) == states[j,13]) {
      if(sum(xor(states[i,c(1:12)],states[j,c(1:12)]))==1) {
        rowcount <- rowcount + 1
        transition[i,j] <- 1/12
      }
    }
    if(i != 256) transition[i,256] <- (12-rowcount)/12
  }
}
# Calculate probabilities --------------------------------------
probs <- matrix(data = c(1,rep(0,255)),nrow = 1, ncol = 256) %*% 
  transition
cum_hits <- c(0,probs[256])
for(i in 1:50) {
  probs <- probs %*% transition
  cum_hits <- c(cum_hits,probs[256])
}
hits <- rep(0,length(cum_hits))
for(i in 2:length(hits)) {
  hits[i] <- cum_hits[i]-cum_hits[i-1]
}
print(round(cum_hits,3))
print(round(hits,3))
print(sum(hits*1:length(hits)))