Welch T-Test ("TVLA") Sliders with Python 3 matplotlib. The thick lines are the distributions of averages of N Gaussian variables, which are also Gaussians! When they don't overlap, then the distributions are clearly separate and T value also crosses the +-4.5 line commonly used in side-channel work. This is close to P value 10^-5 or confidence …

tslider.py

#!/usr/bin/env python3

#   tslider.py
#   2021-08-02  Markku-Juhani O. Saarinen <mjos@pqshield.com>
#   Free to play with!

#   === Welch T test sliders (for teaching TVLA statistics)

#   Welch T-Test ("TVLA") Sliders with Python 3 matplotlib. The thick lines
#   are the distributions of averages of N Gaussian variables, which are
#   also Gaussians! When they don't overlap, then the distributions are
#   clearly separate and T value also crosses the +-4.5 line commonly used
#   in side-channel work. This is close to P value 10^-5 or confidence
#   1-p = 99.999%.

import math
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button
import matplotlib.mathtext as mathtext

#   gaussian height at mu

def pdf_h(sig):
    return 1 / (sig * np.sqrt(2*np.pi));

#   normal distribution scaled at 1 at x==mu

def pdf_1(x, mu, sig):
    return np.exp(-0.5*((x - mu)/sig)**2)

#   welch t value

def welch(mu_a, sig_a, n_a, mu_b, sig_b, n_b):
    return (mu_a - mu_b) / np.sqrt( sig_a**2/n_a + sig_b**2/n_b )

#   this is the two-tail p-value (probability of false positve) for normal
#   distribution and also for the t-test with a large degree of freedom (n).

def pvalue_t(t):
    return math.erfc(abs(t)*2**-0.5)

#   ranges

x_min   = -5.0
x_max   = 5.0
sig_max = 5.0
n_max   = 1000

#   initial parameters
mu_a    = -0.6
sig_a   = 2.0
mu_b    = 0.6
sig_b   = 2.0
n_ab    = 100

color_a = '#88FF88'
color_b = '#CCCCFF'
color_n = '#CC88FF'

#   create the plot

xv = np.linspace(x_min, x_max, 500)
y0 = len(xv)*[0.0]

fig, ax = plt.subplots()
line_a, =   plt.plot(xv, y0, lw = 1, color = color_a)
line_an, =  plt.plot(xv, y0, lw = 2, color = color_a)
line_b, =   plt.plot(xv, y0, lw = 1, color = color_b)
line_bn, =  plt.plot(xv, y0, lw = 2, color = color_b)

ax.margins(x = 0)
plt.subplots_adjust(left = 0.15, right = 0.84,bottom = 0.20)
tval = plt.text(0.5, 1.05,  '0', size = 'x-large',
    transform = ax.transAxes, ha = 'center')


#   horizontal slider to control mu (averages)

sl_mu_a = Slider(
    ax      = plt.axes([0.15, 0.09, 0.69, 0.04]),
    label   = r"$\mu_A$",
    valmin  = x_min,
    valmax  = x_max,
    valfmt  = '%4.2f',
    valinit = mu_a,
    color   = color_a,
    initcolor = 'none',
)

sl_mu_b = Slider(
    ax      = plt.axes([0.15, 0.03, 0.69, 0.04]),
    label   = r"$\mu_B$",
    valmin  = x_min,
    valmax  = x_max,
    valfmt  = '%4.2f',
    valinit = mu_b,
    color   = color_b,
    initcolor = 'none',
)

#   vertical slider to control sigma (std. deviation)

sl_sig_a = Slider(
    ax      = plt.axes([0.87, 0.20, 0.04, 0.68]),
    label   = r"$\sigma_A$",
    valmin  = 0.01,
    valmax  = sig_max,
    valfmt  = '%4.2f',
    valinit = sig_a,
    color   = color_a,
    initcolor = 'none',
    orientation = "vertical"
)

sl_sig_b = Slider(
    ax      = plt.axes([0.93, 0.20, 0.04, 0.68]),
    label   = r"$\sigma_B$",
    valmin  = 0.01,
    valmax  = sig_max,
    valfmt  = '%4.2f',
    valinit = sig_b,
    color   = color_b,
    initcolor = 'none',
    orientation = "vertical"
)

#   just one vertical slider to control n (number of samples) for both a and b

sl_n = Slider(
    ax = plt.axes([0.03, 0.20, 0.04, 0.68]),
    label = r"$N$",
    valmin = 10,
    valmax =  n_max,
    valfmt = '%4.0f',
    valinit = n_ab,
    color = color_n,
    initcolor = 'none',
    orientation = "vertical"
)

#   t-value string

def welch_str():

    n_ab    =   sl_n.val
    mu_a    =   sl_mu_a.val
    mu_b    =   sl_mu_b.val
    sig_a   =   sl_sig_a.val
    sig_b   =   sl_sig_b.val

    t = welch(  mu_a, sig_a, n_ab,
                mu_b, sig_b, n_ab);
    p = pvalue_t(t)

    if p < 1E-5:
        s = f't = {t:.4f}  p = 10$^{{{math.log10(p):.2f}}}$'
    else:
        s = f't = {t:.4f}  p = {p:.6f}'

    tval.set_text(s)

    #   magical t-value 4.5, corresponding roughly to p=10^-5
    if t > -4.5 and t < 4.5:
        tval.set_color('black')
    else:
        tval.set_color('red')


#   update the canvas

def update(val):

    n_ab    =   sl_n.val
    mu_a    =   sl_mu_a.val
    mu_b    =   sl_mu_b.val
    sig_a   =   sl_sig_a.val
    sig_b   =   sl_sig_b.val

    n_sc    =   1/np.sqrt(n_ab)
    nsig_a  =   n_sc * sig_a
    nsig_b  =   n_sc * sig_b

    yv_a    = pdf_h(sig_a) * pdf_1(xv, mu_a, sig_a)
    ny_a    = pdf_h(nsig_a) * pdf_1(xv, mu_a, nsig_a)

    yv_b    = pdf_h(sig_b) * pdf_1(xv, mu_b, sig_b)
    ny_b    = pdf_h(nsig_b) * pdf_1(xv, mu_b, nsig_b)

    ax.set_ylim(0, max(pdf_h(nsig_a), pdf_h(nsig_b)))

    line_a.set_ydata(yv_a)
    line_an.set_ydata(ny_a)
    line_b.set_ydata(yv_b)
    line_bn.set_ydata(ny_b)

    welch_str()
    fig.canvas.draw_idle()

#   register the update function with each slider
sl_mu_a.on_changed(update)
sl_sig_a.on_changed(update)
sl_mu_b.on_changed(update)
sl_sig_b.on_changed(update)
sl_n.on_changed(update)

#   initial
update(0)

plt.show()