# Feller-Pareto distribution in stan

Hello everyone! I am new to stan and trying to model Feller-Pareto distribution X ~ FP(sigma, delta_1, delta_2, gamma), which has the following representation:

X = sigma * (U_1 / U_2)^gamma
U_1 ~ Gamma(delta_1, 1)
U_2 ~ Gamma(delta_2, 1) independently.

My goal is to estimate the parameters (sigma, delta_1, delta_2, gamma) based on independently observed samples X ∼ FP(sigma, delta_1, delta_2, gamma).

My question is where to define the intermediate variable U_1 and U_2. I tried to write a model block like this:

``````model {
real U1[N];
real U2[N];
U1 ~ gamma(delta1, 1);
U2 ~ gamma(delta2, 1);
for(i in 1:N) {
X[i] = sigma * pow((U1[i]/U2[i]), gamma);
}
}
``````

But an error was thrown saying “attempt to assign variable in wrong block. left-hand-side variable origin=data”. I must miss something obvious. Please give me some hints.

Just realized that I need to use the fact that the distribution of X conditioned on U2 is a generalized gamma.

You need to learn the structure of the Stan language before attempting to use it. Both `U1` and `U2` have to be declared in the parameters block, while `X` should be declared and defined in the generated quantities block (and not in the data block).

Thanks to your kindly reply. I tried your suggestions as

``````data {
int<lower=1> N; //the number of observations
}

parameters {
real<lower=0> delta1;
real<lower=0> delta2;
real<lower=0> gamma;
real<lower=0> sigma;
real U1[N];
real U2[N];

}

model {
U1 ~ gamma(delta1, 1);
U2 ~ gamma(delta2, 1);
}

generated quantities {
real<lower=0> X[N];
for(i in 1:N) {
X[i] = sigma * pow((U1[i]/U2[i]), gamma);
}
}

``````

But there is an initialization error like "Initialization between (-2, 2) failed after 100 attempts. ". Any thoughts?

With a lack of a prior distribution for `delta1`, `delta2`, `gamma`, and `sigma`, this posterior distribution may be improper or at least is very difficult to sample from. If your goal is just to sample from the Feller-Pareto distribution, those four quantities would typically be passed as data rather than declared as parameters.

If `U1` has a gamma distribution, you need to declare:

``````real<lower = 0> U1[N];
``````

If the variable is unconstrained, initialization is `uniform(-2, 2)`, so there’s a 50-50 chance each initialization will cause the gamma sampling statement to fail.