You cannot set boundaries for a gamma distribution but you can clamp resulting values to [minimum, maximum]
import numpy as np
from scipy import stats
data = stats.gamma.rvs(2, loc=1.5, scale=2, size=10000)
print(data) # [1.55651126 7.25914918 5.03753265 ... 5.56524255 3.05888237 2.17275238]
# clamp to [minimum, maximum]
minimum = 3
maximum = 5
# assert minimum < maximum
data = np.where(data < minimum, minimum, data)
data = np.where(data > maximum, maximum, data)
print(data) # [3. 5. 5. ... 5. 3.05888237 3. ]
Adjusting for your comment: if you want to sample n
elements in [minimum, maximum] without clamping you could do something like:
from scipy import stats
minimum = 3
maximum = 5
sampled, n = [], 10_000
while len(sampled) != n:
data = stats.gamma.rvs(2, loc=1.5, scale=2, size=n - len(sampled))
sampled.extend(data[(minimum <= data) & (data <= maximum)])