FPCA Procedure

This example demonstrates how to use the FPCA procedure on simulated functional data that have a known low-dimensional structure. The data consist of 50 functional observations; each observation is recorded at 100 equally spaced time points in the interval . Each curve is generated as a linear combination of two orthonormal basis functions, and , and class-specific mean functions and random coefficients are drawn from a standardized Gaussian distribution that is scaled by eigenvalues. The first 25 curves belong to class 0 and use a mean function whose peak is centered at . The remaining 25 curves belong to class 1 and have a mean peak around .

The following IML procedure statements generate the simulated data set func_curves. Each column (y1 to y50) corresponds to a functional observation, and each row represents a time point. The variable _TIMEPOINTS_ stores the common time grid across all observations.

proc iml;
   call randseed(123);
   N = 50;     /* number of curves */
   M = 100;    /* number of time points */
   s = do(0, 10, 10 / (M - 1));

   /* Class-specific mean functions */
   start meanFunct_class0(s);
      return (s + 10 # exp(-(s - 5)##2));
   finish;
   start meanFunct_class1(s);
      return (s + 8 # exp(-(s - 7)##2));
   finish;

   /* Shared eigenfunctions */
   start eigFunct1(s);
      return (cos(2 * s * constant('PI') / 10) / sqrt(5));
   finish;
   start eigFunct2(s);
      return (-sin(2 * s * constant('PI') / 10) / sqrt(5));
   finish;

   /* Generate scores */
   Ksi = j(N, 2);
   call randgen(Ksi, "Normal");
   Ksi = (Ksi - Ksi[:,]) / sqrt(var(Ksi)[,]);
   Ksi = Ksi * diag({5, 2});

   /* Evaluate basis and means */
   eig1 = eigFunct1(s);
   eig2 = eigFunct2(s);
   basis = eig1 // eig2;
   meanVec0 = meanFunct_class0(s);
   meanVec1 = meanFunct_class1(s);

   /* Build curves */
   y = j(N, M, .);
   do i = 1 to N;
      if i <= N/2 then
         y[i,] = Ksi[i,] * basis + meanVec0;
      else
         y[i,] = Ksi[i,] * basis + meanVec1;
   end;

   /* Export wide-format data */
   s_y = s` || y`;
   y_colnames = "y1":"y50";
   colnames = {"_TIMEPOINTS_"} || y_colnames;

   create func_curves from s_y[colname=colnames];
   append from s_y;
   close func_curves;
quit;

After generating the data, you can use the following statements to convert the data to long format and visualize the curves by class. The results are shown in Output 12.1.1.

/* Transpose to long format */
proc transpose data=func_curves out=func_long(rename=(col1=value)) name=id;
   by _TIMEPOINTS_;
run;

/* Create attribute map to color curves by class */
data attrmap;
   length id value $10 linecolor $8;
   do i = 1 to 50;
      value = "y" || strip(put(i, best.));
      id = "idmap";
      if i <= 25 then linecolor = "blue";
      else linecolor = "red";
      output;
   end;
run;

/* Attach attrid to plot data */
data func_long_class;
   set func_long;
   attrid = "idmap";
run;

/* Sort data */
proc sort data=func_long_class;
   by id _TIMEPOINTS_;
run;

/* Plot curves colored by class */
ods graphics / reset width=7in height=5in;

proc sgplot data=func_long_class dattrmap=attrmap;
   series x=_TIMEPOINTS_ y=value / group=id attrid=idmap
          lineattrs=(thickness=1) transparency=0.7;
   xaxis label="t";
   yaxis label="X(t)";
   title "Simulated Functional Curves Colored by Class";
run;

Output 12.1.1: Plot of Functional Curves by Mean Curve Class

Output 12.1.1 displays all 50 functional curves over the common time domain. Each curve is plotted as a line, and the coloring reflects class membership: curves from class 0 (observations y1 to y25) are shown in blue, and curves from class 1 (y26 to y50) are shown in red. This color coding makes it easier to visually assess differences between the two classes. In this example, the two classes exhibit different average behaviors over time, as defined by their respective mean functions. The plot highlights these differences, providing useful information about the functional structure and demonstrating the value of using FPCA for dimension reduction.

You can use the following DATA step to load the generated func_curves data set into your CAS session by using your CAS engine libref:

   data mylib.func_curves;
      set func_curves;
   run;

Then you fit a functional principal component model by using the following PROC FPCA statement:

proc fpca data=mylib.func_curves nBins=50;
   input y1-y50;
   savestate rstore=mylib.fpca_model;
run;

The SAVESTATE statement saves the trained FPCA model to an analytic store, enabling you to efficiently reuse it for future scoring. A downstream example that shows how to score new observations and visualize class separation by using FPCA scores is provided in the documentation of the FPCASCORE procedure.