У вас есть 2 ошибки.
Сначала вы используете 2 для множителя для CI. Более точное значение составляет 1,96. «2» - это просто удобная оценка. Это делает ваш CI, сгенерированный вручную, слишком толстым.
Во-вторых, вы сравниваете нормальное распределение с t-распределением. Это, вероятно, не вызывает разницы больше, чем десятичная пыль, потому что у вас есть 199 степеней свободы для t-dist, что в основном нормально.
Ниже приведен z-показатель 1,96 и вычисление CI с сравнение яблок с яблоками с нормой распределения по сравнению с т.
In [45]: st.norm.cdf(1.96)
Out[45]: 0.9750021048517795
In [46]: print('(',sample_mean-1.96*standard_error,',',sample_mean+1.96*standard_error,')')
( 57558.007862202685 , 61510.37559873406 )
In [47]: st.norm.interval(0.95, loc=sample_mean, scale=se)
Out[47]: (57558.044175045005, 61510.33928589174)