F stat computation for inhouse LMM

lib/toolbox/inhouse_lmm.ml

@@ -30,7 +30,24 @@ let correlations (t : Convergence_tree.u) : (string * string * float) list
 
 let design_matrix ~m ~n ~aa_at_site (al : Alignment.t) =
   L.Mat.init_rows m n (fun i j ->
-      if Char.(al.sequences.(i).[j] = aa_at_site.(j)) then 1. else 0.)
+      if j = 0 then 1.
+      else if Char.(al.sequences.(i).[j - 1] = aa_at_site.(j - 1)) then 1.
+      else 0.)
+
+let predict_y ~_X_ ~theta = L.gemv _X_ theta
+
+let squares_sum ~y_r ~y_q = L.Vec.sub y_r y_q |> L.Vec.sqr_nrm2
+
+let f_stat ~y ~_X_r ~theta_r ~_X_q ~theta_q =
+  let y_r = predict_y ~_X_:_X_r ~theta:theta_r in
+  let y_q = predict_y ~_X_:_X_q ~theta:theta_q in
+  let scm_rq = squares_sum ~y_r ~y_q in
+  let scr_r = squares_sum ~y_r ~y_q:y in
+  (* assuming X has full rank since we are solving using LM equations *)
+  let rank_r = Float.of_int (L.Mat.dim2 _X_r) in
+  let rank_q = 1. in
+  let n = Float.of_int (L.Vec.dim y) in
+  scm_rq /. (rank_r -. rank_q) /. (scr_r /. (n -. rank_r))
 
 let solve ~y ~_X_ ~_T_ =
   let _Xtilde_ = L.Mat.mul _X_ _T_ in