Resampling based on an interpolated CDF, as described in Malik and Pitt. 

    Parameters
    ----------
    W: (N,) array
        weights 
    x: (N,) array 
        particles

    Returns
    -------
    xrs: (N,) array
        the resampled particles

def interpoled_resampling(W, x):
    """Resampling based on an interpolated CDF, as described in Malik and Pitt. 

    Parameters
    ----------
    W: (N,) array
        weights 
    x: (N,) array 
        particles

    Returns
    -------
    xrs: (N,) array
        the resampled particles 

    """
    N = W.shape[0] 
    idx = np.argsort(x)
    xs = x[idx] 
    ws = W[idx]
    cs = np.cumsum(avg_n_nplusone(ws)) 
    u = random.rand(N)
    xrs = np.empty(N)
    where = np.searchsorted(cs, u)
    # costs O(N log(N)) but algorithm has O(N log(N)) complexity anyway
    for n in range(N): 
        m = where[n]
        if m==0:
            xrs[n] = xs[0]
        elif m==N:
            xrs[n] = xs[-1]
        else:
            xrs[n] = interpol(cs[m-1], cs[m], xs[m-1], xs[m], u[n])
    return xrs