Why does my probability plot not agree with my histogram? I am trying to determine the best distribution to characterize a dataset. My dataset visually looks like Cauchy distribution would be a good fit, so I used python to fit a Cauchy distribution to my dataset. Visually, with the fitted Cauchy distribution plotted over the histogram of my data, it looks to fit well. However, the probability plot says otherwise. How could these two plots disagree so much? My coefficient of determination is 0.04. The top plot is sample data (y-axis) against the quantiles of the Cauchy distribution (x-axis), and the bottom plot is the probability density histogram with the x-axis being the values in my dataset.

My "data" is a NumPy array of ~2 million data points imported from a txt file.
out = stats.cauchy.fit(data)
tmp = np.linspace(-3, 3, 1000)
y = stats.cauchy.pdf(tmp, *out)
plt.figure()
ax = plt.subplot(2,1,1);
(_, o2) = stats.probplot(data, dist='cauchy', fit=True, plot=ax);
plt.subplot(2,1,2)
plt.hist(data, 100, density=True)
plt.plot(tmp, y)
print(f"slope: {o2[0]}")
print(f"intercept: {o2[1]}")
print(f"r2: {o2[2]}")


It looks to me as if you have the quantiles and the order statistics interchanged, so that the quantiles are on the vertical axis.
Large i.i.d. samples from the Cauchy distribution typically have a substantial number of wildly extreme outliers. That is what you are seeing here.
