Topic: Utilities

Topic Description:

Basic R utilities used in many packages and functions, such as date and string manipulations, statistical distributions, geometry of lines and distances. The R package date is required for the two data functions.



File: utilities/statistics.r

View File Source Download File No help file available

Function: metrop1step

Function Description: metrop1step

Takes a single metropolis step on a single parameter for any given likelihood function. The arguments start.param and scale.param are atomic (single values), as are adjust and target. The ellipses handle all other arguments to the function. The function func must accept the test parameter as the first argument, plus any additional arguments which come in the ellipses. Note the metropolis rule: if rejected, the old value is returned to be re-used. The return value includes a one if accepted, zero if rejected. The step size, refered to as scale.param, is adjusted following Helene's rule. For every acceptance, scale.param is multiplied by adjust, which is a small number > 1 (1.01, 1.02, 1.1 all seem to work). For every rejection, scale.param is multiplied by (1/adjust); for every acceptance, by adjust^AdjExp, the latter based on the target acceptance rate. When the target acceptance rate is 0.25, which is recommended for any model with > 4 parameters, AdjExp=3. It's easy to see how this system arrives at an equilibrium acceptance rate=target. The program calling metrop1step has to keep track of the scaling parameter: submitting it each time metrop1step is called, and saving the adjusted value for the next call. Given many parameters, a scale must be stored separately for every one. Note the return value is a vector of 6: 1) the new parameter value; 2) the new scale (step size); 3) a zero or a one to keep track of the acceptance rate; 4) the likelihood of original parameter (if rejected) or new parameter (if accepted) 5) the likelihood of original parameter (if accepted) or new parameter (if rejected) 6) the new parameter tested (whether accepted or not)

Function Arguments:

ArgumentDefault Value
func
start.param
scale.param
adjust
target
...

Function Source:

metrop1step=function(func,start.param,scale.param,adjust,target,...)
{
origlike=func(start.param,...)
newval=rnorm(1,mean=start.param,sd=scale.param)
newlike=func(newval,...)

if(is.na(newlike)) browser()

AdjExp=(1-target)/target

## It's possible with shifting parameters to get trapped where start.param is invalid.
# This test allows an escape
if(origlike==(-Inf))
{ likeratio=IfElse(newlike==(-Inf),0,1); browser }
else likeratio=exp(newlike-origlike)

if(is.na(likeratio)) browser()

if(runif(1) {
newscale=scale.param*adjust^AdjExp
return(c(newval,newscale,1,newlike,origlike,newval))
}
else ## Reject
{
newscale=scale.param*(1/adjust)
return(c(start.param,newscale,0,origlike,newlike,newval))
}

}