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.