############################################################################
#
# $RCSfile: psp $
# $Revision: 2.0 $
# $Date: 1997/10/03 14:28:00 $
# $Author: ssolie $
#
############################################################################
#
# Copyright (c) 1997 Software Evolution.  All Rights Reserved.
#
############################################################################
#
# PSP Support Library -- Personal Software Process (PSP) support functions
#
# This set of Maple V Release 3 functions contains many functions which can
# be used to process PSP data.  Each function is designed to perform a
# specific action or produce a specific report.
#
# These functions may or may not be exactly what is required for your
# specific PSP means.  However, the functions do provide a basis from which
# to build upon to improve and enhance your PSP processes.
#
# For more information about the PSP, please refer to:
#	Humphrey, Watts S., "A Discipline for Software Engineering"
#	(Don Mills, Ontario: Addison-Wesley, 1995), ISBN 0-201-54610-8.
#
# WARNING:
#  Please note that not all of the functions properly protect themselves
#  from incorrect parameters and may have unexpected results when used
#  improperly.  Always check the source code to be sure.


###### prediction_interval ######
`help/text/prediction_interval` := TEXT(
	`FUNCTION: prediction_interval - Calculate range and prediction interval`,
	``,
	`CALLING SEQUENCE:`,
	`   prediction_interval(x_data, y_data, estimate, percent)`,
	``,
	`PARAMETERS:`,
	`   x_data   - set of x regression data (estimated values)`,
	`   y_data   - set of y regression data (actual values)`,
	`	estimate - planned or estimated value`,
	`   percent  - prediction interval percentage`,
	``,
	`SYNOPSIS:`,
	`- Data is supplied as arrays of data equal in length.`,
	``,
	`- It is useful to know the prediction error when making estimates.  This`,
	`  function uses your historical data and the linear regression method to`,
	`  predict your likely estimation error.  This will give you the ability to`,
	`  know the likely minimum and maximum estimation based on your previous`,
	`  performance.`,
	``,
	`- The prediction interval percentage means that most of the time the`,
	`  true result will be within the interval.  For example, if you ask for`,
	`  a 90% prediction interval then 90% of the time the true result will`,
	`  be in the interval.  The most common values for the prediction`,
	`  interval used in the PSP are 70% and 90%.`,
	``,
	`- As with all predictions, there is a danger of misinterpretation.  If`,
	`  your process changes such that it invalidates your historical data,`,
	`  the prediction interval will be meaningless.  Exercise caution before`,
	`  committing to the results.`,
	``,
	`- This function requires the stats package, and so can be used only after`,
	`  performing the command with(stats).`,
	``,
	`EXAMPLES:`,
	`> with(stats);`,
	`> plan_loc := [100,327,75,92]: actual_loc := [110,228,90,100];`,
	`> prediction_interval(plan_loc,actual_loc,300,70);`
):

prediction_interval := proc(x_data, y_data, estimated, percent)
	local xk, i, n, beta0, beta1, variance, std_dev, p, t, r, x_avg, y_avg,
		  range, lower_limit, upper_limit, projected;

    ### Parse arguments ###
	if percent < 0 or percent > 100 then
		ERROR(`Percent must be in the range 0..100`);
	fi;

	if estimated < 1 then
		ERROR(`Estimated must be greater than 1`);
	fi;

	### calculate and print out the results ###
	xk		:= estimated;
	n		:= describe[count](x_data);
	r 		:= describe[linearcorrelation](x_data, y_data);
	x_avg   := describe[mean](x_data);
	y_avg   := describe[mean](y_data);

	beta1 := (sum(x_data[i]*y_data[i], i=1..n) - n*x_avg*y_avg) /
			 (sum(x_data[i]^2, i=1..n) - n*x_avg^2);
	beta0 := y_avg - beta1*x_avg;

	projected := beta0 + beta1*xk;
	projected := ceil(projected);

	variance := 1/(n-2) *
		sum((y_data[i] - beta0 - beta1*x_data[i])^2, i=1..n);
	std_dev := sqrt(variance);

	p := (1 + (percent/100))/2;
	t := statevalf[icdf, studentst[n - 2]](p);
	range := t * std_dev * sqrt(1 + 1/n + (xk-x_avg)^2/sum((x_data[i] - x_avg)^2, i=1..n));
	range := ceil(range);

	lower_limit := ceil(projected - range);
	upper_limit := ceil(projected + range);

	print(`N` = n);
	print('r^2' = evalf(r^2));
	print('beta[0]' = evalf(beta0));
	print('beta[1]' = evalf(beta1));
	print('Projected' = projected);
	print(`Range` = range);
	print(`UPI` = upper_limit);
	print(`LPI` = lower_limit);
end;


###### plot_relationship ######
`help/text/plot_relationship` := TEXT(
	`FUNCTION: plot_relationship - Plots data relationship`,
	``,
	`CALLING SEQUENCE:`,
	`   plot_relationship(x_data, y_data)`,
	``,
	`PARAMETERS:`,
	`   x_data - set of x data`,
	`   y_data - set of y data`,
	``,
	`SYNOPSIS:`,
	`- Data is supplied as arrays of numbers equal in length.`,
	``,
	`- Two measurements are used to determine the strength of the linear`,
	`  relationship between the two sets of data: r^2 and p.`,
	`    r^2 - square of the correlation (r^2 >= 0.9 is very good)`,
	`    p   - correlation significance (p <= 0.05 is very good)`,
	``,
	`- A high correlation and significance does not guarantee a linear relation-`,
	`  ship exists between the two data sets.  It only provides evidence that a`,
	`  linear relationship might exist.  You must fully understand the processes`,
	`  involving the data values before making any meaningful conclusions.`,
	``,
	`- This function requires the stats package, and so can be used only after`,
	`  performing the command with(stats)`,
	``,
	`EXAMPLES:`,
	`> with(stats);`,
	`> plan_loc := [100,327,75,92]: actual_loc := [110,228,90,100];`,
	`> plot_relationship(plan_loc,actual_loc);`
):

plot_relationship := proc(x_data, y_data)
	local sigma, r, t, n, p, x, x_avg, y_avg, beta1, beta0, line,
	      graph, points, max_x;

	n 		:= describe[count](x_data);
	r 		:= describe[linearcorrelation](x_data, y_data);
	t 		:= abs(r) * sqrt(n - 2) / sqrt(1 - r^2);
	p 		:= 2 * (1 - statevalf[cdf, studentst[n-2]](t));
	x_avg   := describe[mean](x_data);
	y_avg   := describe[mean](y_data);
	max_x	:= evalf(op(2, describe[range](x_data)));
	beta1	:= (sum(x_data[i]*y_data[i], i=1..n) - n*x_avg*y_avg) /
			   (sum(x_data[i]^2, i=1..n) - n*x_avg^2);
	beta0	:= y_avg - beta1*x_avg;

	points := convert([x_data, y_data], matrix);
	points := linalg[transpose](points);
	points := convert(points, listlist);

	line := [[0, beta0], [max_x, evalf(beta1*max_x + beta0)]];

	graph := PLOT(
		SYMBOL(CROSS),
		STYLE(PATCH),
		POINTS(points),
		CURVES(line)
	);

	print(`N` = n);
	print('beta[0]' = evalf(beta0));
	print('beta[1]' = evalf(beta1));
	print('r^2' = evalf(r^2));
	print('p' = evalf(p));
	print(graph);
end;


save `Maple:PSP/psp.m`;
#quit