MTB > # summarizing a beta density by simulation
MTB > #
MTB > # suppose I want to learn about p=Prob(coin lands heads)
MTB > # I use a uniform prior and toss the coin 100 times, getting
MTB > # 108 heads and 92 tails
MTB > #
MTB > # posterior density is beta(109, 93)
MTB > # I want to compute P(.48<p<.52) and find a 90% interval estimate
MTB > # for p
MTB > #
MTB > # We simulate 1000 values of a beta(109, 93) distribution
MTB > # and summarize the simulated sample.
MTB > 
MTB > name c1 'post'
MTB > rand 1000 'post';
SUBC> beta 109 93.
MTB > dotplot 'post'

Each dot represents 4 points

                                         .
                                       : :
                                    :..::::
                                  ..::::::::: .
                                :.::::::::::: :
                               ::::::::::::::::.. .
                            .:.:::::::::::::::::: :.
                          .::::::::::::::::::::::::::.
                 .  ......::::::::::::::::::::::::::::::.... .
          -+---------+---------+---------+---------+---------+-----post    
       0.400     0.450     0.500     0.550     0.600     0.650


MTB > let k1=sum('post'<.48)/1000
MTB > let k2=sum('post'<.52)/1000
MTB > prin k1 k2
K1       0.0370000
K2       0.269000
MTB > sort 'post' 'post'
MTB > let k3='post'(25)
MTB > let k4='post'(975)
MTB > prin k3 k4
K3       0.475152
K4       0.610683