Convert MATLAB use of Probability Density Function (PDF) to Python

Hi All

After asking in StackOverflow question without getting any answer I'm trying my luck here...

I'm working to convert below MATLAB code to Python:

    data                        = [ 44374 44034 44150 44142 44332 43986 44423 44346 44199 44129 44173 43981 44152 43797 44167 43944 44069 44327 44083 44473 44530 44361 44067 44154 44212 44537 44158 44428 43911 44397];
    RMS                         = sqrt( data );
    MA                          = movmean( RMS , 15 );
    x                           = RMS - MA;
    y                           = linspace( -10 , 10 , 21 );
    f_x                         = pdf( x , y );

So far I have this Python code:

    
    from decimal import *
    import numpy as np
    from scipy.stats import gaussian_kde
    from numpy.lib.stride_tricks import sliding_window_view
    
    def movmean_edge_cases(arr, window_size):   
        # adding NaN for edges
        padding = (window_size - 1) // 2
        arr = np.pad(arr, (padding, padding), 'constant')
        # convert zeros to NaN
        arr[arr == 0] = np.nan
        sliding_windows = sliding_window_view(arr, window_shape=window_size)
        left_edge = np.nanmean(sliding_windows[:padding], axis=1)
        right_edge = np.nanmean(sliding_windows[-padding:], axis=1)
        return left_edge, right_edge
    
    
    def movmean(arr, window_size):
        # Create a windowed filter with equal weights
        weights = np.ones(window_size) / window_size
        
        # Compute the convolution between the signal and the filter
        mean_values = np.convolve(arr, weights, mode='valid')
        # Compute the edge values
        left_edge, right_edge = movmean_edge_cases(arr, window_size)
        # add edge values to the original list
        mean_values = np.insert(mean_values, 0, left_edge)
        mean_values = np.append(mean_values, right_edge)
        return mean_values
        
    data = [44374, 44034, 44150, 44142, 44332, 43986, 44423, 44346, 44199, 44129, 44173, 43981, 44152, 43797, 44167, 43944, 44069, 44327, 44083, 44473, 44530, 44361, 44067, 44154, 44212, 44537, 44158, 44428, 43911, 44397]
    getcontext().prec = 18  # Change the precision
    RMS = [Decimal(x).sqrt() for x in data]
    RMS = [float(x) for x in RMS]
    RMS = np.array(RMS)
    MA = movmean(RMS, 15)
    x = np.subtract(RMS, MA)
    y = np.linspace(-10, 10, 21)
    kde = gaussian_kde(x)
    f_x = kde.pdf(y)

I have implemented the movmean function to be the same as MATLAB. Comparing both code I have ensure that x values and y values are the same for both MATLAB and Python.

Unfortunately the results after running the pdf Probability density function are not equal...

Results in MATLAB:

    0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0666666666666667, 0.800000000000000, 0.133333333333333, 0, 0, 0, 0, 0, 0, 0, 0, 0

Results in Python:

    0, 0, 0, 0, 0, 0, 0, 0, 0.0000005849, 0.1135087105, 0.6936653409, 0.0836530196, 0.0000000072, 0, 0, 0, 0, 0, 0, 0, 0

I tried various combinations of bw_methods in gaussian_kde function without luck.

Thank you!

Answer 1

[Askic V]

I have made a quick test using Matlab's built in function pdf and Python's stats.norm.pdf

x = [-2 -1 0 1 2];
mu = 1;
sigma = 5;
y = pdf('Normal',x,mu,sigma);
vpa(y,7)

ans = 

The following is the Python implementation:

from scipy import stats
import numpy as np

x = np.array([-2, -1, 0, 1, 2])
mu = 1
sigma = 5
pdf_var = stats.norm.pdf(x, loc=mu, scale=sigma)
print(pdf_var)

And Python returns the following:

[0.06664492 0.07365403 0.07820854 0.07978846 0.07820854]

[David Recanati]

Thank you for you answer!

I tried this but my MATLAB code use:

f_x                         = pdf( x , y );

which doesn't specified if it use Noraml distribution or something else...

Do you know what is the defualt when using pdf function without probability distribution name?

[Askic V]

Trying to execute the code in this online version (2023a) gives and error:

data                        = [ 44374 44034 44150 44142 44332 43986 44423 44346 44199 44129 44173 43981 44152 43797 44167 43944 44069 44327 44083 44473 44530 44361 44067 44154 44212 44537 44158 44428 43911 44397];
RMS                         = sqrt( data );
MA                          = movmean( RMS , 15 );
x                           = RMS - MA;
y                           = linspace( -10 , 10 , 21 );
f_x                         = pdf( x , y );

Error using pdf
Requires the first input to be the name of a distribution.

[David Recanati]

I'm using Release R2017b sorry for not mentioned it

You can see this option here: https://www.mathworks.com/help/stats/prob.normaldistribution.pdf.html#:~:text=y%20%3D%20pdf(pd%2Cx)%20returns%20the%20pdf%20of%20the%20probability%20distribution%20object%20pd%2C%20evaluated%20at%20the%20values%20in%20x.