Beta distribution with unknown boundaries not converging

Modeling
1

Dear all,
I want to fit data to a Beta distribution (4 parameters (alpha, beta, a, c) as def. in Sec. 1.2.3. Beta distribution - Wikipedia). In the code below I also treat the parameters a and c (the boundaries) as uncertain parameters.

The code below samples Data with given parameters. The parameters that are returned by stan are far from what is expected.

I am looking forward to your comments and hints.


############################################################
#   Python code
    n_sim = 1000
    alpha = 5
    beta = 7
    y = stats.beta.rvs(alpha, beta, size=n_sim)
    scale_ = 15
    loc_ = 80
    x = y * scale_ + loc_
    cdf = np.linspace(0, 1, len(x))
    plt.figure(1)
    plt.plot(np.sort(x), cdf)

    ############################################################
    #   build stan model
    my_data = {}  # create data dictionary
    my_data["data_n"] = n_sim
    my_data["X"] = x

    print("compile stan model")
    fname = "scr_84_model.stan"
    sm = pystan.StanModel(file=fname)  # build stan model

    ############################################################
    #   fit with stan model
    print("fit with stan")
    fit = sm.sampling( data = my_data, iter = 500
                       )

############################################################
    #   STAN model
//  MODEL 84
//  4 parameter beta distribution

functions{
  real betaself_lpdf( vector y, real alpha_, real beta, real a, real c) {
    vector[rows(y)] lp;
    int N;
    N = rows(y);

    for(n in 1:N){
      if ((y[n]-a)<=0)
        reject("y[n]-a<=0; y[n], a =", y[n], a)
      if ((c-y[n])<=0)
        reject("c-y[n]<=0; c, y[n] =", c, y[n])
      if ((c-a)<=0)
        reject("c-a<=0; c, a =", c, a)

      lp[n]=(alpha_-1)*log(y[n]-a)+(beta-1)*log(c-y[n])-(alpha_+beta-1)*log(c-a)-lbeta(alpha_,beta) ;
   }
    return -sum( lp );
  }

}


data {
    int data_n;             // number of values
    vector[data_n] X;       // data values
}

parameters {
    real <lower= 1e-15, upper = 15 > alpha_;
    real <lower= 1e-15, upper = 15 > beta_;
    real <lower= 1e-15, upper = 100> a;
    real <lower= 1e-15, upper = 150> c;
}

model {
    // priors
    alpha_  ~  normal(5, 3);
    beta_   ~  normal(5, 3);
    a   ~  normal(80, 8);
    c   ~  normal(95, 8);

    //  likelihood
    //    X ~ betaself(alpha_, beta_, a, c);

    target += betaself_lpdf( X | alpha_, beta_, a, c);
}
2

I replaced: return -sum( lp );

with: return sum( lp );

and obtained meaningful values.

Does STAN maximize the loglikelihood and not minimize the neg. loglikelihood? I would appreciate a short comment.

3

Hi,
a Stan program is a function to compute the (positive) log density given the data. I.e. it is not only the log likelihood, but also log prior. Stan can maximize this log density (if you use the optimizing mode), but by default it uses it to run Hamiltonian MC, which works with the positive log density. (I am not sure about the historical roots of the convention to minimize negative log likelihood, but I presume it has to do with what was implemented in early software packages, not any fundamental theoretical reason)

Additionally, you are quite likely to get faster sampling by constraining a to be smaller than the smallest observed value and c to be larger than the largest observed value (as the model already asumes that).

Best of luck with your model!

2 Likes