Hey all,

I have a fairly simple model I’d like to implement in Stan. It’s a mixture of doubly truncated Pareto distributions. K mixture components have their own ymin and alpha, but they share an upper truncation parameter ymax.

```
data {
int<lower=1> K; // number of mixture components
int<lower=0> N; // number of planets
real y[N]; // y values for N planets
}
parameters {
simplex[K] theta; // mixing proportions
real alpha[K]; // power-law indices
real ymin[K]; // ymin for each component
real ymax; // ymax for all K components
}
model {
real ps[K];
real ylower[K];
real yupper[K];
real integral[K];
ylower[1] = ymin[1];
ylower[2] = ymin[2];
yupper[1] = ymax;
yupper[2] = ymax;
alpha ~ uniform(-10.,10.); // prior on alpha
ymin ~ uniform(0., 30.); // prior on ymin
ymax ~ uniform(ymin[2], 30.); // prior on ymax
integral[1] = (yupper[1]^(alpha[1]) - ylower[1]^(alpha[1]))/alpha[1];
integral[2] = (yupper[2]^(alpha[2]) - ylower[2]^(alpha[2]))/alpha[2];
for (n in 1:N) { // loop over planets
for (k in 1:K){ // loop over mixture components
if ((y[n] > ylower[k]) && (y[n] < yupper[k]))
{
ps[k] = theta[k] * (y[n])^(alpha[k]-1.)/integral[k];
}
else
{
ps[k] = 0.;
}
}
target += log(sum(ps));
}
}
```

I realize Stan supports Pareto distributions (`y ~ pareto(ymin, alpha)`

), but I’d like to allow for monotonically increasing distributions, i.e., a Pareto with a positive exponent, alpha<-1.

The above code works, but I doubt it’s an efficient implementation. For N=1000, it takes roughly ~0.1s per model evaluation. Any advice for speeding up this Stan model?