Bootstrap Confidence Intervals for XGBoost regression (Python)

I want to construct Bootstrap Confidence Intervals for XGBoost regression using python. I developed my case based on codes (https://machinelearningmastery.com/calculate-bootstrap-confidence-intervals-machine-learning-results-python/#comment-528118).
Question: I am getting a one bin histogram. I get the single value for the score when we do n_iterations for the bootstrap. This is the problem and it is related to the way I am getting RMSE. Though I tried to find RMSE in different ways. yet, I could not solve the problem How can we solve it?


import numpy
from pandas import read_csv
from sklearn.datasets import load_boston
from sklearn.utils import resample
from matplotlib import pyplot
from xgboost import XGBRegressor
import pandas as pd
import numpy as np
from sklearn.metrics import mean_squared_error
# load dataset
boston_dataset = load_boston()

df = pd.DataFrame(boston_dataset.data, columns=boston_dataset.feature_names)

df['MEDV'] = boston_dataset.target
values1 = df.values
# configure bootstrap
n_iterations = 1000
n_size = int(len(df) * 0.50)
# run bootstrap
stats = list()

# prepare train and test sets

for i in range(n_iterations):
# prepare train and test sets
train = resample(values1, n_samples=n_size)
test = numpy.array([x for x in values1 if x.tolist() not in train.tolist()])

model = XGBRegressor() ## Final for the papers

X_train = train[:,:-1]
y_train = train[:,-1]
X_test = test[:,:-1]
y_test = test[:,-1]

model.fit(X_train,y_train)
predictions = model.predict(X_test)
# make predictions

def rmse_calculator(predicted, actual):

assert len(predicted) == len(actual)
return np.sqrt(
np.mean(
np.power(predicted- actual, 2)))
score = rmse_calculator(y_test , predictions)

#score = mean_squared_error(y_test, predictions) ** 0.5
yt = np.asarray(y_test)
y_pred = np.asarray(predictions)
score = np.sqrt(mean_squared_error(yt,y_pred))
print(score)
stats.append(score)
# plot scores
pyplot.hist(stats)
pyplot.show()
# confidence intervals
alpha = 0.95
p = ((1.0-alpha)/2.0) * 100
lower = max(0.0, numpy.percentile(stats, p))
p = (alpha+((1.0-alpha)/2.0)) * 100
upper = min(1.0, numpy.percentile(stats, p))
print('%.1f confidence interval %.1f%% and %.1f%%' % (alpha*100, lower*100, upper*100))