/*

Script to cross-validate the regression model defined by the model dialog "Model"
in the current data table.  

It currently does split-sample cross-validation (half-half) for 20 random splits. These
are both easy to change.  Change the variables:
     holdBackPct (set to 0.5 for half)
     numberOfIterations (set to 50 now)

The script opens a new data table with two columns, one for the in-sample errors and one
for the out-of-sample errors.

*/

//    pick the percentage to hold back (its set by default to half)
holdBackPct = 0.1;
//    pick the number of iterations 
numberOfIterations = 50;

// ------------- dont make changes below here unless you know what you're doing -------------

// Make sure that you have opened the correct data table and use column names from that table
dt = current data table();

// Read model script from data table
theModel = dt << GetTableProperty("Model");
Insert Into( theModel, Expr( Invisible ) );

// Extract first argument of model script then extract its first argument as a string
yVar = Char(Arg( Arg( theModel, 1 ), 1 ));

//    figure out how many rows to exclude
n = nRows(dt);
k = nCols(dt);
nExclude = floor(n*holdBackPct);

//    new data table to hold results
resultsTable = new table("Cross-Validation Results.jmp", 
							new column("Training_rmse", numeric),
							new column("Testing_rmse", numeric));		

for (iter=0, iter<numberOfIterations, iter++,
	current data table(dt);
	//    exclude rows at random
	dt << clear row states();
	dt << select randomly(nExclude) << exclude;
	//    fit model, save prediction formula
	regr = theModel << Run() << Run Model;
	inestimates = (regr<<report)["Summary of Fit",NumberColBox(1)] << get as matrix;
	DF = (regr<<report)["Analysis of Variance",NumberColBox(1)] << get as matrix;

	// show( inestimates );
	regr << save columns(predicted values) << close window;
	// get the column that holds response and predicted values
	yCol = Column(dt,yVar);
	predCol = Column(dt, N Col( dt ) );
	// compute new column (y-y^)^2, adjusting for DF in the prediction out of sample (RMSE in sample)
	sqrError = Eval(
		Eval Expr(
			dt << New Column( "Squared Error",
				Numeric,
				Formula( (((Expr( yCol )[Row()] - Expr( predCol )[Row()])) ^ 2) )
             )
     ));
    // get properties of insample errors; this would need adjustment for df which is what RMSE does
	// dist = dt << Distribution(invisible, Column(sqrError));	
	// inestimates = (dist<<report)["Moments",NumberColBox(1)] << get as matrix;
	// dist << close window;
	// flip cases used
	dt << exclude << invert row selection << exclude;
	dist = dt << Distribution(invisible, Column(sqrError));
	outestimates = (dist<<report)["Moments",NumberColBox(1)] << get as matrix;
	dist << close window;
	dt << delete columns({sqrError,predCol});
	// put squared errors into the results table
	resultsTable<< add row({Training_rmse=sqrt(1+DF[1]/DF[3])*inestimates[3],Testing_rmse=sqrt(outestimates[1])});
);

dt << clear row states();