The final program:

import statistics


def remove_outliers(vec: list[float]) -> list[float]:
    mu = statistics.fmean(vec)
    sigma = statistics.stdev(vec)
    normed = [(mu - v) / sigma for vin vec]
    return [v for v, n in zip(vec, normed) if -2 <= n <= 2]

Using stdev means that the input data is a smaller sample from a wider population. If the data is the whole population, use pstdev (population standard deviation).