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
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;
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.
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;
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;
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
Are you ready for the spotlight? We're accepting content ideas for SAS Innovate 2025 to be held May 6-9 in Orlando, FL. The call is open until September 25. Read more here about why you should contribute and what is in it for you!
ANOVA, or Analysis Of Variance, is used to compare the averages or means of two or more populations to better understand how they differ. Watch this tutorial for more.
Find more tutorials on the SAS Users YouTube channel.