# Robust mean

## Synopsis

    use oakleaf

    subroutine rmean(data,mean,stderr,stdsig,scale,reltol,robfun,rankinit,flag,verbose,report)

       real, dimension(:), intent(in) :: data
       real, intent(in out) :: mean
       real, intent(out), optional :: stderr, stdsig
       real, intent(in out), optional :: scale
       real, intent(in) :: reltol
       character(len=*), intent(in) :: robfun
       logical, intent(in), optional :: rankinit
       integer, intent(out), optional :: flag
       logical, intent(in), optional :: verbose
       character(len=*), intent(in), optional :: report

    end subroutine rmean

The subroutine is generic, supporting REAL32, REAL64 and REAL128 (available on 64-bit platforms) datatypes. The reccomended, if the nature of a given problem allows it, is REAL64. REAL32 has low accuracy, whilst REAL128 is slow.

## Description

The robust mean subroutine `rmean()` estimates the mean, and related statistical quantities. The data are supposed to be sample from the Normal distribution with possible contamination of a unknown origin.

The parameters for both the location and the dispersion of the Normal distribution N(mean,stdsig) are estimated. The true value X of the sample falls, with 68% probability, into the interval mean - stderr < X < mean + stderr.

## Parameters

### On input:

data
  array of data,
mean
  (optional) an initial estimate of the mean,
scale
  (optional) an initial estimate of the scale (i.e. stdsig),
reltol
  (optional) the relative accuracy of determination of the optimum (default: 2.4%),
robfun
  (optional) select the robust function among: ``hampel'` (default), ``tukey', `huber'` or ``square'` (non-robust least square),
rankinit
  (optional) estimate initial values by the order estimates (default: `.true.`)
verbose
  (optional) print a detailed information about the calulations.
report
  (optional) a filename of a report

`data` are the mandatory input parameter. The procedure initialisation is done internally if `rankinit == .true.` (default); otherwise, both `mean,scale` must be set by the calling procedure. The initial estimates are taken by the median and MAD described in qmean.

### On output:

mean
  the mean
stderr
  (optional) the standard error
stdsig
  (optional) the standard deviation
scale
  (optional) the scale
flag
  (optional) a flag, see the status codes.

The mandatory output parameter is the robust mean. The check of the status is recommended.

## Example

Save the program to the file `example.f08`:

    program example

       use oakleaf
       use iso_fortran_env

       real(REAL64), dimension(5) :: data = [ 1, 2, 3, 4, 5 ]
       real(REAL64) :: mean,stderr,stdsig

       call rmean(data,mean,stderr,stdsig)
       write(*,*) 'rmean:',mean,' stderr:',stderr,'stdsig:',stdsig

    end program example

Than compile and run it (the paths of `-I, -L` are an example):

    $ gfortran -I/usr/local/include example.f08 -L/usr/local/lib -loakleaf -lminpack
    $ ./a.out
    rmean:   3.0001622773505234       stderr:  0.68885212488300640      stdsig:   1.5403201776835767

## Reports

The report file keeps the detailed track of the algorithm for visualisation or debugging purposes. It is not recommended on a common use as it slow-down a run.

The format is (see `src/sqpfmean.f08`) self-descibing and includes: initial estimates, the optimisation steps (current values, trust-region radius, the step, the gradient, the quasi-Newton hessian and the trust-step flag) and final number of iterations, function evaluations, restores, steps out of the trust-region, indefinite and hard steps and the final hessian computed analytically.

## References

Hroch,F.: A contribution to the estimate of the robust mean value, in preparation