後輩たちに正規分布を教えるための正規分布グラフ作成コード

実験TAや後輩の質問対応をしているときに測定と正規分布と3σの関係についてよく聞かれる。教える際には是非確率分布と積分した面積(%換算したやつ)が欲しくなる。しかし3σ以上の確率分布があり、なおかつ確率等がグラフに書かれていて著作フリーなものが見当たらなかったため自作した。後輩の世話をする人は是非使用されたし。
これは以前卒論で確率計算するためのコードだったものを基にして、細かい機能を省き画像を出力できるようにしただけのコードなので細かいところは気にしないでほしい。動くはず。

import numpy as np
import matplotlib.pyplot as plt
import math
import scipy
from scipy import integrate

#######def const######
mu = 10
sigma = 1
# const = 10
fig_range = 4 #(times sigma)

#######cal const######
fig_min = mu - sigma * fig_range
fig_max = mu + sigma * fig_range
const = 1 / math.sqrt(2 * math.pi * sigma **2 )
####################

def gauss_func(x, const, mu, sigma):
    return const * np.exp(- (x - mu) ** 2 / (2 * sigma ** 2))

def int_gaus(rmin, rmax):
    def gf(x):
        return const * np.exp(- (x - mu) ** 2 / (2 * sigma ** 2))
    P = integrate.quad(gf, rmin, rmax)
    return P

def vartical_line(how_times):
    Xs = []
    Ys = []
    Xs.append(mu + how_times * sigma)
    Xs.append(mu + how_times * sigma)
    Ys.append(-0.05 * const)
    Ys.append(gauss_func(mu+how_times*sigma, const, mu, sigma))
    return Xs, Ys


if __name__ == "__main__":
    
    for k in range(1, fig_range+1):
        plt.figure(figsize=(12,9))
        plt.rcParams["font.size"] = 15
        xd = np.arange(fig_min, fig_max, 0.1)
        yd = gauss_func(xd, const, mu, sigma)
        plt.plot(xd, yd, c='b')
        plt.plot([fig_min,fig_max],[0,0],c='k')
        for i in range(1 - fig_range, fig_range):
            Xs, Ys = vartical_line(i)
            plt.plot(Xs, Ys, c='k')
            if i<0:
                plt.text(Xs[0], Ys[0], "μ"+str(i)+"σ", ha='center', va='top')
            if i==0:
                plt.text(Xs[0], Ys[0], "μ", ha='center', va='top')
            if i>0:
                plt.text(Xs[0], Ys[0], "μ+"+str(i)+"σ", ha='center', va='top')
        x_range =  np.arange(mu - k * sigma, mu + k * sigma,0.1)
        x_range = np.append(x_range, mu + k * sigma)
        y_range = gauss_func(x_range, const, mu, sigma)
        
        print(str(k)+"σ = \nP 1-P (1-P)/2")
        P, dP = int_gaus(min(x_range), max(x_range))
        plt.rcParams["font.size"] = 20
        print(P, 1-P,(1-P)/2)
        plt.text(mu,const,str(k)+"σ = "+'{:.5f}'.format(100*P)+"%",ha='center', va='bottom')
        plt.fill_between(x_range, y_range,facecolor='y',alpha=0.5)
        plt.savefig("gauss_"+str(k)+".jpg")
        plt.close()
    print("F!!!!")

お世辞にもきれいなコードではないので改善点を見つけた場合報告いただけるとありがたいです。
以下出力画像

卒論(天文学系)では7σの計算が必要だったためコマンドライン出力は大丈夫だが、画像は小数点以下の出力を増やさないと概数で100%になってしまう。

        plt.text(mu,const,str(k)+"σ = "+'{:.5f}'.format(100*P)+"%",ha='center', va='bottom')

の行を変えて対応すべし。