FPCASCORE Procedure

This example demonstrates how to use the FPCASCORE procedure to score new functional data by using a previously trained FPCA model. The data consist of 50 simulated functional curves that are evaluated at 100 equally spaced time points over the interval . The curves are split evenly between two classes that have different mean functions: class 0 corresponds to curves y1 through y25, and class 1 corresponds to curves y26 through y50. Class 0 curves exhibit a symmetric peak that is centered at , and class 1 curves have a smaller peak that is shifted toward . Random variation is added along two orthonormal eigenfunctions.

The following IML procedure statements generate the simulated data set func_curves. Each column (y1 to y50) corresponds to a curve, and each row corresponds to a time point.

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

   /* Define two 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));   /* Slightly shifted & smaller bump */
   finish;

   /* Define two 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;

   /* Simulate normalized FPCA scores */
   Ksi = j(N, 2);
   call randgen(Ksi, "Normal");
   Ksi = (Ksi - mean(Ksi)) / std(Ksi);   /* standardize */
   Ksi = Ksi * diag({5, 2});

   /* Evaluate eigenfunctions and class-specific means */
   eig1 = eigFunct1(s);
   eig2 = eigFunct2(s);
   meanVec0 = meanFunct_class0(s);  /* Class 0 */
   meanVec1 = meanFunct_class1(s);  /* Class 1 */

   /* Allocate curve matrix */
   y = j(N, M, .);

   /* Fill curves: first half = class 0, second half = class 1 */
   do i = 1 to N;
      if i <= N/2 then
         y[i,] = Ksi[i,] * (eig1 // eig2) + meanVec0;
      else
         y[i,] = Ksi[i,] * (eig1 // eig2) + meanVec1;
   end;

   /* Prepare output: time points + y1 to y50 */
   s_col = s`;
   s_y = s_col || y`;  /* Transpose so each y is a column */
   y_colnames = "y1":"y50";
   colnames = {"_TIMEPOINTS_"} || y_colnames;

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

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;

You can then use the FPCA procedure as follows to train a functional principal component model and save it to an analytic store. Then you apply PROC FPCASCORE as shown to score the same data by using the saved model and extract the first two principal component scores for visualization. The results are displayed in Output 13.1.1.

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

proc fpcascore data=mylib.func_curves model=mylib.fpca_model;
   input y1-y50;
   output out=mylib.fpca_scores npc=2;
run;

/* Assign class label by observation order */
data fpca_scores_labeled;
   set mylib.fpca_scores;
   if _N_ <= 25 then class = 0;
   else class = 1;
run;

/* Plot PC1 vs PC2 colored by class */
ods graphics / reset width=6in height=5in;

proc sgplot data=fpca_scores_labeled;
   scatter x=_PCS_1 y=_PCS_2 / group=class
           markerattrs=(symbol=CircleFilled size=9);
   xaxis label="PC1 Score";
   yaxis label="PC2 Score";
   title "FPCA Score Plot Colored by Class";
run;

Output 13.1.1: First Two FPCA Scores Colored by Class

The plot in Output 13.1.1 shows that the first two functional principal component scores mostly separate the two simulated classes, showing distinct clustering along with some overlap. This demonstrates that FPCA effectively captures the dominant variations that are introduced by the class-specific mean structure in the functional data.