import numpy as np
import matplotlib.pyplot as plt

line = [0,1]
depth = 4

def divide(line, level=0):
    plt.plot(line,[level,level], color="k", lw=5, solid_capstyle="butt")
    if level < depth:
        s = np.linspace(line[0],line[1],4)
        divide(s[:2], level+1)
        divide(s[2:], level+1)

divide(line)
plt.gca().invert_yaxis()
plt.show()

plt.rcParams['figure.figsize'] = [20, 20]

def devilsstaircase(n,W):
    bin = (0.5)**(1+np.array(range(n)))
    tert = (1/3)**(1+np.array(range(n)))
    A=[]
    B=[]
    for i in range(W):
        Y = np.random.binomial(1,0.5,n)*2
        a = np.sum(Y*tert)
        b = np.sum(0.5*Y*bin)
        A = A + [a]
        B = B + [b]
    return A , B

A, B = devilsstaircase(10000,1000)
A.sort()
B.sort()
plt.plot(A,B)
plt.plot(A,np.zeros(1000),',',color='red')
plt.show()

plt.rcParams['figure.figsize'] = [20, 20]

def devilsstaircase_script(n,W):
    bin = (0.5)**(1+np.array(range(n)))
    tert = (1/3)**(1+np.array(range(n)))
    A=[]
    B=[]
    for i in range(W):
        Y = np.random.randint(0,3,n)
        a = np.sum(Y*tert)
        nY = np.min(np.argwhere(Y==1))
        helper = np.concatenate([np.ones(nY),np.zeros(n-nY)])
        b = np.sum(0.5*Y*helper*bin)+1/2**(nY+1)
        A = A + [a]
        B = B + [b]
    return A , B

A, B = devilsstaircase_script(1000,10000)
A.sort()
B.sort()
plt.plot(A,B)
plt.show()

A,B=devilsstaircase(1000,10000)
A.sort()
B.sort()

C,D=devilsstaircase_script(1000,100)
C.sort()
D.sort()

plt.plot(A,B)
plt.plot(C,D,'p',color='red')
plt.show()

y=np.array([n*np.pi%1 for n in range(100000)])
#plt.plot(y,np.zeros(10000),'p')
#plt.show()
a=0.278
b=0.567
y = y[y > a]
y = y[y < b]
len(y)/100000

0.28902

We have learned how to construct measures on products $ {\{0,1\}}^\mathbb{N} $ and can use to create a little panarama on simple nullsets. Take a Bernoulli law with success probability $ p \in ]0,1[ $ and its infinite product $\mu_p$. Then apparently by the law of large numbers (which is just a statement on orthogonality and a Borel Cantelli type wrap up) $$ \frac{1}{n}\sum_{i=1}^n X_i \to p $$ almost surely. Then the measurable set of sequences with average $ q $ has measure $ 1 $ with respect ot $ \mu_p $ if and only if $ p = q $, hence the infinite product measures $ \mu_q $ are not absolutely continuous to each other.

Considering now the laws on $ [0,1] $ by mapping an infinite sequence of $ 0 $ and $ 1 $ to its binary expansion leads to laws which are not absolutely continuous to each other: nota bene that the law for $ p = 0.5 $ corresponds to the Lebesgue measure on $ [0,1] $.

import numpy as np
import matplotlib.pylab as plt

n= 2000
p =2/3
W = 10**4
Z = []
X = np.array(range(W))/W
bin = (0.5)**(1+np.array(range(n)))
for i in range(W):
    Y = np.random.binomial(1,p,n)
    helper = np.sum(Y*bin)
    Z = Z + [helper]

Z.sort()
plt.plot(X,Z)
plt.show()

Some historical remarks on the theory of integration: please take a look at the bigraphies of, e.g., Georg Cantor, Henri Lebesgue, Constantin Caratheodory, Hans Hahn, Johann Radon, Otto Nikodym, Pavel Urysohn and Stefan Banach, Frigyes Riesz, Tatyana Afanasyeva and Ludwig Boltzmann who have been named several times in the lecture or contributed substantially to the development of the field.

It shows a lot about the fascinating development of measure theory in the first part of the twentieth century against the background of its historical events.