I have the following code that works for continuous distributions; I need to modify it to work for discrete distributions (e.g., target becomes Poisson):

ClearAll["Global`*"];

rho = 0.95

r = Sqrt[1 - rho^2];

\[ScriptCapitalD] = (Gamma[
      2*a]/(Gamma[a])^2)*(r*(x*(1 - x))^(r*a))/(x*(1 - 
        x)*(x^r + (1 - x)^r)^(2*a));

a = 7.56158;

\[ScriptCapitalD]1 = 
  ProbabilityDistribution[\[ScriptCapitalD], {x, 0, 1}];

target = UniformDistribution[];

\[ScriptCapitalD]2 = 
  TransformedDistribution[Refine[InverseCDF[target, p], 0 <= p <= 1], 
   p \[Distributed] \[ScriptCapitalD]1];

data = SetPrecision[RandomVariate[\[ScriptCapitalD]2, 100000], 16];

Mean[data]

StandardDeviation[data]

Skewness[data]

Kurtosis[data] - 3.0