Specifying a non-inclusive distribution in Stan: Kumarasuamy Regression

Modeling
,
1

Hi,

I’m new here and new to Stan, I would like to know how I can build a kumarasuamy regression model for my data. How to specify Y ~ kumarasuamy(a,b) in Stan?

modeling proposal:

data {
int<lower=1> N; // sample size
int<lower=1> K; // K predictors
vector<lower=0,upper=1>[N] y; // response
matrix[N,K] X; // predictor matrix
}

parameters {
vector[K] theta; // reg coefficients
real<lower=0> phi; // dispersion parameter
}

transformed parameters{
vector[K] beta;

beta = theta * 5; // same as beta ~ normal(0, 5); fairly diffuse
}

model {
// model calculations
vector[N] LP; // linear predictor
vector[N] mu; // transformed linear predictor

LP = X * beta;

for (i in 1:N) {
mu[i] = inv_logit(LP[i]);
}

// priors
theta ~ normal(0, 1);
phi ~ cauchy(0, 5); // different options for phi
//phi ~ inv_gamma(.001, .001);
//phi ~ uniform(0, 500); // put upper on phi if using this

// likelihood
y ~ kumarasuamy (a, b);
}

from this model I want to build the kumarasuamy regression model.

thanks in advance.

1 Like
2

Welcome to the Stan forums.

You can specify custom distributions in the functions block. See the user’s guide. Here’s what it might look like in your case (though, of course, double-check my work to make sure I’m describing the distribution you have in mind).

functions{
    real kumarasuamy_lpdf(real x, real a, real b){
        // x ~ abx^(a-1) * (1 - x^a)^b-1
        return log(a) + log(b) + (a - 1.0) * log(x) + (b - 1.0) * log(1.0 - x^a);
    }
}
data{
    int N; // Number of observations
    vector<lower = 0, upper = 1>[N] Y; // Outcome
}
parameters{
    real<lower = 0> alpha_k;
    real<lower = 0> beta_k;
}
model{
    alpha_k ~ gamma(0.001, 0.001);
    beta_k ~ gamma(0.001, 0.001);

    for(n in 1:N){
        Y[n] ~ kumarasuamy(alpha_k, beta_k);
    }
}
3 Likes
3

Thanks for the answer. So from this function I can build my kumarasuamy regression model?

4

Yes. You’ll just need to specify how mu and phi in your code relate to a and b in kumarasuamy_lpdf. Assuming a = \mu \phi and b = (1 - \mu) \phi (as in Stan’s built-in beta_proportion distribution, see here), then it might look like this.

model{
   ...
   for(n in 1:N){
      y[n] ~ kumarasuamy(mu[n] * phi, (1 - mu[n]) * phi);
   }
}

Sidenote: you can surround your Stan code with three backticks and “stan” to style it for easier reading on Discourse. e.g.

```stan
   y ~ normal(0, 1)
```

becomes

y ~ normal(0, 1)
3 Likes