Coding with pm.DensityDist() model, random method and sample_posterior_predictive problems

Sorry to bother you, I am quite a newbie here.
I have sampled this data from the mielke distribution:

“”“array([ 9.8, 11.1, 6.8, 13.6, 11.4, 10.5, 9.2, 10. , 9.9, 10.5, 10.8,
8.7, 8.7, 12.2, 10.1, 4.2, 10.3, 11.7, 6.8, 11.7, 11.4, 11.9,
11.3, 8.5, 8.3, 12. , 11.4, 10.2, 12.5])”“”

“”“scipy.stats.mielke.rvs(k=9, s=40, loc= -7, scale=19)”“”

I just would like to use pm.DensitDist() and try to reverse engineering the posterior, in order to understand how to sample posterior predictive for density distribution not in the pool already coded.

I defined logp as:

“”" def logp_mielke(x, k, s, loc, scale):
return k*x**(k-1)/(1+xs)(1+k/s) #Mielke “”"

whic is the definition of mielke distribution

I am having very hard times in understand how to create a random method for this distribution, and I do not find any example that I could follow.

Can you give me an hint?
P.S. I am having hard time in trying to install PyMC4.

Here so far I got:

“”"
data1 = [ 9.8, 11.1, 6.8, 13.6, 11.4, 10.5, 9.2, 10. , 9.9, 10.5, 10.8,
8.7, 8.7, 12.2, 10.1, 4.2, 10.3, 11.7, 6.8, 11.7, 11.4, 11.9,
11.3, 8.5, 8.3, 12. , 11.4, 10.2, 12.5]

def logp_mielke(x, k, s, loc, scale):
return k*x**(k-1)/(1+xs)(1+k/s) #Mielke

with pm.Model() as M1:
loc = pm.Normal(‘loc’, mu=-7, sd=5)
scale = pm.Normal(‘scale’, mu= 19, sd=5)
k = pm.Normal(‘k’, mu=9, sd=5)
s = pm.Normal(‘s’, mu=40, sd=5)

y = pm.DensityDist('y',logp=logp_mielke, observed={'loc':loc, 'scale':scale, 'k':k, 's':s, 'x': data1})

with M1:
trace = pm.sample(draws=1000, tune= 1000, return_inferencedata=True, idata_kwargs={“density_dist_obs”: False}, chains=4)

with M1:
az.plot_trace(trace)
summary = az.summary(trace)
summary

“”"
So far it seems to work:

But without a random method I am not able to sample the
posterior predictive. I sifted the whole on line data without success.

Thank you for your patience.