Clear Filters
Clear Filters

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

13 views (last 30 days)
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!

Answers (1)

Askic V
Askic V on 13 Apr 2023
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]
  3 Comments
Askic V
Askic V on 13 Apr 2023
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
David Recanati on 13 Apr 2023
Edited: David Recanati on 13 Apr 2023
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.

Sign in to comment.

Products


Release

R2017b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!