/*______________________________________________
|                                               |
|                   A1inv                       |
|                   _____                       |
|                                               |
| A1inv macro computes the inverse of the ratio |
| of the first and zeroth order Bessel function |
| of the first kind. It is used to compute the  |
| maximum likelihood estimate of kappa, the     |
| concentration parameter of the von Mises      |
| circular distribution. It returns the value   |
| of a SAS variable.                            |
|_______________________________________________|
|           |           |                       |
| Parameter |  Default  |  Description          |
|___________|___________|_______________________|
|           |           |                       |
| x         |           | numeric SAS variable  |
|           |           | with values between 0 |
|           |           | 1
|           |           |                       |
| result    |           | name of output SAS    |
|           |           | variable              |
|___________|___________|_______________________|
|                                               |
| Written by: Kristin Bishop (04/08/2004)       |
| Modified by: Ulric Lund (08/16/2006)          |
|_____________________________________________*/;

%macro A1inv(x, result);

  if (0 <= &x and &x < 0.53)
    then &result = 2*(&x) + (&x)**3 + 5*((&x)**5)/6;

  else if (0.53 <= &x and &x < 0.85)
    then &result = -0.4 + 1.39*&x + 0.43/(1-&x);

  else if &x >= 0.85
    then &result = 1/(&x**3 - 4*(&x**2) + 3*&x);

%mend A1inv;


/*______________________________________________
|                                               |
|                 circANOVA                     |
|                   _____                       |
|                                               |
| circANOVA runs a one way ANOVA given an       |
| angular response variable and a grouping      |
| variable. It returns an output data set with  |
| the following variables:                      |
|                                               |
| group - the name of each group, numeric or    |
|         character; see below: group is a macro|
|         parameter                             |
| n - number of observations in each group      |
| meanDir - group mean direction of angles      |
| kappaHat - group estimated concentration      |
|            parameter                          |
| r - group resultant length                    |
| F - approximate F-statistic                   |
| testStat - test statistic computed from F;    |
|            fits F-distribution better than F  |
| pvalue - p-value of the test                  |
|_______________________________________________|
|                                               |
| Dependencies: uses macros circSummary         |
|                           estKappa            |
|                           circMean            |
|                           circDisp            |
|                           A1inv               |
|_______________________________________________|
|                                               |
|  Details: The circular ANOVA implemented here |
|  is the high concentration multisample F-test |
|  for a common mean direction assuming a von   |
|  Mises distribution and equal kappas across   |
|  groups.  For details, please see K. Mardia   |
|  and P. Jupp, (2000); Directional Statistics; |
|  Section 7.4; Wiley: Chichester, England.     |
|_______________________________________________|
|           |           |                       |
| Parameter |  Default  |  Description          |
|___________|___________|_______________________|
|           |           |                       |
| dataIn    |           | input SAS data set    |
|           |           |                       |
| dataOut   |           | output SAS data set   |
|           |           |                       |
| y         |           | angular response      |
|           |           | variable in radians   |
|           |           |                       |
| group     |           | variable containing   |
|           |           | group designation     |
|           |           |                       |
| print     |    T      | OPTIONAL: print out   |
|           |           | ANOVA table to output |
|           |           | window. Set to F to   |
|           |           | suppress printing.    |
|___________|___________|_______________________|
|                                               |
| Written by: Blake Sweeney (04/28/2006)        |
| Modified by: Ulric Lund (09/13/2006)          |
|_____________________________________________*/;

%macro circANOVA(dataIn, dataOut, y, group, print=T);

  proc means data=&dataIn noprint nonobs;
    var &y;
    output out=_size (drop=_type_ _freq_) n=totN;
  run;

  data _null_;
    set _size;
    if _n_ = 1 then call symput('n',compress(totN));
  run;

  proc sort data=&dataIn out=_groups nodupkey;
   by &group;
  run;

  proc contents data=_groups out=_groups (keep=nobs);
  run;

  data _null_;
    set _groups;
    if _n_ = 1 then call symput('q',compress(nobs));
  run;

  proc transpose data=&dataIn out=_temp (drop=_name_ &group);
    var &y;
    by &group;
  run;

  proc transpose data=_temp out=_grouped (drop=_name_);
    var _numeric_;
  run;

  %do i=1 %to &q;

    %circSummary(_grouped, _var&i, col&i, print=F)
    %estKappa(_grouped, _kappa&i, col&i, bias=F, print=F)

    data _var&i;
      set _var&i(rename=(meanDir=mu&i n=n&i r=r&i));
    run;

    data _kappa&i;
      set _kappa&i(rename=(kappaHat=kappa&i));
    run;

  %end;

  %let var = _var1;
  %let kappa = _kappa1;

  %do i=1 %to %eval(&q-1);
    %let t = %eval(&i+1);
    %let var = &var _var&t;
    %let kappa = &kappa _kappa&t;
  %end;

  data _kappas;
    merge &kappa;
  run;

  data _stats (keep=n1-n&q mu1-mu&q r1-r&q);
    merge &var;
  run;

  proc transpose data=_stats out=_ns (keep=col1);
    var n1-n&q;
  run;

  proc transpose data=_stats out=_mus (keep=col1);
    var mu1-mu&q;
  run;

  proc transpose data=_stats out=_rs (keep=col1);
    var r1-r&q;
  run;

  proc transpose data=_kappas out=_kappas (keep=col1);
    var kappa1 - kappa&q;
  run;

  %circSummary(&dataIn, _summaryAll, &y, print=F)
  %estKappa(&dataIn, _kappaAll, &y, bias=F, print=F)
     
  data _all;
    merge _summaryAll _kappaAll;
    keep n meanDir kappaHat r;
  run;
 
  data _preCalc;
    merge _ns     (rename=(col1=n_i))
          _mus    (rename=(col1=meanDir_i))
          _kappas (rename=(col1=kappaHat_i))
          _rs     (rename=(col1=r_i));
  run;
  
  data _preCalc;
    if _n_=1 then set _all;
    set _preCalc;
  run;

  proc means data=_preCalc noprint;
    var r_i;
    output out=_sumRs (drop=_type_ _freq_) sum=sumRs;
  run;

  data _calc;
    merge _all (keep=n r kappaHat) _sumRs;
    F = ((sumRs - R)/(&q - 1)) / ((&n - sumRs)/(&n - &q));
    testStat = (1 + 3/(8*kappaHat))*F;
    pvalue = 1 - probf(testStat, &q-1, &n-&q);
  run;
 
  data _anova;
    set _calc (drop=F kappaHat testStat n);
    Source = "Between";
    DF = &q - 1;
    SS = sumRs - R;
    MS = (sumRs - R)/(&q - 1);
    F  = ((sumRs - R)/(&q - 1)) / ((&n - sumRs)/(&n - &q));
    p  = pvalue;
    output;
    
    Source = "Within";
    DF = &n - &q;
    SS = &n - sumRs;
    MS = (&n - sumRs)/(&n - &q);
    F  = .;
    p = .;
    output;

    Source = "Total";
    DF = &n - 1;
    SS = &n - R;
    MS = .;
    F  = .;
    p = .;
    output;
    drop sumRs R pvalue;
  run;

  %if &print=T %then %do;
    proc print data=_anova noobs;
    run;
    title;
  %end;

  proc freq data=&dataIn;
    tables &group / noprint out=_groupNames (drop=count percent);
  run;

  data &dataOut;
    merge _groupNames 
          _ns     (rename=(col1=n)) 
          _mus    (rename=(col1=meanDir))
          _kappas (rename=(col1=kappaHat))
          _rs     (rename=(col1=r)) 
          _calc   (keep=F testStat pvalue);
  run; 

  proc datasets nodetails nolist;
    delete _size _groups _temp _grouped _var1 - _var&q _kappa1 - _kappa&q
           _kappas _stats _ns _mus _rs _sumRs _summaryAll _kappaAll _all _preCalc 
           _groupNames _calc _anova;
  run;quit;

