If you just want to estimate the parameters of a multivariate normal with the Iris data, then the following will work:

m = {m1, m2, m3, m4};
cov = {{v11, v12, v13, v14},
          {v12, v22, v23, v24},
          {v13, v23, v33, v34},
          {v14, v24, v34, v44}};
FindDistributionParameters[Values[data], MultinormalDistribution[m, cov]]
(* {m1 -> 5.84333, m2 -> 3.05733, m3 -> 3.758, m4 -> 1.19933, 
 v12 -> -0.0421511, v13 -> 1.26582, v14 -> 0.512829, v23 -> -0.327459,
  v24 -> -0.120828, v34 -> 1.28697, v11 -> 0.681122, v22 -> 0.188713, 
 v33 -> 3.0955, v44 -> 0.577133} *)