更新,以反映意見移植一個MATLAB代碼到Python
我需要端口MATLAB代碼爲蟒蛇。不幸的是,儘管在調試和谷歌搜尋方面有很多不眠之夜,但我無法讓我的代碼運行。這裏是MATLAB代碼問題:
clc;
serie =0:50;
serie = serie - mean(serie);
y = cumsum(serie);
L = length(y);
%Calculate the detrended fluctuation short term coefficient
npuntos = 10;
f1=zeros(1,npuntos);
for n1=4:16
%Segmentation
seg =(1:n1:L);
%disp(length(seg))
yn = zeros(1,L);
for k = 1:length(seg)-1
yaux = y(seg(k):seg(k+1)-1);
x = 1:length(yaux);
A=[sum(x.^2),sum(x); sum(x),length(x)];
C=[sum(x.*yaux);sum(yaux)];
v=inv(A)*C;
a=v(1); b=v(2);
pen=a;
ord=b;
ytrend = pen*x + ord;
yn(seg(k):seg(k+1)-1) = ytrend';
end
f1(n1) = sqrt((1/seg(end)).*sum((y(1:seg(end)-1)-yn(1:seg(end)-1)).^2));
end
n1=4:16;
f1=f1(4:end);
p1 = polyfit(log10(n1(1:end-2)),log10(f1(1:end-2)),1);
alpha1 = p1(1);
disp(alpha1)
我試圖將代碼轉換成蟒蛇去如下:
import numpy as np
data = np.arange(start=0, stop=51, step=1)
data = data.transpose()
data = data - np.mean(data)
y = np.cumsum(data)
L = len(y)
# Estimate the value of alpha1
npuntos = 12
f1 = [0] * npuntos
for i, n1 in enumerate(np.arange(start=4, stop=16, step=1)):
seg = np.arange(start=0, stop=L, step=n1) # Potential error
yn = [0] * L
for k in np.arange(start=0, stop=len(seg)-1, step=1): # Potential Error
yaux = y[seg[k]:seg[k + 1]-1] # Potential Error
x = np.arange(start=1, stop=len(yaux) + 1, step=1)
A = [[sum(x ** 2), sum(x)], [sum(x), len(x)]]
C = [[sum(x * yaux)], [sum(yaux)]]
v = (np.linalg.inv(A)).dot(C)
pen = v[0]
ord = v[1]
ytrend = pen * x + ord
yn[seg[k]: seg[k + 1] - 1] = ytrend
f1[i] = np.sqrt((1/seg[-1]) * sum((y[1:seg[-1] - 1] - yn[1:seg[-1] - 1]) ** 2))
n1 =np.arange(start=4, stop=16, step=1)
f1 = f1[4:]
xx =np.log10(n1[1: - 2])
yy=np.log10(f1[1: - 2])
print(len(xx))
print(len(yy))
p1 = np.polyfit(xx, yy, 1)
alpha1 = p1[1]
print(alpha1)
不幸的是,我得到類型錯誤在程序執行時,此行
p1 = np.polyfit(xx, yy, 1)
這確實有望爲xx是長度爲9,而xx是剛剛5.使用try/catch塊如評論所說,由輸出是完全錯誤的
try:
f1[n1] = np.sqrt((1/seg[-1]) * sum((y[1:seg[-1] - 1] - yn[1:seg[-1] - 1]) ** 2))
except IndexError:
f1.append(np.sqrt((1/seg[-1]) * sum((y[1:seg[-1] - 1] - yn[1:seg[-1] - 1]) ** 2)))
錯誤是固定的。
我已經通過調試器,但我不能如何錯誤。有人可以幫忙嗎? P.S-如果有人有興趣,上面的代碼應該計算出Detrended Fluctuation Analysis (DFA)。
這我都試過,但它給了我同樣的錯誤。你的意思是這樣的'f1 [n1] = np.sqrt((1/seg [-1])* sum((y [0:seg [-1] - 1] - yn [0:seg [-1 ] - 1])** 2))'?這也沒有用。我得到了相同的'IndexError:列表分配索引超出範圍' –
您正在使用哪個版本的python? – TheBlackCat
這是Python 3.5 –