If I am understanding you correctly, you might like the np.percentile function. To get the 95% credible interval, I often do:
lower = np.percentile(pred_samples['f_pred'], 2.5, axis=0)
upper = np.percentile(pred_samples['f_pred'], 97.5, axis=0)
Like in @cluhmann 's earlier post, the axis keyword will compute the percentile across the samples for each data point, so you should get N lower and upper values, representing N 95% credible intervals. There shouldn’t be a need for a for loop, this should work regardless of the number of data points.