%mend;


/*______________________________________________
|                                               |
|                 circCor                       |
|                   _____                       |
|                                               |
| circCor macro computes a correlation          |
| coefficient of two circular variables.        |
| Optionally, a test of significance can be     |
| requested. It returns an output data set with |
| the following variables:                      |
|                                               |
| n - number of observations in input data set  |
| meanDir1 - mean direction for first variable  |
| meanDir2 - mean direction for second variable |
| r - correlation coefficient                   |
|                                               |
| If a significance test is requested, the      |
| following variables are also included:        |
|                                               |
| testStat - value of the test statistic        |
| pvalue - p-value of the test                  |
|_______________________________________________|
|                                               |
| Dependencies: uses macro circMean             |
|_______________________________________________|
|
|  Details: the correlation coefficient computed|
|  here is due to Jammalamadaka, S. and         |
|  Sarma, Y. (1988); A correlation coefficient  |
|  for angular variables; Statistical Theory and|
|  Data Analysis 2; North Holland: New York.    |
|  See also Jammalamadaka, S. Rao and SenGupta, |
|  A. (2001); Topics in Circular Statistics,    |
|  Section 8.2; World Scientific Press:         |
|  Singapore.                                   |
|_______________________________________________|
|           |           |                       |
| Parameter |  Default  |  Description          |
|___________|___________|_______________________|
|           |           |                       |
| dataIn    |           | input SAS data set    |
|           |           |                       |
| dataOut   |           | output SAS data set   |
|           |           |                       |
| var1      |           | angular variable in   |
|           |           | radians               |
|           |           |                       |
| var2      |           | angular variable in   |
|           |           | radians               |
|           |           |                       |
| test      |    F      | OPTIONAL: set test=T  |
|           |           | to obtain a test of   |
|           |           | significance.         |
|           |           | Default is to not do  |
|           |           | the test.             |
|           |           |                       |
| print     |    T      | OPTIONAL: print out   |
|           |           | results to output     |
|           |           | window. Set to F to   |
|           |           | suppress printing.    |
|           |           |                       |
|___________|___________|_______________________|
|                                               |
| Written by: Kristin Bishop (04/16/2004)       |
| Modified by: Ulric Lund (08/17/2006)          |
|_____________________________________________*/;

%macro circCor(dataIn, dataOut, var1, var2, test=F, print=T);

  %circMean(&dataIn, _meanOut1, &var1, print=F);
  %circMean(&dataIn, _meanOut2, &var2, print=F);

  data _calc;
    if _n_=1 then do;
      set _meanOut1 (rename=(meanDir=meanDir1));
      set _meanOut2 (rename=(meanDir=meanDir2));
    end;
    set &dataIn;

    sinVar1   = sin(&var1-meanDir1);
    sinVar2   = sin(&var2-meanDir2);
    sinVar1Sq = (sin(&var1-meanDir1))**2;
    sinVar2Sq = (sin(&var2-meanDir2))**2;

    Num   + (sinVar1   * sinVar2);
    DenSum1 + sinVar1Sq;
    DenSum2 + sinVar2Sq;
    Denom = sqrt(DenSum1*DenSum2);
    r=Num/Denom;

    %if &test=T %then %do;
      l20 + (sinVar1Sq/n);
      l02 + (sinVar2Sq/n);
      l22 + (sinVar1Sq*sinVar2Sq/n);

      testStat = sqrt((n*l20*l02)/l22)*r;
      pValue = 2 * (1-cdf('normal',(abs(testStat))));
    %end;
  run;

  data &dataout;
    set _calc end=last;
    if last;
    %if &test=F %then %do;
      keep n meanDir1 meanDir2 r;
    %end;
    %if &test=T %then %do;
      keep n meanDir1 meanDir2 r testStat pValue;
    %end;
  run;
  
  proc datasets nodetails nolist;
    delete _meanOut1 _meanOut2 _calc;
  run;quit;
  
  %if &print=T %then %do;
    proc print data=&dataOut noobs;
    run;
  %end;

%mend circCor;


/*______________________________________________
|                                               |
|                 circDisp                      |
|                   _____                       |
|                                               |
| circDisp macro computes measures of dispersion|
| for a circular variable. It returns an output |
| data set with the following variables:        |
|                                               |
| n - number of observations in input data set  |
| r - resultant length                          |
| rbar - mean resultant length                  |
| circVar - circular variance                   |
| circSD - circular standard deviation          |
|_______________________________________________|
|           |           |                       |
| Parameter |  Default  |  Description          |
|___________|___________|_______________________|
|           |           |                       |
| dataIn    |           | input SAS data set    |
|           |           |                       |
| dataOut   |           | output SAS data set   |
|           |           |                       |
| var       |           | angular variable in   |
|           |           | radians               |
|           |           |                       |
| print     |   T       | OPTIONAL: print out   |
|           |           | results to output     |
|           |           | window. Set to F to   |
|           |           | suppress printing.    |
|___________|___________|_______________________|
|                                               |
| Written by: Kristin Bishop (02/09/2004)       |
| Modified by: Ulric Lund (08/13/2006)          |
|_____________________________________________*/;

%macro circDisp(dataIn, dataOut, var, print=T);

  data _use;
    set &dataIn;
    if &var ^= .;
  run;

  data _calc;
    set _use;                             
    n=_n_;
    cosVar=cos(&var);                          
    sinVar=sin(&var);                           

    sumCos + cosVar;                            
    sumSin + sinVar;                            

    r=sqrt((sumCos**2)+(sumSin**2));
    rbar=r/n;

    circVar=1-rbar;
    circSD = sqrt(-2*log(rbar));
  run;

  data &dataOut;
    set _calc end=last;
    if last;
    keep n r rbar circVar circSD;
  run;

  %if &print=T %then %do;
    proc print noobs;
    run;
  %end;

  proc datasets nodetails nolist;
    delete _use _calc;
  run;quit;

%mend circDisp;


/*______________________________________________
|                                               |
|                 circMean                      |
|                   _____                       |
|                                               |
| circMean macro computes the mean direction    |
| of a circular variable. It returns an output  |
| data set with the following variables:        |
|                                               |
| n - number of observations in input data set  |
| aveCos - average of cosines of angles         |
| aveSin - average of sines of angles           |
| meanDir - mean direction of angles            |
|_______________________________________________|
|           |           |                       |
| Parameter |  Default  |  Description          |
|___________|___________|_______________________|
|           |           |                       |
| dataIn    |           | input SAS data set    |
|           |           |                       |
| dataOut   |           | output SAS data set   |
|           |           |                       |
| var       |           | angular variable in   |
|           |           | radians               |
|           |           |                       |
| print     |   T       | OPTIONAL: print out   |
|           |           | results to output     |
|           |           | window. Set to F to   |
|           |           | suppress printing.    |
|___________|___________|_______________________|
|                                               |
| Written by: Kristin Bishop (01/27/2004)       |
| Modified by: Ulric Lund (08/13/2006)          |
|_____________________________________________*/;

%macro circMean(dataIn, dataOut, var, print=T);

  data _use;
    set &dataIn;
    if &var ^= .;
  run;

  data _calc;
    set _use;                             
    n=_n_;
    cosVar=cos(&var);                          
    sinVar=sin(&var);                           

    sumCos + cosVar;                            
    sumSin + sinVar;                            

    aveCos=sumCos/n;
    aveSin=sumSin/n;

    meanDir=atan2(aveSin, aveCos);
  run;

  data &dataOut;
    set _calc end=last;
    if last;
    keep n aveCos aveSin meanDir;
  run;

  %if &print=T %then %do;
    proc print noobs;
    run;
  %end;

  proc datasets nodetails nolist;
    delete _use _calc;
  run;quit;

