|
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
|
| Vorgabe einer Korrelationskoeffizientenfunktion (einseitig) |
dt=1.0
K=50
Kc=20
rho=zeros(K)
for k in range(0,Kc):
rho[k]=(1-k/float(Kc))*cos(1.5*pi*k*dt/float(Kc)) # rho[0]=1
plot(arange(0,K)*dt,rho,'o')
show()
|
dt=1;
K=50;
Kc=20;
rho=zeros(1,K);
for k=1:Kc
rho(k)=(1-(k-1)/Kc)*cos(1.5*pi*(k-1)*dt/Kc); % rho(1)=1
end
plot((0:K-1)*dt,rho,'o')
|
| Entwurf der Koeffizienten eines ARp-Prozesses aus der vorgegebenen Korrelationskoeffizientenfunktion sowie der Varianz |
from numpy.linalg import *
rho/=float(rho[0]) # rho[0]=1
v=1.0 # Varianz
p=len(rho)-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]))))
plot(phi,'o')
show()
|
rho=rho/rho(1); % rho(1)=1
v=1; % Varianz
p=length(rho)-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=linsolve(RHO,eq')';
theta0=sqrt(v*(1-sum(phi.*rho(2:p+1))));
plot(phi,'o')
|
| Kovarianzfunktion des entworfenen Prozesses (Erweiterung der Vorgabe) |
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=200
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]
Cpro=copy(C)
rhopro=C/float(C[0])
plot(tau,C,'.')
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=200;
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
Cpro=C;
rhopro=C/C(1);
plot(tau,C,'.')
|
| Vorbereitung der Realisierungen |
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))
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
v=theta0**2/float(1-sum(phi*rho[1:p+1]))
|
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')'];
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);
v=theta0^2/(1-sum(phi.*rho(2:p+1)));
|
| zwei Realisierungen (inklusive der Initialisierung) |
M=2
N=100
x=zeros([M,N])
for j in range(0,M):
e=standard_normal(p)
u=sqrt(v)*dot(A,e.T)
for i in range(0,N):
u=append(sum(u*phi)+normal(0,theta0),u[0:p-1])
x[j,i]=u[0]
t=arange(0,N)*dt
plot(t,x[0],'o',t,x[1],'o')
show()
|
M=2;
N=100;
x=zeros(M,N);
for j=1:M
e=randn(1,p);
u=(sqrt(v)*A*e')';
for i=1:N
u=[sum(u.*phi)+normrnd(0,theta0),u(1:p-1)];
x(j,i)=u(1);
end
end
t=(1:N)*dt;
plot(t,x(1,:),'o',t,x(2,:),'o')
|
| Realisierungen, Verarbeitung zu einer zeitlich aufgelösten Wahrscheinlichkeitsdichte, Mittelwert, empirische Varianz und Prädiktionen |
M=1000
N=100
t=arange(0,N)*dt
x=zeros([M,N])
Kk=25
xk=0.2*arange(0,Kk)-2.4
f=zeros([Kk-1,N])
ex=zeros([3,N])
mx=zeros([3,N])
for j in range(0,M):
e=standard_normal(p)
u=sqrt(v)*dot(A,e.T)
for i in range(0,N):
u=append(sum(u*phi)+normal(0,theta0),u[0:p-1])
x[j,i]=u[0]
for i in range(0,N):
f[:,i]=histogram(x[0:j+1,i],bins=xk)[0][0:Kk-1]/float(len(x[0:j+1,i]))/(xk[1:Kk]-xk[0:Kk-1])
clf()
imshow(clip(f,0.0,1.0),cmap='gray',extent=[-0.5,N-0.5,-2.4,2.4],aspect='auto',origin='lower')
for i in arange(0,N):
ex[0,i]=0.0 # Erwartungswert
ex[1,i]=ex[0,i]-sqrt(v) # Erwartungswert-Standardabweichung
ex[2,i]=ex[0,i]+sqrt(v) # Erwartungswert+Standardabweichung
mx[0,i]=mean(x[0:j+1,i]) # Mittelwert
if j>0:
vx=var(x[0:j+1,i],ddof=1)
mx[1,i]=mx[0,i]-sqrt(vx) # Mittelwert-Standardabweichung
mx[2,i]=mx[0,i]+sqrt(vx) # Mittelwert+Standardabweichung
plot(t,ex[0],'blue',t,ex[1],'red',t,ex[2],'red',mx[0],'cyan',t,mx[1],'orange',t,mx[2],'orange')
draw()
pause(0.01)
show()
|
M=1000;
N=100;
t=(1:N)*dt;
x=zeros(M,N);
Kk=25;
xk=0.2*(0:Kk-1)-2.4;
f=zeros(Kk-1,N);
ex=zeros(3,N);
mx=zeros(3,N);
for j=1:M
e=randn(1,p);
u=(sqrt(v)*A*e')';
for i=1:N
u=[sum(u.*phi)+normrnd(0,theta0),u(1:p-1)];
x(j,i)=u(1);
end
for i=1:N
fx=histc(x(1:j,i),xk)/length(x(1:j,i));
f(:,i)=fx(1:length(xk)-1);
end
imshow(imresize(mat2gray(min(1,f)),10))
set(gca,'YDir','normal')
hold on
for i=1:N
ex(1,i)=0.0; % Erwartungswert
ex(2,i)=ex(1,i)-sqrt(v); % Erwartungswert-Standardabweichung
ex(3,i)=ex(1,i)+sqrt(v); % Erwartungswert+Standardabweichung
mx(1,i)=mean(x(1:j,i)); % Mittelwert
if j>1
vx=var(x(1:j,i));
mx(2,i)=mx(1,i)-sqrt(vx); % Mittelwert-Standardabweichung
mx(3,i)=mx(1,i)+sqrt(vx); % Mittelwert+Standardabweichung
end
end
ti=10*(t-0.5*dt);
exi=10*(Kk-1)/4.8*(ex+2.4)+0.5;
mxi=10*(Kk-1)/4.8*(mx+2.4)+0.5;
plot(ti,exi(1,:),'blue',ti,exi(2,:),'red',ti,exi(3,:),'red',ti,mxi(1,:),'cyan',ti,mxi(2,:),'yellow',ti,mxi(3,:),'yellow')
hold off
pause(0.01)
end
|
| zwei Realisierungen (inklusive der Initialisierung) mit einer Beobachtung |
M=2
N=100
Ne=40
xe=1.0
x=zeros([M,N])
for j in range(0,M):
x[j,Ne]=-inf
while abs(x[j,Ne]-xe)>0.1:
e=standard_normal(p)
u=sqrt(v)*dot(A,e.T)
for i in range(0,N):
u=append(sum(u*phi)+normal(0,theta0),u[0:p-1])
x[j,i]=u[0]
t=arange(0,N)*dt
plot(t,x[0],'o',t,x[1],'o')
show()
|
M=2;
N=100;
Ne=40;
xe=1;
x=zeros(M,N);
for j=1:M
x(j,Ne)=-inf;
while abs(x(j,Ne)-xe)>0.1
e=randn(1,p);
u=(sqrt(v)*A*e')';
for i=1:N
u=[sum(u.*phi)+normrnd(0,theta0),u(1:p-1)];
x(j,i)=u(1);
end
end
end
t=(1:N)*dt;
plot(t,x(1,:),'o',t,x(2,:),'o')
|
| Realisierungen, Verarbeitung zu einer zeitlich aufgelösten Wahrscheinlichkeitsdichte, Mittelwert, empirische Varianz und Prädiktionen |
M=1000
N=100
Ne=40
xe=1.0
t=arange(0,N)*dt
x=zeros([M,N])
Kk=25
xk=0.2*arange(0,Kk)-2.4
f=zeros([Kk-1,N])
ex=zeros([3,N])
mx=zeros([3,N])
for j in range(0,M):
x[j,Ne]=-inf
while abs(x[j,Ne]-xe)>0.1:
e=standard_normal(p)
u=sqrt(v)*dot(A,e.T)
for i in range(0,N):
u=append(sum(u*phi)+normal(0,theta0),u[0:p-1])
x[j,i]=u[0]
for i in range(0,N):
f[:,i]=histogram(x[0:j+1,i],bins=xk)[0][0:Kk-1]/float(len(x[0:j+1,i]))/(xk[1:Kk]-xk[0:Kk-1])
clf()
imshow(clip(f,0.0,1.0),cmap='gray',extent=[-0.5,N-0.5,-2.4,2.4],aspect='auto',origin='lower')
for i in arange(0,N):
ex[0,i]=xe*rhopro[i-Ne] # Erwartungswert
vi=(1-rhopro[i-Ne]**2)*v
ex[1,i]=ex[0,i]-sqrt(vi) # Erwartungswert-Standardabweichung
ex[2,i]=ex[0,i]+sqrt(vi) # Erwartungswert+Standardabweichung
mx[0,i]=mean(x[0:j+1,i]) # Mittelwert
if j>0:
vx=var(x[0:j+1,i],ddof=1)
mx[1,i]=mx[0,i]-sqrt(vx) # Mittelwert-Standardabweichung
mx[2,i]=mx[0,i]+sqrt(vx) # Mittelwert+Standardabweichung
plot(t,ex[0],'blue',t,ex[1],'red',t,ex[2],'red',mx[0],'cyan',t,mx[1],'orange',t,mx[2],'orange')
draw()
pause(0.01)
show()
|
M=1000;
N=100;
Ne=40;
xe=1;
t=(1:N)*dt;
x=zeros(M,N);
Kk=25;
xk=0.2*(0:Kk-1)-2.4;
f=zeros(Kk-1,N);
ex=zeros(3,N);
mx=zeros(3,N);
for j=1:M
x(j,Ne)=-inf;
while abs(x(j,Ne)-xe)>0.1
e=randn(1,p);
u=(sqrt(v)*A*e')';
for i=1:N
u=[sum(u.*phi)+normrnd(0,theta0),u(1:p-1)];
x(j,i)=u(1);
end
end
for i=1:N
fx=histc(x(1:j,i),xk)/length(x(1:j,i));
f(:,i)=fx(1:length(xk)-1);
end
imshow(imresize(mat2gray(min(1,f)),10))
set(gca,'YDir','normal')
hold on
for i=1:N
ex(1,i)=xe*rhopro(abs(i-Ne)+1); % Erwartungswert
vi=(1-rhopro(abs(i-Ne)+1)^2)*v;
ex(2,i)=ex(1,i)-sqrt(vi); % Erwartungswert-Standardabweichung
ex(3,i)=ex(1,i)+sqrt(vi); % Erwartungswert+Standardabweichung
mx(1,i)=mean(x(1:j,i)); % Mittelwert
if j>1
vx=var(x(1:j,i));
mx(2,i)=mx(1,i)-sqrt(vx); % Mittelwert-Standardabweichung
mx(3,i)=mx(1,i)+sqrt(vx); % Mittelwert+Standardabweichung
end
end
ti=10*(t-0.5*dt);
exi=10*(Kk-1)/4.8*(ex+2.4)+0.5;
mxi=10*(Kk-1)/4.8*(mx+2.4)+0.5;
plot(ti,exi(1,:),'blue',ti,exi(2,:),'red',ti,exi(3,:),'red',ti,mxi(1,:),'cyan',ti,mxi(2,:),'yellow',ti,mxi(3,:),'yellow')
hold off
pause(0.01)
end
|
| zwei Realisierungen (inklusive der Initialisierung) mit zwei Beobachtung |
M=2
N=100
Ne1=40
xe1=1.0
Ne2=50
xe2=0.5
x=zeros([M,N])
for j in range(0,M):
x[j,Ne1]=-inf
x[j,Ne2]=-inf
while (abs(x[j,Ne1]-xe1)>0.1) or (abs(x[j,Ne2]-xe2)>0.1):
e=standard_normal(p)
u=sqrt(v)*dot(A,e.T)
for i in range(0,N):
u=append(sum(u*phi)+normal(0,theta0),u[0:p-1])
x[j,i]=u[0]
t=arange(0,N)*dt
plot(t,x[0],'o',t,x[1],'o')
show()
|
M=2;
N=100;
Ne1=40;
xe1=1;
Ne2=50;
xe2=0.5;
x=zeros(M,N);
for j=1:M
x(j,Ne1)=-inf;
x(j,Ne2)=-inf;
while (abs(x(j,Ne1)-xe1)>0.1) || (abs(x(j,Ne2)-xe2)>0.1)
e=randn(1,p);
u=(sqrt(v)*A*e')';
for i=1:N
u=[sum(u.*phi)+normrnd(0,theta0),u(1:p-1)];
x(j,i)=u(1);
end
end
end
t=(1:N)*dt;
plot(t,x(1,:),'o',t,x(2,:),'o')
|
| Realisierungen, Verarbeitung zu einer zeitlich aufgelösten Wahrscheinlichkeitsdichte, Mittelwert, empirische Varianz und Prädiktionen |
M=1000
N=100
Ne1=40
xe1=1.0
Ne2=50
xe2=0.5
t=arange(0,N)*dt
x=zeros([M,N])
Kk=25
xk=0.2*arange(0,Kk)-2.4
f=zeros([Kk-1,N])
ex=zeros([3,N])
mx=zeros([3,N])
for j in range(0,M):
x[j,Ne1]=-inf
x[j,Ne2]=-inf
while (abs(x[j,Ne1]-xe1)>0.1) or (abs(x[j,Ne2]-xe2)>0.1):
e=standard_normal(p)
u=sqrt(v)*dot(A,e.T)
for i in range(0,N):
u=append(sum(u*phi)+normal(0,theta0),u[0:p-1])
x[j,i]=u[0]
for i in range(0,N):
f[:,i]=histogram(x[0:j+1,i],bins=xk)[0][0:Kk-1]/float(len(x[0:j+1,i]))/(xk[1:Kk]-xk[0:Kk-1])
clf()
imshow(clip(f,0.0,1.0),cmap='gray',extent=[-0.5,N-0.5,-2.4,2.4],aspect='auto',origin='lower')
for i in arange(0,N):
ex[0,i]=xe1*(rhopro[i-Ne1]-rhopro[i-Ne2]*rhopro[Ne2-Ne1])/float(1-rhopro[Ne2-Ne1]**2)+xe2*(rhopro[i-Ne2]-rhopro[i-Ne1]*rhopro[Ne2-Ne1])/float(1-rhopro[Ne2-Ne1]**2) # Erwartungswert
vi=(1-(rhopro[i-Ne1]-rhopro[i-Ne2]*rhopro[Ne2-Ne1])/float(1-rhopro[Ne2-Ne1]**2)*rhopro[i-Ne1]-(rhopro[i-Ne2]-rhopro[i-Ne1]*rhopro[Ne2-Ne1])/float(1-rhopro[Ne2-Ne1]**2)*rhopro[i-Ne2])*v
ex[1,i]=ex[0,i]-sqrt(vi) # Erwartungswert-Standardabweichung
ex[2,i]=ex[0,i]+sqrt(vi) # Erwartungswert+Standardabweichung
mx[0,i]=mean(x[0:j+1,i]) # Mittelwert
if j>0:
vx=var(x[0:j+1,i],ddof=1)
mx[1,i]=mx[0,i]-sqrt(vx) # Mittelwert-Standardabweichung
mx[2,i]=mx[0,i]+sqrt(vx) # Mittelwert+Standardabweichung
plot(t,ex[0],'blue',t,ex[1],'red',t,ex[2],'red',mx[0],'cyan',t,mx[1],'orange',t,mx[2],'orange')
draw()
pause(0.01)
show()
|
M=1000;
N=100;
Ne1=40;
xe1=1;
Ne2=50;
xe2=0.5;
t=(1:N)*dt;
x=zeros(M,N);
Kk=25;
xk=0.2*(0:Kk-1)-2.4;
f=zeros(Kk-1,N);
ex=zeros(3,N);
mx=zeros(3,N);
for j=1:M
x(j,Ne1)=-inf;
x(j,Ne2)=-inf;
while (abs(x(j,Ne1)-xe1)>0.1) || (abs(x(j,Ne2)-xe2)>0.1)
e=randn(1,p);
u=(sqrt(v)*A*e')';
for i=1:N
u=[sum(u.*phi)+normrnd(0,theta0),u(1:p-1)];
x(j,i)=u(1);
end
end
for i=1:N
fx=histc(x(1:j,i),xk)/length(x(1:j,i));
f(:,i)=fx(1:length(xk)-1);
end
imshow(imresize(mat2gray(min(1,f)),10))
set(gca,'YDir','normal')
hold on
for i=1:N
ex(1,i)=xe1*(rhopro(abs(i-Ne1)+1)-rhopro(abs(i-Ne2)+1)*rhopro(abs(Ne2-Ne1)+1))/(1-rhopro(abs(Ne2-Ne1)+1)^2)+xe2*(rhopro(abs(i-Ne2)+1)-rhopro(abs(i-Ne1)+1)*rhopro(abs(Ne2-Ne1)+1))/(1-rhopro(abs(Ne2-Ne1)+1)^2); % Erwartungswert
vi=(1-(rhopro(abs(i-Ne1)+1)-rhopro(abs(i-Ne2)+1)*rhopro(abs(Ne2-Ne1)+1))/(1-rhopro(abs(Ne2-Ne1)+1)^2)*rhopro(abs(i-Ne1)+1)-(rhopro(abs(i-Ne2)+1)-rhopro(abs(i-Ne1)+1)*rhopro(abs(Ne2-Ne1)+1))/(1-rhopro(abs(Ne2-Ne1)+1)^2)*rhopro(abs(i-Ne2)+1))*v;
ex(2,i)=ex(1,i)-sqrt(vi); % Erwartungswert-Standardabweichung
ex(3,i)=ex(1,i)+sqrt(vi); % Erwartungswert+Standardabweichung
mx(1,i)=mean(x(1:j,i)); % Mittelwert
if j>1
vx=var(x(1:j,i));
mx(2,i)=mx(1,i)-sqrt(vx); % Mittelwert-Standardabweichung
mx(3,i)=mx(1,i)+sqrt(vx); % Mittelwert+Standardabweichung
end
end
ti=10*(t-0.5*dt);
exi=10*(Kk-1)/4.8*(ex+2.4)+0.5;
mxi=10*(Kk-1)/4.8*(mx+2.4)+0.5;
plot(ti,exi(1,:),'blue',ti,exi(2,:),'red',ti,exi(3,:),'red',ti,mxi(1,:),'cyan',ti,mxi(2,:),'yellow',ti,mxi(3,:),'yellow')
hold off
pause(0.01)
end
|