Esta nueva sesión de OCaml voy a implementar la receta del libro de recetas de OCaml (cookbook recipes) para el problema de calcular los vecinos próximos (nearest neighbour) de un conjunto de datos.
Vecinos más próximos
Es un problema que se puede formular dentro de la estadística aplicada; así pues, imaginemos N elementos en un espacio (qué podré representar por vectores) que son muestras de una variable aleatoria multidimensional. Imaginemos un elemento de referencia (a menudo es un elemento de este conjunto). ¿Cuáles son los elementos del conjunto más cercanos a esta referencia? La respuesta inicial es que depende de como definamos la distancia entre dos elementos (que en matemáticas se llama métrica). Una distancia es una función que dados dos elementos de este espacio devuelve la distancia entre ellos, pero debe tener una serie de propiedades para poder calificarse como distancia.
Vamos a considerar que los elementos son vectores con tres elementos (con valores reales). En nuestro caso hacemos uso de que los elementos serán parte un espacio vectorial real tridimensional y la distancia más usual se llama distancia euclídea que es la que implementaremos, pero nuestra definición de vecinos próximos (como función) intentará tener la definición de distancia entre sus argumentos y no predeterminada.
Operaciones
Como vamos a utilizar una distancia euclídea entre dos vectores (que serán dos arrays cuyos elementos son floats) entonces su definición será:
$d(\bo x,\bo y)=√(∑↙{i=0}↖3 (x_i-y_i)^2)$
pero esto es análogo a decir que
$d(\bo x, \bo y)=√[(\bo x- \bo y)∙(\bo x-\bo y)]$
donde $∙$ es el producto escalar (dot product en inglés), de dos vectores (que representamos como $\bo x$ e $\bo y$). Por tanto, veamos como quedan estas definiciones en OCaml usando el módulo Array de la librería standard:
(* Operator `*@` is the typical dot product between vectors, in this case with floats elements.*)
let ( *@ ) u v = Array.map2 ( *. ) u v |> Array.fold_left ( +. ) 0.
(* Operator `*!` defines multiplication of a scalar (float) by a vector *)
let ( *! ) a v = Array.map ( fun x -> a *. x ) v
(* Operator `+:` define the sum of two vectors *)
let ( +: ) u v = Array.map2 ( +. ) u v
(* Operator `-:` is introduced for convenience/simplicity in the `metric_p2`, we could also define
it directly as `let ( -: ) u v = Array.map2 ( +. ) u ( -1.0 *! v )` *)
let ( -: ) u v = u +: (-1.0 *! v)
(* `metric_p2` is the Euclidean distance based on the p-norm with p=2, it is the most typical distance. *)
let metric_p2 u v = sqrt (( u -: v ) *@ ( u -: v) )
Distancias de un grupo de puntos
Como decíamos tenemos que calcular las distancias entre un elemento de referencia y un conjunto de muestras de una variable aleatoria, definimos, pues, una función distance_in_set:
(* `distance_in_set` calculates all the distances in a set of elements `s` with respect to a reference
element `ref` using a `metric` function. It returns a tuple `int * float` for the index of the element
in the set and the distance of that element with respect to the reference. *)
let distance_in_set ref s metric =
let n = Array.length s in
Array.init n (fun i -> (i, metric ref (Array.get s i)) )
Como vemos ref es nuestro elemento de referencia, s la muestra que posee el conjunto de elementos a evaluar, y distance es la función que calcula las distancias entre elementos (vectores). Hemos elegido que esta función devuelva una pareja (tuple) de datos: el índice del elemento en nuestro conjunto de datos (el conjunto de datos será un array de un array), y la distancia de ese elemento a nuestra referencia.
La segunda parte es ordenar este conjunto basándonos en el valor de las distancias:
(* `sort_by_distance` will sort the output of the `distance_in_set` with in-place replacement, so the
indices of the k nearest neighbours will be the k first elements of the array returned by this function *)
let sort_by_distance (d:(int * float) array) =
let cmd x y =
if (snd x -. snd y) = 0.0 then 0
else if (snd x -. snd y) > 0.0 then 1
else -1
in
Array.sort cmd d
De este modo, obtener el vecino más próximo a una referencia será:
(* `get_nn s d` returns the nearest neighbour (first element of the sorted by distance set *)
let get_nn s d =
let nn_idx = fst (Array.get d 0) in
(nn_idx, Array.get s nn_idx)
Programa completo
El programa completo tiene algunas funciones de soporte para definir elementos de nuestro espacio muestral (en nuestro caso usando una distribución uniforme) y otras funciones auxiliares para imprimir los resultados.
(* Operator `*@` is the typical dot product between vectors, in this case with floats elements.*)
let ( *@ ) u v = Array.map2 ( *. ) u v |> Array.fold_left ( +. ) 0.
(* Operator `*!` defines multiplication of a scalar (float) by a vector *)
let ( *! ) a v = Array.map ( fun x -> a *. x ) v
(* Operator `+:` define the sum of two vectors *)
let ( +: ) u v = Array.map2 ( +. ) u v
(* Operator `-:` is introduced for convenience/simplicity in the `distance_p2`, we could also define
it directly as `let ( -: ) u v = Array.map2 ( +. ) u ( -1.0 *! v )` *)
let ( -: ) u v = u +: (-1.0 *! v)
(* `metric_p2` is the Euclidean distance based on the p-norm with p=2, it is the most typical distance. *)
let metric_p2 u v = sqrt (( u -: v ) *@ ( u -: v) )
(* `distance_in_set` calculates all the distances in a set of elements `s` with respect to a reference
element `ref` using a `metric` function. It returns a tuple `int * float` for the index of the element
in the set and the distance of that element with respect to the reference. *)
let distance_in_set ref s metric =
let n = Array.length s in
Array.init n (fun i -> (i, metric ref (Array.get s i)) )
(* `sort_by_distance` will sort the output of the `distance_in_set` with in-place replacement, so the
indices of the k nearest neighbours will be the k first elements of the array returned by this function *)
let sort_by_distance (d:(int * float) array) =
let cmd x y =
if (snd x -. snd y) = 0.0 then 0
else if (snd x -. snd y) > 0.0 then 1
else -1
in
Array.sort cmd d
(* `get_nn s d` returns the nearest neighbour (first element of the sorted by distance set *)
let get_nn s d = let nn_idx = fst (Array.get d 0) in
(nn_idx, Array.get s nn_idx)
(* Supporting functions: we create a set of functions to sample in space points. The first is `sample_in_line`
which draws a point in a segment [x,y] for a uniform distribution (all points in a segment have the same
probability). We use the standard library, but there is a complete set of functions for that in Owl.Stats,
so in the case of `sample_in_line` we could use for example `Owl.Stats.uniform_rvs ~a:x ~b:y` *)
let sample_in_interval x y = (Random.float 1.0 ) *. (y -. x)
(* `sample_in_cube` draws an element in a cube (the position of the vector is within a cube), the element
(vector) will have 3 coordinates, so an array of length 3. *)
let sample_in_cube a b = [| (sample_in_interval a b) ; (sample_in_interval a b) ; (sample_in_interval a b) |]
(* `sample_set_in_cube` create n random elements of the space within a cube *)
let sample_set_in_cube n a b = Array.init n (fun _ -> (sample_in_cube a b))
(* Auxiliary functions: three functions to print our results. *)
let print_point v =
let s = Array.map (Printf.sprintf "%6f") v |> Array.to_list in
Printf.printf "[ %s ]\n" (String.concat ", " s)
let print_set s = Array.iter print_point s
let print_distance x =
let ind = Array.map (fst) x in
let dis = Array.map (snd) x in
let s = Array.map2 (Printf.sprintf "(index %d -> distance %6f)") ind dis |> Array.to_list in
Printf.printf "[ %s ]\n" (String.concat ", " s)
(* Testing the functions *)
let () = let ref = sample_in_cube (-10.0) 10.0 in
let _ = Printf.printf "Our reference point is:" in
let _ = print_point ref in
let cc = sample_set_in_cube 5 (-10.0) 10.0 in
let _ = Printf.printf "The set of points to evaluate is:\n" in
let _ = print_set cc in
let _ = Printf.printf "Their distances with respect to reference are:\n" in
let _ = print_distance (distance_in_set ref cc distance_p2) in
let ss = distance_in_set ref cc distance_p2 in
let _ = sort_by_distance ss in
let _ = Printf.printf "The sorted points by distance to reference are:\n" in
let _ = print_distance ss in
let _ = Printf.printf "Therefore the nearest neighbour is:\n" in
let nn = get_nn cc ss in
let _ = Printf.printf "Index: %d, Coordinates: " (fst nn) in
print_point (snd nn);;
Discusión
Es posible que la tarea más difícil sea entender las funciones Array.map, Array.map2 y Array.fold_left para definir nuestros operadores entre vectores, y la función Array.sort donde es clave poder crear una función que compare los segundos elementos de una tupla (tuple). Por lo demás, es sencillamente aplicar definiciones matemáticas de un modo directo. Hago hincapié en que hemos definido operadores inflix lo que posible y creo que sencillo en OCaml.