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