Log Likelihood for BETA distribution in Proc NLMIXED

[LadyIDO]

Hello,

I'm trying to estimate the beta distribution parameters (alpha, beta) using Proc NLMIXED. Because beta is not one of the distributions with built-in log likelihood, I used the general statement as follows:

proc nlmixed data=test;
  parms a=5 b=5;
  *s=beta(a,b);
  s=gamma(a*b)/(gamma(a)*gamma(b));
  ll=(1-a)*xbeta+(1-b)*(1-xbeta)-log(s);
  model xbeta~general(ll);
run;

I tried using the beta function directly and with gamma function. Both gave me the 'QUANEW optimization cannot be completed' error.

The test data is just 1000 beta variates with alpha=beta=5. I even set the initial values to 5, but still got this error. 

Anne

[StatDave]

If you want to fit a model to a response variable that is distributed as beta, then you can use the FMM or GLIMMIX procedure. Specify DIST=BETA in the MODEL statement. If your goal is to estimate the parameters of the beta distribution for one population, you can use the HISTOGRAM statement in PROC UNIVARIATE as mentioned in this note.

[LadyIDO]

I was trying to understand how to use nlmixed to estimate parameters of 
less common distribution like the following for beta-binomial. I didn't know about
FMM before. I'll check it out. Thanks.

proc nlmixed data=&data fconv=1E-14 df=&rdf alpha=0.05 ecov cov;
  parms mu=0.5 gamma=0.5;
  bounds mu>=0, mu<=1, gamma>0, gamma<1;
  ll=0;
  theta=gamma/(1-gamma);
  do i=1 to &ntrials;
    if i<=&nsucc
      then ll=ll+log(mu+(i-1)*theta);
        else ll=ll+log((1-mu)+(i-&nsucc-1)*theta);
    ll=ll-log(1+(i-1)*theta);
  end;
  alpha=mu/theta;
  beta=(1-mu)/theta;
  model &nsucc~general(ll);
run;

[Rick_SAS · Accepted solution]

First, the beta function is 

s=gamma(a+b)/(gamma(a)*gamma(b));

not

s=gamma(a*b)/(gamma(a)*gamma(b));

Second, the beta distribution when a=b is somewhat degenerate, in the sense that the mean is exactly 1/2. This might be relevant for your nonconvergence issues. I'm not sure.

But if you want to use PROC NLMIXED, the easiest way is to specify the LOGPDF of the beta distribution directly:

data test;
call streaminit(54321);
do i = 1 to 1000;
   xbeta = rand("Beta", 5, 5);
   output;
end;
run;

proc nlmixed data=test;
  parms a=4 b=5.5;
  ll= logpdf("Beta", xbeta, a, b);;
  model xbeta~general(ll);
run;

[LadyIDO]

Thank you so much Rick! Using logpdf did it! I didn't know about logpdf function. This would really help with future work needing log likelihood.

Thanks again!

Anne