|
Python |
Matlab/Octave |
| Vorbereitung (Laden von Modulen/Paketen) |
#python -m pip install --upgrade numpy
#python -m pip install --upgrade matplotlib
#pip install --upgrade numpy
#pip install --upgrade matplotlib
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
|
%Matlab: erfordert die Statistics and Machine Learning Toolbox
%Octave: pkg load statistics
if exist('OCTAVE_VERSION','builtin')~=0
pkg load statistics
end
|
| Parametriesierung eines AR4-Model |
p=4 #Prozessordnung (AR4)
phi=array([-0.7,0.2,0.4,0.3]) #4 AR-Koeffizienten
theta0=0.5 #1 MA-Koeffizient für Eingangsrauschen
dt=1.0
|
p=4; %Prozessordnung (AR4)
phi=[-0.7,0.2,0.4,0.3]; %4 AR-Koeffizienten
theta0=0.5; %1 MA-Koeffizient für Eingangsrauschen
dt=1;
|
| Kovarianzfunktion des simulierten Prozesses |
from numpy.linalg import *
p=len(phi)
PHI=zeros([p,p])
for i in range(0,p):
for j in range(0,p):
PHI[i,j]=0
if (i+j+1>=0) and (i+j+1<p):
PHI[i,j]+=phi[i+j+1]
if (i-j-1>=0) and (i-j-1<p):
PHI[i,j]+=phi[i-j-1]
if i==j:
PHI[i,j]-=1
eq=zeros(p)
for i in range(0,p):
eq[i]=-phi[i]
rho=append([1],linalg.solve(PHI,eq))
v=theta0**2/float(1-sum(phi*rho[1:p+1]))
K=100
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
C=zeros(K)
for k in range(maximum(-(K//2),-p),minimum((K+1)//2,p+1)):
C[k]=v*rho[abs(k)]
for k in range(minimum((K+1)//2,p+1),(K+1)//2):
C[k]=0
for j in range(0,p):
C[k]+=phi[j]*C[k-(j+1)]
for k in range(maximum(-(K//2),-p-1),-(K//2)-1,-1):
C[k]=0
for j in range(0,p):
C[k]+=phi[j]*C[k+j+1]
Csim=copy(C)
plot(tau,Csim,'.')
show()
|
p=length(phi);
PHI=zeros(p,p);
for i=1:p
for j=1:p
PHI(i,j)=0;
if (i+j>0) && (i+j<=p)
PHI(i,j)=PHI(i,j)+phi(i+j);
end
if (i-j>0) && (i-j<=p)
PHI(i,j)=PHI(i,j)+phi(i-j);
end
if i==j
PHI(i,j)=PHI(i,j)-1;
end
end
end
eq=zeros(1,p);
for i=1:p
eq(i)=-phi(i);
end
rho=[1,linsolve(PHI,eq')'];
v=theta0^2/(1-sum(phi.*rho(2:p+1)));
K=100;
tau=circshift((-fix(K/2):K-1-fix(K/2))*dt,[0;fix((K+1)/2)]);
C=zeros(1,K);
for k=max(-fix(K/2),-p):min(K-1-fix(K/2),p)
C(mod(k,K)+1)=v*rho(abs(k)+1);
end
for k=min(K-fix(K/2),p+1):K-1-fix(K/2)
C(mod(k,K)+1)=0;
for j=1:p
C(mod(k,K)+1)=C(mod(k,K)+1)+phi(j)*C(mod(k-j,K)+1);
end
end
for k=max(-fix(K/2),-p)-1:-1:-fix(K/2)
C(mod(k,K)+1)=0;
for j=1:p
C(mod(k,K)+1)=C(mod(k,K)+1)+phi(j)*C(mod(k+j,K)+1);
end
end
Csim=C;
plot(tau,Csim,'.')
|
| Leistungsdichte des simulierten Prozesses |
K=1000
p=len(phi)
df=1/float(K*dt)
f=roll(arange(-(K//2),(K+1)//2)*df,(K+1)//2)
S=zeros(K)
for i in range(0,len(f)):
S[i]=dt*theta0**2/abs(1-sum(phi*exp(-2j*pi*arange(1,p+1)*f[i]*dt)))**2
Ssim=copy(S)
plot(f,Ssim,'.')
show()
|
K=1000;
p=length(phi);
df=1/(K*dt);
f=circshift((-fix(K/2):fix((K-1)/2))*df,[0;fix((K+1)/2)]);
S=zeros(1,K);
for i=(1:length(f))
S(i)=dt*theta0^2/abs(1-sum(phi.*exp(-2j*pi*(1:p)*f(i)*dt)))^2;
end
Ssim=S;
plot(f,Ssim,'.')
|
| Generieren eines Signals (inklusive der Initialisierung) |
from numpy.linalg import *
PHI=zeros([p,p])
for i in range(0,p):
for j in range(0,p):
PHI[i,j]=0
if (i+j+1>=0) and (i+j+1<p):
PHI[i,j]+=phi[i+j+1]
if (i-j-1>=0) and (i-j-1<p):
PHI[i,j]+=phi[i-j-1]
if i==j:
PHI[i,j]-=1
eq=zeros(p)
for i in range(0,p):
eq[i]=-phi[i]
rho=append([1],linalg.solve(PHI,eq))
RHO=zeros([p,p])
for i in range(0,p):
for j in range(0,p):
RHO[i,j]=rho[abs(i-j)];
A=zeros([p,p])
for i in range(0,p):
A[i,i]=1.0
bl=inf
B=(dot(inv(A.T),RHO)-A)/2.0
b=amax(abs(B))
A=A+B
while b<bl:
bl=b
B=(dot(inv(A.T),RHO)-A)/2.0
b=amax(abs(B))
A=A+B
e=standard_normal(p)
v=theta0**2/float(1-sum(phi*rho[1:p+1]))
u=sqrt(v)*dot(A,e.T)
N=1000
x=zeros(N)
for i in range(0,N):
u=append(sum(u*phi)+normal(0,theta0),u[0:p-1])
x[i]=u[0]
t=arange(0,N)*dt
plot(t,x,'.')
show()
|
PHI=zeros(p,p);
for i=1:p
for j=1:p
PHI(i,j)=0;
if (i+j>0) && (i+j<=p)
PHI(i,j)=PHI(i,j)+phi(i+j);
end
if (i-j>0) && (i-j<=p)
PHI(i,j)=PHI(i,j)+phi(i-j);
end
if i==j
PHI(i,j)=PHI(i,j)-1;
end
end
end
eq=zeros(1,p);
for i=1:p
eq(i)=-phi(i);
end
rho=[1,linsolve(PHI,eq')'];
RHO=zeros(p,p);
for i=1:p
for j=1:p
RHO(i,j)=rho(abs(i-j)+1);
end
end
A=sqrtm(RHO);
e=randn(1,p);
v=theta0^2/(1-sum(phi.*rho(2:p+1)));
u=(sqrt(v)*A*e')';
N=1000;
x=zeros(1,N);
for i=1:N
u=[sum(u.*phi)+normrnd(0,theta0),u(1:p-1)];
x(i)=u(1);
end
t=(0:N-1)*dt;
plot(t,x,'.')
|
| Schätzung der Korrelationskoeffizientenfunktion und der Varianz |
from numpy.fft import *
v=var(x)
X=fft(append(x-mean(x),zeros(N))) # Mittelwert abgezogen, da Ergebnis sehr empfindlich
E=dt**2*abs(X)**2
CE=real(ifft(E))/float(dt)
K=2*p+1
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
rho=zeros(K)
for k in range(-(K//2),(K+1)//2):
rho[k]=CE[k]/float(dt)/float(N-abs(k))
rho/=rho[0]
|
v=var(x,1);
X=fft([x-mean(x),zeros(1,N)]); % Mittelwert abgezogen, da Ergebnis sehr empfindlich
E=dt^2*abs(X).^2;
CE=real(ifft(E))/dt;
K=2*p+1;
tau=circshift((-fix(K/2):fix((K-1)/2))*dt,fix((K+1)/2));
rho=zeros(1,K);
for k=-fix(K/2):fix(K-1)/2
rho(mod(k,K)+1)=CE(mod(k,2*N)+1)/dt/(N-abs(k));
end
rho=rho/rho(1);
|
| Bestimmung der Parameter des AR-Prozesses mit dem Yule-Walker-Verfahren |
from numpy.linalg import *
v*=rho[0] # falls rho[0]!=1
rho/=float(rho[0]) # rho[0]=1
RHO=zeros([p,p])
for i in range(0,p):
for j in range(0,p):
RHO[i,j]=rho[abs(i-j)]
eq=zeros(p)
for i in range(0,p):
eq[i]=rho[i+1]
phi=real(linalg.solve(RHO,eq))
theta0=real(sqrt(v*(1-sum(phi*rho[1:p+1]))))
print(phi)
print(theta0)
|
v=v*rho(1); % falls rho(1)~=1
rho=rho/rho(1); % rho(1)=1
RHO=zeros(p,p);
for i=1:p
for j=1:p
RHO(i,j)=rho(abs(i-j)+1);
end
end
eq=zeros(1,p);
for i=1:p
eq(i)=rho(i+1);
end
phi=real(linsolve(RHO,eq')')
theta0=real(sqrt(v*(1-sum(phi.*rho(2:p+1)))))
|
| Kovarianzfunktion des geschätzten Prozesses |
from numpy.linalg import *
p=len(phi)
PHI=zeros(p,p)
for i in range(0,p):
for j in range(0,p):
PHI[i,j]=0
if (i+j+1>=0) and (i+j+1<p):
PHI[i,j]+=phi[i+j+1]
if (i-j-1>=0) and (i-j-1<p):
PHI[i,j]+=phi[i-j-1]
if i==j:
PHI[i,j]-=1
eq=zeros(p)
for i in range(0,p):
eq[i]=-phi[i]
rho=append([1],linalg.solve(PHI,eq))
v=theta0**2/float(1-sum(phi*rho[1:p+1]))
K=100
tau=roll(arange(-(K//2),(K+1)//2)*dt,(K+1)//2)
C=zeros(K)
for k in range(maximum(-(K//2),-p),minimum((K+1)//2,p+1)):
C[k]=v*rho[abs(k)]
for k in range(minimum((K+1)//2,p+1),(K+1)//2):
C[k]=0
for j in range(0,p):
C[k]+=phi[j]*C[k-(j+1)]
for k in range(maximum(-(K//2),-p-1),-(K//2)-1,-1):
C[k]=0
for j in range(0,p):
C[k]+=phi[j]*C[k+j+1]
Cest=copy(C)
plot(tau,Cest,'.',tau,Csim,'.')
show()
|
p=length(phi);
PHI=zeros(p,p);
for i=1:p
for j=1:p
PHI(i,j)=0;
if (i+j>0) && (i+j<=p)
PHI(i,j)=PHI(i,j)+phi(i+j);
end
if (i-j>0) && (i-j<=p)
PHI(i,j)=PHI(i,j)+phi(i-j);
end
if i==j
PHI(i,j)=PHI(i,j)-1;
end
end
end
eq=zeros(1,p);
for i=1:p
eq(i)=-phi(i);
end
rho=[1,linsolve(PHI,eq')'];
v=theta0^2/(1-sum(phi.*rho(2:p+1)));
K=100;
tau=circshift((-fix(K/2):K-1-fix(K/2))*dt,[0;fix((K+1)/2)]);
C=zeros(1,K);
for k=max(-fix(K/2),-p):min(K-1-fix(K/2),p)
C(mod(k,K)+1)=v*rho(abs(k)+1);
end
for k=min(K-fix(K/2),p+1):K-1-fix(K/2)
C(mod(k,K)+1)=0;
for j=1:p
C(mod(k,K)+1)=C(mod(k,K)+1)+phi(j)*C(mod(k-j,K)+1);
end
end
for k=max(-fix(K/2),-p)-1:-1:-fix(K/2)
C(mod(k,K)+1)=0;
for j=1:p
C(mod(k,K)+1)=C(mod(k,K)+1)+phi(j)*C(mod(k+j,K)+1);
end
end
Cest=C;
plot(tau,Cest,'.',tau,Csim,'.')
|
| Leistungsdichte des geschätzten Prozesses |
K=1000
p=len(phi)
df=1/float(K*dt)
f=roll(arange(-(K//2),(K+1)//2)*df,(K+1)//2)
S=zeros(K)
for i in range(0,len(f)):
S[i]=dt*theta0**2/abs(1-sum(phi*exp(-2j*pi*arange(1,p+1)*f[i]*dt)))**2
Sest=copy(S)
plot(f,Sest,'.',f,Ssim,'.')
show()
|
K=1000;
p=length(phi);
df=1/(K*dt);
f=circshift((-fix(K/2):fix((K-1)/2))*df,[0;fix((K+1)/2)]);
S=zeros(1,K);
for i=(1:length(f))
S(i)=dt*theta0^2/abs(1-sum(phi.*exp(-2j*pi*(1:p)*f(i)*dt)))^2;
end
Sest=S;
plot(f,Sest,'.',f,Ssim,'.')
|
| vorgefertigtes Yule-Walker-Verfahren |
(aus der Spectrum Toolbox)
#python -m pip install --upgrade spectrum
#pip install --upgrade spectrum
import spectrum
a,ve,refl=spectrum.aryule(x,p) #a-Koeffizienten von gebrochenrationaler Darstellung, Varianz des Eingangsrauschens
phi=-a
theta0=sqrt(ve)
print(phi)
print(theta0)
|
%Matlab: erfordert die Statistics and Machine Learning Toolbox
%Octave: pkg load statistics
% pkg load signal
if exist('OCTAVE_VERSION','builtin')~=0
pkg load statistics
pkg load signal
end
[a,ve,k]=aryule(x,p); %a-Koeffizienten von gebrochenrationaler Darstellung, Varianz des Eingangsrauschens
phi=-a(2:end)
theta0=sqrt(ve)
|
| vorgefertigtes Burg-Verfahren |
(aus der Spectrum Toolbox)
#python -m pip install --upgrade spectrum
#pip install --upgrade spectrum
import spectrum
a,ve,refl=spectrum.arburg(x,p) #a-Koeffizienten von gebrochenrationaler Darstellung, Varianz des Eingangsrauschens
phi=-real(a)
theta0=sqrt(ve)
print(phi)
print(theta0)
|
%Matlab: erfordert die Statistics and Machine Learning Toolbox
%Octave: pkg load statistics
% pkg load signal
if exist('OCTAVE_VERSION','builtin')~=0
pkg load statistics
pkg load signal
end
[a,ve,k]=arburg(x,p); %a-Koeffizienten von gebrochenrationaler Darstellung, Varianz des Eingangsrauschens
phi=-a(2:end)
theta0=sqrt(ve)
|
| vorgefertigte Kovarianzmethode |
(aus der Spectrum Toolbox)
#python -m pip install --upgrade spectrum
#pip install --upgrade spectrum
import spectrum
a,err=spectrum.arcovar(x,p) #a-Koeffizienten von gebrochenrationaler Darstellung, summierter Prädiktionsfehler
phi=-a
theta0=sqrt(err/float(len(x)))
print(phi)
print(theta0)
|
Matlab:
[a,ve]=arcov(x,p); %a-Koeffizienten von gebrochenrationaler Darstellung, Varianz des Eingangsrauschens
phi=-a(2:end)
theta0=sqrt(ve)
Octave:
fehlt noch (Funktion arcov befindet sich im Review) |
| Vergleich der Prädiktion der Kovarianzmethode mit dem tatsächlichen Datensatz |
xc=zeros(len(x))
for i in range(1,len(x)):
for j in range(0,minimum(p,i)):
xc[i]+=x[i-j-1]*phi[j]
print(sum((x-xc)**2)/float(len(x)))
plot(t,x,'.',t,xc,'.')
show()
|
xc=zeros(1,length(x));
for i=2:length(x)
for j=1:min(p,i-1)
xc(i)=xc(i)+x(i-j)*phi(j);
end
end
sum((x-xc).^2)/length(x)
plot(t,x,'.',t,xc,'.')
|
| vorgefertigte, modifizierte Kovarianzmethode |
(aus der Spectrum Toolbox)
#python -m pip install --upgrade spectrum
#pip install --upgrade spectrum
import spectrum
a,err=spectrum.modcovar(x,p) #a-Koeffizienten von gebrochenrationaler Darstellung, summierter Prädiktionsfehler
phi=-a
theta0=sqrt(err/float(2*len(x)))
print(phi)
print(theta0)
|
Matlab:
[a,ve]=armcov(x,p); %a-Koeffizienten von gebrochenrationaler Darstellung, Varianz des Eingangsrauschens
phi=-a(2:end)
theta0=sqrt(ve)
Octave:
fehlt noch (Funktion armcov befindet sich im Review) |
| Vergleich der Prädiktion der modifizierten Kovarianzmethode mit dem tatsächlichen Datensatz |
xmc=zeros(len(x))
for i in range(1,len(x)):
for j in range(0,minimum(p,i)):
xmc[i]+=x[i-j-1]*phi[j]
print(sum((x-xmc)**2)/float(len(x)))
plot(t,x,'.',t,xc,'.',t,xmc,'.')
show()
|
xmc=zeros(1,length(x));
for i=2:length(x)
for j=1:min(p,i-1)
xmc(i)=xmc(i)+x(i-j)*phi(j);
end
end
sum((x-xmc).^2)/length(x)
plot(t,x,'.',t,xc,'.',t,xmc,'.')
|
| Akaike und Bayesian Information Criterion für unterschiedliche AR-Modellordnungen mit der modifizierten Kovarianzmethode zur Prozessidentifikation |
#python -m pip install --upgrade spectrum
#pip install --upgrade spectrum
import spectrum
pest=arange(1,11)
aics=zeros(len(pest))
bics=zeros(len(pest))
for k in range(len(pest)):
a,err=spectrum.modcovar(x,pest[k]) #a-Koeffizienten von gebrochenrationaler Darstellung, summierter Praediktionsfehler
phi=-a
theta0=sqrt(err/float(2*len(x)))
xmc=zeros(len(x))
for i in range(1,len(x)):
for j in range(0,minimum(pest[k],i)):
xmc[i]+=x[i-j-1]*phi[j]
sig2=sum((x-xmc)**2)/float(len(x))
aics[k]=2*pest[k]+N*log(sig2)
bics[k]=log(N)*pest[k]+N*log(sig2)
plot(pest,aics,'o',pest,bics,'o')
show()
|
Matlab:
pest=[1:10];
aics=zeros(1,length(pest));
bics=zeros(1,length(pest));
for k=1:length(pest)
%[a,ve]=armcov(x,pest(k)); %a-Koeffizienten von gebrochenrationaler Darstellung, Varianz des Eingangsrauschens
[a,ve,r]=aryule(x,pest(k)); %a-Koeffizienten von gebrochenrationaler Darstellung, Varianz des Eingangsrauschens
phi=-a(2:end)
theta0=sqrt(ve)
xmc=zeros(1,length(x));
for i=2:length(x)
for j=1:min(pest(k),i-1)
xmc(i)=xmc(i)+x(i-j)*phi(j);
end
end
sig2=sum((x-xmc).^2)/length(x)
aics(k)=2*pest(k)+N*log(sig2);
bics(k)=log(N)*pest(k)+N*log(sig2);
end
plot(pest,aics,'o',pest,bics,'o')
Octave:
fehlt noch (Funktion armcov befindet sich im Review) |