#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh 'README' <<'END_OF_FILE'
X================================================================================
XSUMMARY
X
X
Xprop.ss is a S-PLUS program for calculating sample size requirements using
Xthree different Bayesian criteria in the context of designing an experiment
Xto estimate the difference between two binomial proportions.
X
XThis program was written by Patrick Belisle, Roxane du Berger and
XLawrence Joseph and developed in a SunOS 4.1.3 operating environment
Xon a SPARC station ELC.
X
XThis program contains the complete set of functions necessary to
Xcompute sample sizes for binomial data satisfying any of
Xthe criteria presented in the reference listed below. In particular,
XS-PLUS functions that calculate the required sample sizes for the Average
XLength Criterion, the Average Coverage Criterion and the Worst Outcome
XCriterion in the context of binomial observations are provided.
XIn all cases, estimation of the difference between two binomial
Xproportions is considered.
XFunctions for both the fully Bayesian and the mixed Bayesian/likelihood
Xapproaches are provided. In addition, equally sized groups or unequally
Xsized groups can be specified. In the latter case, optimal sample sizes
Xare selected, in the sense that the variance of the posterior density for
Xthe difference in proportions is minimized.
X
X
XThe results obtained from prop.ss are based on an approximation to
Xthe exact posterior distribution of the difference between two
Xbinomial proportions; however, a FORTRAN 77 program (by Lawrence Joseph
Xand Roxane du Berger) providing sample sizes based on the exact posterior
Xis available from the authors. prop.ss has the advantage of being much
Xfaster to run and still excellent at approximating sample size requirements.
X
X
X------------------------------------------------------------------------
XOther sample size algorithms
X
X
Xbhpd1: A Fortran 77 program to calculate sample sizes based on highest
Xposterior density intervals in the context of a single binomial experiment
Xusing the same approaches as in prop.ss; it is available in the directory
X'general' of StatLib (submitted by Lawrence Joseph, David Wolfson and
XRoxane du Berger).
X
Xmean.ss: A S-PLUS program to calculate sample sizes based on highest posterior
Xdensity intervals for normal means or difference between two normal means using
Xthe same Bayesian criteria as above is available in the directory 'S' of StatLib
X(submitted by Lawrence Joseph and Patrick Belisle).
X
X
X================================================================================
XDOCUMENTATION
X
XDetailed information about the implementation of the Bayesian
Xapproaches to sample size requirements in the context of designing
Xan experiment to estimate the difference between two binomial
Xproportions are found in our manuscript:
X
X
X ``Bayesian and mixed Bayesian/likelihood criteria for
X sample size determination'',
X Lawrence Joseph, Roxane du Berger and Patrick Belisle
X
XQuestions or comments regarding the program prop.ss may be e-mailed to
XLawrence Joseph at
X
X joseph@binky.epi.mcgill.ca
X
XRegular mail requests may be sent to:
X
X Lawrence Joseph
X Division of Clinical Epidemiology
X Montreal General Hospital
X 1650 Cedar Avenue
X Montreal, Quebec, CANADA
X H3G 1A4
X
X
X
X================================================================================
XEXAMPLE OF EXECUTION
X
XThe program can be loaded by typing source('prop.ss') at the S-PLUS prompt;
Xa list of functions along with their respective arguments can always be obtained
Xby typing list.of.functions.prop (without brackets).
X
XResults for example 1 in our paper are obtained by typing:
X
X> acc(len=0.05, level=0.95, c1=3, d1=11, c2=11, d2=54)
X[1] 1809 1809
X> alc(0.05, 0.95, 3, 11, 11, 54)
X[1] 1751 1751
X> woc(0.05, 0.95, 3, 11, 11, 54)
X[1] 3033 3033
X> modwoc(0.05,0.95,3,11,11,54,worst.level=0.95)
X[1] 2582 2582
X> modwoc(0.05, 0.95, 3, 11, 11, 54, 0.99)
X[1] 2687 2687
X> freq(0.05, 0.95, 1, 1, 1, 1)
X[1] 3074
X> freq(0.05, 0.95, 3, 11, 11, 54)
X[1] 1899
X
XNote that sample size requirements obtained above may not be exactly those
Xpublished, due to Monte Carlo error.
X
XThe sample size provided by acc ensures that the probability coverage of the CI
Xof length 0.05 will be on average 95%, using prior information equivalent to
Xobserving 3 successes and 11 failures for the first treatment/drug and 11
Xsuccesses and 54 failures for the second. ALC sample size ensures that CI
Xof probability coverage of 95% will be on average of length 0.05, while WOC
Xensures that the reported length will always be smaller than 0.05. MODWOC
X(with worst-levels equal to 0.95 and 0.99) sample sizes ensure that the
Xreported lengths will be smaller than 0.05 with probability 95% and 99%.
X
XThe two last sample sizes given above are the frequentist sample sizes when
Xfilling in 1) the worst values for pi1 and pi2 in the standard sample size
Xformula and 2) likely values for pi1 and pi2, the probabilities of
Xsuccess for the two different treatments/drugs.
X
X
X
X================================================================================
XCOPYRIGHT
X
X(c) Copyright Lawrence Joseph and Patrick Belisle 1995
X
Xprop.ss is a program written by Lawrence Joseph, Roxane du Berger and
XPatrick Belisle at the Division of Clinical Epidemiology,
XDepartment of Medicine, Montreal General Hospital. This program
Xis an implementation of the submitted manuscript ``Bayesian and mixed
XBayesian/likelihood criteria for sample size determination''.
X
XYou are free to use this program, for non-commercial purposes only,
Xunder two conditions:
X
X(1) This note is not to be removed;
X(2) Publications using prop.ss results should reference the manuscript
X mentioned above.
X
X
XPlease e-mail all comments or questions to
X
X joseph@binky.epi.mcgill.ca
X
X--------------------------------------------------------------------------------
END_OF_FILE
if test 5896 -ne `wc -c <'README'`; then
echo shar: \"'README'\" unpacked with wrong size!
fi
# end of 'README'
fi
if test -f 'prop.ss' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'prop.ss'\"
else
echo shar: Extracting \"'prop.ss'\" \(44593 characters\)
sed "s/^X//" >'prop.ss' <<'END_OF_FILE'
X#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X#
X# prop.ss is a S-PLUS program for calculating sample size requirements using
X# three different Bayesian criteria in the context of designing an experiment
X# to estimate the difference between two binomial proportions.
X#
X# This program was written by Patrick Belisle, Roxane du Berger
X# and Lawrence Joseph and developed in a SunOS 4.1.3 operating environment
X# on a SPARC station ELC.
X#
X# This program contains the complete set of functions necessary to
X# compute sample sizes for binomial data satisfying any of
X# the criteria presented in the reference listed below. In particular,
X# S-PLUS functions that calculate the required sample sizes for the Average
X# Length Criterion, the Average Coverage Criterion and the Worst Outcome
X# Criterion in the context of binomial observations are provided.
X# In all cases, estimation of the difference between two binomial
X# proportions is considered.
X# Functions for both the fully Bayesian and the mixed Bayesian/likelihood
X# approaches are provided. In addition, equally sized groups or unequally
X# sized groups can be specified. In the latter case, optimal sample sizes
X# are selected, in the sense that the variance of the posterior density for
X# the difference in proportions is minimized.
X#
X#
X# The results obtained from prop.ss are based on an approximation
X# of the exact posterior distribution of the difference between two
X# binomial proportions; however, a Fortran 77 program (by Lawrence Joseph
X# and Roxane du Berger) providing sample sizes based on the exact posterior
X# is available from the authors. prop.ss has the advantage of being much
X# faster to run and still excellent at approximating sample size requirements.
X#
X#
X# Detailed information about the implementation of the Bayesian
X# approaches to sample size requirements in the context of designing
X# an experiment to estimate the difference between two binomial
X# proportions are found in our manuscript:
X#
X#
X# ``Bayesian and mixed Bayesian/likelihood criteria for
X# sample size determination'',
X# Lawrence Joseph, Roxane du Berger and Patrick Belisle
X#
X#
X# (c) Copyright Lawrence Joseph and Patrick Belisle 1995
X#
X# prop.ss is a program written by Lawrence Joseph, Roxane du Berger and
X# Patrick Belisle at the Division of Clinical Epidemiology,
X# Department of Medicine, Montreal General Hospital. This program
X# is an implementation of the submitted manuscript ``Bayesian and mixed
X# Bayesian/likelihood criteria for sample size determination''.
X#
X# You are free to use this program, for non-commercial purposes only,
X# under two conditions:
X#
X# (1) This note is not to be removed;
X# (2) Publications using prop.ss results should reference the manuscript
X# mentioned above.
X#
X#
X# Please e-mail all comments or questions to
X#
X# joseph@binky.epi.mcgill.ca
X#
X#
X#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X
X
X
X#################################################################################
X# List of functions to be called by user, with corresponding parameters
X#
X
Xlist.of.functions.prop_function()
X{
X# Below are listed all functions of interest for sample size
X# calculations when comparing two binomial proportions; parameters with
X# values in brackets (the default values) are optional.
X#
X# All parameters are presented in our paper.
X# For further information, feel free to contact authors.
X#
X#
X# ====================================================================
X#
X# D I F F E R E N C E B E T W E E N 2 P R O P O R T I O N S
X#
X# ====================================================================
X#
X#
X# Pure bayesian approach
X# ______________________
X#
X#
X# Average coverage criterion
X# acc *
X# length, level, c1, d1, c2, d2, opt (F), m (1000), mcs (3)
X#
X#
X# Average length criterion
X# alc
X# length, level, c1, d1, c2, d2, opt (F), m (1000), mcs (3)
X#
X#
X# Worst outcome criterion
X# woc
X# length, level, c1, d1, c2, d2, opt (F)
X#
X#
X# Modified worst outcome criterion
X# modwoc
X# length, level, c1, d1, c2, d2, worst.level (0.95), opt (F)
X#
X#
X#
X# Mixed Bayesian / likelihood
X# ___________________________
X#
X#
X# Average coverage criterion
X# mblacc
X# length, level, c1, d1, c2, d2, m (1000), mcs (3)
X#
X#
X# Average length criterion
X# mblalc
X# length, level, c1, d1, c2, d2, m (1000), mcs (3)
X#
X#
X# Worst outcome criterion
X# mblwoc
X# length, level, c1, d1, c2, d2
X#
X#
X# Modified worst outcome criterion
X# mblmodwoc
X# length, level, c1, d1, c2, d2, worst.level (0.95)
X#
X#
X#
X# Frequentist [ with point estimates c1/(c1+d1) and c2/(c2+d2) ]
X# ___________
X#
X# freq
X# length, level, c1, d1, c2, d2
X#
X#
X#
X# * Note:
X# if opt = F, then n1 = n2
X# if opt = T, then (n1,n2) minimize the expected posterior
X# variance given a total of n1+n2 observations
X}
X
X
X#******************************************************************************
X
X
Xfreq_function(len,level,c1,d1,c2,d2)
X{
X ceiling(4*qnorm((1+level)/2)^2*(c1*d1/(c1+d1)^2+c2*d2/(c2+d2)^2)/len^2)
X}
X
X
X
X###############################################################################
X#
X# A V E R A G E C O V E R A G E
X#
X#
X
X
Xacc_function(len,level,c1,d1,c2,d2,opt=F,m=1000,mcs=3)
X{
X # mcs is standing for 'maximal consecutive steps'(in same direction)
X
X min.for.possible.return_2^ceiling(1.5*mcs)
X
X # If we always allow a return, there is a risk of making bad steps
X # when we are close to the answer.
X # Thus, we should not allow any return once some arbitrary 'step' (which is
X # 'min.for.possible.return') is reached.
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #**************************************************************************
X # Define initial step
X # (as a function of any frequentist sample size estimate)
X
X step <- ceiling(log(freq(len,level,c1,d1,c2,d2))/log(2))
X
X # Also define the threshold to cross for the quantity under study (the
X # length or the coverage probability of an HPD region)
X
X threshold_level
X
X # and define a factor, which is +/- 1, depending on if the quantity under
X # study is (almost) surely too large or too small when making no
X # observations [-1 if the quantity to measure is DEcreasing with n
X # +1 if the quantity to measure is INcreasing with n]
X #
X # [ -1 if threshold_len, +1 if thresold_level ]
X
X factor_+1
X
X #**************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X quantity.to.measure_ifelse(factor == +1,0,2*threshold)
X
X step <- 2^step
X
X history.ns_0
X history.steps_0
X
X n1_0
X
X max.cons.steps.same.dir_mcs
X found.upper.bound_F
X possible.to.move.back_T
X cons.steps.same.dir_0
X direction_+1
X
X while(step>=1)
X {
X while(sign(factor*(threshold-quantity.to.measure)) == direction && step >= 1)
X {
X step[found.upper.bound]_max(1,step/2)
X possible.to.move.back[step < min.for.possible.return &&
X found.upper.bound]_F
X
X
X
X n1_n1+direction*step
X
X cons.steps.same.dir_cons.steps.same.dir+1
X
X history.ns_c(n1,history.ns)
X history.steps_c(step*direction,history.steps)
X
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #************************************************************************
X # Define n2 from n1
X
X n2_ifelse(!opt,n1,ceiling(sqrt(c2*d2/c1/d1*(c1+d1)*(c1+d1+1)/(c2+d2)/
X (c2+d2+1))*(c1+d1+n1)-c2-d2))
X
X n2[n2<0]_0
X
X #**********************************
X # Here compute the average coverage
X
X pi1_rbeta(m,c1,d1)
X pi2_rbeta(m,c2,d2)
X x1_rbinom(m,n1,pi1)
X x2_rbinom(m,n2,pi2)
X
X # Precaution: if n1 or n2 is 0, then the correponding x given
X # by rbinom is NA. We correct the situation by setting it to 0.
X
X x1[is.na(x1)]_0
X x2[is.na(x2)]_0
X
X # Posterior variance of the difference between the two proportions
X
X posterior.mean_(c2+x2)/(c2+d2+n2)-(c1+x1)/(c1+d1+n1)
X posterior.var_(x1+c1)*(n1-x1+d1)/(n1+c1+d1)^2/(n1+c1+d1+1)+
X (x2+c2)*(n2-x2+d2)/(n2+c2+d2)^2/(n2+c2+d2+1)
X
X # We make the approximation of the difference between two betas by
X # another beta, with parameters a and b given by
X
X a_-1/2 * (posterior.mean + 1) * (posterior.mean^2 - 1 + posterior.var)/
X posterior.var
X b_1/2*(1-posterior.mean)*(-1*posterior.mean^2+1-posterior.var)/
X posterior.var
X
X # and approximate the probability coverage of the HPD region by the
X # coverage obtained by a symetric region around the posterior mean
X
X pi.min_posterior.mean-len/2
X pi.max_posterior.mean+len/2
X pi.max[pi.min< -1]_-1+len ; pi.min[pi.min< -1]_-1
X pi.min[pi.max>1]_1-len ; pi.max[pi.max>1]_1
X coverage.around.mean_pbeta((pi.max+1)/2,a,b)-pbeta((pi.min+1)/2,a,b)
X
X quantity.to.measure_mean(coverage.around.mean)
X
X #************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X# 2
X if(found.upper.bound &&
X cons.steps.same.dir == max.cons.steps.same.dir+1 &&
X possible.to.move.back)
X {
X
X if(sign(factor*(threshold-quantity.to.measure)) == direction)
X {
X # There was (most likely) a mistake, look for n's in the other direction
X at.n_seq(along=history.ns)[history.ns == n1 ]
X hs.an_history.steps[at.n]
X
X step_abs(hs.an[sign(hs.an) != direction][1])
X
X cons.steps.same.dir_0
X
X # and if there has never been a step coming from the other direction...
X step[is.na(step)]_max(abs(hs.an))
X }
X# 3
X else
X {
X # There was (most likely) no mistake; keep looking around the same n's
X
X direction_-direction
X cons.steps.same.dir_0
X }
X }
X
X# 1
X if(found.upper.bound &&
X cons.steps.same.dir==max.cons.steps.same.dir &&
X sign(factor*(threshold-quantity.to.measure))==direction &&
X possible.to.move.back)
X {
X step_2*step
X }
X
X }
X
X found.upper.bound_T
X
X direction_-direction
X cons.steps.same.dir_0
X step[step==1]_0
X
X
X }
X
X direction[n1==0]_0
X
X n1[direction==+1]_n1+1
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #****************************************************************************
X # Once again, define n2
X
X n2_ifelse(!opt,n1,ceiling(sqrt(c2*d2/c1/d1*(c1+d1)*(c1+d1+1)/(c2+d2)/
X (c2+d2+1))*(c1+d1+n1)-c2-d2))
X
X n2[n2<0]_0
X
X #****************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X # Return
X
X c(n1,n2)
X
X}
X
X
X###############################################################################
X#
X# A V E R A G E L E N G T H
X#
X#
X
Xalc_function(len,level,c1,d1,c2,d2,opt=F,m=1000,mcs=3)
X{
X # mcs is standing for 'maximal consecutive steps'(in same direction)
X
X min.for.possible.return_2^ceiling(1.5*mcs)
X
X # If we always allow a return, there is a risk of making bad steps
X # when we are close to the answer.
X # Thus, we should not allow any return once some arbitrary 'step' (which is
X # 'min.for.possible.return') is reached.
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #**************************************************************************
X # Define initial step
X # (as a function of any frequentist sample size estimate)
X
X step <- ceiling(log(freq(len,level,c1,d1,c2,d2))/log(2))
X
X # Also define the threshold to cross for the quantity under study (the
X # length or the coverage probability of an HPD region)
X
X threshold_len
X
X # and define a factor, which is +/- 1, depending on if the quantity under
X # study is (almost) surely too large or too small when making no
X # observations [-1 if the quantity to measure is DEcreasing with n
X # +1 if the quantity to measure is INcreasing with n]
X #
X # [ -1 if threshold_len, +1 if thresold_level ]
X
X factor_-1
X
X #**************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X quantity.to.measure_ifelse(factor == +1,0,2*threshold)
X
X step <- 2^step
X
X history.ns_0
X history.steps_0
X
X n1_0
X
X max.cons.steps.same.dir_mcs
X found.upper.bound_F
X possible.to.move.back_T
X cons.steps.same.dir_0
X direction_+1
X
X while(step>=1)
X {
X while(sign(factor*(threshold-quantity.to.measure)) == direction && step >= 1)
X {
X step[found.upper.bound]_max(1,step/2)
X possible.to.move.back[step < min.for.possible.return &&
X found.upper.bound]_F
X
X n1_n1+direction*step
X
X cons.steps.same.dir_cons.steps.same.dir+1
X
X history.ns_c(n1,history.ns)
X history.steps_c(step*direction,history.steps)
X
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #************************************************************************
X # Define n2 from n1
X
X n2_ifelse(!opt,n1,ceiling(sqrt(c2*d2/c1/d1*(c1+d1)*(c1+d1+1)/(c2+d2)/
X (c2+d2+1))*(c1+d1+n1)-c2-d2))
X
X n2[n2<0]_0
X
X #*********************************
X
X # Here compute the average length
X
X pi1_rbeta(m,c1,d1)
X pi2_rbeta(m,c2,d2)
X x1_rbinom(m,n1,pi1)
X x2_rbinom(m,n2,pi2)
X
X # Precaution: if n1 or n2 is 0, then the correponding x given
X # by rbinom is NA. We correct the situation by setting it to 0.
X
X x1[is.na(x1)]_0
X x2[is.na(x2)]_0
X
X # Posterior variance of the difference between the two proportions
X
X posterior.mean_(c2+x2)/(c2+d2+n2)-(c1+x1)/(c1+d1+n1)
X posterior.var_(x1+c1)*(n1-x1+d1)/(n1+c1+d1)^2/(n1+c1+d1+1)+
X (x2+c2)*(n2-x2+d2)/(n2+c2+d2)^2/(n2+c2+d2+1)
X
X # We make the approximation of the difference between two betas by
X # another beta, with parameters a and b given by
X
X a_-1/2 * (posterior.mean + 1) * (posterior.mean^2 - 1 + posterior.var)/
X posterior.var
X b_1/2*(1-posterior.mean)*(-1*posterior.mean^2+1-posterior.var)/
X posterior.var
X
X # To estimate the length of the HPD region of level 1-alpha, we
X # compute the length of the region bounded by lower and upper
X # quantiles of order alpha/2. In extreme cases (when both proportions
X # are much different) this approximation will be conservative, but
X # these circumstances are unlikely to be encounteded.
X
X quantity.to.measure_2*mean(qbeta((1+level)/2,a,b)-qbeta((1-level)/2,a,b))
X
X #************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X# 2
X if(found.upper.bound &&
X cons.steps.same.dir == max.cons.steps.same.dir+1 &&
X possible.to.move.back)
X {
X
X if(sign(factor*(threshold-quantity.to.measure)) == direction)
X {
X # There was (most likely) a mistake, look for n's in the other direction
X at.n_seq(along=history.ns)[history.ns == n1 ]
X hs.an_history.steps[at.n]
X
X step_abs(hs.an[sign(hs.an) != direction][1])
X
X cons.steps.same.dir_0
X
X # and if there has never been a step coming from the other direction...
X step[is.na(step)]_max(abs(hs.an))
X }
X# 3
X else
X {
X # There was (most likely) no mistake; keep looking around the same n's
X
X direction_-direction
X cons.steps.same.dir_0
X }
X }
X
X# 1
X if(found.upper.bound &&
X cons.steps.same.dir==max.cons.steps.same.dir &&
X sign(factor*(threshold-quantity.to.measure))==direction &&
X possible.to.move.back)
X {
X step_2*step
X }
X
X }
X
X found.upper.bound_T
X
X direction_-direction
X cons.steps.same.dir_0
X step[step==1]_0
X
X
X }
X
X direction[n1==0]_0
X
X n1[direction==+1]_n1+1
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #****************************************************************************
X # Once again, define n2
X
X n2_ifelse(!opt,n1,ceiling(sqrt(c2*d2/c1/d1*(c1+d1)*(c1+d1+1)/(c2+d2)/
X (c2+d2+1))*(c1+d1+n1)-c2-d2))
X
X n2[n2<0]_0
X
X #****************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X
X # Return
X
X c(n1,n2)
X}
X
X
X###############################################################################
X#
X# W O R S T O U T C O M E
X#
X#
X
X
Xwoc_function(len,level,c1,d1,c2,d2,opt=F)
X{
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #**************************************************************************
X # Define initial step
X # (as a function of any frequentist sample size estimate)
X
X step <- ceiling(log(freq(len,level,c1,d1,c2,d2))/log(2))
X
X # Also define the threshold to cross for the quantity under study (the
X # length or the coverage probability of an HPD region)
X
X threshold_level
X
X # and define a factor, which is +/- 1, depending on if the quantity under
X # study is (almost) surely too large or too small when making no
X # observations [-1 if the quantity to measure is DEcreasing with n
X # +1 if the quantity to measure is INcreasing with n]
X #
X # [ -1 if threshold_len, +1 if thresold_level ]
X
X factor_+1
X
X #**************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X quantity.to.measure_ifelse(factor == +1,0,2*threshold)
X
X step <- 2^step
X
X n1_0
X
X found.upper.bound_F
X direction_+1
X
X while(step>=1)
X {
X while(sign(factor*(threshold-quantity.to.measure)) == direction && step >= 1)
X {
X step[found.upper.bound]_max(1,step/2)
X
X n1_n1+direction*step
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #************************************************************************
X # Define n2 from n1
X
X n2_ifelse(!opt,n1,ceiling(n1+c1+d1-c2-d2))
X
X n2[n2<0]_0
X
X #*********************************
X # Compute the worst coverage
X
X # worst x's are defined by (17) in
X # Joseph, du Berger, Belisle
X
X x1_ifelse(n1 >= abs(c1-d1),ceiling((n1-c1+d1)/2),n1)
X x2_ifelse(n2 >= abs(c2-d2),ceiling((n2-c2+d2)/2),n2)
X
X posterior.mean_(c2+x2)/(c2+d2+n2)-(c1+x1)/(c1+d1+n1)
X posterior.var_(x1+c1)*(n1-x1+d1)/(n1+c1+d1)^2/(n1+c1+d1+1)+
X (x2+c2)*(n2-x2+d2)/(n2+c2+d2)^2/(n2+c2+d2+1)
X
X # We make the approximation of the difference between two betas by
X # another beta, with parameters a and b given by
X
X a_-1/2 * (posterior.mean + 1) * (posterior.mean^2 - 1 + posterior.var)/
X posterior.var
X b_1/2*(1-posterior.mean)*(-1*posterior.mean^2+1-posterior.var)/
X posterior.var
X
X pi.min_posterior.mean-len/2
X pi.max_posterior.mean+len/2
X pi.max[pi.min< -1]_-1+len ; pi.min[pi.min< -1]_-1
X pi.min[pi.max>1]_1-len ; pi.max[pi.max>1]_1
X
X quantity.to.measure_pbeta((pi.max+1)/2,a,b)-pbeta((pi.min+1)/2,a,b)
X
X #************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X }
X
X found.upper.bound_T
X
X direction_-direction
X step[step==1]_0
X
X
X }
X
X direction[n1==0]_0
X
X n1[direction==+1]_n1+1
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #****************************************************************************
X # Once again, define n2
X
X n2_ifelse(!opt,n1,ceiling(n1+c1+d1-c2-d2))
X
X n2[n2<0]_0
X
X #****************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X
X # Return
X
X c(n1,n2)
X}
X
X
X###############################################################################
X#
X# W O R S T O U T C O M E (modified) using Normal Approximation
X# for preposteriors of x1 and x2
X#
X
X
X
Xmodwoc_function(len,level,c1,d1,c2,d2,worst.level=0.95,opt=F)
X{
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #**************************************************************************
X # Define initial step
X # (as a function of any frequentist sample size estimate)
X
X step <- ceiling(log(freq(len,level,c1,d1,c2,d2))/log(2))
X
X # Also define the threshold to cross for the quantity under study (the
X # length or the coverage probability of an HPD region)
X
X threshold_level
X
X # and define a factor, which is +/- 1, depending on if the quantity under
X # study is (almost) surely too large or too small when making no
X # observations [-1 if the quantity to measure is DEcreasing with n
X # +1 if the quantity to measure is INcreasing with n]
X #
X # [ -1 if threshold_len, +1 if thresold_level ]
X
X factor_+1
X
X #**************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X quantity.to.measure_ifelse(factor == +1,0,2*threshold)
X
X step <- 2^step
X
X n1_0
X
X found.upper.bound_F
X direction_+1
X quant_qchisq(worst.level,2)
X
X
X while(step>=1)
X {
X while(sign(factor*(threshold-quantity.to.measure)) == direction && step >= 1)
X {
X step[found.upper.bound]_max(1,step/2)
X
X n1_n1+direction*step
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #************************************************************************
X # Define n2 from n1
X
X n2_ifelse(!opt,n1,ceiling(n1+c1+d1-c2-d2))
X
X n2[n2<0]_0
X
X #*****************************************
X # Compute the worst (worst-level) coverage
X
X # worst x's are defined by (17) in
X # Joseph, du Berger, Belisle
X
X x1_ifelse(n1 >= abs(c1-d1),ceiling((n1-c1+d1)/2),n1)
X x2_ifelse(n2 >= abs(c2-d2),ceiling((n2-c2+d2)/2),n2)
X
X # Preposterior moments of x1 and x2
X
X m1_n1*c1/(c1+d1)
X v1_n1*c1*d1*(n1+c1+d1)/(c1+d1)^2/(c1+d1+1)
X m2_n2*c2/(c2+d2)
X v2_n2*c2*d2*(n2+c2+d2)/(c2+d2)^2/(c2+d2+1)
X
X # (x1,x2) is the worst outcome possible;
X # see if it is in the HPD region of level
X # worst-level (approximately ellipsoidal)
X
X distance_(x1-m1)^2/v1+(x2-m2)^2/v2
X if (distance > quant)
X {
X # If the worst possible outcome is not in HPD region of indicated
X # level, then the worst outcome of that level is on the border
X # on the HPD region (on the border of the ellipse)
X
X # Normal approximation for distributions of X1 and X2
X
X x11_seq(floor(m1-sqrt(v1*quant)),ceiling(m1+sqrt(v1*quant)))
X # Eliminate impossible values and doubles
X x11[x11<0]_0
X x11[x11>n1]_n1
X x11[x11[-length(x11)]!=x11[-1]]
X # Prepare the lower part of the ellipse
X x12_x11[length(x11):1]
X
X # Higher part of the ellipse
X s1_v2*(quant-(x11-m1)^2/v1)
X s1[s1<0]_0
X x21_ceiling(m2+sqrt(s1))
X # Lower part of the ellipse
X s2_v2*(quant-(x12-m1)^2/v1)
X s2[s2<0]_0
X x22_floor(m2-sqrt(s2))
X
X x1_c(x11,x12,x11[1])
X x2_c(x21,x22,x21[1])
X x2[x2<0]_0
X x2[x2>n2]_n2
X }
X
X # Compute the posterior mean & variance
X # (for points on the border of the ellipse if worst possible outcome
X # is not in the HPD region)
X
X posterior.mean_(c2+x2)/(c2+d2+n2)-(c1+x1)/(c1+d1+n1)
X posterior.var_(x1+c1)*(n1-x1+d1)/(n1+c1+d1)^2/(n1+c1+d1+1)+
X (x2+c2)*(n2-x2+d2)/(n2+c2+d2)^2/(n2+c2+d2+1)
X
X # Compute the maximal posterior mean & variance
X # (necessary if considering points on the border of
X # the ellipse)
X
X posterior.mean_posterior.mean[posterior.var==max(posterior.var)][1]
X posterior.var_max(posterior.var)[1]
X
X pi.min_posterior.mean-len/2
X pi.max_posterior.mean+len/2
X pi.max[pi.min< -1]_-1+len ; pi.min[pi.min< -1]_-1
X pi.min[pi.max>1]_1-len ; pi.max[pi.max>1]_1
X
X # We make the approximation of the difference between two betas by
X # another beta, with parameters a and b given by
X
X a_-1/2 * (posterior.mean + 1) * (posterior.mean^2 - 1 + posterior.var)/
X posterior.var
X b_1/2*(1-posterior.mean)*(-1*posterior.mean^2+1-posterior.var)/
X posterior.var
X
X quantity.to.measure_pbeta((1+pi.max)/2,a,b)-pbeta((1+pi.min)/2,a,b)
X
X #************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X }
X
X found.upper.bound_T
X
X direction_-direction
X step[step==1]_0
X
X }
X
X direction[n1==0]_0
X
X n1[direction==+1]_n1+1
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #****************************************************************************
X # Once again, define n2
X
X n2_ifelse(!opt,n1,ceiling(n1+c1+d1-c2-d2))
X
X n2[n2<0]_0
X
X #****************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X
X # Return
X
X c(n1,n2)
X
X}
X
X
X
X
X#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X
X# M I X E D B A Y E S I A N / L I K E L I H O O D
X
X
X#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X
X
X
X###############################################################################
X#
X# M B L - A V E R A G E C O V E R A G E
X#
X#
X
X
Xmblacc_function(len,level,c1,d1,c2,d2,m=1000,mcs=3)
X{
X # mcs is standing for 'maximal consecutive steps'(in same direction)
X
X min.for.possible.return_2^ceiling(1.5*mcs)
X
X # If we always allow a return, there is a risk of making bad steps
X # when we are close to the answer.
X # Thus, we should not allow any return once some arbitrary 'step' (which is
X # 'min.for.possible.return') is reached.
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #**************************************************************************
X # Define initial step
X # (as a function of any frequentist sample size estimate)
X
X step <- ceiling(log(freq(len,level,c1,d1,c2,d2))/log(2))
X
X # Also define the threshold to cross for the quantity under study (the
X # length or the coverage probability of an HPD region)
X
X threshold_level
X
X # and define a factor, which is +/- 1, depending on if the quantity under
X # study is (almost) surely too large or too small when making no
X # observations [-1 if the quantity to measure is DEcreasing with n
X # +1 if the quantity to measure is INcreasing with n]
X #
X # [ -1 if threshold_len, +1 if thresold_level ]
X
X factor_+1
X
X #**************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X quantity.to.measure_ifelse(factor == +1,0,2*threshold)
X
X step <- 2^step
X
X history.ns_0
X history.steps_0
X
X n1_0
X
X max.cons.steps.same.dir_mcs
X found.upper.bound_F
X possible.to.move.back_T
X cons.steps.same.dir_0
X direction_+1
X
X while(step>=1)
X {
X while(sign(factor*(threshold-quantity.to.measure)) == direction && step >= 1)
X {
X step[found.upper.bound]_max(1,step/2)
X possible.to.move.back[step < min.for.possible.return &&
X found.upper.bound]_F
X
X n1_n1+direction*step
X
X cons.steps.same.dir_cons.steps.same.dir+1
X
X history.ns_c(n1,history.ns)
X history.steps_c(step*direction,history.steps)
X
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #************************************************************************
X # Define n2 from n1
X
X n2_n1
X
X #*********************************
X
X # Compute the average coverage
X
X pi1_rbeta(m,c1,d1)
X pi2_rbeta(m,c2,d2)
X x1_rbinom(m,n1,pi1)
X x2_rbinom(m,n2,pi2)
X
X # Precaution: if n1 or n2 is 0, then the correponding x given
X # by rbinom is NA. We correct the situation by setting it to 0.
X
X x1[is.na(x1)]_0
X x2[is.na(x2)]_0
X
X # Posterior variance of the difference between the two proportions
X # (same formulas as for 'acc', replacing c's and d's by 1's.
X
X posterior.mean_(1+x2)/(2+n2)-(1+x1)/(2+n1)
X posterior.var_(x1+1)*(n1-x1+1)/(n1+2)^2/(n1+3)+
X (x2+1)*(n2-x2+1)/(n2+2)^2/(n2+3)
X
X # We make the approximation of the difference between two betas by
X # another beta, with parameters a and b given by
X
X a_-1/2 * (posterior.mean + 1) * (posterior.mean^2 - 1 + posterior.var)/
X posterior.var
X b_1/2*(1-posterior.mean)*(-1*posterior.mean^2+1-posterior.var)/
X posterior.var
X
X # and approximate the probability coverage of the HPD region by the
X # coverage obtained by a symetric region around the posterior mean
X
X pi.min_posterior.mean-len/2
X pi.max_posterior.mean+len/2
X pi.max[pi.min< -1]_-1+len ; pi.min[pi.min< -1]_-1
X pi.min[pi.max>1]_1-len ; pi.max[pi.max>1]_1
X coverage.around.mean_pbeta((pi.max+1)/2,a,b)-pbeta((pi.min+1)/2,a,b)
X
X quantity.to.measure_mean(coverage.around.mean)
X
X #************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X# 2
X if(found.upper.bound &&
X cons.steps.same.dir == max.cons.steps.same.dir+1 &&
X possible.to.move.back)
X {
X
X if(sign(factor*(threshold-quantity.to.measure)) == direction)
X {
X # There was (most likely) a mistake, look for n's in the other direction
X at.n_seq(along=history.ns)[history.ns == n1 ]
X hs.an_history.steps[at.n]
X
X step_abs(hs.an[sign(hs.an) != direction][1])
X
X cons.steps.same.dir_0
X
X # and if there has never been a step coming from the other direction...
X step[is.na(step)]_max(abs(hs.an))
X }
X# 3
X else
X {
X # There was (most likely) no mistake; keep looking around the same n's
X
X direction_-direction
X cons.steps.same.dir_0
X }
X }
X
X# 1
X if(found.upper.bound &&
X cons.steps.same.dir==max.cons.steps.same.dir &&
X sign(factor*(threshold-quantity.to.measure))==direction &&
X possible.to.move.back)
X {
X step_2*step
X }
X
X }
X
X found.upper.bound_T
X
X direction_-direction
X cons.steps.same.dir_0
X step[step==1]_0
X
X
X }
X
X direction[n1==0]_0
X
X n1[direction==+1]_n1+1
X
X # Return
X
X c(n1,n1)
X
X}
X
X
X
X###############################################################################
X#
X# M B L - A V E R A G E L E N G T H
X#
X#
X
X
Xmblalc_function(len,level,c1,d1,c2,d2,m=1000,mcs=3)
X{
X # mcs is standing for 'maximal consecutive steps'(in same direction)
X
X min.for.possible.return_2^ceiling(1.5*mcs)
X
X # If we always allow a return, there is a risk of making bad steps
X # when we are close to the answer.
X # Thus, we should not allow any return once some arbitrary 'step' (which is
X # 'min.for.possible.return') is reached.
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #**************************************************************************
X # Define initial step
X # (as a function of any frequentist sample size estimate)
X
X
X step <- ceiling(log(freq(len,level,c1,d1,c2,d2))/log(2))
X
X # Also define the threshold to cross for the quantity under study (the
X # length or the coverage probability of an HPD region)
X
X threshold_len
X
X # and define a factor, which is +/- 1, depending on if the quantity under
X # study is (almost) surely too large or too small when making no
X # observations [-1 if the quantity to measure is DEcreasing with n
X # +1 if the quantity to measure is INcreasing with n]
X #
X # [ -1 if threshold_len, +1 if thresold_level ]
X
X factor_-1
X
X #**************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X quantity.to.measure_ifelse(factor == +1,0,2*threshold)
X
X step <- 2^step
X
X history.ns_0
X history.steps_0
X
X n1_0
X
X max.cons.steps.same.dir_mcs
X found.upper.bound_F
X possible.to.move.back_T
X cons.steps.same.dir_0
X direction_+1
X
X while(step>=1)
X {
X while(sign(factor*(threshold-quantity.to.measure)) == direction && step >= 1)
X {
X step[found.upper.bound]_max(1,step/2)
X possible.to.move.back[step < min.for.possible.return &&
X found.upper.bound]_F
X
X n1_n1+direction*step
X
X cons.steps.same.dir_cons.steps.same.dir+1
X
X history.ns_c(n1,history.ns)
X history.steps_c(step*direction,history.steps)
X
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #************************************************************************
X # Define n2 from n1
X
X n2_n1
X
X #*********************************
X # Compute the average length
X
X pi1_rbeta(m,c1,d1)
X pi2_rbeta(m,c2,d2)
X x1_rbinom(m,n1,pi1)
X x2_rbinom(m,n2,pi2)
X
X # Precaution: if n1 or n2 is 0, then the correponding x given
X # by rbinom is NA. We correct the situation by setting it to 0.
X
X x1[is.na(x1)]_0
X x2[is.na(x2)]_0
X
X # Posterior variance of the difference between the two proportions
X # (same formulas as for 'alc', replacing c's and d's by 1's.
X
X posterior.mean_(1+x2)/(2+n2)-(1+x1)/(2+n1)
X posterior.var_(x1+1)*(n1-x1+1)/(n1+2)^2/(n1+3)+
X (x2+1)*(n2-x2+1)/(n2+2)^2/(n2+3)
X
X # We make the approximation of the difference between two betas by
X # another beta, with parameters a and b given by
X
X a_-1/2 * (posterior.mean + 1) * (posterior.mean^2 - 1 + posterior.var)/
X posterior.var
X b_1/2*(1-posterior.mean)*(-1*posterior.mean^2+1-posterior.var)/
X posterior.var
X
X # To estimate the length of the HPD region of level 1-alpha, we
X # compute the length of the region bounded by lower and upper
X # quantiles of order alpha/2. In extreme cases (when both proportions
X # are much different) this approximation will be conservative, but
X # these circumstances are unlikely to be encounteded.
X
X quantity.to.measure_2*mean(qbeta((1+level)/2,a,b)-qbeta((1-level)/2,a,b))
X
X #************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X# 2
X if(found.upper.bound &&
X cons.steps.same.dir == max.cons.steps.same.dir+1 &&
X possible.to.move.back)
X {
X
X if(sign(factor*(threshold-quantity.to.measure)) == direction)
X {
X # There was (most likely) a mistake, look for n's in the other direction
X at.n_seq(along=history.ns)[history.ns == n1 ]
X hs.an_history.steps[at.n]
X
X step_abs(hs.an[sign(hs.an) != direction][1])
X
X cons.steps.same.dir_0
X
X # and if there has never been a step coming from the other direction...
X step[is.na(step)]_max(abs(hs.an))
X }
X# 3
X else
X {
X # There was (most likely) no mistake; keep looking around the same n's
X
X direction_-direction
X cons.steps.same.dir_0
X }
X }
X
X# 1
X if(found.upper.bound &&
X cons.steps.same.dir==max.cons.steps.same.dir &&
X sign(factor*(threshold-quantity.to.measure))==direction &&
X possible.to.move.back)
X {
X step_2*step
X }
X
X }
X
X found.upper.bound_T
X
X direction_-direction
X cons.steps.same.dir_0
X step[step==1]_0
X
X }
X
X direction[n1==0]_0
X
X n1[direction==+1]_n1+1
X
X # Return
X
X c(n1,n1)
X
X}
X
X
X
X###############################################################################
X#
X# M B L - W O R S T O U T C O M E
X#
X#
X
Xmblwoc_function(len,level,c1,d1,c2,d2)
X{
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #**************************************************************************
X # Define initial step
X # (as a function of any frequentist sample size estimate)
X
X
X step <- ceiling(log(freq(len,level,c1,d1,c2,d2))/log(2))
X
X # Also define the threshold to cross for the quantity under study (the
X # length or the coverage probability of an HPD region)
X
X threshold_level
X
X # and define a factor, which is +/- 1, depending on if the quantity under
X # study is (almost) surely too large or too small when making no
X # observations [-1 if the quantity to measure is DEcreasing with n
X # +1 if the quantity to measure is INcreasing with n]
X #
X # [ -1 if threshold_len, +1 if thresold_level ]
X
X factor_+1
X
X #**************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X quantity.to.measure_ifelse(factor == +1,0,2*threshold)
X
X step <- 2^step
X
X n1_0
X
X found.upper.bound_F
X direction_+1
X
X while(step>=1)
X {
X while(sign(factor*(threshold-quantity.to.measure)) == direction && step >= 1)
X {
X step[found.upper.bound]_max(1,step/2)
X
X n1_n1+direction*step
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #************************************************************************
X # Define n2 from n1
X
X n2_n1
X
X #*********************************
X
X # Here compute the worst coverage
X
X # worst x's are defined by (17) in
X # Joseph, du Berger, Belisle, for general c's and d's, here set to 1.
X
X x1_ceiling(n1/2)
X x2_ceiling(n2/2)
X
X posterior.mean_(1+x2)/(2+n2)-(1+x1)/(2+n1)
X posterior.var_(x1+1)*(n1-x1+1)/(n1+2)^2/(n1+3)+
X (x2+1)*(n2-x2+1)/(n2+2)^2/(n2+3)
X
X # We make the approximation of the difference between two betas by
X # another beta, with parameters a and b given by
X
X a_-1/2 * (posterior.mean + 1) * (posterior.mean^2 - 1 + posterior.var)/
X posterior.var
X b_1/2*(1-posterior.mean)*(-1*posterior.mean^2+1-posterior.var)/
X posterior.var
X
X # and approximate the probability coverage of the HPD region by the
X # coverage obtained by a symetric region around the posterior mean
X
X pi.min_posterior.mean-len/2
X pi.max_posterior.mean+len/2
X pi.max[pi.min< -1]_-1+len ; pi.min[pi.min< -1]_-1
X pi.min[pi.max>1]_1-len ; pi.max[pi.max>1]_1
X coverage.around.mean_pbeta((pi.max+1)/2,a,b)-pbeta((pi.min+1)/2,a,b)
X
X quantity.to.measure_mean(coverage.around.mean)
X
X #************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X }
X
X found.upper.bound_T
X
X direction_-direction
X step[step==1]_0
X }
X
X direction[n1==0]_0
X
X n1[direction==+1]_n1+1
X
X # Return
X
X c(n1,n1)
X
X}
X
X
X
X###############################################################################
X#
X# M B L - W O R S T O U T C O M E (modified)
X#
X#
X
X
Xmblmodwoc_function(len,level,c1,d1,c2,d2,worst.level=.95)
X{
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #**************************************************************************
X # Define initial step
X # (as a function of any frequentist sample size estimate)
X
X step <- ceiling(log(freq(len,level,c1,d1,c2,d2))/log(2))
X
X # Also define the threshold to cross for the quantity under study (the
X # length or the coverage probability of an HPD region)
X
X threshold_level
X
X # and define a factor, which is +/- 1, depending on if the quantity under
X # study is (almost) surely too large or too small when making no
X # observations [-1 if the quantity to measure is DEcreasing with n
X # +1 if the quantity to measure is INcreasing with n]
X #
X # [ -1 if threshold_len, +1 if thresold_level ]
X
X factor_+1
X
X #**************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X quantity.to.measure_ifelse(factor == +1,0,2*threshold)
X
X step <- 2^step
X
X n1_0
X
X found.upper.bound_F
X direction_+1
X quant_qchisq(worst.level,2)
X
X while(step>=1)
X {
X while(sign(factor*(threshold-quantity.to.measure)) == direction && step >= 1)
X {
X step[found.upper.bound]_max(1,step/2)
X
X n1_n1+direction*step
X
X #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
X #************************************************************************
X # Define n2 from n1
X
X n2_n1
X
X #*****************************************
X # Compute the worst (worst-level) coverage
X
X # worst x's are defined by (17) in
X # Joseph, du Berger, Belisle
X
X x1_ifelse(n1 >= abs(c1-d1),ceiling((n1-c1+d1)/2),n1)
X x2_ifelse(n2 >= abs(c2-d2),ceiling((n2-c2+d2)/2),n2)
X
X # Preposterior moments of x1 and x2
X
X m1_n1*c1/(c1+d1)
X v1_n1*c1*d1*(n1+c1+d1)/(c1+d1)^2/(c1+d1+1)
X m2_n2*c2/(c2+d2)
X v2_n2*c2*d2*(n2+c2+d2)/(c2+d2)^2/(c2+d2+1)
X
X # (x1,x2) is the worst outcome possible;
X # see if it is in the HPD region of level
X # worst-level (approximately ellipsoidal)
X
X distance_(x1-m1)^2/v1+(x2-m2)^2/v2
X if (distance > quant)
X {
X # If the worst possible outcome is not in HPD region of indicated
X # level, then the worst outcome of that level is on the border
X # on the HPD region (on the border of the ellipse)
X
X # Normal approximation for distributions of X1 and X2
X
X x11_seq(floor(m1-sqrt(v1*quant)),ceiling(m1+sqrt(v1*quant)))
X # Eliminate impossible values and doubles
X x11[x11<0]_0
X x11[x11>n1]_n1
X x11[x11[-length(x11)]!=x11[-1]]
X # Prepare the lower part of the ellipse
X x12_x11[length(x11):1]
X
X # Higher part of the ellipse
X s1_v2*(quant-(x11-m1)^2/v1)
X s1[s1<0]_0
X x21_ceiling(m2+sqrt(s1))
X # Lower part of the ellipse
X s2_v2*(quant-(x12-m1)^2/v1)
X s2[s2<0]_0
X x22_floor(m2-sqrt(s2))
X
X x1_c(x11,x12,x11[1])
X x2_c(x21,x22,x21[1])
X x2[x2<0]_0
X x2[x2>n2]_n2
X }
X
X # Compute the posterior mean & variance
X # (for points on the border of the ellipse if worst possible outcome
X # is not in the HPD region)
X # setting c's and d's to 1's.
X
X posterior.mean_(1+x2)/(2+n2)-(1+x1)/(2+n1)
X posterior.var_(x1+1)*(n1-x1+1)/(n1+2)^2/(n1+3)+
X (x2+1)*(n2-x2+1)/(n2+2)^2/(n2+3)
X
X # Compute the maximal posterior mean & variance
X # (necessary if considering points on the border of
X # the ellipse)
X
X posterior.mean_posterior.mean[posterior.var==max(posterior.var)][1]
X posterior.var_max(posterior.var)[1]
X
X pi.min_posterior.mean-len/2
X pi.max_posterior.mean+len/2
X pi.max[pi.min< -1]_-1+len ; pi.min[pi.min< -1]_-1
X pi.min[pi.max>1]_1-len ; pi.max[pi.max>1]_1
X
X # We make the approximation of the difference between two betas by
X # another beta, with parameters a and b given by
X
X a_-1/2 * (posterior.mean + 1) * (posterior.mean^2 - 1 + posterior.var)/
X posterior.var
X b_1/2*(1-posterior.mean)*(-1*posterior.mean^2+1-posterior.var)/
X posterior.var
X
X quantity.to.measure_pbeta((1+pi.max)/2,a,b)-pbeta((1+pi.min)/2,a,b)
X
X #************************************************************************
X #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
X
X }
X
X found.upper.bound_T
X
X direction_-direction
X step[step==1]_0
X
X }
X
X direction[n1==0]_0
X
X n1[direction==+1]_n1+1
X
X # Return
X
X c(n1,n1)
X
X}
X
END_OF_FILE
if test 44593 -ne `wc -c <'prop.ss'`; then
echo shar: \"'prop.ss'\" unpacked with wrong size!
fi
# end of 'prop.ss'
fi
echo shar: End of shell archive.
exit 0