======  ESTIMATION OF WEIBULL PARAMETERS (BASIC, C, PASCAL) ======

NOTES:

1.  Two-parameter Weibull, or three-parameter with location parameter
    (origin) assumed known.

2.  The C and Pascal programs use the parametrization F(x) = 1 -
    exp[-a(x-c)^b], instead of the alternative F(x) = 1 -
    exp{-[(x-c)/a]^b}.  See the BASIC version for the relationship
    between the two.

3.  Method of moments, i.e., fitted distribution agrees with the data
    in mean and variance.  We normally fit the distribution to the
    tree basal areas or diameters squared, instead of to the
    diameters (note that with two parameters both are Weibulls). 
    Then the mean is the mean basal area, or the usual (quadratic)
    mean dbh.

4.  The implementation avoids iterations by using an approximation in
    Garcia, O.  "Simplified method-of-moments estimation for the
    Weibull distribution", New Zealand Journal of Forestry Science
    11:304-306, 1981.  This gives an accuracy of at least 5 or 6
    significant figures (if I remember correctly, I do not have the
    paper at hand right now).  Limited to coefficients of variation
    less than 120%, which covers the range where the Weibull is unimodal.


===========================  MS BASIC  ===============================

(Note lines wrapped-around, followed by blank lines)

100 PRINT: PRINT TAB(15);"MOMENT ESTIMATES FOR WEIBULL PARAMETERS"
110 PRINT: PRINT TAB(13);"(Garcia,O., N.Z.J.For.Sci.11,304-306,1981.)"
120 PRINT: PRINT "If the location parameter, c, is not zero, the mean and coef
ficient"

130 PRINT "of variation must be for x-c."
140 PRINT: INPUT "Mean, coefficient of variation (%)";M,Z: Z=Z/100
150 IF Z>1.2 THEN PRINT "OUT OF RANGE": GOTO 140
160 Z=(((((7.454537E-03*Z*Z-8.354348E-02)*Z+.153109251#)*Z-1.946641E-03)*Z-.22
000991#)*(1-Z)*(1-Z)+1)*Z

170 B=1/Z
180 G=1
190 IF Z > 1 THEN G=Z: Z=Z-1
200 G=((((((((.035868343#*Z-.193527818#)*Z+.482199394#)*Z-.756704078#)*Z+.9182
06857#)*Z-.897056937#)*Z+.988205891#)*Z-.577191652#)*Z+1)*G

210 PRINT "     F(x) = 1 - exp[-";(G/M)^B;"(x-c)^";B;"]"
220 PRINT "or"
230 PRINT "     F(x) = 1 - exp[-{(x-c)/";M/G;"}^";B;"]."
240 GOTO 140

===============================  C  =============================

#include <stdio.h>
#include <math.h>

/*-----------------------------------------------------------------------*/
/* Moment estimates for Weibull parameters.
    Garcia, O., N.Z.J.For.Sci. 11,304-306,1981.
*/
#define OUTofRANGE 0
/* returns OUTofRANGE if arguments are out of range, !OUTofRANGE if OK */

int Weibull(double mean, double stdev, double orign,        /* input */
                 double *a, double *b, double *c)           /* output */
{
    double x;

    *c = orign;
    x = stdev / (mean - *c);
    if (x > 1.2 || x <= 0)
        return OUTofRANGE;
    x = x*(((((0.007454537*x*x-0.08354348)*x+0.153109251)*x-0.001946641)
        *x-0.22000991)*(1-x)*(1-x)+1);
    *b = 1 / x;
    if (x <= 1)
        *a = 1;
    else
    {   *a = x;
        x -= 1;
    }
    *a = ((((((((0.035868343*x-0.193527818)*x+0.482199394)*x-0.756704078)
            *x+0.918206857)*x-0.897056937)*x+0.988205891)*x-0.577191652)
            *x+1)* *a;
    *a = exp(*b * log(*a / (mean - *c)));
    return !OUTofRANGE;
}

/*-----------------------------------------------------------------------*/
/* Test, example */

void main()
{
    double m, s, o, a, b, c;
    int n;

    do
    {   printf("\nEnter the mean, standard deviation, and origin\n\
(separated by blanks, exit with 0 0 0): ");
        if ((n = scanf("%lf %lf %lf", &m, &s, &o)) != 3)
        {   printf("\nInvalid input (%d values read).\n", n);
            scanf("%*s");   /* read and ignore entry */
        }
        else if (m > 0)
        {
            printf("Mean: %g, std.dev.: %g, origin: %g\n", m, s, o);
            if (Weibull (m, s, o, &a, &b, &c))
            {   printf("OK:\n");
                printf(
  "\tF(x) = 1 - exp[-%g (x-%g)^%g]\nor\n\tF(x) = 1 - exp[-{(x-%g)/%g}^%g]\n",
                    a, c, b, c, pow(a, -1/b), b);
            }
            else
                printf("*** WEIBULL: Parameters out of range ***\n");
        }
    } while (m > 0);
}

=========================  Pascal  ================================

PROGRAM TryWeibull (INPUT,OUTPUT);

VAR m, s, o, a, b, c : REAL;

{----------------------------------------------------------------------------}
PROCEDURE Weibull (mean, stdev, orign : REAL; VAR a, b, c : REAL);

    {Moment estimates for Weibull parameters.
     Garcia,O., N.Z.J.For.Sci.11,304-306,1981.}
    
    VAR x : REAL;
    
    BEGIN
    c := orign;
    x := stdev / (mean - c);
    IF (x > 1.2) OR (x <= 0.0) THEN
	Writeln ('*** WEIBULL: Parameters out of range')
    ELSE
	BEGIN
	x := x*(((((0.007454537*x*x-0.08354348)*x+0.153109251)*x-0.001946641)
	     *x-0.22000991)*(1.0-x)*(1.0-x)+1.0);
	b := 1.0 / x;
	IF x <= 1.0 THEN
	    a := 1.0
	ELSE
	    BEGIN
	    a := x;
	    x := x - 1.0
	    END;
	a := ((((((((0.035868343*x-0.193527818)*x+0.482199394)*x-0.756704078)
	     *x+0.918206857)*x-0.897056937)*x+0.988205891)*x-0.577191652)
	     *x+1.0)*a;
	a := Exp(b * Ln(a / (mean - c)))
	END
    END; {Weibull}
{----------------------------------------------------------------------------}

BEGIN {Main}
REPEAT
    Write ('Mean, standard deviation, origin? ');
    Readln (m, s, o);
    IF m > 0 THEN
	BEGIN
	Weibull (m, s, o, a, b, c);
	Writeln ('a = ', a:11, ',  b = ', b:11, ',  c = ', c:11)
	END
UNTIL m < 0
    
END.

======================================================================