initM[list_List, k_Integer, distFunc_Symbol] :=
  Module[{m = {RandomChoice[list]}, n, d},
   While[Length[m] < k,
    n = RandomChoice@Nearest[m, #] & /@ list;
    d = Apply[distFunc, Transpose[{n, list}], {1}];
    m = Append[m, RandomChoice[d -> list]]
    ];
   m
   ];
kmeanspp[list_, k_,
   opts : OptionsPattern[{DistanceFunction ->
       SquaredEuclideanDistance, "RandomSeed" -> {}}]] :=
  BlockRandom[SeedRandom[OptionValue["RandomSeed"]];
   Module[{m = initM[list, k, OptionValue[DistanceFunction]], update,
     partition, clusters}, update[] := m = Mean /@ clusters;
    partition[_] := (clusters =
       GatherBy[list,
        RandomChoice@
          Nearest[m, #, (# -> OptionValue[#] &@DistanceFunction)] &];
      update[]);
    FixedPoint[partition, list];
    {clusters, m}
    ]
   ];