If you just want a distribution that looks approximately normal and satisfies your descriptive stats, here is one possible approach. Start with a normally distributed sample of 148 numbers and appliesapply a series of transformations to (approximately) satisfy the descriptive stats. Of course, there are many distributions that could satisfy the problem...
# function for descriptive stats
stats = function(x) c(min(x),max(x),median(x),mean(x),sd(x))
# simple power transformation (hold min and max constant)
pow = function(x,lam) {
t = (x-min(x))^lam
(t/max(t))*(max(x)-min(x))+min(x)
}
# power transform of upper and lower halves of data (hold min,max,median constant)
pow2 = function(par, x) {
m = median(x)
t1 = pow(m-x[1:74], par[1])
t2 = pow(x[75:148]-m, par[2])
c(m-t1, t2+m)
}
# transformation to fit minimum and maximum
t1 = function(x) {
x = ((x-min(x))/diff(range(x)) *110) + 50
}
# optimise power transformation match median
t2 = function(x) {
l = optimise(function(l) { (median(pow(x,l))-97.7)^2 }, c(-5,5))$min
pow(x,l)
}
# optimise power transformation of upper and lower halves to fit mean and sd
t3 = function(x) {
l2 = optim(c(1,1), function(par) {
r = pow2(par,x); (mean(r)-101.73)^2 + (sd(r)-20.45)^2 })$par
pow2(l2, x)
}
d = t1(sort(rnorm(148)))
stats(d)
d = t2(d)
stats(d)
d = t3(d)
stats(d) # result should match your descriptive stats
hist(d) # looks normal-ish
# repeat and plot many distributions that satisfy requirements
plot(d,cumsum(d), type="l")
for(n in 1:500) {
d = t3(t2(t1(sort(rnorm(148)))))
lines(d,cumsum(d), col=rgb(1,0,0,0.05))
}