%mend circMean;


/*__________________________________________________
|                                                   |                                                  
|                     circPlot                      |
|                     ________                      |
|                                                   |
| circPlot creates a circular plot of data measured |
| in radians. If the stack option is selected can   |
| return an output data set with the following      |
| variables:                                        |
|                                                   |
| bin - bin number                                  |
| low - lower bound of bin                          |
| high - upper bound of bin                         |
| count - number of observations in each bin        |
| prop - proportion of observations in each bin     |
|                                                   |
| Graph can be optionally saved to an external file.|
|___________________________________________________|
|           |               |                       |
| Parameter | Default       | Description           |
|___________|_______________|_______________________|
|           |               |                       |
| dataIn    |               | the input data set    |
|           |               | measured in radians   |
|           |               |                       |
| var       |               | angular variable in   |
|           |               | in radians            |
|           |               |                       |
| stack     |      F        | OPTIONAL: if F macro  |
|           |               | plots each point on   |
|           |               | the circumfrenc of    |
|           |               | the circle.  If T,    |
|           |               | macro creates bins and|
|           |               | stacks points at the  |
|           |               | midpoint of each bin. |
|           |               |                       |
| dataOut   |     NA        | OPTIONAL: name of     |
|           |               | output data set with  |
|           |               | variables noted above |
|           |               |                       |
| bins      |     36        | OPTIONAL: numeric     |
|           |               | value above zero      |
|           |               | spec number of bins   |
|           |               |                       |
| radius    |      1        | OPTIONAL: positive    |
|           |               | number that increases |
|           |               | or decreases the      |
|           |               | radius of the circle. |
|           |               | Values below 1 shrink,|
|           |               | and values above 1    |
|           |               | expand the circle.
|           |               |                       |
| ptsize    |      1        | OPTIONAL: numeric     |
|           |               | value above zero that |
|           |               | increases/decreases   |
|           |               | size of points        |
|           |               |                       |
| ptspace   |      0.1      | OPTIONAL: numeric     |
|           |               | value above zero that |
|           |               | increases/decreases   |
|           |               | space between points  |
|           |               |                       |
| txtsize   |      1        | OPTIONAL: numeric     |
|           |               | value above zero that |
|           |               | increases/decreases   |
|           |               | the text size         |
|           |               |                       |
| txtspace  |      0.05     | OPTIONAL: numeric     |
|           |               | value above zero that |
|           |               | increases/decreases   |
|           |               | space between text    |
|           |               | and plotted circle    |
|           |               |                       |
| axisco    |     white     | OPTIONAL: color of the|
|           |               | axis                  |
|           |               |                       |
| axis      |      1        | OPTIONAL: positive    |
|           |               | value indicating the  |
|           |               | +/- endpoints of axes |
|           |               |                       |
| hshift    |      0        | OPTIONAL: shifts the  |
|           |               | the center of the     |   
|           |               | circle along the      |
|           |               | x-axis                |
|           |               |                       |
| vshift    |      0        | OPTIONAL: shifts the  |
|           |               | the center of the     |   
|           |               | circle along the      |
|           |               | y-axis                |
|           |               |                       |
| length    |     4in       | OPTIONAL: horizontal &|
|           |               | vertical length of the|
|           |               | graphics output area  |
|___________|_______________|_______________________|
|                                                   |
| Written by: Michael Cummings (April 2004)         |
| Modified by: Ulric Lund (12/30/2006)              |
|_________________________________________________*/;

%macro circPlot(dataIn, var, stack=F, dataOut=NA, bins=36, radius=1, 
                ptsize=1, ptspace=0.1, txtsize=1, txtspace=0.05, axisco=white,
                axis=1, hshift=0, vshift=0, length=4in);
  
  data _use;  /* Take variable modulo 2 pi */
    set &dataIn;
    &var = mod(&var+2*constant('pi'), 2*constant('pi'));
  run;

  %if &stack=F %then %do;
    data _points;
      set _use;
      x = &radius * cos(&var) + &hshift;
      y = &radius * sin(&var) + &vshift;
    run;
  %end;

  %if &stack=T %then %do;
    data _binAssign (keep=&var bin);
      set _use;
      angle = 2*constant('pi')/&bins;        
      cut1 = 0;                              
      cut2 = angle;                           
      bin = 1; 
      flag = 0; 

      do while(flag=0);  /* Find bin number for each point */
        if &var >= cut1 and &var < cut2 then do;
          flag = 1;
          bin = bin;                      
        end;
        else do;
          cut1 = cut2; 
          cut2 = cut2 + angle;
          bin = bin + 1;
        end;
      end; 
    run;

    proc freq data=_binAssign noprint;
      tables bin / out=_binFreq;
    run;

    data _points;
      set _binFreq;
      n = 1;
      r = &radius;

      do while(n <= count);  /* Get x and y coordinate of stacked points for each bin */
        x = r * cos((constant('pi')/&bins)*(2*bin - 1)) + &hshift;
        y = r * sin((constant('pi')/&bins)*(2*bin - 1)) + &vshift;
        output;
        n = n + 1;            **increment n until = freq for specific bin;
        r = r + &ptspace;     **increment radius;
      end;
    run;

  
    data _dummy;  /* Create dummy data set to fill in values for bins with no observations */
      angle = 2*constant('pi')/&bins;         
 
      do i=1 to &bins;
        bin = i;
        low = (i-1)*angle;
        high = i*angle;
        output;
      end;
    run;

    data _tempOut;  /* Create output data set */
      merge _binfreq _dummy;  
      by bin;
      if count = . then count = 0;      
      if percent = . then percent = 0;
      drop i angle;
    run;

    %if &dataOut ne NA %then %do;
      data &dataOut;
        set _tempOut;
      run;
    %end;

  %end;

  data _text;  /* Data set with text to annotate the graph */
    length text $ 3;
    xsys='2';
    ysys='2';
    color = 'black';
    size=&txtsize;

    x = &radius - &txtspace + &hshift;
    y = 0 + &vshift;
    text = '0';
    position = '4';
    output;

    x = -1*&radius + &txtspace + &hshift;
    y = 0 + &vshift;
    text = '180';
    position = '6';
    output;

    x = 0 + &hshift;
    y = &radius - &txtspace + &vshift;
    text = '90';
    position = '8';
    output;

    x = 0 + &hshift;
    y = -1*&radius + &txtspace + &vshift;
    text = '270';
    position = '2';
    output;

    x = 0 + &hshift;
    y = 0 + &vshift;
    text = '+';
    position = '5';
    output;
  run;

  data _circle;  /* Create data set to annotate graph with a circle */
    function = 'pie';
    line=0;
    xsys = '2';
    ysys = '2';
    hsys = '2';
    x = 0 + &hshift;
    y = 0 + &vshift;
    angle = 0;
    rotate = 360;
    size = &radius;
  run;

  data _annoData;
    set _text _circle;
  run;

  /* Set options for the axis and symbol */
  axis1 length=&length  order=(-&axis &axis) color=&axisco;
  axis2 length=&length  order=(-&axis &axis) color=&axisco;
  symbol1 color=black value=dot height=&ptsize;

  proc gplot data=_points annotate=_annoData;
    plot y*x / haxis=axis1 vaxis=axis2;
  run;quit;

  /* Restore default options for axes and symbol */
  axis1;
  axis2;
  symbol1;

  %if &stack=F %then %do;
    proc datasets;
      delete _use _points _text _circle _annoData;
    run;quit;
  %end;

  %if &stack=T %then %do;
    proc datasets;
      delete _use _points _text _circle _annoData _binassign _binfreq _dummy _tempOut;
    run;quit;
  %end;

%mend circPlot;


