Tweedie distribution in Stan

model.stan

data {
  int N;
  int M;
  real<lower=0> Y[N];
}

parameters {
  real<lower=0> mu;
  real<lower=0> phi;
  real<lower=1, upper=2> theta;
}

transformed parameters {
  real lambda = 1/phi*mu^(2-theta)/(2-theta);
  real alpha = (2-theta)/(theta-1);
  real beta = 1/phi*mu^(1-theta)/(theta-1);
}

model {
  mu ~ cauchy(0, 5);
  phi ~ cauchy(0, 5);

  for (n in 1:N) {
    if (Y[n] == 0) {
      target += -lambda;
    } else {
      vector[M] ps;
      for (m in 1:M)
        ps[m] = poisson_lpmf(m | lambda) + gamma_lpdf(Y[n] | m*alpha, beta);
      target += log_sum_exp(ps);
    }
  }
}

run-model.R

library(rstan)
library(tweedie)

stanmodel <- stan_model(file='model/model.stan')

N1 <- 200
M1 <- 30
Y1 <- rtweedie(N1, power=1.3, mu=1, phi=1)
data1 <- list(N=N1, M=M1, Y=Y1)
fit1 <- sampling(stanmodel, data=data1)

N2 <- 1000
M2 <- 30
Y2 <- rtweedie(N2, power=1.01, mu=3, phi=1)
data2 <- list(N=N2, M=M2, Y=Y2)
fit2 <- sampling(stanmodel, data=data2)

I got it, many thanks.