/****************************************************************************
*                                                                           *
*   STD.SPL     Copyright 2007, 2023 (C) DSP Development Corporation        *
*                                                                           *
*   Author:     Randy Race                                                  *
*                                                                           *
*   Synopsis:   Column based standard deviation                             *
*                                                                           *
*   Revisions:  13 Aug 2007     RRR     Creation                            *
*                8 Feb 2023     RRR     optional args                       *
*                                                                           *
****************************************************************************/

#if @HELP_STD

    STD

    Purpose: Calculates the sample or population standard deviation of
             a series.

    Syntax:  STD(s, dof, dim, "naflag")

             (sd, mn) = STD(s, dof, dim, "naflag")

                     s - A series.

                  dof - Optional. An integer, the normalization mode.
                         For N = length(s)

                          0: sample standard deviation, normalize
                             by 1/(N-1) (default)

                          1: population standard deviation, normalize
                             by 1/N

                   dim - Optional. An integer or string, the computation
                         dimension.

                              0 : calculate over full array

                              1 : calculate column-wise (default)

                              2 : calculate row-wise

                          "all" : calculate over full array


              "naflag" - Optional. A string, the NA handling method.

                             "omitnan" : ignore NA values (default)

                          "includenan" : include NA values

    Returns: A real scalar for a one column series or computation over the
             entire array else a 1xM real table for an M column series.

             (sd, mn) = STD(s, dof, dim, "naflag") returns the standard
             deviation and mean value as separate variables.

    Example:
             W1: 1..10
             std(w1)

             Returns 3.0277, the sample standard deviation.

    Example:
             W1: 1..10
             std(w1+1e7)

             Returns 3.0277, same as above since the true standard deviation
             is independent of the mean value.

    Example:
             W1: 1..10
             std(w1, 1)

             Returns 2.8723, the population standard deviation.

    Example:
             W1: ravel(1..9, 3)^2
             W2: std(w1)

             W2 == {{4.0415, 10.0167, 16.0104}}, the standard deviation of
             each column.

    Example:
             W1: ravel(1..9, 3)^2
             W2: std(w1, 0, 2)

             W2 == {24.5561, 30.4467, 36.3731}, the standard deviation of
             each row.

    Example:
             W1: ravel(1..9, 3)^2
             W2: {std(w1, 0, "all")}

             W2 == {28.0802}, the standard deviation of the entire array.

    Example:
             W1: ravel(1..9, 3)^2
             (sd, mn) = std(w1, 0, "all");

             sd == 28.0802
             mn == 31.6667

             Returns the standard deviation and mean of the entire array.

    Example:
             W1: {1, 10, nan, 3, 4, nan, 6}
             W2: {std(w1, 0)}
             W3: {std(w1, 0, "omitnan")}
             W4: {std(w1, 0, "includenan")}

             W2 == W3 == {3.4205}, the standard deviation of the series with
             NA values ignored.

             W4 == nan, since NA values are included.

    Remarks:
             If s is the input series and N the length of the series,
             the basic form of the sample standard deviation (dof == 0)
             is:

                std(s) = sqrt(sum((s - mean(s))^2/(N-1)))

             The population standard deviation (dof == 1) is defined as:

                std(s) = sqrt(sum((s - mean(s))^2/N))


             STD uses a fast, highly accurate corrected two-pass Neumaier
             sum algorithm that exhibits insensitivity to round-off errors.

             As shown in the second example, the standard deviation is
             independent of an additive mean value.

             STD returns the standard deviation of each column or row for
             a multi-column series. Use DIM = 0 or DIM = "all" or use
             STDEV to compute the standard deviation of the entire array.

    See Also:
             Mean
             Stdev
#endif


/* standard deviation of each column */
std(s, dof, dim, naflag, method)
{
        local N, sd, mn;

        if (argc < 1 || not(isarray(s)))
        {
                /* default to current series */
                (dof, dim, naflag, method) = stdvar_parse_args(FALSE, s, dof, dim, naflag);

                s = refseries(w0);
        }
        else
        {
                (dof, dim, naflag, method) = stdvar_parse_args(TRUE, dof, dim, naflag, method);
        }

        if (dim > 1)
        {
                s = s';
        }

        if (dim == 0)
        {
                /* full array */
                (sd, mn) = stdev(s, -1, -1, dof, naflag, method);
        }
        else
        {
                /* stdev of each column or row */
                (sd, mn) = colstdev(s, -1, -1, dof, naflag, method);

                if (numcols(sd) == 1)
                {
                        /* return scalars */
                        sd = sd[1];
                        mn = mn[1];
                }
                else if (dim > 1)
                {
                        sd = sd';
                        mn = mn';
                }
        }
        
        return(sd, mn);
}