/*______________________________________________
|                                               |
|                 circRange                     |
|                   _____                       |
|                                               |
| circRange macro computes the range of a       |
| circular variable. It returns an output data  |
| set with the following variables:             |
|                                               |
| n - number of observations in input data set  |
| range - circular range of angles; smallest    |
|         arc encompassing all angles.          |
|_______________________________________________|
|           |           |                       |
| Parameter |  Default  |  Description          |
|___________|___________|_______________________|
|           |           |                       |
| dataIn    |           | input SAS data set    |
|           |           |                       |
| dataOut   |           | output SAS data set   |
|           |           |                       |
| var       |           | angular variable in   |
|           |           | radians               |
|           |           |                       |
| print     |   T       | OPTIONAL: print out   |
|           |           | results to output     |
|           |           | window. Set to F to   |
|           |           | suppress printing.    |
|___________|___________|_______________________|
|                                               |
| Written by: Kristin Bishop (01/27/2004)       |
| Modified by: Ulric Lund (08/19/2006)          |
|_____________________________________________*/;

%macro circRange(dataIn, dataOut, var, print=T);

  data _use;
    set &dataIn;
    if &var ^= .;
  run;

  data _mod;
    set _use;
    varMod=mod(&var, 2*constant('PI'));
  run;

  proc sort data=_mod;
    by varMod;
  run;

  proc means data=_mod noprint;
    var varMod;
    output out=_minMax (drop=_type_ _freq_) min=low max=high;
  run;

  data _calc;
    if _n_=1 then do;
      set _minMax;
      space = low - high + (2*constant('PI'));
      output;
    end;
    set _mod;
    space = dif(varMod);
    output;
  run;

  proc means data=_calc noprint;
    var space;
    output out=_range (drop=_type_ _freq_) n=n max=maxSpace;
  run;

  data &dataout; 
    set _range;
    range = 2 * constant('PI') - maxSpace;
    keep n range;
  run;

  %if &print=T %then %do;
    proc print noobs;
    run;
  %end;

  proc datasets nodetails nolist;
    delete _use _mod _minMax _range _calc;
  run;quit;

%mend circRange;


/*______________________________________________
|                                               |
|                 circReg                       |
|               ___________                     |
|                                               |
| circReg fits a simple regression predicting   |
| one angular response variable from another    |
| angular predictor variable. Two output data   |
| sets can be generated, with the following     |
| variables:                                    |
|                                               |
| dataOut1:                                     |
| x - the predictor variable in radians         |
| y - the response variable in radians          |
| y_hat - the predicted values for y in radians |
|                                               |
| datatOut2:
| type - the type of statistic computed         |
| response - the response variable (sine or     |
|            cosine of y                        |
| intercept - estimate, t-value, or p-value of  |
|             the intercept in the estimated    |
|             regression models                 |
| cosx1 - cosx# - estimate, t-value, or p-value |
| sinx1 - sinx#   of the coefficients in the    |
|                 estimated regression models   |
| F - partial F-tests for the cosine and sine   |
|     estimated regression models, testing      |
|     whether the current models are sufficient |
|     relative to the next higher order models  |
| pvalues - p-values of the F-tests             |
|_______________________________________________|
|                                               |
| Details: The simple circular regression model |
| implemented here predicts one angular         |
| variable from one other.  This is done by     |
| fitting models to both the sine and cosine    |
| of the response variable using trigonometric  |
| polynomials of the predictor variable. The    |
| user can either specify the order of the      |
| trigonometric polynomials to use, or have the |
| macro add terms until partial F-tests         |
| determine that additional terms are not       |
| necessary in both sine and cosine models.     |
| For more details see Jammalamadaka, S. Rao and|
| SenGupta, A. (2001); Topics in Circular       |
| Statistics, Section 8.6; World Scientific     |
| Press: Singapore.                             |
|_______________________________________________|
|           |           |                       |
| Parameter |  Default  |  Description          |
|___________|___________|_______________________|
|           |           |                       |
| dataIn    |           | input SAS data set    |
|           |           |                       |
| dataOut1  |           | output SAS data set   |
|           |           | to contain the        |
|           |           | observed x & y values |
|           |           | along with the        |
|           |           | predicted y values    |
|           |           |                       |
| dataOut2  |           | output SAS data set   |
|           |           | to contain the table  |
|           |           | of estimated          |
|           |           | regression coeffs,    |
|           |           | F-statistics, and     |
|           |           | p-values              |
|           |           |                       |
| x         |           | angular predictor     |
|           |           | variable in radians   |
|           |           |                       |
| y         |           | angular response      |
|           |           | variable in radians   |
|           |           |                       |
| m         |    0      | OPTIONAL: order of    |
|           |           | trig polynomial to be |
|           |           | fit. Default value of |
|           |           | 0, allows macro to    |
|           |           | increase levels of    |
|           |           | models until terms    |
|           |           | are no longer needed. |
|           |           |                       |
| alpha     |    0.05   | OPTIONAL: level of    |
|           |           | signifigance to use   |
|           |           | for partial F-tests.  |
|           |           |                       |
| print     |    T      | OPTIONAL: print out   |
|           |           | regression output.    |
|           |           | Set to F to suppress  |
|           |           | printing.             |
|___________|___________|_______________________|
|                                               |
| Written by: Blake Sweeney (04/17/2006)        |
| Modified by: Ulric Lund (12/28/2006)          |
|_____________________________________________*/;


