stats.ml

open Eval

exception Illegal
type funct = float list
type dataset = (float * float) list

(* POST: [rsq dset] is a float in the range 0..1. that represents the r squared
 * value for a linear regression performed on [dset].
 * PRE: [dset] is a (float * float) list of finite length.
 *)
let rsq dset =
	let abr = lin_reg_stats dset in
	if List.length abr < 3 then raise Illegal
	else List.nth abr 2

(* POST: [residuals ?func dset] is a float list that represents the "residual"
 * for each element in [dset]; that is, the difference between the y-value and the
 * expected y-value given by the regression model [?func]. [residuals dset]
 * will be the list of residuals for a linear regression on [dset]. [residuals
 * func dset] will be the list of residuals for the user-input function and [dset].
 * PRE: [dset] is a (float * float) list of finite length. [func] if provided is
 * a float list of length 2, 3, or 4.
 *)
let residuals ?func dset =
	let xylst = List.split dset in
	let xlst = fst xylst in
	let ylst = snd xylst in
	let func = match func with None -> lin_reg dset | Some x -> x in
	if List.length dset < 2 then raise Illegal else
	if List.length func < 2 then raise Illegal
	else 
	let yhatlst = generate_y_lst func xlst in
	let ypairlst = List.combine ylst yhatlst in
	List.map (fun (y1,y2) -> y1 -. y2) ypairlst

(* POST: [residuals_plot ?func dset] is a (float * float) list with elements of
 * the form (x-coordinate, residual) where the residuals are from the float list
 * generated by [residuals ?func dset].
 * PRE: [dset] is a (float * float) list of finite length. [func] if provided is
 * a float list of length 2, 3, or 4.
 *)
let residuals_plot ?func dset =
	let xlst = List.map (fun (x,y) -> x) dset in
	let reslst = residuals ?func dset in
	List.combine xlst reslst

(* POST: [residualssq ?func dset] is a float list with elements that represent
 * the square of the residuals generated by [residuals ?func dset].
 * PRE: [dset] is a (float * float) list of finite length. [func] if provided is
 * a float list of length 2, 3, or 4.
 *)
 let residualssq ?func dset =
	let reslst = residuals ?func dset in
	List.map (fun x -> x**2.) reslst


(* POST: [variance ?func dset] is a float that represents the variance of
 * a regression model [func] applied to a dataset [dset]. If the optional argument
 * [func] is provided, it will return the variance for the dataset modeled on
 * the user-input [func]; if [func] is not provided, it will return the variance
 * for the dataset modeled on a linear regression on [dset].
 * PRE: [dset] is a (float * float) list of finite length. [func] if provided is
 * a float list of length 2, 3, or 4.
 *)
let variance ?func dset =
	let ressq = residualssq ?func dset in
	let n = float_of_int (List.length dset) in
	sum(ressq)/.(n -. 1.)

(* POST: [stderror ?func dset] is a float that represents the standard error of
 * a regression model [func] applied to a dataset [dset]. If the optional argument
 * [func] is provided, it will return the standard error for the dataset modeled on
 * the user-input [func]; if [func] is not provided, it will return the standard
 * error for the dataset modeled on a linear regression on [dset].
 * PRE: [dset] is a (float * float) list of finite length. [func] if provided is
 * a float list of length 2, 3, or 4.
 *)
 let stderror ?func dset =
	let n = float_of_int (List.length dset) in
	let ressq = residualssq ?func dset in
	sqrt(sum(ressq)/.(n))

(* POST: [stderror ?func dset] is a float that represents the standard deviation of
 * the elements in [dset] applied to the regression model [func]. If the optional
 * argument [func] is provided, it will return the standard deviation for the
 * dataset modeled on the user-input [func]; if [func] is not provided, it will
 * return the standard deviation of the dataset modeled on a linear regression
 * on [dset].
 * PRE: [dset] is a (float * float) list of finite length. [func] if provided is
 * a float list of length 2, 3, or 4.
 *)
 let stddev ?func dset = sqrt(variance ?func dset)

(* POST: [chisq ?func dset] is a float that represents the chi squared test on the
 * regression model [func] with the dataset [dset]. If the optional
 * argument [func] is provided, it will perform the chi squared test on the
 * user-input [func]; if [func] is not provided, it will perform the chi squared
 * test on a linear regression on [dset].
 * PRE: [dset] is a (float * float) list of finite length where observed values,
 * i.e. y-values, are >=1.. [func] if provided is a float list of length 2, 3, or 4.
 *)
 let chisq ?func dset =
	let ressq = residualssq ?func dset in
	let ylst = List.map (fun (x,y) -> y) dset in
	let resylst = List.combine ressq ylst in
	sum(List.map (fun (e,y) -> e/.y) resylst)


let rec fact x =
	match x with
	| 0. -> 1.
	| y -> x *. fact (x-.1.)

(*POST: [poisson mean x] is a float calculating the Poisson distribution
* from the given values
* PRE: mean is a float greater than 0 and x is a nonnegative int
*)
let poisson mean x =
	if mean <= 0. || x<0 then raise Illegal else
	(2.71828**(-.mean))*.(mean**(float x))/.(fact (float x))

(*POST: [poissonc mean x] is a float calculating the cumulative Poisson 
* distribution from the given values
* PRE: mean is a float greater than 0 and x is a nonnegative int
*)
let rec poissonc mean x =
	if mean <= 0. || x<0 then raise Illegal else
	match x with
	| 0 -> poisson mean 0
	| y -> (poisson mean y) +. (poissonc mean (y-1))

(*POST: [binomial n p x] is a float calculating the binomial distribution
* from the given values
* PRE: n is an int greater than 0, p is a float between 0 and 1, and x is
* an int less than or equal to n
*)
let binomial n p x =
	if n<=0 || p<0. || p>1. || x>n then raise Illegal else
	(((fact (float n))/.((fact (float x))*.
		(fact (float (n-x))))))*.(p**(float x))*.((1.-.p)**(float (n-x)))

(*POST: [binomialc n p x] is a float calculating the cumulative binomial 
* distribution from the given values 
* PRE: n is an int greater than 0, p is a float between 0 and 1, and x is
* an int less than or equal to n
*)
let rec binomialc n p x =
	if n<=0 || p<0. || p>1. || x>n then raise Illegal else
	match x with
	| 0 -> binomial n p x
	| y -> (binomial n p y)+. (binomialc n p (y-1))