필터 지우기
필터 지우기

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

조회 수: 17 (최근 30일)
David Recanati
David Recanati 2023년 4월 13일
편집: David Recanati 2023년 4월 13일
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!

답변 (1개)

Askic V
Askic V 2023년 4월 13일
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
Askic V
Askic V 2023년 4월 13일
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 2023년 4월 13일
편집: David Recanati 2023년 4월 13일
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.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Call Python from MATLAB에 대해 자세히 알아보기

제품


릴리스

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by