random.gamma incorrect results #13138

In addition to parameterization (see: #3928 ), numpy.random.gamma has an incorrect result, not aligned with scipy.stats.gamma . Please check example below:
Numpy version: 1.15.2
Scipy version: 1.1.0

Code example:

import numpy as np
import scipy.stats as st
from matplotlib import pyplot as plt
alpha = 7.
lamdba = 0.0025
theta = 1/lamdba
# Option 1:
mean = alpha / lamdba
var = alpha / (lamdba ** 2)
std = np.sqrt(var)
print(mean, std)
# Option 2 (scipy.stats):
print(st.gamma.mean(a=alpha, scale=1/lamdba), st.gamma.std(a=alpha, scale=1/lamdba))
# Option 3 & 4:
distribution = st.gamma(a=alpha, scale=1/lamdba)
print(distribution.mean(), distribution.std())
m, v = distribution.stats('mv')
print(m, np.sqrt(v))
print




    
("percentile:")
print("-"*10)
for i in [0.5, 0.75, 0.95, 0.995]:
	print(round(st.gamma.ppf(i, a=alpha, scale=1/lamdba),2))
print("10.000 simulations:")
print("-"*10)
items = [1, lamdba, std, var, 1/lamdba, 1/std, 1/var, theta**0.5, theta**0.5/alpha]
for item in items:
	scale = item
	obs = []
	for i in range(10000):
		s = np.random.gamma(shape=mean, scale=scale)
		obs.append(s)
	print(np.mean(obs), np.percentile(obs,99.5)) # ~2800, ~6263.87
# Plot:
x = np.linspace (0, 9000, 200) 
y = st.gamma.pdf(x, a=alpha, scale=1/lamdba)
plt.plot(x, y)
plt.show()

2800.0 1058.3005244258363
2800.0 1058.3005244258363
2800.0 1058.3005244258363
2800.0 1058.3005244258363
percentile:
----------
2667.85
3423.39
4736.96
6263.87
10.000 simulations:
----------
2799.897006963713 2939.4130691768314
6.999216070936152 7.354750837444137
2962502.3180644396 3107142.964045034
3135519955.3881283 3286186790.337828
1120179.56533464 1174490.7069547528
2.6462321958833925 2.7776389287956964
0.00249983363715886 0.002619975990984695
55999.8440596136 58808.98014137368
7999.3703145983245 8386.056327629416