%macro circReg(dataIn, dataOut1, dataOut2, x, y, m=0, alpha=0.05, print=T);

  %if &print=F %then %do; ods listing close; %end;
  %if &m ne 0 %then %do; /* Begin: order specified by user */

    data _reg;
      set &dataIn;
      cos_y = cos(&y);
      sin_y = sin(&y);

      %do i=1 %to %eval(&m + 1);
        cosx&i = cos(&i*&x);
        sinx&i = sin(&i*&x);
      %end;
    run;

    /*******************************************************************/
    /* Extracting the ANOVA tables for trigonometric regression models */
    /*******************************************************************/
    
    proc reg data=_reg alpha=&alpha outest=_cosCoef&m tableout; 
      model cos_y = cosx1 - cosx&m sinx1 - sinx&m;
      title "Circular Regression w/ Order=&m, Alpha=&alpha";
      output out=_fitg2 p=g2;
      ods output anova=_cosAnova&m;
    run;quit;

    proc reg data=_reg alpha=&alpha outest=_sinCoef&m tableout;
      model sin_y = cosx1 - cosx&m sinx1 - sinx&m;
      title "Circular Regression w/ Order=&m, Alpha=&alpha";
      output out=_fitg1 p=g1;
      ods output anova=_sinAnova&m;
    run;quit;


    %let n = %eval(&m + 1);  /* Get next higher order model to test for */
                             /* significance thereof.                   */

    proc reg data=_reg alpha=&alpha outest=_cosCoef&n tableout;
      model cos_y = cosx1 - cosx&n sinx1 - sinx&n;
      title "Circular Regression w/ Order=&n, Alpha=&alpha";
      ods output anova=_cosAnova&n;
    run;quit;

    proc reg data=_reg alpha=&alpha outest=_sinCoef&n tableout; 
      model sin_y = cosx1 - cosx&n sinx1 - sinx&n;
      title "Circular Regression w/ Order=&n, Alpha=&alpha";
      ods output anova=_sinAnova&n;
    run;quit;


    /*******************************************/
    /* Compute partial F-test for cosine model */
    /*******************************************/

    data _mergeCos (drop=dependent source ms fvalue probf);
      merge _cosAnova&m (rename=(ss=ssReduced df=dfReduced))
            _cosAnova&n (rename=(ss=ssFull    df=dfFull));
    run;

    data _modelCos (drop=model dfReduced ssReduced dfFull ssFull);
      set _mergeCos;
      if _n_=1 then do;
        ssrFull    = ssFull;
        ssrReduced = ssReduced;
        dfrReduced = dfReduced;
        dfrFull    = dfFull;
      end;
      else delete;
    run;

    data _errorCos (drop=model dfReduced ssReduced dfFull ssFull);
      set _mergeCos;
      if _n_=2 then do;
        sseFull = ssFull;
        sseReduced = ssReduced;
        dfeReduced = dfReduced;
        dfeFull = dfFull;
      end;
      else delete;
    run;

    data _ftestCos;   /* Calculate F-stat and p-value for the cosine model */
      merge _modelCos _errorCos;
      numSS = (sseReduced-sseFull);
      denSS = sseFull;
      numDf = dfeReduced-dfeFull;
      denDf = dfeFull;
      F = (numSS/numDf) / (denSS/denDf);
      FProb = 1-probf(F,numDf,denDf);
    run;


    /*****************************************/
    /* Compute partial F-test for sine model */
    /*****************************************/

    data _mergeSin (drop=dependent source ms fvalue probf);
      merge _sinAnova&m (rename=(ss=ssReduced df=dfReduced))
            _sinAnova&n (rename=(ss=ssFull    df=dfFull));
    run;

    data _modelSin (drop=model dfReduced ssReduced dfFull ssFull);
      set _mergeSin;
      if _n_=1 then do;
        ssrFull    = ssFull;
        ssrReduced = ssReduced;
        dfrReduced = dfReduced;
        dfrFull    = dfFull;
      end;
      else delete;
    run;

    data _errorSin (drop=model dfReduced ssReduced dfFull ssFull);
      set _mergeSin;
      if _n_=2 then do;
        sseFull = ssFull;
        sseReduced = ssReduced;
        dfeReduced = dfReduced;
        dfeFull = dfFull;
      end;
      else delete;
    run;

    data _ftestSin;   /* Calculate F-stat and p-value for the sine model */
      merge _modelSin _errorSin;
      numSS = (sseReduced-sseFull);
      denSS = sseFull;
      numDf = dfeReduced-dfeFull;
      denDf = dfeFull;
      F = (numSS/numDf) / (denSS/denDf);
      FProb = 1-probf(F,numDf,denDf);
    run;

    data _pvalues (keep=pSin pCos);
      merge _ftestCos (rename=(FProb=pCos))
            _ftestSin (rename=(FProb=pSin));
    run;


    /***************************************************************/
    /* Create output data set with fitted values and original data */
    /***************************************************************/

    data &dataOut1 (keep= &x &y y_hat);
      merge _fitg1 
            _fitg2;
       y_hat = atan2(g1,g2);
    run;


    /**************************************************************************/
    /* Create output data set with estimated coefficients and test statistics */
    /**************************************************************************/

    data _fstats;
      set _ftestCos _ftestSin;
      keep f;
    run;

    data _coeffp;
      set _sinCoef&m (where=(_type_='PVALUE'))
          _cosCoef&m (where=(_type_='PVALUE'));
      keep _type_ _depvar_ intercept cosx1-cosx&m sinx1-sinx&m;
    run;

    data _coefft;
      set _sinCoef&m (where=(_type_='T'))
          _cosCoef&m (where=(_type_='T'));
      keep _type_ _depvar_ intercept cosx1-cosx&m sinx1-sinx&m;
    run;

    data _coeff (drop= _model_ _rmse_ cos_y sin_y);
      merge _cosCoef&m _sinCoef&m;
      where _type_ = "PARMS";
      by _depvar_;
    run;

    proc transpose data=_pvalues out=_pvalues;
    run;

    data _preOut;
      merge _coeff _fstats _pvalues (rename=(col1=pvalues));
      drop _name_;
    run;

    data &dataOut2;
      set _preOut _coefft _coeffp;
      if _type_ = "PARMS" then _type_="ESTIMATE";
    run;

    proc sort data=&dataOut2 out=&dataOut2 (rename=(_type_=type _depvar_=response));
      by _depvar_;
    run;

    %if &print=T %then %do;
      title;
      proc print data=&dataOut2 noobs;
      run;
    %end;

    /* Delete data sets created */
    proc datasets;
      delete _coeff _coeffp _coefft 
             _cosAnova&m _cosAnova&n _cosCoef&m _cosCoef&n
             _sinAnova&m _sinAnova&n _sinCoef&m _sinCoef&n
             _errorCos _errorSin _fitg1 _fitg2 _fstats _ftestCos _ftestSin
             _mergeCos _mergeSin _modelCos _modelSin _preOut _pvalues _reg; 
    run;quit;

    %if &print=F %then %do; ods listing; %end;

  %end; /* End: order specified by user */


  /****************************************/
  /* Let the macro find the optimal order */
  /****************************************/

  %if &m=0 %then %do; /* Begin: order not specified by user */

    %let flag=0;

    %do %while(&flag=0);

      %let m=%eval(&m + 1);
      %let n=%eval(&m + 1);

      data _reg;
        set &dataIn;
        cos_y = cos(&y);
        sin_y = sin(&y);

        %do i=1 %to &n;
          cosx&i = cos(&i*&x);
          sinx&i = sin(&i*&x);
        %end;
      run;


      /*******************************************************************/
      /* Extracting the ANOVA tables for trigonometric regression models */
      /*******************************************************************/
    
      proc reg data=_reg alpha=&alpha outest=_cosCoef&m tableout; 
        model cos_y = cosx1 - cosx&m sinx1 - sinx&m;
        title "Circular Regression w/ Order=&m, Alpha=&alpha";
        output out=_fitg2 p=g2;
        ods output anova=_cosAnova&m;
      run;quit;

      proc reg data=_reg alpha=&alpha outest=_sinCoef&m tableout;
        model sin_y = cosx1 - cosx&m sinx1 - sinx&m;
        title "Circular Regression w/ Order=&m, Alpha=&alpha";
        output out=_fitg1 p=g1;
        ods output anova=_sinAnova&m;
      run;quit;

      proc reg data=_reg alpha=&alpha outest=_cosCoef&n tableout;
        model cos_y = cosx1 - cosx&n sinx1 - sinx&n;
        title "Circular Regression w/ Order=&n, Alpha=&alpha";
        ods output anova=_cosAnova&n;
      run;quit;

      proc reg data=_reg alpha=&alpha outest=_sinCoef&n tableout; 
        model sin_y = cosx1 - cosx&n sinx1 - sinx&n;
        title "Circular Regression w/ Order=&n, Alpha=&alpha";
        ods output anova=_sinAnova&n;
      run;quit;


      /*******************************************/
      /* Compute partial F-test for cosine model */
      /*******************************************/

      data _mergeCos (drop=dependent source ms fvalue probf);
        merge _cosAnova&m (rename=(ss=ssReduced df=dfReduced))
              _cosAnova&n (rename=(ss=ssFull    df=dfFull));
      run;

      data _modelCos (drop=model dfReduced ssReduced dfFull ssFull);
        set _mergeCos;
        if _n_=1 then do;
          ssrFull    = ssFull;
          ssrReduced = ssReduced;
          dfrReduced = dfReduced;
          dfrFull    = dfFull;
        end;
        else delete;
      run;

      data _errorCos (drop=model dfReduced ssReduced dfFull ssFull);
        set _mergeCos;
        if _n_=2 then do;
          sseFull = ssFull;
          sseReduced = ssReduced;
          dfeReduced = dfReduced;
          dfeFull = dfFull;
        end;
        else delete;
      run;

      data _ftestCos;   /* Calculate F-stat and p-value for the cosine model */
        merge _modelCos _errorCos;
        numSS = (sseReduced-sseFull);
        denSS = sseFull;
        numDf = dfeReduced-dfeFull;
        denDf = dfeFull;
        F = (numSS/numDf) / (denSS/denDf);
        FProb = 1-probf(F,numDf,denDf);
      run;


      /*****************************************/
      /* Compute partial F-test for sine model */
      /*****************************************/

      data _mergeSin (drop=dependent source ms fvalue probf);
        merge _sinAnova&m (rename=(ss=ssReduced df=dfReduced))
              _sinAnova&n (rename=(ss=ssFull    df=dfFull));
      run;

      data _modelSin (drop=model dfReduced ssReduced dfFull ssFull);
        set _mergeSin;
        if _n_=1 then do;
          ssrFull    = ssFull;
          ssrReduced = ssReduced;
          dfrReduced = dfReduced;
          dfrFull    = dfFull;
        end;
        else delete;
      run;

      data _errorSin (drop=model dfReduced ssReduced dfFull ssFull);
        set _mergeSin;
        if _n_=2 then do;
          sseFull = ssFull;
          sseReduced = ssReduced;
          dfeReduced = dfReduced;
          dfeFull = dfFull;
        end;
        else delete;
      run;

      data _ftestSin;   /* Calculate F-stat and p-value for the sine model */
        merge _modelSin _errorSin;
        numSS = (sseReduced-sseFull);
        denSS = sseFull;
        numDf = dfeReduced-dfeFull;
        denDf = dfeFull;
        F = (numSS/numDf) / (denSS/denDf);
        FProb = 1-probf(F,numDf,denDf);
      run;

      data _pvalues (keep=pSin pCos);
        merge _ftestCos (rename=(FProb=pCos))
              _ftestSin (rename=(FProb=pSin));
        if pSin > &alpha & pCos > &alpha then call symput("flag",'1');
      run;

    %end;   /* End while loop */


    /***************************************************************/
    /* Create output data set with fitted values and original data */
    /***************************************************************/

    data &dataOut1 (keep= &x &y y_hat);
      merge _fitg1 
            _fitg2;
       y_hat = atan2(g1,g2);
    run;


    /**************************************************************************/
    /* Create output data set with estimated coefficients and test statistics */
    /**************************************************************************/

    data _fstats;
      set _ftestCos _ftestSin;
      keep f;
    run;

    data _coeffp;
      set _sinCoef&m (where=(_type_='PVALUE'))
          _cosCoef&m (where=(_type_='PVALUE'));
      keep _type_ _depvar_ intercept cosx1-cosx&m sinx1-sinx&m;
    run;

    data _coefft;
      set _sinCoef&m (where=(_type_='T'))
          _cosCoef&m (where=(_type_='T'));
      keep _type_ _depvar_ intercept cosx1-cosx&m sinx1-sinx&m;
    run;

    data _coeff (drop= _model_ _rmse_ cos_y sin_y);
      merge _cosCoef&m _sinCoef&m;
      where _type_ = "PARMS";
      by _depvar_;
    run;

    proc transpose data=_pvalues out=_pvalues;
    run;

    data _preOut;
      merge _coeff _fstats _pvalues (rename=(col1=pvalues));
      drop _name_;
    run;

    data &dataOut2;
      set _preOut _coefft _coeffp;
      if _type_ = "PARMS" then _type_="ESTIMATE";
    run;

    proc sort data=&dataOut2 out=&dataOut2 (rename=(_type_=type _depvar_=response));
      by _depvar_;
    run;

    %if &print=T %then %do;
      title;
      proc print data=&dataOut2 noobs;
      run;
    %end;

    /* Delete data sets created */
    proc datasets;
      delete _coeff _coeffp _coefft 
             _cosAnova1 - _cosAnova&n _cosCoef1 - _cosCoef&n
             _sinAnova1 - _sinAnova&n _sinCoef1 - _sinCoef&n
             _errorCos _errorSin _fitg1 _fitg2 _fstats _ftestCos _ftestSin
             _mergeCos _mergeSin _modelCos _modelSin _preOut _pvalues _reg; 
    run;quit;

    %if &print=F %then %do; ods listing; %end;

  %end; /* End: order not specified by user */

