Esta nueva sesión de OCaml va a continuar un ejercicio de aplicacion similar a los que he encontrado en el cookbook de OCaml, en este caso va a ser el cálculo de la regresión lineal simple. Es además un buen ejemplo para proponer como receta.

Regresion Lineal Simple

La regresion lineal simple es un método muy usado en técnicas experimentales, para relacionar dos variables. El contexto es un experimento donde tomamos muestras de dos variable aleatorias $X$ e $Y$, donde consideramos $X$ la variable aleatoria independiente y $Y$ la variable aleatoria dependiente o de respuesta. Además introducimos como premisa que poseen una relacion lineal de modo que los valores siguen la relacion $y~α+βx$. Esto ya nos situa en un marco de análisis de datos supervisado (proponemos nosotros la relación funcional).

Así pues cada muestra posee dos valores $(x_i, y_i)$ de modo que nuestra primera hipótesis es que $y_i = α + βx_i +e_i$ donde $e_i$ son los errores de muestreo o experimentales. La segunda hipótesis es la distribución de la variable aleatoria $e$ que se suele considerar normal con media cero $N(0,σ)$. Sin embargo, en los cursos iniciales de técnicas experimentales se habla de método de mínimos cuadrados, ya que para simplificar la explicación e introdución del método se considera sencillamente que el objetivo es minimizar la suma de los cuadrados $(e_i)^2$, es decir de los errores (tomados siempre como positivos). Este método da lugar a dos ecuaciones para estimar los valores de $\hat α $ y $\hat β$ que puedes ver en wikipedia y que son las que vamos a implementar en OCaml.

Código usando Array

Programar esto es ya sencillo desde estás expresiones y podemos hacerlo de modo sistemático definiendo funciones que son útiles más alla de este cálculo. Todo ello usando la libreria Array

Definimos una funcion sum que suma los elementos de un Array unidimensional (o vector), una funcion media o avg que calcula la media de los valores de un array. Y por ultimo una funcion bias que construye un array donde a cada elemento se le sustrae la media.

let sum  v   = Array.fold_left (+.) 0. v 
let avg  v   = (sum v) /. float_of_int (Array.length v)
let bias v   = Array.map (fun x -> x -. (avg v)) v

También si definimos el producto escalar de dos vectores como:

let dot_vec a b = Array.map2 ( *. ) a b |> Array.fold_left ( +. ) 0.

Entonces las ecuaciones se convierten en:

let slope x y      = (dot_vec (bias x) (bias y) ) /. ( dot_vec (bias x) (bias x) )
let inter x y beta = avg y -. beta *. (avg x)
let linear_regression x y = 
        let beta  = slope x y in 
        let alpha = inter x y beta  in 
        (alpha, beta)

Podemos comprobar el algoritmo:

let () = 
    let x = [|1.0;2.0;3.0;4.0|]  in 
    let y = [|2.0001;4.009;6.000;8.005|]  in 
    let a,b = linear_regression x y in 
    Printf.printf "y = a + bx ; a = %f ; b=%f" a b

Conclusión

Vemos como implementar esto en OCaml usando la libreria Array es bastante directo y rápido, es posible que no sea la implementación más rápida posible dentro de OCaml ya que nuestro objetivo ha sido usar la expresividad del lenguage más que buscar un cálculo optimizado.