2

According to the original paper by Huang

https://arxiv.org/pdf/1401.4211.pdf

The marginal Hibert spectrum is given by:

enter image description here

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?

4

1 回答 1

1

您可以从 2D 直方图中提取概率并将其用于积分：

# Added some numbers to have something to run
import numpy as np
import matplotlib.pyplot as plt
IA = np.random.rand(100,100)
IF = np.random.rand(100,100)
bins_F = np.linspace(0,1,20)
bins_A = np.linspace(0,1,100)
min_f = 0
fs = 1.0


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)

f_values = (f_edges[1:]+f_edges[:-1])/2 
A_values = (A_edges[1:]+A_edges[:-1])/2 
dA = A_values[1]-A_values[0] # for the integral

#Pjoint.shape (19,99)

h = np.zeros(f_values.shape)
for i in range(len(f_values)):
    f = f_values[i]

    # column of the histogram with frequency f, probability
    p = Pjoint[i]  

    # summatory equivalent to the integral 
    integral_result = np.sum(p*A_values**2*dA ) 
    h[i] = integral_result

plt.figure()    
plt.plot(f_values,h)
于 2021-10-06T16:54:37.647 回答