{"action":"created","issue":{"url":"https://api.github.com/repos/tatome/dip_test/issues/1","labels_url":"https://api.github.com/repos/tatome/dip_test/issues/1/labels{/name}","comments_url":"https://api.github.com/repos/tatome/dip_test/issues/1/comments","events_url":"https://api.github.com/repos/tatome/dip_test/issues/1/events","html_url":"https://github.com/tatome/dip_test/issues/1","id":88901091,"number":1,"title":"Errors when data has duplicate values","user":{"login":"deto","id":551879,"avatar_url":"https://avatars.githubusercontent.com/u/551879?v=3","gravatar_id":"","url":"https://api.github.com/users/deto","html_url":"https://github.com/deto","followers_url":"https://api.github.com/users/deto/followers","following_url":"https://api.github.com/users/deto/following{/other_user}","gists_url":"https://api.github.com/users/deto/gists{/gist_id}","starred_url":"https://api.github.com/users/deto/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/deto/subscriptions","organizations_url":"https://api.github.com/users/deto/orgs","repos_url":"https://api.github.com/users/deto/repos","events_url":"https://api.github.com/users/deto/events{/privacy}","received_events_url":"https://api.github.com/users/deto/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":2,"created_at":"2015-06-17T03:29:02Z","updated_at":"2015-06-17T08:25:18Z","closed_at":null,"body":"I was excited to use this for a project I'm working on - it's much faster than the python implementation that was floating around.  I was hoping that someone would make a real python implementation like this - I was debating trying myself, but I found the original paper a bit difficult!\r\n\r\nHowever, in playing around with this, I noticed an error between this version and the other Python and R versions.  When the input data has repeating values, the dip statistic is not correct (or, at least, it is not the same as what the other Python and R versions of this code report).\r\n\r\nContrived example:  If your data vector is [1, 1, 1, 1, 2, 2, 3, 3, 3, 3]\r\n\r\nThe \"old\" Python version and the R diptest package report a statistic of 0.2, but this version reports just zero.\r\n\r\nUnfortunately, the data I\"m working with has many distributions with repeated values (usually repeated zeros).  I ran a comparison between this dip test and the old python version and found that they agreed on only around 700 out of 7000 sample sets.  However, for data vectors where the samples are all unique, the dip statistics always matched between tests, and for vectors with repeated samples, the dip statistics usually disagreed.  In your verification, you used random data, so I'm guessing you were only testing with unique sample sets.\r\n\r\nAny ideas what might be happening?  "},"comment":{"url":"https://api.github.com/repos/tatome/dip_test/issues/comments/112715359","html_url":"https://github.com/tatome/dip_test/issues/1#issuecomment-112715359","issue_url":"https://api.github.com/repos/tatome/dip_test/issues/1","id":112715359,"user":{"login":"tatome","id":6285646,"avatar_url":"https://avatars.githubusercontent.com/u/6285646?v=3","gravatar_id":"","url":"https://api.github.com/users/tatome","html_url":"https://github.com/tatome","followers_url":"https://api.github.com/users/tatome/followers","following_url":"https://api.github.com/users/tatome/following{/other_user}","gists_url":"https://api.github.com/users/tatome/gists{/gist_id}","starred_url":"https://api.github.com/users/tatome/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/tatome/subscriptions","organizations_url":"https://api.github.com/users/tatome/orgs","repos_url":"https://api.github.com/users/tatome/repos","events_url":"https://api.github.com/users/tatome/events{/privacy}","received_events_url":"https://api.github.com/users/tatome/received_events","type":"User","site_admin":false},"created_at":"2015-06-17T08:25:18Z","updated_at":"2015-06-17T08:25:18Z","body":"Hi David,\n\nI haven't really had time to look into the reason for the error too \ndeeply, but I've found that the small change in the attached version of \nthe dip implementation seems to fix it.  I haven't tested it rigorously, \nbut on a reasonably-sized random data set with repeated values it seems \nto compute the same values as the R implementation.\n\nThanks for letting me know about this apparent bug!  I'll do some more \ntesting and then update the repo.\n\nCheers,\nJohannes\n\nOn , David DeTomaso wrote:\n> I was excited to use this for a project I'm working on - it's much\n> faster than the python implementation that was floating around.  I was\n> hoping that someone would make a real python implementation like this\n> - I was debating trying myself, but I found the original paper a bit\n> difficult!\n> \n> However, in playing around with this, I noticed an error between this\n> version and the other Python and R versions.  When the input data has\n> repeating values, the dip statistic is not correct (or, at least, it\n> is not the same as what the other Python and R versions of this code\n> report).\n> \n> Contrived example:  If your data vector is [1, 1, 1, 1, 2, 2, 3, 3, 3, \n> 3]\n> \n> The \"old\" Python version and the R diptest package report a statistic\n> of 0.2, but this version reports just zero.\n> \n> Unfortunately, the data I\"m working with has many distributions with\n> repeated values (usually repeated zeros).  I ran a comparison between\n> this dip test and the old python version and found that they agreed on\n> only around 700 out of 7000 sample sets.  However, for data vectors\n> where the samples are all unique, the dip statistics always matched\n> between tests, and for vectors with repeated samples, the dip\n> statistics usually disagreed.  In your verification, you used random\n> data, so I'm guessing you were only testing with unique sample sets.\n> \n> Any ideas what might be happening?\n> \n> ---\n> Reply to this email directly or view it on GitHub:\n> https://github.com/tatome/dip_test/issues/1\n\"\"\"\n    Module for computing The Hartigans' dip statistic\n\n    The dip statistic measures unimodality of a sample from a random process.\n\n    See: \n    Hartigan, J. A.; Hartigan, P. M. The Dip Test of Unimodality. The Annals \n    of Statistics 13 (1985), no. 1, 70--84. doi:10.1214/aos/1176346577. \n    http://projecteuclid.org/euclid.aos/1176346577.\n\"\"\"\n\nimport numpy as np\nimport collections\n\ndef _gcm_(cdf, idxs):\n    work_cdf = cdf\n    work_idxs = idxs\n    gcm = [work_cdf[0]]\n    touchpoints = [0]\n    while len(work_cdf) > 1:\n        distances = work_idxs[1:] - work_idxs[0]\n        slopes = (work_cdf[1:] - work_cdf[0]) / distances\n        minslope = slopes.min()\n        minslope_idx = np.where(slopes == minslope)[0][0] + 1\n        gcm.extend(work_cdf[0] + distances[:minslope_idx] * minslope)\n        touchpoints.append(touchpoints[-1] + minslope_idx)\n        work_cdf = work_cdf[minslope_idx:]\n        work_idxs = work_idxs[minslope_idx:]\n    return np.array(np.array(gcm)),np.array(touchpoints)\n\ndef _lcm_(cdf, idxs):\n    g,t = _gcm_(1-cdf[::-1], idxs.max() - idxs[::-1])\n    return 1-g[::-1], len(cdf) - 1 - t[::-1]\n\ndef _touch_diffs_(part1, part2, touchpoints):\n    diff = np.abs((part2[touchpoints] - part1[touchpoints]))\n    return diff.max(), diff\n\ndef dip(histogram=None, idxs=None):\n    \"\"\"\n        Compute the Hartigans' dip statistic either for a histogram of\n        samples (with equidistant bins) or for a set of samples.\n    \"\"\"\n    if idxs is None:\n        idxs = np.arange(len(histogram))\n    elif histogram is None:\n        h = collections.Counter(idxs)\n        idxs = np.msort(h.keys())\n        histogram = np.array([h[i] for i in idxs])\n    else:\n        if len(histogram) != len(idxs):\n            raise ValueError(\"Need exactly as many indices as histogram bins.\")\n        if len(idxs) != len(set(idxs)):\n            raise ValueError(\"idxs must be unique if histogram is given.\")\n        if not np.array_equal(np.msort(idxs), idxs):\n            idxs_s = np.argsort(idxs)\n            idx = np.asarray(idxs)[idxs_s]\n            histogram = np.asarray(histogram)[idxs_s]\n\n    cdf = np.cumsum(histogram, dtype=float)\n    cdf /= cdf[-1]\n\n    work_idxs = idxs\n    work_histogram = np.asarray(histogram, dtype=float) / np.sum(histogram)\n    work_cdf = cdf\n\n    D = 0\n    left = [0]\n    right = [1]\n\n    while True:\n        left_part, left_touchpoints   = _gcm_(work_cdf - work_histogram, work_idxs)\n        right_part, right_touchpoints = _lcm_(work_cdf, work_idxs)\n\n        d_left, left_diffs   = _touch_diffs_(left_part, right_part, left_touchpoints)\n        d_right, right_diffs = _touch_diffs_(left_part, right_part, right_touchpoints)\n\n        if d_right > d_left:\n            xr = right_touchpoints[d_right == right_diffs][-1]\n            xl = left_touchpoints[left_touchpoints <= xr][-1]\n            d  = d_right\n        else:\n            xl = left_touchpoints[d_left == left_diffs][0]\n            xr = right_touchpoints[right_touchpoints >= xl][0]\n            d  = d_left\n\n        left_diff  = np.abs(left_part[:xl+1] - work_cdf[:xl+1]).max()\n        right_diff = np.abs(right_part[xr:]  - work_cdf[xr:] + work_histogram[xr:]).max()\n\n        if d <= D or xr == 0 or xl == len(work_cdf):\n            the_dip = max(np.abs(cdf[:len(left)] - left).max(), np.abs(cdf[-len(right)-1:-1] - right).max())\n            return the_dip/2, (cdf, idxs, left, left_part, right, right_part)\n        else:\n            D = max(D, left_diff, right_diff)\n\n        work_cdf = work_cdf[xl:xr+1]\n        work_idxs = work_idxs[xl:xr+1]\n        work_histogram = work_histogram[xl:xr+1]\n\n        left[len(left):] = left_part[1:xl+1]\n        right[:0] = right_part[xr:-1]\n\ndef plot_dip(histogram=None, idxs=None):\n    from matplotlib import pyplot as plt\n\n    d,(cdf,idxs,left,left_part,right,right_part) = dip(histogram,idxs)\n\n\n    plt.plot(idxs[:len(left)], left, color='red')\n    plt.plot(idxs[len(left)-1:len(left)+len(left_part) - 1], left_part, color='gray')\n    plt.plot(idxs[-len(right):], right, color='blue')\n    plt.plot(idxs[len(cdf) - len(right) + 1 - len(right_part):len(cdf) - len(right) + 1], right_part, color='gray')\n\n    the_dip = max(np.abs(cdf[:len(left)] - left).max(), np.abs(cdf[-len(right)-1:-1] - right).max())\n    l_dip_idxs = np.abs(cdf[:len(left)] - left) == the_dip\n    r_dip_idxs = np.abs(cdf[-len(right)-1:-1] - right) == the_dip\n    print(the_dip/2, d)\n\n    plt.vlines(x=idxs[:len(left)][l_dip_idxs], ymin=cdf[:len(left)][l_dip_idxs], ymax = cdf[:len(left)][l_dip_idxs] - the_dip)\n    plt.vlines(x=idxs[-len(right):][r_dip_idxs], ymin=cdf[-len(right)-1:-1][r_dip_idxs], ymax = cdf[-len(right)-1:][r_dip_idxs] + the_dip)\n\n    plt.plot(np.repeat(idxs,2)[1:], np.repeat(cdf,2)[:-1], color='black')\n    plt.scatter(idxs, cdf)\n\n    plt.show()\n\n\ndef crit_points(random_function, quantiles, sample_size, n_samples):\n    \"\"\"\n        Compute the quantiles for the dip statistic for n_samples \n        samples of size sample_size from the random process given by \n        random_function.\n\n        Parameters:\n        random_function : a paramter-free function which returns random values.\n        quantiles : a sequence of values between 0 and 1\n        sample_size : the size of the samples to draw from random_function\n        n_samples : the number of samples for which to compute dips\n\n        Returns: a list such that the i'th value is the greatest dip observed\n        such that the fraction of dips less than or equal to that value is less\n        than the i'th value from quantiles.\n    \"\"\"\n    data = [[random_function() for _ in range(sample_size)] for __ in range(n_samples)]\n    dips = np.array([dip(idxs=samples)[0] for samples in data])\n    \n    return np.percentile(dips, [p * 100 for p in quantiles])\n"}}