Marginal Hilbert Spectrum, Does my integration scheme match the analytical expression? According to the original paper by Huang
https://arxiv.org/abs/1401.4211
The marginal Hilbert spectrum is given by:
$$h(\omega)=\int_0^\infty p(\omega,\mathcal A)\mathcal A^2\mathrm d\mathcal A$$
where $A = A(w,t)$ (i.e., a function time and frequency) and $p(w,A)$
the joint probability density function of $P(ω, A)$ of the frequency $ωi$ and amplitude $Ai$.
I am trying to estimate 1) The joint probability density using the plt.hist2d 2) the integral shown below using a sum.
The code I am using is the following:
IA_flat1  = np.ravel(IA)       ### Turn matrix to 1 D array
IF_flat1  = np.ravel(IF)       ### Here IA corresponds to A

IF_flat  = IF_flat1[(IF_flat1>min_f) & (IF_flat1<fs)]   ### Keep only desired frequencies
IA_flat  = IA_flat1[(IF_flat1>min_f) & (IF_flat1<fs)]   ### Keep IA that correspond to desired frequencies

### return the Joint probability density
Pjoint,f_edges, A_edges,_ = plt.hist2d(IF_flat,IA_flat,bins=[bins_F,bins_A], density=True)
plt.close()

n1 = np.digitize(IA_flat, A_edges).astype(int)          ### Return the indices of the bins to which  
n2 = np.digitize(IF_flat, f_edges).astype(int)          ### each value in input array belongs.

### define integration function
from numba import jit, prange   ### Numba is added for speed
@jit(nopython=True, parallel= True)
def get_int(A_edges, Pjoint ,IA_flat, n1, n2):
    
    dA = np.diff(A_edges)[0]                     ### Find dx for integration
    sum_h = np.zeros(np.shape(Pjoint)[0])        ### Intitalize array
    
    for j in prange(np.shape(Pjoint)[0]):
        h = np.zeros(np.shape(Pjoint)[1])        ### Intitalize array
        for k in prange(np.shape(Pjoint)[1]):
            needed = IA_flat[(n1==k) & (n2==j)]         ### Keep only the elements of arrat that 
                                                        ### are related to PJoint[j,k]   
                
            h[k]   = Pjoint[j,k]*np.nanmean(needed**2)*dA  ### Pjoint*A^2*dA
        sum_h[j]   = np.nansum(h)                       ### Sum_{i=0}^{N}(Pjoint*A^2*dA)
    return sum_h

### Now run previously defined function
sum_h = get_int(A_edges, Pjoint ,IA_flat, n1, n2)

1) I am not sure that everything is correct though. Any suggestions or comments on what I might be doing wrong?
2) Is there a way to do the same using a scipy integration scheme?
 A: 

*

*I am not sure that everything is correct though. Any suggestions or comments on what I might be doing wrong?


For numerical integration, merely summing up the integrand times a discretized dx value for each integrand element usually does not provide very accurate outputs.  The better thing to do is to either find or write your own numerical integration scheme.  I have found that Simpson's 1/3 rule works very well.  The wiki link above even has an example Python implementation (not sure how efficient it is so make sure to optimize as necessary).
I make this point as I do not know to what A_edges or f_edges or Pjoint refers (I can guess but I doubt it will change my answer if I knew).  If this is just an array of start/stop locations of a histogram, then my comment above is relevant.



*Is there a way to do the same using a scipy integration scheme?


I am not familiar with all of the built-in functions of scipy but if it has a Simpson's or trapezoid rule integration functions I would use those.
As an aside, I have vectorized examples of the Simpson's rule algorithm written in IDL if you want to see an example.