%mend circReg;

/*______________________________________________
|                                               |
|                 circSummary                   |
|                   _____                       |
|                                               |
| circSummary macro provides a summary of       |
| information about a circular variable.        |
| It returns an output data set with the        |
| following variables:                          |
|                                               |
| n - number of observations in input data set  |
| aveCos - average of cosines of angles         |
| aveSin - average of sines of angles           |
| meanDir - mean direction of angles            |
| r - resultant length                          |
| rbar - mean resultant length                  |
| circVar - circular variance                   |
| circSD - circular standard deviation          |
|_______________________________________________|
|                                               |
| Dependencies: uses macros circMean            |
|                           circDisp            |
|_______________________________________________|
|           |           |                       |
| Parameter |  Default  |  Description          |
|___________|___________|_______________________|
|           |           |                       |
| dataIn    |           | the input data set    |
|           |           | measured in radians   |
|           |           |                       |
| dataOut   |           | the output data set   |
|           |           | returned by the macro |
|           |           |                       |
| var       |           | the variable in the   |
|           |           | input data set that   |
|           |           | you want to find      |
|           |           | the summary of        |
|           |           |                       |
| print     |   T       | OPTIONAL: print out   |
|           |           | results to output     |
|           |           | window. Set to F to   |
|           |           | suppress printing.    |
|___________|___________|_______________________|
|                                               |
| Written by: Kristin Bishop (02/19/2004)       |
| Modified by: Ulric Lund (08/17/2006)          |
|_____________________________________________*/;

%macro circSummary(dataIn, dataout, var, print=T);

  %circMean(&dataIn, _meanOut, &var, print=F);
  %circDisp(&dataIn, _dispOut, &var, print=F);

  data &dataout;
    merge _meanOut _dispOut;
  run;
  
  %if &print=T %then %do;
    proc print noobs;
    run;
  %end;

  proc datasets nodetails nolist;
    delete _meanOut _dispOut;
  run;quit;

%mend circSummary;


/*______________________________________________
|                                               |
|                    deg                        |
|                   _____                       |
|                                               |
| rad macro converts an angular variable        |
| measured in radians into degrees.             |
| It returns the value of a SAS variable.       |
|_______________________________________________|
|           |           |                       |
| Parameter |  Default  |  Description          |
|___________|___________|_______________________|
|           |           |                       |
| radVar    |           | angular input variable|
|           |           | measured in radians   |
|           |           |                       |
| degVar    |           | angular output        |
|           |           | variable measured in  |
|           |           | degrees               |
|___________|___________|_______________________|
|                                               |
| Written by: Kristin Bishop (01/27/2004)       |
|_____________________________________________*/;

%macro deg(radVar, degVar);

  &degVar = &radVar*180/constant('PI');

%mend deg;


/*__________________________________________________
|                                                   |
|                       dvm                         |
|                     ________                      |
|                                                   |
| dvm returns the probability density function of   |
| the von Mises distribution evaluated at a         |
| specified value.  It returns the value of a SAS   |
| variable.                                         |
|___________________________________________________|
|                                                   |
| Note: this macro uses a DROP statement.  It is    |
|       advisable not to use a KEEP statement in    |
|       the DATA step that uses this macro.         |
|___________________________________________________|
|           |               |                       |
| Parameter | Default       | Description           |
|___________|_______________|_______________________|
|           |               |                       |
| theta     |               | angular input variable|
|           |               | measured in radians   |
|           |               |                       |
| mu        |               | mean direction of the |
|           |               | von Mises distribution|
|           |               |                       |
| kappa     |               | concentration         |
|           |               | parameter of the von  |
|           |               | Mises distribution    |
|           |               |                       |
| density   |               | name of output SAS    |
|           |               | variable              |
|___________|_______________|_______________________|
|                                                   |
| Written by: Michael A. Cummings (04/15/2004)      |
| Modified by: Ulric Lund (08/17/2006)              |
|_________________________________________________*/;

