How can I specify similar to the code below so that now considering the kumarasuamy distribution ? I believe it is necessary to implement the kumarasuamy_rng function, how can I proceed?

Or rather, how to generalize to any distribution that does not have an accumulated distribution function with a closed expression?

```
generated quantities {
vector[N] y_rep;
for (i in 1:N) {
real mu;
real A;
real B;
mu = inv_logit(X[i] * beta);
A = mu * phi;
B = (1.0 - mu) * phi;
y_rep[i] = beta_rng(A, B);
}
}
```