betaHPD                 package:pscl                 R Documentation

_c_o_m_p_u_t_e _a_n_d _o_p_t_i_o_n_a_l_l_y _p_l_o_t _b_e_t_a _H_D_R_s

_D_e_s_c_r_i_p_t_i_o_n:

     Compute and optionally plot highest density regions for the Beta
     distribution.

_U_s_a_g_e:

        betaHPD(alpha,beta,p=.95,plot=FALSE,xlim=NULL,debug=FALSE)

_A_r_g_u_m_e_n_t_s:

   alpha: scalar, first shape parameter of the Beta density.  Must be
          greater than 1, see details

    beta: scalar, second shape parameter of the Beta density.  Must be
          greater than 1, see details

       p: scalar, content of HPD, must lie between 0 and 1

    plot: logical flag, if 'TRUE' then plot the density and show the
          HDR

    xlim: numeric vector of length 2, the limits of the density's
          support to show when plotting; the default is 'NULL', in
          which case the function will confine plotting to where the
          density is non-neglible

   debug: logical flag, if 'TRUE' produce messages to the console

_D_e_t_a_i_l_s:

     The Beta density arises frequently in Bayesian models of binary
     events, rates, and proportions, which take on values in the open
     unit interval.  For instance, the Beta density is a conjugate
     prior for the unknown success probability in binomial trials. 
     With shape parameters alpha > 1 and beta > 1, the Beta density is
     unimodal.

     In general, suppose theta in Theta subseteq R^k is a random
     variable with density f(theta).  A highest density region (HDR) of
     f(theta) with content p in (0,1] is a set mathcal{Q} subseteq
     Theta with the following properties:

                  int_mathcal{Q} f(theta) dtheta = p

     and

 f(theta) > f(theta^*) , forall theta in mathcal{Q}, theta^* notin mathcal{Q}.

     For a unimodal Beta density (the class of Beta densities handled
     by this function), a HDR of content 0 < p < 1 is simply an
     interval mathcal{Q} in (0,1).

     This function uses numerical methods to solve for the end points
     of a HDR for a Beta density with user-specified shape parameters,
     via repeated calls to the functions 'dbeta', 'pbeta' and 'qbeta'.
     The function 'optimize' is used to find points v and w such that 

                             f(v) = f(w)

     subject to the constraint

              int_v^w f(theta; alpha, beta) dtheta = p,

     where f(theta; alpha, beta) is a Beta density with shape
     parameters alpha and beta. 

     In the special case of alpha = beta > 1, the end points of a HDR 
     with content p are given by the (1 pm p)/2 quantiles of the Beta
     density, and are computed with the 'qbeta' function. 

     Again note that the function will only compute a HDR for a
     unimodal Beta density, and exit with an error if 'alpha<=1 | beta
     <=1'. Note that the uniform density results with alpha = beta = 1,
     which does not have a unique HDR with content 0 < p < 1.  With
     shape parameters alpha<1 and beta>1 (or vice-versa, respectively),
     the Beta density is infinite at 0 (or 1, respectively), but still
     integrates to one, and so a HDR is still well-defined (but not
     implemented here, at least not yet). Similarly, with 0 < alpha,
     beta < 1 the Beta density is infinite at both 0 and 1, but
     integrates to one, and again a HDR of content p<1 is well-defined
     in this case, but will be a set of two disjoint intervals (again,
     at present, this function does not cover this case).

_V_a_l_u_e:

     If the numerical optimization is successful an vector of length 2,
     containing v and w, defined above.    If the optimization fails
     for whatever reason, a vector of 'NAs' is returned.

     The function will also produce a plot of the density with area
     under the density supported by the HDR shaded, if the user calls
     the function with 'plot=TRUE'; the plot will appear on the current
     graphics device.

     Debugging messages are printed to the console if the 'debug'
     logical flag is set to 'TRUE'.

_A_u_t_h_o_r(_s):

     Simon Jackman jackman@stanford.edu.  Thanks to John Bullock who
     discovered a bug in an earlier version.

_S_e_e _A_l_s_o:

     'pbeta', 'qbeta', 'dbeta', 'uniroot'

_E_x_a_m_p_l_e_s:

     betaHPD(4,5)
     betaHPD(2,120)
     betaHPD(120,45,p=.75,xlim=c(0,1))

