0
$\begingroup$

I want to program the chi-square between two distribution Xu and Vu with length N. This function is defined as below: Xu and Vu are real and positive.

chiDistance[U_, V_] := Module[{vecs, vselect}, (
(*delete the {0,0} occurrence to avoid the division by zero*)
   vecs = Join[Transpose[{U}], Transpose[{V}], 2];
   vselect = Select[vecs, # != {0., 0.} &];
   (1/2) Total[(vselect[[All, 1]] - 
         vselect[[All, 2]])^2/(vselect[[All, 1]] + 
        vselect[[All, 2]])]

   )]

This function take

0.047625 when N=40000

When I replaced Select[] by DeleteCases[], I got the same complexity. enter image description here

$\endgroup$
4
$\begingroup$

I am not sure whether the intention is vectors with all positive entries. If not then potential pairs (a,-a) will also be division by zero.

Here is another implementation of formula, removing zero denominators:

cd[u_, v_] := Module[{pos, us, vs},
  pos = Position[u + v, _?(# != 0 &)];
  us = Extract[u, pos];
  vs = Extract[v, pos];
  Total[(us - vs)^2/(us + vs)]/2]

or exploiting SparseArray properties:

cdsa[u_, v_] := Module[{pos, us, vs},
  pos = SparseArray[u + v]["NonzeroPositions"];
  us = Extract[u, pos];
  vs = Extract[v, pos];
  Total[(us - vs)^2/(us + vs)]/2]

Appears to confer only small advantage:

Needs["GeneralUtilities`"]
BenchmarkPlot[{cdsa @@ # &, 
  chiDistance @@ # &}, {RandomInteger[{0, 10}, #], 
   RandomInteger[{0, 10}, #]} &, {100, 1000, 10000, 40000}, 
 "IncludeFits" -> True]

enter image description here

$\endgroup$
  • $\begingroup$ using Total[N@(us - vs)^2/(us + vs)]/2 in the last line of cdsa seems to make it much faster (+1). $\endgroup$ – kglr Nov 8 '14 at 9:17
  • $\begingroup$ @kguler thank you, I should have tested on reals but done in haste... $\endgroup$ – ubpdqn Nov 8 '14 at 9:22
  • $\begingroup$ @ubpdqn, I edited my question. the two vectors are real and positive. So, I think that the use of Total[N@(us - vs)^2/(us + vs)]/2 does not change the complexity. $\endgroup$ – BetterEnglish Nov 8 '14 at 15:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.