/*****************************************************************************
*                                                                            *
*   POWFIT.SPL    Copyright (C) 2000 DSP Development Corporation             *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:       Randy Race                                                 *
*                                                                            *
*   Synopsis:     Fits y(x) = A * x^B                                        *
*                                                                            *
*   Revisions:    20 Jun 2000  RRR  Creation - from series.mac               *
*                                                                            *
*****************************************************************************/

#include 

#if @HELP_POWFIT

    POWFIT

    Purpose: Fits y(x) = A * x^B using linearization

    Syntax:  POWFIT(s)

             s   - input series or array


    Returns: A series, the fitted power curve

             (fit, coef) = powfit(s)

             returns both the fit and the coefficients as a series

             (fit, A, B) = powfit(s)

             returns the fit as a series and the coefficients as
             separate scalars


    Example:
             W1: 10 * (1..100)^0.5

             W2: powfit(w1);overp(w1, lred)

             Overplots the original data with the calculated power fit.


    Example:
             (fit, coef) = powfit(w1)

             fit is the same series as in W2

             coef == {10.0, 0.5}

    Remarks:
             POWFIT fits a power curve of the form y = Ax^b. The fit is
             accomplished by fitting a line to the following equation:

                            ln(y) = ln(A) + b*ln(x)

             Note that both x and y must be positive.

    See Also:
             Expfit
             Polyfit
             Trend
#endif


/* fit y(x) = A * x ^ B */
ITERATE powfit(s)
{
        local coef, result;

        if (argc < 1) error("powfit - input series required");

        /* transform and fit */
        coef    = polyfit(ln(s), ln(xvals(s)), 1, -1);
        coef[1] = exp(coef[1]);

        result = coef[1] * xvals(s) ^ coef[2];

        if (isxy(s))
        {
                result = xy(xvals(s), result);
        }

        if (outargc > 1)
        {
                if (outargc > 2)
                {
                        return(result, coef[1], coef[2]);
                }

                return(result, coef);
        }

        return(result);
}