%macro dvm(theta, mu, kappa, density);

  _bessel = ibessel(0,&kappa,0);
  &density = (1/(2*constant('PI')*_bessel)) * exp(&kappa*cos(&theta - &mu));
  drop _bessel;

%mend dvm;


/*______________________________________________
|                                               |
|                 estKappa                      |
|                   _____                       |
|                                               |
| estKappa macro computes the maximum likelihood|
| estimate of kappa, the concentration parameter|
| of the von Mises distribution, given a        |
| circular data set. It returns an output data  |
| set with the following variables:             |
|                                               |
| n - number of observations in input data set  |
| rbar - mean resultant length                  |
| kappaHat - estimate of kappa                  |
|_______________________________________________|
|                                               |
| Dependencies: uses macros circDisp            |
|                           A1inv               |
|_______________________________________________|
|           |           |                       |
| Parameter |  Default  |  Description          |
|___________|___________|_______________________|
|           |           |                       |
| dataIn    |           | input SAS data set    |
|           |           |                       |
| dataOut   |           | output SAS data set   |
|           |           |                       |
| var       |           | angular variable in   |
|           |           | radians               |
|           |           |                       |
| bias      |   F       | OPTIONAL: set bias=T  |
|           |           | if you want use the   |
|           |           | bias correction in the|
|           |           | computation of the MLE|
|           |           | Default is to NOT use |
|           |           | the bias correction.  |
|           |           |                       |
| print     |   T       | OPTIONAL: print out   |
|           |           | results to output     |
|           |           | window. Set to F to   |
|           |           | suppress printing.    |
|___________|___________|_______________________|
|                                               |
| Written by: Kristin Bishop (04/13/2004)       |
| Modified by: Ulric Lund (08/16/2006)          |
|_____________________________________________*/;

%macro estKappa(dataIn, dataOut, var, bias=F, print=T);

  %circDisp(&dataIn, _getR, &var, print=F);

  data &dataOut;
    set _getR;
    %A1inv(rbar, kappaHat);

    %if &bias=T %then %do;
      kappaML = kappaHat;
      if kappaML <  2 then kappaHat = max(kappaML-2*(n*kappaML)**(-1),0);
      if kappaML >= 2 then kappaHat = ((n-1)**3 * kappaML)/(n**3+n);
      drop kappaML;
    %end;
    keep n rbar kappaHat;
  run;

  %if &print=T %then %do;
    proc print noobs;
    run;
  %end;

  proc datasets nodetails nolist;
    delete _getR;
  run;quit;

%mend estKappa;


/*______________________________________________
|                                               |
|                    rad                        |
|                   _____                       |
|                                               |
| rad macro converts an angular variable        |
| measured in degrees into radians. It returns  |
| the value of a SAS variable.                  |
|_______________________________________________|
|           |           |                       |
| Parameter |  Default  |  Description          |
|___________|___________|_______________________|
|           |           |                       |
| degVar    |           | angular input variable|
|           |           | measured in degrees   |
|           |           |                       |
| radVar    |           | angular output        |
|           |           | variable measured in  |
|           |           | radians               |
|___________|___________|_______________________|
|                                               |
| Written by: Kristin Bishop (01/27/2004)       |
| Modified by: Ulric Lund (08/17/2006)          |
|_____________________________________________*/;

%macro rad(degVar, radVar);

  &radVar=(&degVar*constant('PI')/180);

%mend rad;


/*______________________________________________
|                                               |
|                 rayleigh                      |
|                   _____                       |
|                                               |
| rayleigh macro computes a Rayleigh test       |
| of uniformity for an angular variable. It     |
| returns an output data set with the following |
| variables:                                    |
|                                               |
| n - number of observations in input data set  |
| meanDir - mean direction of angles            |
| rbar - mean resultant length                  |
| testStat - value of the test statistic        |
| pvalue - p-value of the test                  |
|                                               |
| If a mean direction is specified the following|
| variable is also included:                    |
|                                               |
| muNot - the hypothesized mean direction       |
|_______________________________________________|
|                                               |
| Note: if a mean direction is not hypothesized |
|       then the testStat will coincide with    |
|       mean resultant length.                  |
|_______________________________________________|
|                                               |
| Dependencies: uses macro circSummary          |
|_______________________________________________|
|           |           |                       |
| Parameter |  Default  |  Description          |
|___________|___________|_______________________|
|           |           |                       |
| dataIn    |           | input SAS data set    |
|           |           |                       |
| dataOut   |           | output SAS data set   |
|           |           |                       |
| var       |           | angular variable in   |
|           |           | radians               |
|           |           |                       |
| mu        |   NA      | OPTIONAL: hypothesized|
|           |           | mean direction;       |
|           |           | the default is to     |
|           |           | perform Rayleigh's    |
|           |           | test without a mean   |
|           |           | direction hypothesized|
|           |           |                       |
| print     |   T       | OPTIONAL: print out   |
|           |           | results to output     |
|           |           | window. Set to F to   |
|           |           | suppress printing.    |
|___________|___________|_______________________|
|                                               |
| Written by: Blake Sweeney (02/02/2006)        |
| Modified by: Ulric Lund (08/19/2006)          |
|_____________________________________________*/;

%macro rayleigh(dataIn, dataOut, var, mu=NA, print=T);

  %circSummary(&dataIn, _summaryOut, &var, print=F);

  %if &mu ^= NA %then %do;
    data &dataOut;
      set _summaryOut;
      muNot = &mu;
      rNotBar = rbar*cos(meanDir-muNot);
      testStat = rNotBar;
      zNot = sqrt(2*n)*rNotBar;
      fz = exp(-.5*(zNot**2))/sqrt(2*constant('PI'));
      pz = probnorm(zNot);
      pvalue = 1-pz+fz*(((3*zNot-(zNot**3))/(16*n)) + ((15*zNot)+(305*zNot**3)-(125*zNot**5)+(9*zNot**7))/(4608*n**2));
      testStat = rNotBar;
      keep n meanDir rbar muNot testStat pvalue;
    run;
  %end;

  %else %do;
    data &dataOut;
      set _summaryOut;
      z=(n*rbar**2);
      testStat=rbar;
      pvalue = exp(-z)*(1+(((2*z)-(z**2))/(4*n)) - ((24*z) - (132*z**2) + (76*z**3) - (9*z**4))/(288*n**2));
      keep n meanDir rbar testStat pvalue;
    run;
  %end;

  %if &print=T %then %do;
    proc print noobs;
    run;
  %end;

  proc datasets nodetails nolist;
    delete _summaryOut;
  run;quit;

%mend Rayleigh;


/*__________________________________________________
|                                                   |
|                       rvm                         |
|                     ________                      |
|                                                   |
| rvm generates pseudo random numbers from the von  |
| Mises distribution.  It returns an output data set|
| with the following variable:                      |
|                                                   |
| rvm - contains random numbers generated           |
|___________________________________________________|
|           |               |                       |
| Parameter | Default       | Description           |
|___________|_______________|_______________________|
|           |               |                       |
| n         |               | number of observations|
|           |               | to be generated       |
|           |               |                       |
| mu        |               | mean direction of     |
|           |               | von Mises distribution|
|           |               |                       |
| kappa     |               | concentration         |
|           |               | parameter of von Mises|
|           |               | distribution; number  |
|           |               | greater than 0        |
|           |               |                       |
| dataOut   |               | output SAS data set   |
|           |               |                       |
| seed      |      0        | seed value; default   |
|           |               | value of zero leads   |
|           |               | to non-replicable     |
|           |               | results               |
|___________|_______________|_______________________|
|                                                   |
| Written by: Michael A. Cummings (04/15/2004)      |
| Modified by: Ulric Lund (08/17/2006)              |
|_________________________________________________*/;

