*! Date : 2 Mar 2011 *! Version : 1.01 *! Author : Adrian Mander *! Email : adrian.mander@mrc-bsu.cam.ac.uk *! For teaching how to find the Beta prior of choice *! *! betaprior, mean() var() /* 28Feb2011 v1.00 The command is born 2Mar2011 v1.01 bug fix */ program define betaprior, rclass version 10.0 preserve syntax [, Mean(real 0.2) Variance(real 0.1) Graph] /******************************************** * Mean of beta is a/(a+b) * Var of beta(a,b) is ab/((a+b)^2*(a+b+1)) ********************************************/ if `mean'>=1 | `mean'<=0 { di "{err}ERROR: you can not have a mean outside the interval (0,1)" exit(196) } di "{txt}The prior mean is {res}`mean'" di "{txt} and variance is {res}`variance'" local a = ((`mean')^2-(`mean')^3-`mean'*`variance')/`variance' if `a'<0 { di "{err}ERROR: the variance is too large for this mean, {res} `mean'" di "{err}Select a smaller variance" exit(196) } local b = `a'/`mean'-`a' local a3 : di %5.3f `a' local b3 : di %5.3f `b' local mn = `a'/(`a'+`b') local va = `a'*`b'/((`a'+`b')^2*(`a'+`b'+1)) di "{txt}Your prior is {res}Beta(`a3',`b3') twoway function y=betaden(`a3',`b3',x), title(Beta(`a3',`b3')) yscale(off) restore end