This page presents SAS code for implementing Bayesian conjugate models using standard prior–likelihood formulations and their associated posterior distributions.
- Binomial Beta model
- Poisson Gamma model
- Normal Normal model
- Normal Gamma model
- Normal Inverse gamma model
For preliminary theoretical part, refer
The choice of parameters should always follow the native argument convention of the corresponding R function. Several probability distributions, such as the Gamma distribution, permit alternative but equivalent parameterizations (e.g., shape–scale or shape–rate). Consequently, we should exercise appropriate care while specifying parameters to ensure consistency with the underlying R implementation.
The implementations use closed-form analytical expressions for posterior summaries wherever available, thereby avoiding numerical approximation except where necessary (e.g., posterior quantiles).
Numerical output includes posterior summaries such as the mean, median, and mode (where defined and computable), with user-selectable posterior quantiles and credible intervals of any desired probability level. Appropriate safeguards have been incorporated for cases in which the posterior mode does not exist or cannot be determined analytically.
<div class="copy-content">
<p id="copyText">* ============================================================================= Bayesian Conjugate Models – Closed Form Analytical Solutions SAS Macro Port – exact parallel to R and Python versions Models : Binomial-Beta | Poisson-Gamma | Normal-Normal Normal-Gamma | Normal-Inverse-Gamma Input : SAS dataset + var + group_var OR scalar / inline values Output : ODS HTML tables (proc print) + sgplot – group-wise complete output All moments : closed form Quantiles : SAS quantile() function – native SAS parameterisation quantiles= accepts integer percentages e.g. %str(10 25 75) ============================================================================= */ /* ============================================================================= NOTES ON NATIVE SAS PARAMETERISATION —————————————————————————– Beta : pdf(‘BETA’, x, alpha, beta) quantile(‘BETA’, p, alpha, beta) Gamma : pdf(‘GAMMA’, x, alpha, beta) alpha=shape, beta=scale quantile(‘GAMMA’, p, alpha, beta) Mean = alpha*beta, Var = alpha*beta^2, Mode = (alpha-1)*beta Normal : pdf(‘NORMAL’, x, mu, sigma) sigma = SD (not variance) quantile(‘NORMAL’, p, mu, sigma) InvGamma : No direct SAS distribution – derived via: If X ~ Gamma(alpha, beta_scale) then 1/X ~ InvGamma(alpha, 1/beta_scale) quantile InvGamma(p, a, b_scale) = 1 / quantile(‘GAMMA’, 1-p, a, 1/b_scale) pdf InvGamma(x, a, b_scale) = pdf(‘GAMMA’, 1/x, a, 1/b_scale) / x^2 Mean = b_scale/(a-1) for a>1, Mode = b_scale/(a+1) Users pass parameters native to each model as documented below. ============================================================================= */ /* =============================================================================
</p>
</div>
SAS Code is
Bionimal thoery
Show Code
/* =============================================================================
Bayesian Conjugate Models - Closed Form Analytical Solutions
SAS Macro Port - exact parallel to R and Python versions
Models : Binomial-Beta | Poisson-Gamma | Normal-Normal
Normal-Gamma | Normal-Inverse-Gamma
Input : SAS dataset + var + group_var OR scalar / inline values
Output : ODS HTML tables (proc print) + sgplot - group-wise complete output
All moments : closed form
Quantiles : SAS quantile() function - native SAS parameterisation
quantiles= accepts integer percentages e.g. %str(10 25 75)
============================================================================= */
/* =============================================================================
NOTES ON NATIVE SAS PARAMETERISATION
-----------------------------------------------------------------------------
Beta : pdf('BETA', x, alpha, beta)
quantile('BETA', p, alpha, beta)
Gamma : pdf('GAMMA', x, alpha, beta) alpha=shape, beta=scale
quantile('GAMMA', p, alpha, beta)
Mean = alpha*beta, Var = alpha*beta^2, Mode = (alpha-1)*beta
Normal : pdf('NORMAL', x, mu, sigma) sigma = SD (not variance)
quantile('NORMAL', p, mu, sigma)
InvGamma : No direct SAS distribution - derived via:
If X ~ Gamma(alpha, beta_scale) then 1/X ~ InvGamma(alpha, 1/beta_scale)
quantile InvGamma(p, a, b_scale) = 1 / quantile('GAMMA', 1-p, a, 1/b_scale)
pdf InvGamma(x, a, b_scale) = pdf('GAMMA', 1/x, a, 1/b_scale) / x^2
Mean = b_scale/(a-1) for a>1, Mode = b_scale/(a+1)
Users pass parameters native to each model as documented below.
============================================================================= */
/* =============================================================================
1. BINOMIAL - BETA
=============================================================================
Likelihood : X | theta ~ Binomial(n, theta)
Prior : theta ~ Beta(alpha, beta)
Posterior : theta | x ~ Beta(alpha + x, beta + n - x)
Parameters:
data = input SAS dataset name (or leave blank for inline x/n_trials)
var = binary variable column name (character or 0/1 numeric)
group_var = grouping variable (optional, blank for none)
success_level = level that counts as success (default: first sorted level)
x = number of successes (scalar, used when data= is blank)
n_trials = number of trials (scalar, used when data= is blank)
alpha = prior Beta shape1
beta = prior Beta shape2
cri_levels = space-separated CrI levels e.g. %str(90 95 99)
quantiles = space-separated integer percentiles e.g. %str(10 25 75)
cri_type = label for CrI column (default: CrI)
============================================================================= */
%macro conjugate_binomial_beta(
data = ,
var = ,
group_var = ,
success_level = ,
x = ,
n_trials = ,
alpha = ,
beta = ,
cri_levels = 90 95,
quantiles = ,
cri_type = CrI
);
%local i g lev lo hi q nlevs nq ngrps grp_val label_val;
/* --- Input path: dataset or scalar --- */
%if %length(&data) > 0 %then %do;
%let dsid = %sysfunc(open(&data));
%if &dsid = 0 %then %do;
%put ERROR: Dataset &data not found.;
%return;
%end;
%if %sysfunc(varnum(&dsid, &var)) = 0 %then %do;
%put ERROR: Variable &var not found in dataset &data.;
%let rc = %sysfunc(close(&dsid));
%return;
%end;
%let rc = %sysfunc(close(&dsid));
data _bb_raw_;
set &data;
where &var ne '';
%if %length(&success_level) > 0 %then %do;
_success_ = (&var = "&success_level");
%end;
%else %do;
_success_ = .;
%end;
run;
%if %length(&success_level) = 0 %then %do;
proc sql noprint;
select min(&var) into :success_level trimmed
from _bb_raw_;
quit;
data _bb_raw_;
set _bb_raw_;
_success_ = (&var = "&success_level");
run;
%end;
%if %length(&group_var) > 0 %then %do;
proc sql noprint;
create table _bb_summ_ as
select catx(' = ', "&group_var", &group_var) as _label_ length=200,
sum(_success_) as _x_,
count(*) as _nobs_
from _bb_raw_
group by &group_var
order by &group_var;
quit;
%end;
%else %do;
proc sql noprint;
create table _bb_summ_ as
select 'All' as _label_ length=200,
sum(_success_) as _x_,
count(*) as _nobs_
from _bb_raw_;
quit;
%end;
%end;
%else %do;
data _bb_summ_;
_label_ = "1";
_x_ = &x;
_nobs_ = &n_trials;
run;
%end;
/* --- Compute posterior parameters and moments for all groups --- */
data _bb_results_;
set _bb_summ_;
a_prior = α
b_prior = β
a_post = &alpha + _x_;
b_post = &beta + (_nobs_ - _x_);
pr_mean = &alpha / (&alpha + &beta);
pr_var = (&alpha * &beta) / ((&alpha + &beta)**2 * (&alpha + &beta + 1));
pr_sd = sqrt(pr_var);
pr_median = quantile('BETA', 0.5, &alpha, &beta);
if &alpha > 1 and &beta > 1 then
pr_mode = (&alpha - 1) / (&alpha + &beta - 2);
else pr_mode = .;
po_mean = a_post / (a_post + b_post);
po_var = (a_post * b_post) / ((a_post + b_post)**2 * (a_post + b_post + 1));
po_sd = sqrt(po_var);
po_median = quantile('BETA', 0.5, a_post, b_post);
if a_post > 1 and b_post > 1 then
po_mode = (a_post - 1) / (a_post + b_post - 2);
else po_mode = .;
/* =============================================================================
Bayesian Conjugate Models - Closed Form Analytical Solutions
SAS Macro Port - exact parallel to R and Python versions
Models : Binomial-Beta | Poisson-Gamma | Normal-Normal
Normal-Gamma | Normal-Inverse-Gamma
Input : SAS dataset + var + group_var OR scalar / inline values
Output : ODS HTML tables (proc print) + sgplot - group-wise complete output
All moments : closed form
Quantiles : SAS quantile() function - native SAS parameterisation
quantiles= accepts integer percentages e.g. %str(10 25 75)
============================================================================= */
/* =============================================================================
NOTES ON NATIVE SAS PARAMETERISATION
-----------------------------------------------------------------------------
Beta : pdf('BETA', x, alpha, beta)
quantile('BETA', p, alpha, beta)
Gamma : pdf('GAMMA', x, alpha, beta) alpha=shape, beta=scale
quantile('GAMMA', p, alpha, beta)
Mean = alpha*beta, Var = alpha*beta^2, Mode = (alpha-1)*beta
Normal : pdf('NORMAL', x, mu, sigma) sigma = SD (not variance)
quantile('NORMAL', p, mu, sigma)
InvGamma : No direct SAS distribution - derived via:
If X ~ Gamma(alpha, beta_scale) then 1/X ~ InvGamma(alpha, 1/beta_scale)
quantile InvGamma(p, a, b_scale) = 1 / quantile('GAMMA', 1-p, a, 1/b_scale)
pdf InvGamma(x, a, b_scale) = pdf('GAMMA', 1/x, a, 1/b_scale) / x^2
Mean = b_scale/(a-1) for a>1, Mode = b_scale/(a+1)
Users pass parameters native to each model as documented below.
============================================================================= */
/* =============================================================================
1. BINOMIAL - BETA
=============================================================================
Likelihood : X | theta ~ Binomial(n, theta)
Prior : theta ~ Beta(alpha, beta)
Posterior : theta | x ~ Beta(alpha + x, beta + n - x)
Parameters:
data = input SAS dataset name (or leave blank for inline x/n_trials)
var = binary variable column name (character or 0/1 numeric)
group_var = grouping variable (optional, blank for none)
success_level = level that counts as success (default: first sorted level)
x = number of successes (scalar, used when data= is blank)
n_trials = number of trials (scalar, used when data= is blank)
alpha = prior Beta shape1
beta = prior Beta shape2
cri_levels = space-separated CrI levels e.g. %str(90 95 99)
quantiles = space-separated integer percentiles e.g. %str(10 25 75)
cri_type = label for CrI column (default: CrI)
============================================================================= */
%macro conjugate_binomial_beta(
data = ,
var = ,
group_var = ,
success_level = ,
x = ,
n_trials = ,
alpha = ,
beta = ,
cri_levels = 90 95,
quantiles = ,
cri_type = CrI
);
%local i g lev lo hi q nlevs nq ngrps grp_val label_val;
/* --- Input path: dataset or scalar --- */
%if %length(&data) > 0 %then %do;
%let dsid = %sysfunc(open(&data));
%if &dsid = 0 %then %do;
%put ERROR: Dataset &data not found.;
%return;
%end;
%if %sysfunc(varnum(&dsid, &var)) = 0 %then %do;
%put ERROR: Variable &var not found in dataset &data.;
%let rc = %sysfunc(close(&dsid));
%return;
%end;
%let rc = %sysfunc(close(&dsid));
data _bb_raw_;
set &data;
where &var ne '';
%if %length(&success_level) > 0 %then %do;
_success_ = (&var = "&success_level");
%end;
%else %do;
_success_ = .;
%end;
run;
%if %length(&success_level) = 0 %then %do;
proc sql noprint;
select min(&var) into :success_level trimmed
from _bb_raw_;
quit;
data _bb_raw_;
set _bb_raw_;
_success_ = (&var = "&success_level");
run;
%end;
%if %length(&group_var) > 0 %then %do;
proc sql noprint;
create table _bb_summ_ as
select catx(' = ', "&group_var", &group_var) as _label_ length=200,
sum(_success_) as _x_,
count(*) as _nobs_
from _bb_raw_
group by &group_var
order by &group_var;
quit;
%end;
%else %do;
proc sql noprint;
create table _bb_summ_ as
select 'All' as _label_ length=200,
sum(_success_) as _x_,
count(*) as _nobs_
from _bb_raw_;
quit;
%end;
%end;
%else %do;
data _bb_summ_;
_label_ = "1";
_x_ = &x;
_nobs_ = &n_trials;
run;
%end;
/* --- Compute posterior parameters and moments for all groups --- */
data _bb_results_;
set _bb_summ_;
a_prior = α
b_prior = β
a_post = &alpha + _x_;
b_post = &beta + (_nobs_ - _x_);
pr_mean = &alpha / (&alpha + &beta);
pr_var = (&alpha * &beta) / ((&alpha + &beta)**2 * (&alpha + &beta + 1));
pr_sd = sqrt(pr_var);
pr_median = quantile('BETA', 0.5, &alpha, &beta);
if &alpha > 1 and &beta > 1 then
pr_mode = (&alpha - 1) / (&alpha + &beta - 2);
else pr_mode = .;
po_mean = a_post / (a_post + b_post);
po_var = (a_post * b_post) / ((a_post + b_post)**2 * (a_post + b_post + 1));
po_sd = sqrt(po_var);
po_median = quantile('BETA', 0.5, a_post, b_post);
if a_post > 1 and b_post > 1 then
po_mode = (a_post - 1) / (a_post + b_post - 2);
else po_mode = .;
obs_prop = _x_ / _nobs_;
run;
/* --- Get number of groups for loop --- */
proc sql noprint;
select count(*) into :ngrps trimmed from _bb_results_;
quit;
/* === GROUP LOOP === */
%do g = 1 %to &ngrps;
/* Extract this group's row */
data _bb_grp_;
set _bb_results_(firstobs=&g obs=&g);
run;
proc sql noprint;
select _label_ into :label_val trimmed from _bb_grp_;
quit;
/* --- Data Summary --- */
title "Binomial-Beta | var: &var | &label_val | Data Summary";
proc print data=_bb_grp_ noobs label;
var _x_ _nobs_ obs_prop;
label _x_ = 'Successes'
_nobs_ = 'Trials'
obs_prop= 'Observed Proportion';
format obs_prop 8.6;
run;
title;
/* --- Parameters & Moments --- */
data _bb_mom_;
set _bb_grp_;
length Quantity $40 Prior $15 Posterior $15;
Quantity = 'alpha';
Prior = put(a_prior, 12.6); Posterior = put(a_post, 12.6); output;
Quantity = 'beta';
Prior = put(b_prior, 12.6); Posterior = put(b_post, 12.6); output;
Quantity = 'Mean';
Prior = put(pr_mean, 12.6); Posterior = put(po_mean, 12.6); output;
Quantity = 'Median';
Prior = put(pr_median, 12.6); Posterior = put(po_median, 12.6); output;
Quantity = 'Mode';
if pr_mode = . then Prior = 'undefined'; else Prior = put(pr_mode, 12.6);
if po_mode = . then Posterior = 'undefined'; else Posterior = put(po_mode, 12.6);
output;
Quantity = 'Variance';
Prior = put(pr_var, 12.6); Posterior = put(po_var, 12.6); output;
Quantity = 'SD';
Prior = put(pr_sd, 12.6); Posterior = put(po_sd, 12.6); output;
keep Quantity Prior Posterior;
run;
title "Binomial-Beta | var: &var | &label_val | Parameters and Moments";
proc print data=_bb_mom_ noobs; run;
title;
/* --- Credible Intervals --- */
%if %length(&cri_levels) > 0 %then %do;
data _bb_cri_;
set _bb_grp_;
length Level $20 Lower $15 Upper $15;
%let nlevs = %sysfunc(countw(&cri_levels));
%do i = 1 %to &nlevs;
%let lev = %scan(&cri_levels, &i);
%let lo = %sysevalf((1 - &lev/100) / 2);
%let hi = %sysevalf(1 - (1 - &lev/100) / 2);
Level = "&lev% &cri_type";
Lower = put(quantile('BETA', &lo, a_post, b_post), 12.6);
Upper = put(quantile('BETA', &hi, a_post, b_post), 12.6);
output;
%end;
keep Level Lower Upper;
run;
title "Binomial-Beta | var: &var | &label_val | Credible Intervals (&cri_type)";
proc print data=_bb_cri_ noobs; run;
title;
%end;
/* --- Posterior Quantiles --- */
%if %length(&quantiles) > 0 %then %do;
data _bb_qtile_;
set _bb_grp_;
length Quantile $15 Value $15;
%let nq = %sysfunc(countw(&quantiles));
%do i = 1 %to &nq;
%let q = %sysevalf(%scan(&quantiles, &i, %str( )) / 100);
Quantile = "Q(&q)";
Value = put(quantile('BETA', &q, a_post, b_post), 12.6);
output;
%end;
keep Quantile Value;
run;
title "Binomial-Beta | var: &var | &label_val | Posterior Quantiles";
proc print data=_bb_qtile_ noobs; run;
title;
%end;
/* --- Plots --- */
data _bb_plot_;
set _bb_grp_;
if min(a_prior, b_prior) < 1 then do; pr_lo=0.0001; pr_hi=0.9999; end;
else do; pr_lo=0; pr_hi=1; end;
if min(a_post, b_post) < 1 then do; po_lo=0.0001; po_hi=0.9999; end;
else do; po_lo=0; po_hi=1; end;
do _i_ = 1 to 1000;
theta_pr = pr_lo + (_i_-1)*(pr_hi-pr_lo)/999;
theta_po = po_lo + (_i_-1)*(po_hi-po_lo)/999;
prior_density = pdf('BETA', theta_pr, a_prior, b_prior);
post_density = pdf('BETA', theta_po, a_post, b_post);
output;
end;
keep theta_pr theta_po prior_density post_density;
run;
ods graphics on;
title "Binomial-Beta | var: &var | &label_val | Prior [Beta(&alpha, &beta)]";
proc sgplot data=_bb_plot_;
series x=theta_pr y=prior_density /
lineattrs=(color="#2C7BB6" thickness=1.5);
xaxis label='theta';
yaxis label='Density';
run;
title;
title "Binomial-Beta | var: &var | &label_val | Posterior";
proc sgplot data=_bb_plot_;
series x=theta_po y=post_density /
lineattrs=(color="#2C7BB6" thickness=1.5);
xaxis label='theta';
yaxis label='Density';
run;
title;
ods graphics off;
%end; /* end group loop */
/* Cleanup */
proc datasets library=work nolist;
delete _bb_raw_ _bb_summ_ _bb_results_ _bb_grp_ _bb_mom_
_bb_cri_ _bb_qtile_ _bb_plot_;
quit;
%mend conjugate_binomial_beta;
/* =============================================================================
2. POISSON - GAMMA
=============================================================================
Likelihood : X_i ~ Poisson(lambda)
Prior : lambda ~ Gamma(alpha, beta) alpha=shape, beta=scale (SAS native)
Posterior : lambda | x ~ Gamma(alpha + sum(x), beta_post_scale)
beta_post_scale = 1 / (1/beta + n)
Parameters:
data = input SAS dataset
var = count variable (non-negative integers)
group_var = grouping variable (optional)
x = space-separated counts (used when data= is blank)
alpha = prior Gamma shape
beta = prior Gamma scale (SAS native: mean = alpha*beta)
cri_levels, quantiles, cri_type as above
============================================================================= */
%macro conjugate_poisson_gamma(
data = ,
var = ,
group_var = ,
x = ,
alpha = ,
beta = ,
cri_levels= 90 95,
quantiles = ,
cri_type = CrI
);
%local i g lev lo hi q nlevs nq ngrps label_val;
%if %length(&data) > 0 %then %do;
%let dsid = %sysfunc(open(&data));
%if &dsid = 0 %then %do;
%put ERROR: Dataset &data not found.;
%return;
%end;
%if %sysfunc(varnum(&dsid, &var)) = 0 %then %do;
%put ERROR: Variable &var not found in dataset &data.;
%let rc = %sysfunc(close(&dsid));
%return;
%end;
%let vartype = %sysfunc(vartype(&dsid, %sysfunc(varnum(&dsid,&var))));
%let rc = %sysfunc(close(&dsid));
%if &vartype ne N %then %do;
%put ERROR: Variable &var must be numeric for Poisson-Gamma model.;
%return;
%end;
%if %length(&group_var) > 0 %then %do;
proc sql noprint;
create table _pg_summ_ as
select catx(' = ', "&group_var", &group_var) as _label_ length=200,
sum(&var) as _sumx_,
count(*) as _nobs_,
mean(&var) as _meanx_
from &data
where &var >= 0
group by &group_var
order by &group_var;
quit;
%end;
%else %do;
proc sql noprint;
create table _pg_summ_ as
select 'All' as _label_ length=200,
sum(&var) as _sumx_,
count(*) as _nobs_,
mean(&var) as _meanx_
from &data
where &var >= 0;
quit;
%end;
%end;
%else %do;
data _pg_summ_;
length _label_ $200;
_label_ = 'All';
_sumx_ = 0;
_nobs_ = %sysfunc(countw(&x));
%let _nxvals_ = %sysfunc(countw(&x));
%do i = 1 %to &_nxvals_;
_sumx_ = _sumx_ + %scan(&x, &i);
%end;
_meanx_ = _sumx_ / _nobs_;
run;
%end;
data _pg_results_;
set _pg_summ_;
a_prior = α
b_prior = β
a_post = &alpha + _sumx_;
b_post_scale = 1 / (1/&beta + _nobs_);
pr_mean = &alpha * β
pr_var = &alpha * &beta**2;
pr_sd = sqrt(pr_var);
pr_median = quantile('GAMMA', 0.5, &alpha, &beta);
if &alpha >= 1 then pr_mode = (&alpha - 1) * β
else pr_mode = .;
po_mean = a_post * b_post_scale;
po_var = a_post * b_post_scale**2;
po_sd = sqrt(po_var);
po_median = quantile('GAMMA', 0.5, a_post, b_post_scale);
if a_post >= 1 then po_mode = (a_post - 1) * b_post_scale;
else po_mode = .;
run;
proc sql noprint;
select count(*) into :ngrps trimmed from _pg_results_;
quit;
%do g = 1 %to &ngrps;
data _pg_grp_;
set _pg_results_(firstobs=&g obs=&g);
run;
proc sql noprint;
select _label_ into :label_val trimmed from _pg_grp_;
quit;
title "Poisson-Gamma | var: &var | &label_val | Data Summary";
proc print data=_pg_grp_ noobs label;
var _nobs_ _sumx_ _meanx_;
label _nobs_ = 'n'
_sumx_ = 'sum(x)'
_meanx_ = 'mean(x)';
format _meanx_ 12.6;
run;
title;
data _pg_mom_;
set _pg_grp_;
length Quantity $40 Prior $15 Posterior $15;
Quantity = 'alpha (shape)';
Prior = put(a_prior, 12.6); Posterior = put(a_post, 12.6); output;
Quantity = 'beta (scale)';
Prior = put(b_prior, 12.6); Posterior = put(b_post_scale, 12.6); output;
Quantity = 'Mean';
Prior = put(pr_mean, 12.6); Posterior = put(po_mean, 12.6); output;
Quantity = 'Median';
Prior = put(pr_median, 12.6); Posterior = put(po_median, 12.6); output;
Quantity = 'Mode';
if pr_mode = . then Prior = 'undefined'; else Prior = put(pr_mode, 12.6);
if po_mode = . then Posterior = 'undefined'; else Posterior = put(po_mode, 12.6);
output;
Quantity = 'Variance';
Prior = put(pr_var, 12.6); Posterior = put(po_var, 12.6); output;
Quantity = 'SD';
Prior = put(pr_sd, 12.6); Posterior = put(po_sd, 12.6); output;
keep Quantity Prior Posterior;
run;
title "Poisson-Gamma | var: &var | &label_val | Parameters and Moments";
proc print data=_pg_mom_ noobs; run;
title;
%if %length(&cri_levels) > 0 %then %do;
data _pg_cri_;
set _pg_grp_;
length Level $20 Lower $15 Upper $15;
%let nlevs = %sysfunc(countw(&cri_levels));
%do i = 1 %to &nlevs;
%let lev = %scan(&cri_levels, &i);
%let lo = %sysevalf((1 - &lev/100) / 2);
%let hi = %sysevalf(1 - (1 - &lev/100) / 2);
Level = "&lev% &cri_type";
Lower = put(quantile('GAMMA', &lo, a_post, b_post_scale), 12.6);
Upper = put(quantile('GAMMA', &hi, a_post, b_post_scale), 12.6);
output;
%end;
keep Level Lower Upper;
run;
title "Poisson-Gamma | var: &var | &label_val | Credible Intervals (&cri_type)";
proc print data=_pg_cri_ noobs; run;
title;
%end;
%if %length(&quantiles) > 0 %then %do;
data _pg_qtile_;
set _pg_grp_;
length Quantile $15 Value $15;
%let nq = %sysfunc(countw(&quantiles));
%do i = 1 %to &nq;
%let q = %sysevalf(%scan(&quantiles, &i, %str( )) / 100);
Quantile = "Q(&q)";
Value = put(quantile('GAMMA', &q, a_post, b_post_scale), 12.6);
output;
%end;
keep Quantile Value;
run;
title "Poisson-Gamma | var: &var | &label_val | Posterior Quantiles";
proc print data=_pg_qtile_ noobs; run;
title;
%end;
data _pg_plot_;
set _pg_grp_;
pr_lo = max(1e-6, quantile('GAMMA', 0.001, a_prior, b_prior));
pr_hi = quantile('GAMMA', 0.999, a_prior, b_prior);
po_lo = max(1e-6, quantile('GAMMA', 0.001, a_post, b_post_scale));
po_hi = quantile('GAMMA', 0.999, a_post, b_post_scale);
do _i_ = 1 to 1000;
lambda_pr = pr_lo + (_i_-1)*(pr_hi-pr_lo)/999;
lambda_po = po_lo + (_i_-1)*(po_hi-po_lo)/999;
prior_density = pdf('GAMMA', lambda_pr, a_prior, b_prior);
post_density = pdf('GAMMA', lambda_po, a_post, b_post_scale);
output;
end;
keep lambda_pr lambda_po prior_density post_density;
run;
ods graphics on;
title "Poisson-Gamma | var: &var | &label_val | Prior [Gamma(&alpha, &beta)]";
proc sgplot data=_pg_plot_;
series x=lambda_pr y=prior_density /
lineattrs=(color="#2C7BB6" thickness=1.5);
xaxis label='lambda';
yaxis label='Density';
run;
title;
title "Poisson-Gamma | var: &var | &label_val | Posterior";
proc sgplot data=_pg_plot_;
series x=lambda_po y=post_density /
lineattrs=(color="#2C7BB6" thickness=1.5);
xaxis label='lambda';
yaxis label='Density';
run;
title;
ods graphics off;
%end;
proc datasets library=work nolist;
delete _pg_summ_ _pg_results_ _pg_grp_ _pg_mom_
_pg_cri_ _pg_qtile_ _pg_plot_;
quit;
%mend conjugate_poisson_gamma;
/* =============================================================================
3. NORMAL - NORMAL (estimating MEAN | variance known)
=============================================================================
Likelihood : X_i ~ N(mu, sigma2) sigma2 known
Prior : mu ~ N(mu0, tau2)
Posterior : mu | x ~ N(mu_post, tau2_post)
tau2_post = 1 / (1/tau2 + n/sigma2)
mu_post = tau2_post * (mu0/tau2 + n*xbar/sigma2)
SAS Normal: quantile('NORMAL', p, mu, sigma) sigma = SD
Parameters:
mu0 = prior mean
tau2 = prior variance
sigma2 = known likelihood variance (fixed)
============================================================================= */
%macro conjugate_normal_normal(
data = ,
var = ,
group_var = ,
x = ,
mu0 = ,
tau2 = ,
sigma2 = ,
cri_levels= 90 95,
quantiles = ,
cri_type = CrI
);
%local i g lev lo hi q nlevs nq ngrps label_val;
%if %length(&data) > 0 %then %do;
%let dsid = %sysfunc(open(&data));
%if &dsid = 0 %then %do;
%put ERROR: Dataset &data not found.;
%return;
%end;
%if %sysfunc(varnum(&dsid, &var)) = 0 %then %do;
%put ERROR: Variable &var not found in dataset &data.;
%let rc = %sysfunc(close(&dsid));
%return;
%end;
%let rc = %sysfunc(close(&dsid));
%if %length(&group_var) > 0 %then %do;
proc sql noprint;
create table _nn_summ_ as
select catx(' = ', "&group_var", &group_var) as _label_ length=200,
mean(&var) as _xbar_,
count(*) as _nobs_
from &data
where &var is not missing
group by &group_var
order by &group_var;
quit;
%end;
%else %do;
proc sql noprint;
create table _nn_summ_ as
select 'All' as _label_ length=200,
mean(&var) as _xbar_,
count(*) as _nobs_
from &data
where &var is not missing;
quit;
%end;
%end;
%else %do;
data _nn_summ_;
length _label_ $200;
_label_ = 'All';
_nobs_ = %sysfunc(countw(&x));
_xbar_ = 0;
%let _nxvals_ = %sysfunc(countw(&x));
%do i = 1 %to &_nxvals_;
_xbar_ = _xbar_ + %scan(&x, &i);
%end;
_xbar_ = _xbar_ / _nobs_;
run;
%end;
data _nn_results_;
set _nn_summ_;
tau2_post = 1 / (1/&tau2 + _nobs_/&sigma2);
mu_post = tau2_post * (&mu0/&tau2 + _nobs_*_xbar_/&sigma2);
pr_mean = &mu0;
pr_var = &tau2;
pr_sd = sqrt(&tau2);
po_mean = mu_post;
po_var = tau2_post;
po_sd = sqrt(tau2_post);
run;
proc sql noprint;
select count(*) into :ngrps trimmed from _nn_results_;
quit;
%do g = 1 %to &ngrps;
data _nn_grp_;
set _nn_results_(firstobs=&g obs=&g);
run;
proc sql noprint;
select _label_ into :label_val trimmed from _nn_grp_;
quit;
title "Normal-Normal | var: &var | &label_val | Data Summary";
proc print data=_nn_grp_ noobs label;
var _nobs_ _xbar_;
label _nobs_ = 'n'
_xbar_ = 'xbar';
format _xbar_ 12.6;
run;
title;
data _nn_mom_;
set _nn_grp_;
length Quantity $40 Prior $15 Posterior $15;
Quantity = 'mu (mean)';
Prior = put(&mu0, 12.6); Posterior = put(mu_post, 12.6); output;
Quantity = 'tau2 (variance)';
Prior = put(&tau2, 12.6); Posterior = put(tau2_post, 12.6); output;
Quantity = 'SD';
Prior = put(pr_sd, 12.6); Posterior = put(po_sd, 12.6); output;
Quantity = 'sigma2 (known)';
Prior = put(&sigma2, 12.6); Posterior = put(&sigma2, 12.6); output;
keep Quantity Prior Posterior;
run;
title "Normal-Normal | var: &var | &label_val | Parameters and Moments";
proc print data=_nn_mom_ noobs; run;
title;
%if %length(&cri_levels) > 0 %then %do;
data _nn_cri_;
set _nn_grp_;
length Level $20 Lower $15 Upper $15;
%let nlevs = %sysfunc(countw(&cri_levels));
%do i = 1 %to &nlevs;
%let lev = %scan(&cri_levels, &i);
%let lo = %sysevalf((1 - &lev/100) / 2);
%let hi = %sysevalf(1 - (1 - &lev/100) / 2);
Level = "&lev% &cri_type";
Lower = put(quantile('NORMAL', &lo, mu_post, po_sd), 12.6);
Upper = put(quantile('NORMAL', &hi, mu_post, po_sd), 12.6);
output;
%end;
keep Level Lower Upper;
run;
title "Normal-Normal | var: &var | &label_val | Credible Intervals (&cri_type)";
proc print data=_nn_cri_ noobs; run;
title;
%end;
%if %length(&quantiles) > 0 %then %do;
data _nn_qtile_;
set _nn_grp_;
length Quantile $15 Value $15;
%let nq = %sysfunc(countw(&quantiles));
%do i = 1 %to &nq;
%let q = %sysevalf(%scan(&quantiles, &i, %str( )) / 100);
Quantile = "Q(&q)";
Value = put(quantile('NORMAL', &q, mu_post, po_sd), 12.6);
output;
%end;
keep Quantile Value;
run;
title "Normal-Normal | var: &var | &label_val | Posterior Quantiles";
proc print data=_nn_qtile_ noobs; run;
title;
%end;
data _nn_plot_;
set _nn_grp_;
pr_lo = quantile('NORMAL', 0.001, &mu0, pr_sd);
pr_hi = quantile('NORMAL', 0.999, &mu0, pr_sd);
po_lo = quantile('NORMAL', 0.001, mu_post, po_sd);
po_hi = quantile('NORMAL', 0.999, mu_post, po_sd);
do _i_ = 1 to 1000;
mu_pr = pr_lo + (_i_-1)*(pr_hi-pr_lo)/999;
mu_po = po_lo + (_i_-1)*(po_hi-po_lo)/999;
prior_density = pdf('NORMAL', mu_pr, &mu0, pr_sd);
post_density = pdf('NORMAL', mu_po, mu_post, po_sd);
output;
end;
keep mu_pr mu_po prior_density post_density;
run;
ods graphics on;
title "Normal-Normal | var: &var | &label_val | Prior [N(&mu0, &tau2)]";
proc sgplot data=_nn_plot_;
series x=mu_pr y=prior_density /
lineattrs=(color="#2C7BB6" thickness=1.5);
xaxis label='mu';
yaxis label='Density';
run;
title;
title "Normal-Normal | var: &var | &label_val | Posterior";
proc sgplot data=_nn_plot_;
series x=mu_po y=post_density /
lineattrs=(color="#2C7BB6" thickness=1.5);
xaxis label='mu';
yaxis label='Density';
run;
title;
ods graphics off;
%end;
proc datasets library=work nolist;
delete _nn_summ_ _nn_results_ _nn_grp_ _nn_mom_
_nn_cri_ _nn_qtile_ _nn_plot_;
quit;
%mend conjugate_normal_normal;
/* =============================================================================
4. NORMAL - GAMMA (estimating PRECISION | mean known)
=============================================================================
Likelihood : X_i ~ N(mu_known, 1/phi) phi = precision
Prior : phi ~ Gamma(alpha, beta) alpha=shape, beta=scale (SAS native)
Posterior : phi | x ~ Gamma(alpha + n/2, beta_post_scale)
SS = sum((x - mu_known)^2)
beta_post_scale = 1 / (1/beta + SS/2)
============================================================================= */
%macro conjugate_normal_gamma(
data = ,
var = ,
group_var = ,
x = ,
mu_known = ,
alpha = ,
beta = ,
cri_levels= 90 95,
quantiles = ,
cri_type = CrI
);
%local i g lev lo hi q nlevs nq ngrps label_val;
%if %length(&data) > 0 %then %do;
%let dsid = %sysfunc(open(&data));
%if &dsid = 0 %then %do;
%put ERROR: Dataset &data not found.;
%return;
%end;
%if %sysfunc(varnum(&dsid, &var)) = 0 %then %do;
%put ERROR: Variable &var not found in dataset &data.;
%let rc = %sysfunc(close(&dsid));
%return;
%end;
%let rc = %sysfunc(close(&dsid));
%if %length(&group_var) > 0 %then %do;
proc sql noprint;
create table _ng_summ_ as
select catx(' = ', "&group_var", &group_var) as _label_ length=200,
count(*) as _nobs_,
sum((&var - &mu_known)**2) as _ss_
from &data
where &var is not missing
group by &group_var
order by &group_var;
quit;
%end;
%else %do;
proc sql noprint;
create table _ng_summ_ as
select 'All' as _label_ length=200,
count(*) as _nobs_,
sum((&var - &mu_known)**2) as _ss_
from &data
where &var is not missing;
quit;
%end;
%end;
%else %do;
data _ng_summ_;
length _label_ $200;
_label_ = 'All';
_nobs_ = %sysfunc(countw(&x));
_ss_ = 0;
%let _nxvals_ = %sysfunc(countw(&x));
%do i = 1 %to &_nxvals_;
_ss_ = _ss_ + (%scan(&x, &i) - &mu_known)**2;
%end;
run;
%end;
data _ng_results_;
set _ng_summ_;
a_prior = α
b_prior = β
a_post = &alpha + _nobs_ / 2;
b_post_scale = 1 / (1/&beta + _ss_/2);
pr_mean = &alpha * β
pr_var = &alpha * &beta**2;
pr_sd = sqrt(pr_var);
pr_median = quantile('GAMMA', 0.5, &alpha, &beta);
if &alpha >= 1 then pr_mode = (&alpha - 1) * β
else pr_mode = .;
po_mean = a_post * b_post_scale;
po_var = a_post * b_post_scale**2;
po_sd = sqrt(po_var);
po_median = quantile('GAMMA', 0.5, a_post, b_post_scale);
if a_post >= 1 then po_mode = (a_post - 1) * b_post_scale;
else po_mode = .;
run;
proc sql noprint;
select count(*) into :ngrps trimmed from _ng_results_;
quit;
%do g = 1 %to &ngrps;
data _ng_grp_;
set _ng_results_(firstobs=&g obs=&g);
run;
proc sql noprint;
select _label_ into :label_val trimmed from _ng_grp_;
quit;
title "Normal-Gamma | var: &var | &label_val | Data Summary";
proc print data=_ng_grp_ noobs label;
var _nobs_ _ss_;
label _nobs_ = 'n'
_ss_ = 'SS = sum((x-mu)^2)';
format _ss_ 12.6;
run;
title;
data _ng_mom_;
set _ng_grp_;
length Quantity $40 Prior $15 Posterior $15;
Quantity = 'alpha (shape)';
Prior = put(a_prior, 12.6); Posterior = put(a_post, 12.6); output;
Quantity = 'beta (scale)';
Prior = put(b_prior, 12.6); Posterior = put(b_post_scale, 12.6); output;
Quantity = 'Mean';
Prior = put(pr_mean, 12.6); Posterior = put(po_mean, 12.6); output;
Quantity = 'Median';
Prior = put(pr_median, 12.6); Posterior = put(po_median, 12.6); output;
Quantity = 'Mode';
if pr_mode = . then Prior = 'undefined'; else Prior = put(pr_mode, 12.6);
if po_mode = . then Posterior = 'undefined'; else Posterior = put(po_mode, 12.6);
output;
Quantity = 'Variance';
Prior = put(pr_var, 12.6); Posterior = put(po_var, 12.6); output;
Quantity = 'SD';
Prior = put(pr_sd, 12.6); Posterior = put(po_sd, 12.6); output;
keep Quantity Prior Posterior;
run;
title "Normal-Gamma | var: &var | &label_val | Parameters and Moments (Precision)";
proc print data=_ng_mom_ noobs; run;
title;
%if %length(&cri_levels) > 0 %then %do;
data _ng_cri_;
set _ng_grp_;
length Level $20 Lower $15 Upper $15;
%let nlevs = %sysfunc(countw(&cri_levels));
%do i = 1 %to &nlevs;
%let lev = %scan(&cri_levels, &i);
%let lo = %sysevalf((1 - &lev/100) / 2);
%let hi = %sysevalf(1 - (1 - &lev/100) / 2);
Level = "&lev% &cri_type";
Lower = put(quantile('GAMMA', &lo, a_post, b_post_scale), 12.6);
Upper = put(quantile('GAMMA', &hi, a_post, b_post_scale), 12.6);
output;
%end;
keep Level Lower Upper;
run;
title "Normal-Gamma | var: &var | &label_val | Credible Intervals (&cri_type) [precision scale]";
proc print data=_ng_cri_ noobs; run;
title;
%end;
%if %length(&quantiles) > 0 %then %do;
data _ng_qtile_;
set _ng_grp_;
length Quantile $15 Value $15;
%let nq = %sysfunc(countw(&quantiles));
%do i = 1 %to &nq;
%let q = %sysevalf(%scan(&quantiles, &i, %str( )) / 100);
Quantile = "Q(&q)";
Value = put(quantile('GAMMA', &q, a_post, b_post_scale), 12.6);
output;
%end;
keep Quantile Value;
run;
title "Normal-Gamma | var: &var | &label_val | Posterior Quantiles [precision scale]";
proc print data=_ng_qtile_ noobs; run;
title;
%end;
data _ng_plot_;
set _ng_grp_;
pr_lo = max(1e-6, quantile('GAMMA', 0.001, a_prior, b_prior));
pr_hi = quantile('GAMMA', 0.999, a_prior, b_prior);
po_lo = max(1e-6, quantile('GAMMA', 0.001, a_post, b_post_scale));
po_hi = quantile('GAMMA', 0.999, a_post, b_post_scale);
do _i_ = 1 to 1000;
phi_pr = pr_lo + (_i_-1)*(pr_hi-pr_lo)/999;
phi_po = po_lo + (_i_-1)*(po_hi-po_lo)/999;
prior_density = pdf('GAMMA', phi_pr, a_prior, b_prior);
post_density = pdf('GAMMA', phi_po, a_post, b_post_scale);
output;
end;
keep phi_pr phi_po prior_density post_density;
run;
ods graphics on;
title "Normal-Gamma | var: &var | &label_val | Prior [Gamma(&alpha, &beta)] (Precision)";
proc sgplot data=_ng_plot_;
series x=phi_pr y=prior_density /
lineattrs=(color="#2C7BB6" thickness=1.5);
xaxis label='phi (precision)';
yaxis label='Density';
run;
title;
title "Normal-Gamma | var: &var | &label_val | Posterior (Precision)";
proc sgplot data=_ng_plot_;
series x=phi_po y=post_density /
lineattrs=(color="#2C7BB6" thickness=1.5);
xaxis label='phi (precision)';
yaxis label='Density';
run;
title;
ods graphics off;
%end;
proc datasets library=work nolist;
delete _ng_summ_ _ng_results_ _ng_grp_ _ng_mom_
_ng_cri_ _ng_qtile_ _ng_plot_;
quit;
%mend conjugate_normal_gamma;
/* =============================================================================
5. NORMAL - INVERSE GAMMA (estimating VARIANCE | mean known)
=============================================================================
Likelihood : X_i ~ N(mu_known, sigma2)
Prior : sigma2 ~ InvGamma(alpha, beta_scale)
Posterior : sigma2 | x ~ InvGamma(alpha + n/2, beta_post_scale)
SS = sum((x - mu_known)^2)
beta_post_scale = beta_scale + SS/2
SAS has no direct InvGamma distribution - derived via Gamma relationship.
quantile InvGamma(p, a, b_scale) = 1 / quantile('GAMMA', 1-p, a, 1/b_scale)
pdf InvGamma(x, a, b_scale) = pdf('GAMMA', 1/x, a, 1/b_scale) / x^2
Parameters:
alpha = prior InvGamma shape
beta_scale = prior InvGamma scale
============================================================================= */
%macro conjugate_normal_igamma(
data = ,
var = ,
group_var = ,
x = ,
mu_known = ,
alpha = ,
beta_scale = ,
cri_levels = 90 95,
quantiles = ,
cri_type = CrI
);
%local i g lev lo hi q nlevs nq ngrps label_val;
%if %length(&data) > 0 %then %do;
%let dsid = %sysfunc(open(&data));
%if &dsid = 0 %then %do;
%put ERROR: Dataset &data not found.;
%return;
%end;
%if %sysfunc(varnum(&dsid, &var)) = 0 %then %do;
%put ERROR: Variable &var not found in dataset &data.;
%let rc = %sysfunc(close(&dsid));
%return;
%end;
%let rc = %sysfunc(close(&dsid));
%if %length(&group_var) > 0 %then %do;
proc sql noprint;
create table _ig_summ_ as
select catx(' = ', "&group_var", &group_var) as _label_ length=200,
count(*) as _nobs_,
sum((&var - &mu_known)**2) as _ss_
from &data
where &var is not missing
group by &group_var
order by &group_var;
quit;
%end;
%else %do;
proc sql noprint;
create table _ig_summ_ as
select 'All' as _label_ length=200,
count(*) as _nobs_,
sum((&var - &mu_known)**2) as _ss_
from &data
where &var is not missing;
quit;
%end;
%end;
%else %do;
data _ig_summ_;
length _label_ $200;
_label_ = 'All';
_nobs_ = %sysfunc(countw(&x));
_ss_ = 0;
%let _nxvals_ = %sysfunc(countw(&x));
%do i = 1 %to &_nxvals_;
_ss_ = _ss_ + (%scan(&x, &i) - &mu_known)**2;
%end;
run;
%end;
data _ig_results_;
set _ig_summ_;
a_prior = α
bs_prior = &beta_scale;
a_post = &alpha + _nobs_ / 2;
bs_post = &beta_scale + _ss_ / 2;
if &alpha > 1 then pr_mean = &beta_scale / (&alpha - 1);
else pr_mean = .;
if &alpha > 2 then pr_var = &beta_scale**2 / ((&alpha-1)**2 * (&alpha-2));
else pr_var = .;
if pr_var ne . then pr_sd = sqrt(pr_var); else pr_sd = .;
pr_median = 1 / quantile('GAMMA', 0.5, &alpha, 1/&beta_scale);
pr_mode = &beta_scale / (&alpha + 1);
if a_post > 1 then po_mean = bs_post / (a_post - 1);
else po_mean = .;
if a_post > 2 then po_var = bs_post**2 / ((a_post-1)**2 * (a_post-2));
else po_var = .;
if po_var ne . then po_sd = sqrt(po_var); else po_sd = .;
po_median = 1 / quantile('GAMMA', 0.5, a_post, 1/bs_post);
po_mode = bs_post / (a_post + 1);
run;
proc sql noprint;
select count(*) into :ngrps trimmed from _ig_results_;
quit;
%do g = 1 %to &ngrps;
data _ig_grp_;
set _ig_results_(firstobs=&g obs=&g);
run;
proc sql noprint;
select _label_ into :label_val trimmed from _ig_grp_;
quit;
title "Normal-IGamma | var: &var | &label_val | Data Summary";
proc print data=_ig_grp_ noobs label;
var _nobs_ _ss_;
label _nobs_ = 'n'
_ss_ = 'SS = sum((x-mu)^2)';
format _ss_ 12.6;
run;
title;
data _ig_mom_;
set _ig_grp_;
length Quantity $40 Prior $15 Posterior $15;
Quantity = 'alpha (shape)';
Prior = put(a_prior, 12.6); Posterior = put(a_post, 12.6); output;
Quantity = 'beta (scale)';
Prior = put(bs_prior, 12.6); Posterior = put(bs_post, 12.6); output;
Quantity = 'Mean';
if pr_mean = . then Prior = 'undefined'; else Prior = put(pr_mean, 12.6);
if po_mean = . then Posterior = 'undefined'; else Posterior = put(po_mean, 12.6);
output;
Quantity = 'Median';
Prior = put(pr_median, 12.6); Posterior = put(po_median, 12.6); output;
Quantity = 'Mode';
Prior = put(pr_mode, 12.6); Posterior = put(po_mode, 12.6); output;
Quantity = 'Variance';
if pr_var = . then Prior = 'undefined'; else Prior = put(pr_var, 12.6);
if po_var = . then Posterior = 'undefined'; else Posterior = put(po_var, 12.6);
output;
Quantity = 'SD';
if pr_sd = . then Prior = 'undefined'; else Prior = put(pr_sd, 12.6);
if po_sd = . then Posterior = 'undefined'; else Posterior = put(po_sd, 12.6);
output;
keep Quantity Prior Posterior;
run;
title "Normal-IGamma | var: &var | &label_val | Parameters and Moments (Variance)";
proc print data=_ig_mom_ noobs; run;
title;
%if %length(&cri_levels) > 0 %then %do;
data _ig_cri_;
set _ig_grp_;
length Level $20 Lower $15 Upper $15;
%let nlevs = %sysfunc(countw(&cri_levels));
%do i = 1 %to &nlevs;
%let lev = %scan(&cri_levels, &i);
%let lo = %sysevalf((1 - &lev/100) / 2);
%let hi = %sysevalf(1 - (1 - &lev/100) / 2);
Level = "&lev% &cri_type";
Lower = put(1 / quantile('GAMMA', 1-&lo, a_post, 1/bs_post), 12.6);
Upper = put(1 / quantile('GAMMA', 1-&hi, a_post, 1/bs_post), 12.6);
output;
%end;
keep Level Lower Upper;
run;
title "Normal-IGamma | var: &var | &label_val | Credible Intervals (&cri_type) [variance scale]";
proc print data=_ig_cri_ noobs; run;
title;
%end;
%if %length(&quantiles) > 0 %then %do;
data _ig_qtile_;
set _ig_grp_;
length Quantile $15 Value $15;
%let nq = %sysfunc(countw(&quantiles));
%do i = 1 %to &nq;
%let q = %sysevalf(%scan(&quantiles, &i, %str( )) / 100);
Quantile = "Q(&q)";
Value = put(1 / quantile('GAMMA', 1-&q, a_post, 1/bs_post), 12.6);
output;
%end;
keep Quantile Value;
run;
title "Normal-IGamma | var: &var | &label_val | Posterior Quantiles [variance scale]";
proc print data=_ig_qtile_ noobs; run;
title;
%end;
data _ig_plot_;
set _ig_grp_;
pr_lo = max(1e-6, (1 / quantile('GAMMA', 1-0.001, a_prior, 1/bs_prior)) * 0.5);
pr_hi = 1 / quantile('GAMMA', 1-0.999, a_prior, 1/bs_prior);
po_lo = max(1e-6, (1 / quantile('GAMMA', 1-0.001, a_post, 1/bs_post)) * 0.5);
po_hi = 1 / quantile('GAMMA', 1-0.999, a_post, 1/bs_post);
do _i_ = 1 to 1000;
s2_pr = pr_lo + (_i_-1)*(pr_hi-pr_lo)/999;
s2_po = po_lo + (_i_-1)*(po_hi-po_lo)/999;
prior_density = pdf('GAMMA', 1/s2_pr, a_prior, 1/bs_prior) / s2_pr**2;
post_density = pdf('GAMMA', 1/s2_po, a_post, 1/bs_post) / s2_po**2;
output;
end;
keep s2_pr s2_po prior_density post_density;
run;
ods graphics on;
title "Normal-IGamma | var: &var | &label_val | Prior [InvGamma(&alpha, &beta_scale)] (Variance)";
proc sgplot data=_ig_plot_;
series x=s2_pr y=prior_density /
lineattrs=(color="#2C7BB6" thickness=1.5);
xaxis label='sigma2 (variance)';
yaxis label='Density';
run;
title;
title "Normal-IGamma | var: &var | &label_val | Posterior (Variance)";
proc sgplot data=_ig_plot_;
series x=s2_po y=post_density /
lineattrs=(color="#2C7BB6" thickness=1.5);
xaxis label='sigma2 (variance)';
yaxis label='Density';
run;
title;
ods graphics off;
%end;
proc datasets library=work nolist;
delete _ig_summ_ _ig_results_ _ig_grp_ _ig_mom_
_ig_cri_ _ig_qtile_ _ig_plot_;
quit;
%mend conjugate_normal_igamma;
/* =============================================================================
USAGE EXAMPLES
=============================================================================*/
/* --- 1. Binomial-Beta : scalar ---*/
%conjugate_binomial_beta(x=14, n_trials=20, alpha=2, beta=2,
cri_levels=%str(80 90 95 99));
/*
/*--- 1b. Binomial-Beta : dataset ---
%conjugate_binomial_beta(data=df, var=var1,
success_level=yes, alpha=2, beta=2, cri_levels=%str(90 95),
quantiles=%str(10 20 30));*/
%conjugate_binomial_beta(data=df, var=var1, group_var=var2,
success_level=yes, alpha=2, beta=2, cri_levels=%str(90 95),quantiles=%str(10 20 30));
/* --- 2. Poisson-Gamma : inline data ---*/
%conjugate_poisson_gamma(x=%str(3 5 2 4 6), alpha=2, beta=1,
cri_levels=%str(80 90 95 99));
/* --- 2b. Poisson-Gamma : dataset ---*/
%conjugate_poisson_gamma(data=mylib.mydata, var=counts, group_var=site,
alpha=2, beta=1, cri_levels=%str(90 95));
/* --- 3. Normal-Normal : inline data ---*/
%conjugate_normal_normal(x=%str(2.1 1.9 2.3 2.0), mu0=2, tau2=1,
sigma2=0.5, cri_levels=%str(90 95));
/* --- 3b. Normal-Normal : dataset + group ---*/
%conjugate_normal_normal(data=df, var=var1,
mu0=2, tau2=1, sigma2=0.5,
cri_levels=%str(90 95), quantiles=%str(10 75));
%conjugate_normal_normal(data=df, var=var1, group_var=var2,
mu0=2, tau2=1, sigma2=0.5,
cri_levels=%str(90 95), quantiles=%str(10 75));
/* --- 4. Normal-Gamma : inline data ---*/
%conjugate_normal_gamma(x=%str(2.1 1.9 2.3 2.0), mu_known=2,
alpha=2, beta=1, cri_levels=%str(90 95));
/* --- 4b. Normal-Gamma : dataset ---*/
%conjugate_normal_gamma(data=df, var=var1, mu_known=2,
alpha=2, beta=1, cri_levels=%str(90 95));
/* --- 5. Normal-IGamma : inline data ---*/
%conjugate_normal_igamma(x=%str(2.1 1.9 2.3 2.0), mu_known=2,
alpha=3, beta_scale=1, cri_levels=%str(90 95));
/* --- 5b. Normal-IGamma : dataset ---*/
%conjugate_normal_igamma(data=df, var=var1, mu_known=2,
alpha=3, beta_scale=1, cri_levels=%str(90 95),
quantiles=%str(1 10));
/*=============================================================================*/
...--------- This is the code help we create a copy button ----------- ...
<div class="copy-content-box">
<button type="button" class="copy-btn" onclick="copyHiddenText(this)">
Copy
</button>
<p class="hidden-copy-text">
This is the paragraph content that will be copied when the user clicks the Copy button. You can put your complete content here.
</p>
</div>
<script>
function copyHiddenText(button) {
const text = button.parentElement.querySelector('.hidden-copy-text').innerText;
navigator.clipboard.writeText(text).then(function() {
const originalText = button.innerText;
button.innerText = "Copied!";
setTimeout(function() {
button.innerText = originalText;
}, 2000);
});
}
</script>
<style>
.hidden-copy-text {
display: none;
}
</style>Binomial – Beta
Show Code
This cide is important …. !
%macro conjugate_binomial_beta(
data = ,
var = ,
group_var = ,
success_level = ,
x = ,
n_trials = ,
alpha = ,
beta = ,
cri_levels = 90 95,
quantiles = ,
cri_type = CrI
);
%local i g lev lo hi q nlevs nq ngrps grp_val label_val;
%if %length(&data) > 0 %then %do;
%let dsid = %sysfunc(open(&data));
%if &dsid = 0 %then %do;
%put ERROR: Dataset &data not found.;
%return;
%end;
%if %sysfunc(varnum(&dsid, &var)) = 0 %then %do;
%put ERROR: Variable &var not found in dataset &data.;
%let rc = %sysfunc(close(&dsid));
%return;
%end;
%let rc = %sysfunc(close(&dsid));
data _bb_raw_;
set &data;
where &var ne '';
%if %length(&success_level) > 0 %then %do;
_success_ = (&var = "&success_level");
%end;
%else %do;
_success_ = .;
%end;
run;
%if %length(&success_level) = 0 %then %do;
proc sql noprint;
select min(&var) into :success_level trimmed
from _bb_raw_;
quit;
data _bb_raw_;
set _bb_raw_;
_success_ = (&var = "&success_level");
run;
%end;
%if %length(&group_var) > 0 %then %do;
proc sql noprint;
create table _bb_summ_ as
select catx(' = ', "&group_var", &group_var) as _label_ length=200,
sum(_success_) as _x_,
count(*) as _nobs_
from _bb_raw_
group by &group_var
order by &group_var;
quit;
%end;
%else %do;
proc sql noprint;
create table _bb_summ_ as
select 'All' as _label_ length=200,
sum(_success_) as _x_,
count(*) as _nobs_
from _bb_raw_;
quit;
%end;
%end;
%else %do;
data _bb_summ_;
_label_ = "1";
_x_ = &x;
_nobs_ = &n_trials;
run;
%end;
data _bb_results_;
set _bb_summ_;
a_prior = α
b_prior = β
a_post = &alpha + _x_;
b_post = &beta + (_nobs_ - _x_);
pr_mean = &alpha / (&alpha + &beta);
pr_var = (&alpha * &beta) / ((&alpha + &beta)**2 * (&alpha + &beta + 1));
pr_sd = sqrt(pr_var);
pr_median = quantile('BETA', 0.5, &alpha, &beta);
if &alpha > 1 and &beta > 1 then
pr_mode = (&alpha - 1) / (&alpha + &beta - 2);
else pr_mode = .;
po_mean = a_post / (a_post + b_post);
po_var = (a_post * b_post) / ((a_post + b_post)**2 * (a_post + b_post + 1));
po_sd = sqrt(po_var);
po_median = quantile('BETA', 0.5, a_post, b_post);
if a_post > 1 and b_post > 1 then
po_mode = (a_post - 1) / (a_post + b_post - 2);
else po_mode = .;
obs_prop = _x_ / _nobs_;
run;
proc sql noprint;
select count(*) into :ngrps trimmed from _bb_results_;
quit;
%do g = 1 %to &ngrps;
data _bb_grp_;
set _bb_results_(firstobs=&g obs=&g);
run;
proc sql noprint;
select _label_ into :label_val trimmed from _bb_grp_;
quit;
title "Binomial-Beta | var: &var | &label_val | Data Summary";
proc print data=_bb_grp_ noobs label;
var _x_ _nobs_ obs_prop;
label _x_ = 'Successes'
_nobs_ = 'Trials'
obs_prop= 'Observed Proportion';
format obs_prop 8.6;
run;
title;
data _bb_mom_;
set _bb_grp_;
length Quantity $40 Prior $15 Posterior $15;
Quantity = 'alpha';
Prior = put(a_prior, 12.6); Posterior = put(a_post, 12.6); output;
Quantity = 'beta';
Prior = put(b_prior, 12.6); Posterior = put(b_post, 12.6); output;
Quantity = 'Mean';
Prior = put(pr_mean, 12.6); Posterior = put(po_mean, 12.6); output;
Quantity = 'Median';
Prior = put(pr_median, 12.6); Posterior = put(po_median, 12.6); output;
Quantity = 'Mode';
if pr_mode = . then Prior = 'undefined'; else Prior = put(pr_mode, 12.6);
if po_mode = . then Posterior = 'undefined'; else Posterior = put(po_mode, 12.6);
output;
Quantity = 'Variance';
Prior = put(pr_var, 12.6); Posterior = put(po_var, 12.6); output;
Quantity = 'SD';
Prior = put(pr_sd, 12.6); Posterior = put(po_sd, 12.6); output;
keep Quantity Prior Posterior;
run;
title "Binomial-Beta | var: &var | &label_val | Parameters and Moments";
proc print data=_bb_mom_ noobs; run;
title;
%if %length(&cri_levels) > 0 %then %do;
data _bb_cri_;
set _bb_grp_;
length Level $20 Lower $15 Upper $15;
%let nlevs = %sysfunc(countw(&cri_levels));
%do i = 1 %to &nlevs;
%let lev = %scan(&cri_levels, &i);
%let lo = %sysevalf((1 - &lev/100) / 2);
%let hi = %sysevalf(1 - (1 - &lev/100) / 2);
Level = "&lev% &cri_type";
Lower = put(quantile('BETA', &lo, a_post, b_post), 12.6);
Upper = put(quantile('BETA', &hi, a_post, b_post), 12.6);
output;
%end;
keep Level Lower Upper;
run;
title "Binomial-Beta | var: &var | &label_val | Credible Intervals (&cri_type)";
proc print data=_bb_cri_ noobs; run;
title;
%end;
%if %length(&quantiles) > 0 %then %do;
data _bb_qtile_;
set _bb_grp_;
length Quantile $15 Value $15;
%let nq = %sysfunc(countw(&quantiles));
%do i = 1 %to &nq;
%let q = %sysevalf(%scan(&quantiles, &i, %str( )) / 100);
Quantile = "Q(&q)";
Value = put(quantile('BETA', &q, a_post, b_post), 12.6);
output;
%end;
keep Quantile Value;
run;
title "Binomial-Beta | var: &var | &label_val | Posterior Quantiles";
proc print data=_bb_qtile_ noobs; run;
title;
%end;
data _bb_plot_;
set _bb_grp_;
if min(a_prior, b_prior) < 1 then do; pr_lo=0.0001; pr_hi=0.9999; end;
else do; pr_lo=0; pr_hi=1; end;
if min(a_post, b_post) < 1 then do; po_lo=0.0001; po_hi=0.9999; end;
else do; po_lo=0; po_hi=1; end;
do _i_ = 1 to 1000;
theta_pr = pr_lo + (_i_-1)*(pr_hi-pr_lo)/999;
theta_po = po_lo + (_i_-1)*(po_hi-po_lo)/999;
prior_density = pdf('BETA', theta_pr, a_prior, b_prior);
post_density = pdf('BETA', theta_po, a_post, b_post);
output;
end;
keep theta_pr theta_po prior_density post_density;
run;
ods graphics on;
title "Binomial-Beta | var: &var | &label_val | Prior [Beta(&alpha, &beta)]";
proc sgplot data=_bb_plot_;
series x=theta_pr y=prior_density /
lineattrs=(color="#2C7BB6" thickness=1.5);
xaxis label='theta';
yaxis label='Density';
run;
title;
title "Binomial-Beta | var: &var | &label_val | Posterior";
proc sgplot data=_bb_plot_;
series x=theta_po y=post_density /
lineattrs=(color="#2C7BB6" thickness=1.5);
xaxis label='theta';
yaxis label='Density';
run;
title;
ods graphics off;
%end;
proc datasets library=work nolist;
delete _bb_raw_ _bb_summ_ _bb_results_ _bb_grp_ _bb_mom_
_bb_cri_ _bb_qtile_ _bb_plot_;
quit;
%mend conjugate_binomial_beta;
Bayesian Conjugate Models — SAS Macros
1
Binomial-Beta
2
Poisson-Gamma
3
Normal-Normal
4
Normal-Gamma
5
Normal-Inverse-Gamma