%macro rvm(n, mu, kappa, dataOut, seed=0);

  data &dataOut (keep = rvm);
    a = 1 + sqrt(1 + 4*&kappa**2);
    b = (a - sqrt(2*a))/(2*&kappa);
    r = (1 + b**2)/(2*b);
    obs = 0;

    do while (obs < &n);
      U1 = ranuni(&seed);
      z = cos(U1*constant('PI'));
      f = (1 + r*z)/(r + z);
      c = &kappa*(r - f);
      U2 = ranuni(&seed);

      if c*(2-c)-U2 > 0 then do;
        U3 = ranuni(&seed);
        rvm = sign(U3-0.5)*arcos(f) + &mu;
        rvm = mod(rvm,2*constant('PI'));
        obs = obs + 1;
        output;
      end;

      else do;
        if log(c/U2) + 1 - c >= 0 then do;
          U3 = ranuni(&seed);
          rvm = sign(U3-0.5)*arcos(f) + &mu;
          rvm = mod(rvm,2*constant('PI'));
          obs = obs + 1;
          output;
        end;
      end;
    end;
  run;

%mend rvm;


/*__________________________________________________
|                                                   |
|                       rwc                         |
|                     ________                      |
|                                                   |
| rwc generates pseudo random numbers from the      |
| wrapped Cauchy distribution.  It returns an output|
| data set with the following variable:             |
|                                                   |
| rwc - contains random numbers generated           |
|___________________________________________________|
|           |               |                       |
| Parameter | Default       | Description           |
|___________|_______________|_______________________|
|           |               |                       |
| n         |               | number of observations|
|           |               | to be generated       |
|           |               |                       |
| mu        |               | mean direction of     |
|           |               | wrapped Cauchy        |
|           |               | distribution          |
|           |               |                       |
| rho       |               | concentration         |
|           |               | parameter of wrapped  |
|           |               | Cauchy distribution;  |
|           |               | number between 0 and 1|
|           |               |                       |
| dataOut   |               | output SAS data set   |
|           |               |                       |
| seed      |      0        | seed value; default   |
|           |               | value of zero leads   |
|           |               | to non-replicable     |
|           |               | results               |
|___________|_______________|_______________________|
|                                                   |
| Written by: Michael A. Cummings (04/15/2004)      |
| Modified by: Ulric Lund (08/17/2006)              |
|_________________________________________________*/;

%macro rwc(n, mu, rho, dataOut, seed=0);

  data _calc;
    retain seed &seed;
    do i=1 to &n;
      call rancau (seed,U);
      output;
    end;
  run;

  data &dataOut (keep = rwc);
    set _calc;
    scale = -log(&rho);
    rwc = mod(U*scale + &mu, 2*constant('PI'));
  run;

  proc datasets nodetails nolist;
    delete _calc;
  run;quit;

%mend rwc;


/*__________________________________________________
|                                                   |
|                       rwn                         |
|                     ________                      |
|                                                   |
| rwn generates pseudo random numbers from the      |
| wrapped normal distribution.  It returns an output|
| data set with the following variable:             |
|                                                   |
| rwn - contains random numbers generated           |
|___________________________________________________|
|           |               |                       |
| Parameter | Default       | Description           |
|___________|_______________|_______________________|
|           |               |                       |
| n         |               | number of observations|
|           |               | to be generated       |
|           |               |                       |
| mu        |               | mean direction of     |
|           |               | wrapped normal        |
|           |               | distribution          |
|           |               |                       |
| rho       |               | concentration         |
|           |               | parameter of wrapped  |
|           |               | normal distribution;  |
|           |               | number between 0 and 1|
|           |               |                       |
| dataOut   |               | output SAS data set   |
|           |               |                       |
| seed      |      0        | seed value; default   |
|           |               | value of zero leads   |
|           |               | to non-replicable     |
|           |               | results               |
|___________|_______________|_______________________|
|                                                   |
| Written by: Michael A. Cummings (04/15/2004)      |
| Modified by: Ulric Lund (08/17/2006)              |
|_________________________________________________*/;

%macro rwn(n, mu, rho, dataOut, seed=0);

  data _calc;
    retain seed &seed;
    do i=1 to &n;
      call rannor (seed,U);
      output;
    end;
  run;

  data &dataOut (keep = rwn);
    set _calc;
    sigma = sqrt(-2*log(&rho));
    rwn = mod(U*sigma + &mu, 2*constant('PI'));
  run;

  proc datasets nodetails nolist;
    delete _calc;
  run;quit;

%mend rwn;


/*______________________________________________
|                                               |
|                 trigMoment                    |
|                   _____                       |
|                                               |
| trigMoment macro computes the specified order |
| trigonometric moment for a circular variable. |
| It returns an output data set with the        |
| following variables:                          |
|                                               |
| n - number of observations in input data set  |
| rhoHatp - amplitude of the pth trigonometric  |
|           moment, where the p at the end of   |
|           the variable name will be the       |
|           numeric value of the specified      |
|           order                               |
| muHatp - direction of the pth trigonometric   |
|          moment                               |
|_______________________________________________|
|                                               |
| Dependencies: uses macros circMean            |
|                           circDisp            |
|_______________________________________________|
|           |           |                       |
| Parameter |  Default  |  Description          |
|___________|___________|_______________________|
|           |           |                       |
| dataIn    |           | the input data set    |
|           |           | measured in radians   |
|           |           |                       |
| dataOut   |           | the output data set   |
|           |           | returned by the macro |
|           |           |                       |
| var       |           | the variable in the   |
|           |           | input data set that   |
|           |           | you want to find      |
|           |           | the summary of        |
|           |           |                       |
| p         |   1       | OPTIONAL: the order of|
|           |           | the trigonometric     |
|           |           | moment to be computed.|
|           |           | Default is for a first|
|           |           | order moment.         |
|           |           |                       |
| center    |   F       | OPTIONAL: logical flag|
|           |           | whether to compute    |
|           |           | centered moments or   |
|           |           | not. Default is to    |
|           |           | compute an uncentered |
|           |           | moment.               |
| print     |   T       | OPTIONAL: print out   |
|           |           | results to output     |
|           |           | window. Set to F to   |
|           |           | suppress printing.    |
|___________|___________|_______________________|
|                                               |
| Written by: Kristin Bishop (03/02/2004)       |
| Modified by: Ulric Lund (08/17/2006)          |
|_____________________________________________*/;

%macro trigMoment(dataIn, dataOut, var, p=1, center=F, print=T);

  %circMean(&dataIn, _meanOut , &var, print=F);

  %if &p=1 and &center=F %then %do;
    %circDisp(&dataIn, _dispOut , &var, print=F);
    data &dataOut;
      merge _meanOut (rename=(meanDir=muHat1))
            _dispOut (rename=(rbar=rhoHat1));
      keep n muHat1 rhoHat1;
    run;
    
    proc datasets nodetails nolist;
      delete _meanOut _dispOut;
    run;quit;
  %end;
  
  %else %do;
    data _calc;
      if _n_=1 then set _meanOut;
      set &dataIn;
      %if &center=T %then %do;
        center=1;
      %end;
      %else %if &center=F %then %do;
        center=0;
      %end; 

      sinP=sin(&p*(&var-meanDir*center));
      cosP=cos(&p*(&var-meanDir*center));

      sumSin + sinP;
      sumCos + cosP;

      Sp=sumSin/n;
      Cp=sumCos/n;

      rhoHat&p = sqrt(Cp**2 + Sp**2);
      muHat&p  = atan2(Sp,Cp);
    run;

    data &dataout;
      set _calc end=last;
      if last;
      keep n muHat&p rhoHat&p;
    run;
    
    proc datasets nodetails nolist;
      delete _meanOut _calc;
    run;quit;
  %end;

  %if &print=T %then %do;
    proc print data=&dataOut noobs;
    run;
  %end;

%mend trigMoment;