For calculating the pmf for displaying expected frequencies in the histogram, you assume that the number of Bernoulli trials AND the parameters a and b are known quantities. Of course, when presented with data believed to come from a beta-binomial distribution, the parameters a and b will typically need to be estimated. Even though specified here as a parameter of the beta-binomial distribution, the number of Bernoulli trials will generally be known. When the number of trials are known, the parameters a and b can be estimated. The Wikipedia entry on the beta-binomial distribution mentions both moment estimators of a and b and also states that one can employ likelihood maximization to estimate the unknown parameters. Note that estimation through likelihood maximization allows the parameters themselves to depend on predictor variables . We do not have predictor variables in the example provided. But should predictor variables exist, the ability to allow the parameters to depend on covariates is beneficial.

Below, I present NLMIXED code that produces ML estimates of the parameters a and b. The NLMIXED procedure also computes the expected frequency under the beta-binomial distribution. There is no need for a data step to produce the expected frequencies. The parameter estimates are read in a data step and written to macro variables that are subsequently referenced as &a. and &b. in the SGPLOT procedure.

Note that the default NLMIXED likelihood maximization procedure (QUANEW) produced estimates of the parameters but did not produce standard errors for the parameter estimates. For these data, the DBLDOG maximization technique produced estimates with confidence intervals that contain the true parameters. Other data may require a different likelihood maximization technique.

Note, too, that the parameters a and b must be positive. Consequently, the NLMIXED procedure is parameterized in terms of log_a and log_b. Those parameters are exponentiated to produce a and b which are passed into the log-likelihood computation. Should a and b depend on covariates, the code expressed as "a = exp(log_a)" could be rewritten as "a = exp(b0_a + b1_a*X1 + b2_a*X2 ...;" with similar code for parameter b.

ods output AdditionalEstimates=betaBin_ests;
proc nlmixed data=BetaBin tech=dbldog;
parms log_a log_b 1; /* Parameterize in terms of log_a and log_b */
a = exp(log_a); /* a > 0 ensured */
b = exp(log_b); /* b > 0 ensured */
nTrials = 10;
ll = lcomb(nTrials, x) +
logbeta(x+a, nTrials - x + b) -
logbeta(a, b);
pmf = exp(ll);
model ll ~ general(ll);
estimate "a" exp(log_a);
estimate "b" exp(log_b);
predict &SampleSize. * pmf
out=PDFBetaBinom_est(keep=x pred); /* expected frequency computation */
run;

/* Keep just a single expected value record for each unique value of X. All other */
/* records are duplicates. Also, the predict statement names the predicted value */
/* as "pred". Rename to EX_est for consistency with Rick Wicklin presentation. */
proc sort data=PDFBetaBinom_est(rename=(pred=EX_est)) nodupkey;
by x;
run;

/* Write parameter estimates a and b to macro variables */
data _null_;
set betaBin_ests;
call symput(label, put(estimate,6.4));
run;

/* Merge simulated data and PMF. Overlay PMF on data distribution. */
data All_est;
merge BetaBin PDFBetaBinom_est(rename=(x=t));
run;

title "Estimated Beta-Binomial Distribution";
title2 "Sample Size = &SampleSize";
proc sgplot data=All_est;
vbarbasic x / barwidth=1 legendlabel='Simulated Sample'; /* requires SAS 9.4M3 */
scatter x=t y=EX_est / legendlabel='Expected Value'
markerattrs=GraphDataDefault(symbol=CIRCLEFILLED size=10);
inset "nTrials = 10" "a = &a." "b = &b." / position=topleft border;
yaxis grid;
xaxis label="x" integer type=linear; /* force TYPE=LINEAR */
run;