wave.band              package:waveband              R Documentation

_P_o_s_t_e_r_i_o_r _c_r_e_d_i_b_l_e _i_n_t_e_r_v_a_l_s _f_o_r _w_a_v_e_l_e_t _r_e_g_r_e_s_s_i_o_n

_D_e_s_c_r_i_p_t_i_o_n:

     Computes posterior credible intervals for an unknown regression
     curve.

_U_s_a_g_e:

     wave.band(data = 0, alpha = 0.5, beta = 1., filter.number = 8, family = 
             "DaubLeAsymm", bc = "periodic", dev = var, j0 = 3., plotfn = TRUE, 
             retvalue = TRUE, n = 128, type = "data", rsnr = 3)

_A_r_g_u_m_e_n_t_s:

     Either data or a value of type other than "data" must be supplied. 

    data: If type="data", then data should be a vector of data. The
          length of the vector should be a power of two not greater
          than 1024.

    type: Either type="data", in which case a vector of data should be
          supplied, or type should specify a standard test function and
          wave.band will generate a test data set via a call to
          test.data. Permissible values for type are "blocks", "bumps",
          "doppler", "heavi", or "ppoly"; see the documentation for
          test.data for more details.

alpha, beta: Hyperparameters which determine the priors placed on the
          wavelet coefficients. Both alpha and beta take positive
          values; see Abramovich, Sapatinas, & Silverman (1998) or
          Chipman & Wolfson (1999) for more details on selecting alpha
          and beta.

filter.number: A parameter relating to the smoothness of wavelet that
          you want to use in the decomposition.

 family : Specifies the family of wavelets to be used. Two popular
          options are "DaubExPhase" and "DaubLeAsymm" but see the help
          for filter.select for more possibilities.

      bc: Specifies the boundary handling. If bc="periodic" the
          default, then the function you decompose is assumed to be
          periodic on it's interval of definition. Other boundary
          options exist, but are currently unsupported for this
          function.

     dev: This argument supplies the function to be used to compute the
          spread of the absolute values coefficients. The function
          supplied must return a value of spread on the variance scale
          (i.e. not standard deviation) such as the var() function. A
          popular, useful and robust alternative is the madmad
          function.

      j0: The primary resolution level; used in assessing the universal
          threshold which is used in the empirical Bayes estimation of
          hyperparameters.

  plotfn: If 'plotfn=TRUE', wave.band draws the noisy data, the
          BayesThresh function estimate, and pointwise 99 percent
          credible intervals for the regression function. If the value
          of 'type' is not "data", then the true function will also be
          plotted.

retvalue: If 'retvalue=TRUE', then a lengthy list of results will be
          returned. Note that if both 'plotfn' and 'retvalue' are set
          to 'FALSE', then 'wave.band' will return no results
          whatsoever.

       n: If 'type' is not "data", then a data vector of length 'n'
          will be generated; note that 'n' should be a power of two not
          greater than 1024.

    rsnr: If type is not "data", then the data vector generated will
          have root signal-to-noise ratio as specified by 'rsnr'. 

_D_e_t_a_i_l_s:

     This function implements the WaveBand method of Barber, Nason, &
     Silverman (2001) to compute posterior credible intervals for a
     regression function. The credible intervals are found by
     approximating the posterior distribution of the estimated
     regression curve at each design point. A mixture prior with two
     components (a zero-mean normal and a point mass at zero) is placed
     on each wavelet coefficient and updated by the data to give the
     posteriors for the wavelet coefficients. This is the same prior
     used by Abramovich, Sapatinas, & Silverman (1998) in their
     BayesThresh method, implemented in the function 'BAYES.THR'. 

     The cumulants of these posteriors are computed and stored in the
     'wd.objects' returned by 'wave.band' as 'Kr.wd'. These are summed
     to give the posterior cumulants of the regression curve, which are
     used to fit a Johnson distribution (Johnson, 1949), using the
     algorithm of Hill, Hill, & Holder (1976). Percentage points of
     these distributions are computed by the algorithm of Hill (1976)
     and give the credible intervals themselves. 

     Code to implement the algorithms by Hill (1976) and Hill, Hill, &
     Holder (1976) was obtained from the StatLib archive.

_V_a_l_u_e:

     If 'retvalue=FALSE', the value returned by 'wave.band' is 'NULL'.
     Otherwise, 'wave.band' returns a list with the following
     components: 

    data: The data vector which has been analysed. 

     nts: A list containing four vectors named one, two, three, and
          four. Vector one contains the first cumulants of the
          regression function estimate, vector to the second cumulants
          and so on.

   Kr.wd: A list of four wd objects. These contain the first to fourth
          cumulants of the wavelet coefficients, as well as recording
          the wavelet used in the decomposition. 

   bands: A list containing pointwise upper and lower credible limits
          for the regression function estimate for nominal coverage
          rates 80, 90, 95 and 99 percent. The widths of the credible
          intervals is also included. The vectors are named with "l",
          "u", and "w" indicating lower limits, upper limits, and
          intervals widths, while "80", "90", "95", and "99" refer to
          the nominal coverage rate. 

   param: A record of parameters in the call to 'wave.band'.

_S_I_D_E _E_F_F_E_C_T_S:

     If 'plotfn=TRUE', results are plotted on the current graphics
     device.

_S_e_e _A_l_s_o:

     'BAYES.THR', 'plot.wb', 'power.sum', 'test.data'

_E_x_a_m_p_l_e_s:

     #
     # First, look at the piecewise polynomial example.
     #
     # This plot and the plots for the smooth example below show
     # the data as points, the BayesThresh estimate (thick line), 
     # pointwise 99 percent credible intervals (thin lines), and
     # the true function (dotted thin line).
     #
     #ppoly.wb <- wave.band(type = "ppoly", n = 1024, rsnr=4)
     #
     # Plotting the cumulants shows that there are significant 
     # third and fourth cumulants in some places.
     #
     #t <- (1:1024)/1024
     #plot(t, ppoly.wb$cumulants$one, type="l", xlab="t", ylab = "one")
     #plot(t, ppoly.wb$cumulants$two, type="l", xlab="t", ylab = "two")
     #plot(t, ppoly.wb$cumulants$three, type="l", xlab="t", ylab = "three")
     #plot(t, ppoly.wb$cumulants$four, type="l", xlab="t", ylab = "four")
     #
     # Now consider how much difference the prior can make.
     # Consider a smooth example, first using the default prior,
     # and then using a smoother prior.
     #

     #gs <- sin(2*pi*t) + 2*(t - 0.5)^2
     #gs.noisy <- gs + rnorm(n=1024, sd=sqrt(var(gs))/2)
     #gs.wb1 <- wave.band(data=gs.noisy)

     #gs.wb2 <- wave.band(data=gs.noisy, alpha=4, beta=1)

