|
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 sowie die Signal Processing Toolbox
%Octave: pkg load statistics
if exist('OCTAVE_VERSION','builtin')~=0
pkg load statistics
end
|
| Mehrere Mittelwert- und Varianzschätzungen mit bekanntem Erwartungswert und mit geschätztem Mittelwert, ohne und mit Besselsche Korrektur für unabhängige Werte |
M=10
N=10
mu=5
sigma=1
x=normal(mu,sigma,[M,N])
plot(arange(0,M),mean(x,axis=1),'o',arange(0,M),sum((x-mu)**2,axis=1)/float(N),'+',arange(0,M),var(x,axis=1),'o',arange(0,M),var(x,ddof=1,axis=1),'x')
show()
oder
M=10
mx=zeros(M)
s2a=zeros(M)
s2b=zeros(M)
s2c=zeros(M)
N=10
mu=5
sigma=1
for j in range(0,M):
x=normal(mu,sigma,N)
mx[j]=mean(x)
s2a[j]=sum((x-mu)**2)/float(N) # mit bekanntem Erwartungswert
s2b[j]=var(x) # ohne Bessel-Korrektur
s2c[j]=var(x,ddof=1) # mit Bessel-Korrektur
plot(arange(0,M),mx,'o',arange(0,M),s2a,'+',arange(0,M),s2b,'o',arange(0,M),s2c,'x')
show()
|
M=10;
N=10;
mu=5;
sigma=1;
x=normrnd(mu,sigma,[M,N]);
plot((1:M),mean(x,2),'o',(1:M),sum((x-mu).^2,2)/N,'+',(1:M),var(x,1,2),'o',(1:M),var(x,0,2),'x')
oder
M=10;
mx=zeros(1,M);
s2a=zeros(1,M);
s2b=zeros(1,M);
s2c=zeros(1,M);
N=10;
mu=5;
sigma=1;
for j=1:M
x=normrnd(mu,sigma,[1,N]);
mx(j)=mean(x);
s2a(j)=sum((x-mu).^2)/N; % mit bekanntem Erwartungswert
s2b(j)=var(x,1); % ohne Bessel-Korrektur
s2c(j)=var(x); % mit Bessel-Korrektur
end
plot((1:M),mx,'o',(1:M),s2a,'+',(1:M),s2b,'o',(1:M),s2c,'x')
|
| Eigenschaften des Mittelwertschätzers: Mittelung über 10000 Realisierungen |
M=10000
mx=zeros(M)
N=10
mu=5
sigma=1
for m in range(0,M):
x=normal(mu,sigma,N)
mx[m]=mean(x)
print("empir. syst. Fehler: ",mean(mx)-mu)
print("empir. Var. zuf. Fehler: ",var(mx,ddof=1))
print("empir. mittl. quadrat. Fehler: ",sum((mx-mu)**2)/float(M))
|
M=10000;
mx=zeros(1,M);
N=10;
mu=5;
sigma=1;
for m=1:M
x=normrnd(mu,sigma,[1,N]);
mx(m)=mean(x);
end
"empir. syst. Fehler: ",mean(mx)-mu
"empir. Var. zuf. Fehler: ",var(mx)
"empir. mittl. quadrat. Fehler: ",sum((mx-mu).^2)/M
|
| Vergleich der Varianzschätzer mit bekanntem Erwartungswert und mit geschätztem Mittelwert, ohne und mit Besselsche Korrektur für unabhängige Werte: Mittelung über 10000 Realisierungen |
M=10000
s2a=zeros(M)
s2b=zeros(M)
s2c=zeros(M)
N=10
mu=5
sigma=1
for m in range(0,M):
x=normal(mu,sigma,N)
s2a[m]=sum((x-mu)**2)/float(N) # mit bekanntem Erwartungswert
s2b[m]=var(x) # ohne Bessel-Korrektur
s2c[m]=var(x,ddof=1) # mit Bessel-Korrektur
print("empir. syst. Fehler: ",mean(s2a)-sigma**2,mean(s2b)-sigma**2,mean(s2c)-sigma**2)
print("empir. Var. zuf. Fehler: ",var(s2a,ddof=1),var(s2b,ddof=1),var(s2c,ddof=1))
print("empir. mittl. quadrat. Fehler: ",sum((s2a-sigma**2)**2)/float(M),sum((s2b-sigma**2)**2)/float(M),sum((s2c-sigma**2)**2)/float(M))
|
M=10000;
s2a=zeros(1,M);
s2b=zeros(1,M);
s2c=zeros(1,M);
N=10;
mu=5;
sigma=1;
for m=1:M
x=normrnd(mu,sigma,[1,N]);
s2a(m)=sum((x-mu).^2)/N; % mit bekanntem Erwartungswert
s2b(m)=var(x,1); % ohne Bessel-Korrektur
s2c(m)=var(x); % mit Bessel-Korrektur
end
"empir. syst. Fehler: ",mean(s2a)-sigma^2,mean(s2b)-sigma^2,mean(s2c)-sigma^2
"empir. Var. zuf. Fehler: ",var(s2a),var(s2b),var(s2c)
"empir. mittl. quadrat. Fehler: ",sum((s2a-sigma^2).^2)/M,sum((s2b-sigma^2).^2)/M,sum((s2c-sigma^2).^2)/M
|
| Mittlerer linearer Fehler und mittlerer quadratischer Fehler der Varianzschätzer mit bekanntem Erwartungswert und mit geschätztem Mittelwert, ohne und mit Besselsche Korrektur für verschiedene Stichprobenumfänge unabhängiger Werte |
M=10000
K=10
N=zeros(K)
MLEa=zeros(K)
MLEb=zeros(K)
MLEc=zeros(K)
MSEa=zeros(K)
MSEb=zeros(K)
MSEc=zeros(K)
mu=5
sigma=1
for k in range(0,K):
N[k]=maximum(2,int(round(10**((k+1)/4.0))))
s2a=zeros(M)
s2b=zeros(M)
s2c=zeros(M)
for m in range(0,M):
x=normal(mu,sigma,int(N[k]))
s2a[m]=sum((x-mu)**2)/float(len(x)) # mit bekanntem Erwartungswert
s2b[m]=var(x) # ohne Bessel-Korrektur
s2c[m]=var(x,ddof=1) # mit Bessel-Korrektur
MLEa[k]=sum(s2a-sigma**2)/float(M)
MLEb[k]=sum(s2b-sigma**2)/float(M)
MLEc[k]=sum(s2c-sigma**2)/float(M)
MSEa[k]=sum((s2a-sigma**2)**2)/float(M)
MSEb[k]=sum((s2b-sigma**2)**2)/float(M)
MSEc[k]=sum((s2c-sigma**2)**2)/float(M)
Ke=100
Ne=zeros(Ke)
for ke in range(0,Ke):
Ne[ke]=10**((ke+5)/40.0)
# Mittlerer linearer Fehler
semilogx(Ne,zeros(Ke),N,MLEa,'+',Ne,-1/Ne,N,MLEb,'o',Ne,zeros(Ke),N,MLEc,'x')
show()
# Mittlerer quadratischer Fehler
loglog(Ne,2/Ne,N,MSEa,'+',Ne,(2*Ne-1)/Ne**2,N,MSEb,'o',Ne,2/(Ne-1),N,MSEc,'x')
show()
|
M=10000;
K=10;
N=zeros(1,K);
MLEa=zeros(1,K);
MLEb=zeros(1,K);
MLEc=zeros(1,K);
MSEa=zeros(1,K);
MSEb=zeros(1,K);
MSEc=zeros(1,K);
mu=5;
sigma=1;
for k=1:K
N(k)=max(2,round(10^(k/4)));
s2a=zeros(1,M);
s2b=zeros(1,M);
s2c=zeros(1,M);
for m=1:M
x=normrnd(mu,sigma,[1,N(k)]);
s2a(m)=sum((x-mu).^2)/length(x); % mit bekanntem Erwartungswert
s2b(m)=var(x,1); % ohne Bessel-Korrektur
s2c(m)=var(x); % mit Bessel-Korrektur
end
MLEa(k)=sum(s2a-sigma^2)/M;
MLEb(k)=sum(s2b-sigma^2)/M;
MLEc(k)=sum(s2c-sigma^2)/M;
MSEa(k)=sum((s2a-sigma^2).^2)/M;
MSEb(k)=sum((s2b-sigma^2).^2)/M;
MSEc(k)=sum((s2c-sigma^2).^2)/M;
end
Ke=100;
Ne=zeros(1,Ke);
for ke=1:Ke
Ne(ke)=10^((ke+5)/40);
end
% Mittlerer linearer Fehler
semilogx(Ne,zeros(1,Ke),N,MLEa,'+',Ne,-1./Ne,N,MLEb,'o',Ne,zeros(1,Ke),N,MLEc,'x')
waitforbuttonpress
% Mittlerer quadratischer Fehler
loglog(Ne,2./Ne,N,MSEa,'+',Ne,(2*Ne-1)./Ne.^2,N,MSEb,'o',Ne,2./(Ne-1),N,MSEc,'x')
|
| Vergleich der Varianzschätzer mit bekanntem linearen Trend und bekanntem Erwartungswert und mit geschätztem linearen Trend und geschätztem Mittelwert, ohne und mit Besselsche Korrektur für unabhängige Werte: Mittelung über 10000 Realisierungen |
from scipy.signal import *
M=10000
s2a=zeros(M)
s2b=zeros(M)
s2c=zeros(M)
N=10
mu=5
nu=0.2
sigma=1
for m in range(0,M):
x=normal(mu,sigma,N)+nu*(arange(0,N)-(N-1)/2.0)
s2a[m]=sum((x-mu-nu*(arange(0,N)-(N-1)/2.0))**2)/float(N) # mit bekanntem linearen Trend und bekanntem Erwartungswert
x=detrend(x)
s2b[m]=var(x) # ohne Bessel-Korrektur
s2c[m]=var(x,ddof=2) # mit Bessel-Korrektur
print("empir. syst. Fehler: ",mean(s2a)-sigma**2,mean(s2b)-sigma**2,mean(s2c)-sigma**2)
print("empir. Var. zuf. Fehler: ",var(s2a,ddof=1),var(s2b,ddof=1),var(s2c,ddof=1))
print("empir. mittl. quadrat. Fehler: ",sum((s2a-sigma**2)**2)/float(M),sum((s2b-sigma**2)**2)/float(M),sum((s2c-sigma**2)**2)/float(M))
|
M=10000;
s2a=zeros(1,M);
s2b=zeros(1,M);
s2c=zeros(1,M);
N=10;
mu=5;
nu=0.2;
sigma=1;
for m=1:M
x=normrnd(mu,sigma,[1,N])+nu*((0:N-1)-(N-1)/2);
s2a(m)=sum((x-mu-nu*((0:N-1)-(N-1)/2)).^2)/N; % mit bekanntem linearen Trend und bekanntem Erwartungswert
x=detrend(x);
s2b(m)=var(x,1); % ohne Bessel-Korrektur
s2c(m)=N/(N-2)*var(x,1); % mit Bessel-Korrektur
end
"empir. syst. Fehler: ",mean(s2a)-sigma^2,mean(s2b)-sigma^2,mean(s2c)-sigma^2
"empir. Var. zuf. Fehler: ",var(s2a),var(s2b),var(s2c)
"empir. mittl. quadrat. Fehler: ",sum((s2a-sigma^2).^2)/M,sum((s2b-sigma^2).^2)/M,sum((s2c-sigma^2).^2)/M
|
| Mittlerer linearer und mittlerer quadratischer Fehler der Varianzschätzer mit bekanntem linearen Trend und bekanntem Erwartungswert und mit geschätztem linearen Trend und geschätztem Mittelwert, ohne und mit Besselsche Korrektur für verschiedene Stichprobenumfänge unabhängiger Werte |
from scipy.signal import *
M=10000
K=10
N=zeros(K)
MLEa=zeros(K)
MLEb=zeros(K)
MLEc=zeros(K)
MSEa=zeros(K)
MSEb=zeros(K)
MSEc=zeros(K)
mu=5
nu=0.2
sigma=1
for k in range(0,K):
N[k]=maximum(3,int(round(10**((k+2)/4.0))))
s2a=zeros(M)
s2b=zeros(M)
s2c=zeros(M)
for m in range(0,M):
x=normal(mu,sigma,int(N[k]))+nu*(arange(0,int(N[k]))-(N[k]-1)/2.0)
s2a[m]=sum((x-mu-nu*(arange(0,int(N[k]))-(N[k]-1)/2.0))**2)/float(len(x)) # mit bekanntem linearen Trend und bekanntem Erwartungswert
x=detrend(x)
s2b[m]=var(x) # ohne Bessel-Korrektur
s2c[m]=var(x,ddof=2) # mit Bessel-Korrektur
MLEa[k]=sum(s2a-sigma**2)/float(M)
MLEb[k]=sum(s2b-sigma**2)/float(M)
MLEc[k]=sum(s2c-sigma**2)/float(M)
MSEa[k]=sum((s2a-sigma**2)**2)/float(M)
MSEb[k]=sum((s2b-sigma**2)**2)/float(M)
MSEc[k]=sum((s2c-sigma**2)**2)/float(M)
Ke=100
Ne=zeros(Ke)
for ke in range(0,Ke):
Ne[ke]=10**((ke+13)/40.0)
# Mittlerer linearer Fehler
semilogx(Ne,zeros(Ke),N,MLEa,'+',Ne,-2/Ne,N,MLEb,'o',Ne,zeros(Ke),N,MLEc,'x')
show()
# Mittlerer quadratischer Fehler
loglog(Ne,2/Ne,N,MSEb,'o',Ne,2/Ne,N,MSEa,'+',Ne,2/(Ne-2),N,MSEc,'x')
show()
|
M=10000;
K=10;
N=zeros(1,K);
MLEa=zeros(1,K);
MLEb=zeros(1,K);
MLEc=zeros(1,K);
MSEa=zeros(1,K);
MSEb=zeros(1,K);
MSEc=zeros(1,K);
mu=5;
sigma=1;
for k=1:K
N(k)=max(3,round(10^((k+1)/4)));
s2a=zeros(1,M);
s2b=zeros(1,M);
s2c=zeros(1,M);
for m=1:M
x=normrnd(mu,sigma,[1,N(k)])+nu*((0:N(k)-1)-(N(k)-1)/2);
s2a(m)=sum((x-mu-nu*((0:N(k)-1)-(N(k)-1)/2)).^2)/length(x); % mit bekanntem linearen Trend und bekanntem Erwartungswert
x=detrend(x);
s2b(m)=var(x,1); % ohne Bessel-Korrektur
s2c(m)=N(k)/(N(k)-2)*var(x,1); % mit Bessel-Korrektur
end
MLEa(k)=sum(s2a-sigma^2)/M;
MLEb(k)=sum(s2b-sigma^2)/M;
MLEc(k)=sum(s2c-sigma^2)/M;
MSEa(k)=sum((s2a-sigma^2).^2)/M;
MSEb(k)=sum((s2b-sigma^2).^2)/M;
MSEc(k)=sum((s2c-sigma^2).^2)/M;
end
Ke=100;
Ne=zeros(1,Ke);
for ke=1:Ke
Ne(ke)=10^((ke+13)/40);
end
% Mittlerer linearer Fehler
semilogx(Ne,zeros(1,Ke),N,MLEa,'+',Ne,-2./Ne,N,MLEb,'o',Ne,zeros(1,Ke),N,MLEc,'x')
waitforbuttonpress
% Mittlerer quadratischer Fehler
loglog(Ne,2./Ne,N,MSEb,'o',Ne,2./Ne,N,MSEa,'+',Ne,2./(Ne-2),N,MSEc,'x')
|