import numpy as np
import matplotlib.pyplot as plt
import scipy as sp

X = sp.stats.expon(scale = 1_000)

N = 20
data = X.rvs(size = N)
data

array([ 611.94994999,   32.73060543,  260.90819498,  973.65795288,
       1632.10218268, 4381.7841653 ,  835.46812813,  528.23626679,
        254.15655304,  441.6069858 , 2919.16986341, 2251.14528326,
        235.38532241,  162.70802574,  444.60070433,  557.57739079,
         76.06937215,  154.1181732 ,   17.14743599,  876.32466646])

mhat = np.mean(data)

mhat

np.float64(882.3423611380082)

T = sp.stats.t(df = N - 1) # df = degrees of freedom

t = T.ppf(0.975)

s = np.std(data)

mhat + t * s / np.sqrt(N)

np.float64(1393.9459800414802)

mhat - t * s / np.sqrt(N)

np.float64(370.73874223453623)

R = 300
b = np.zeros(R)
for r in range(R):
    data = X.rvs(size = N)
    mhat = np.mean(data)
    s = np.std(data)
    ub = mhat + t * s / np.sqrt(N)
    lb = mhat - t * s / np.sqrt(N)
    if lb <= X.mean() <= ub:
        b[r] = 1
    else:
        print(lb, ub)

453.9537556994002 850.4374755559503
472.78925179802263 969.2202032215766
401.3036839198024 922.9629909214336
497.56431925865024 968.0345750398307
345.2210743336254 795.058582518935
524.2389366909345 927.9887092515028
419.0184087085173 946.1539418458832
358.94499283345345 939.7740992441109
439.35808400766615 939.9633577953923
436.5919885338835 988.4433657061343
425.4287665070773 979.3029945668578
416.87828016667794 999.6568330626903
1008.0904235724335 1617.9305335847307
390.55419146578095 932.9463894164827
312.918882779649 971.5254379461342
354.1858152185465 710.0067431749425
283.3343873631666 811.1640287979503
257.1346110475754 795.2412205752319
395.5975749528006 913.6656017393395
447.5933094835603 859.0398166196821
342.4168714980951 810.1722189496563
381.0829847389331 958.4878554860232
430.53656333175303 953.1713122605302
431.3415387888875 992.1820341479252
1022.2790818797441 1684.2586352978108
301.1602963763309 988.2192654477269
1084.9189160480055 2325.7594694117774
390.74720430112006 935.6306460156641
437.11023062346624 919.504289684404

np.mean(b)

np.float64(0.9129333333333334)