This is a recursive implementation of evaluating PDF of a Cumulative Density Function (CDF) at a point `u`, in arbitrary dimension, using finite difference methods.

multiderivative.py

def cdf2pdf( f, u, delta=0.001, kwargs={} ):
    """
    Computes every partial derivative of the CDF to find the underlying PDF.

    f: a numpy vectorized python/numpy function
    u: the point to evaluate the CDF at.
    delta: the precision
    kwargs: additional keywords args to f
    """

    def _wrapper(*args):
        u = np.array(args)
        return f(u, **kwargs)
    n = u.shape[0]
    p= _pdf( _wrapper, u, delta)
    return np.exp( np.log(p) - n*np.log( delta ) )
    #return p / delta**n
    

def _pdf(f, u, delta = 0.001 ):
    n = u.shape[0]
    if n==1:
        t= f(u[0]+delta/2) - f(u[0]-delta/2)
        return t
    else:
        f_plus = lambda *x: f( u[0] + delta/2, *x)
        f_minus = lambda *x: f( u[0] - delta/2, *x)
        return _pdf(f_plus, u[1:], delta ) - _pdf(f_minus, u[1:], delta ) 


if __name__=="__main__":

    def F( u ):
      return (u**2).sum()**2
  
    u = np.array( [0.4, 0.4] )
    print cdf2pdf( F, u )