Signal- und Messdatenverarbeitung


Testcodes für Kapitel 3: Verteilungsfunktionen

weitere Codes zu Verteilungsfunktionen ohne bzw. mit Gewichtung sowie zu Zufallszahlen und Rauschgeneratoren

systematische Zusammenstellung von kombinierbaren Programmabschnitten zur Signal- und Messdatenverarbeitung

Python Matlab/Octave
Vorbereitung (Laden von Modulen/Paketen)
#python -m pip install --upgrade numpy
#python -m pip install --upgrade matplotlib
#python -m pip install --upgrade scipy
#pip install --upgrade numpy
#pip install --upgrade matplotlib
#pip install --upgrade scipy
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
from scipy.special import * #für die Fehler- und die Falkultätsfunktion
%Matlab: erfordert die Statistics and Machine Learning Toolbox
%Octave: pkg load statistics
if exist('OCTAVE_VERSION','builtin')~=0
pkg load statistics
end
Generieren von normalverteilten Zufallszahlen
N=10000
mu=3
sigma=1
x=normal(mu,sigma,N) #Hier muss die Standardabweichung angegeben werden.
plot(x,'.')
show()
N=10000;
mu=3;
sigma=1;
x=normrnd(mu,sigma,[1,N]); %Hier muss die Standardabweichung angegeben werden.
plot(x,'.')
gleichmäßig verteilte Abtastpunkte …
K=20
xi=0.5*arange(0,K)-2
plot(xi,'o')
show()
K=20;
xi=0.5*(0:K-1)-2;
plot(xi,'o')
… oder ungleichmäßig verteilte Abtastpunkte
K=20
xi=0.5*(1+arange(0,K)/float(5))**2-3
plot(xi,'o')
show()
K=20;
xi=0.5*(1+(0:K-1)/5).^2-3;
plot(xi,'o')
kumulative Verteilungsfunktion s. Anhang D.5
F(x)=1σ2πxe12(ξμσ)2dξ=12[1+erf(xμσ2)]
=12πσ2xe(ξμ)22σ2dξ=12[1+erf(xμ2σ2)]
linksseitige Stetigkeit:
#Bei histogram() werden jeweils Werte gleich der unteren Intervallgrenze mitgezählt.
#Im letzten Intervall werden dann aber Werte beider Intervallgrenzen gezählt.
#Unterhalb des ersten Abfragewertes werden keine Werte gezählt.
#Deshalb werden ein Intervall von -∞ bis zur ersten Intervallgrenze und ein Intervall ab der letzten Intervallgrenze bis +∞ ergänzt.
F=cumsum(histogram(x,bins=append(append(-inf,xi),inf))[0][0:len(xi)])/float(len(x))
plot(0.05*arange(-60,180),0.5*(1+erf((0.05*arange(-60,180)-mu)/sqrt(2*sigma**2))),xi,F,'o')
show()
rechtsseitige Stetigkeit:
#Für rechtsseitige Stetigkeit werden alle Werte und Intervallgrenzen invertiert und ihre Reihenfolge umgekehrt.
#Hier hilft dann, dass von histogram() Werte gleich der unteren Intervallgrenzen mitgezählt werden.
#Es muss dann nur noch ein Intervall zur Unterscheidung der Werte bis zur letzten Intervallgrenze und oberhalb der letzten Intervallgrenze ergänzt werden.
F=cumsum(flip(histogram(-x,bins=flip(-append(-inf,xi)))[0]))/float(len(x))
plot(0.05*arange(-60,180),0.5*(1+erf((0.05*arange(-60,180)-mu)/sqrt(2*sigma**2))),xi,F,'o')
show()
linksseitige Stetigkeit:
%Bei histc() werden jeweils Werte gleich der unteren Intervallgrenze mitgezählt.
%Unterhalb des ersten Abfragewertes werden keine Werte gezählt.
%Deshalb wird ein Intervall von -∞ bis zur ersten Intervallgrenze ergänzt.
F=cumsum(histc(x,[-inf,xi]))/length(x);
F=F(1:end-1);
plot((-3:0.05:9),normcdf((-3:0.05:9),mu,sigma),xi,F,'o')
rechtsseitige Stetigkeit:
%Für rechtsseitige Stetigkeit werden alle Werte und Intervallgrenzen invertiert und ihre Reihenfolge umgekehrt.
%Auch hier muss ein Intervall ergänzt werden, um alle Werte vor der ersten Intervallgrenze mitzuzählen.
F=cumsum(flip(histc(-x,flip(-[-inf,xi]))))/length(x);
F=F(2:end);
plot((-3:0.05:9),normcdf((-3:0.05:9),mu,sigma),xi,F,'o')
Wahrscheinlichkeitsdichtefunktion f(x)=1σ2πe12(xμσ)2=12πσ2e(xμ)22σ2
K=len(xi)
f=histogram(x,bins=xi)[0][0:K-1]/float(len(x))/(xi[1:K]-xi[0:K-1])
plot(0.05*arange(-100,100)+mu,exp(-(0.05*arange(-100,100))**2/(2*sigma**2))/(sqrt(2*pi*sigma**2)))
step(xi,append(0,f))
show()
f=histc(x,xi)/length(x);
f=f(1:length(xi)-1);
f=f./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(f));
plot((-5:0.05:5)+mu,exp(-(-5:0.05:5).^2/(2*sigma^2))/sqrt(2*pi*sigma^2))
hold on
stairs(xi,[f,0])
hold off
Generieren von Poisson-verteilten Zufallszahlen
N=10000
Lambda=3
x=poisson(Lambda,N)
plot(x,'.')
show()
N=10000;
lambda=3;
x=poissrnd(lambda,[1,N]);
plot(x,'.')
Abtastpunkte (mit Zwischenwerten)
K=60
xi=0.2*arange(0,K)-1
plot(xi,'o')
show()
K=60;
xi=0.2*(0:K-1)-1;
plot(xi,'o')
kumulative Verteilungsfunktion (Programm unverändert) s. Anhang D.4
linksseitige Stetigkeit: F(x)=eλk=0,k<xλkk!
rechtsseitige Stetigkeit: F(x)=eλk=0,kxλkk!
linksseitige Stetigkeit:
#Bei histogram() werden jeweils Werte gleich der unteren Intervallgrenze mitgezählt.
#Im letzten Intervall werden dann aber Werte beider Intervallgrenzen gezählt.
#Unterhalb des ersten Abfragewertes werden keine Werte gezählt.
#Deshalb werden ein Intervall von -∞ bis zur ersten Intervallgrenze und ein Intervall ab der letzten Intervallgrenze bis +∞ ergänzt.
F=cumsum(histogram(x,bins=append(append(-inf,xi),inf))[0][0:len(xi)])/float(len(x))
Fp=zeros(len(xi))
for i in range(0,len(xi)):
  Fp[i]=sum(exp(-Lambda)*Lambda**(arange(0,int(ceil(xi[i]))))/factorial(arange(0,int(ceil(xi[i])))))

plot(xi,Fp,'o',xi,F,'.')
show()
rechtsseitige Stetigkeit:
#Für rechtsseitige Stetigkeit werden alle Werte und Intervallgrenzen invertiert und ihre Reihenfolge umgekehrt.
#Hier hilft dann, dass von histogram() Werte gleich der unteren Intervallgrenzen mitgezählt werden.
#Es muss dann nur noch ein Intervall zur Unterscheidung der Werte bis zur letzten Intervallgrenze und oberhalb der letzten Intervallgrenze ergänzt werden.
F=cumsum(flip(histogram(-x,bins=flip(-append(-inf,xi)))[0]))/float(len(x))
Fp=zeros(len(xi))
for i in range(0,len(xi)):
  Fp[i]=sum(exp(-Lambda)*Lambda**(arange(0,int(floor(xi[i]))+1))/factorial(arange(0,int(floor(xi[i]))+1)))

plot(xi,Fp,'o',xi,F,'.')
show()
linksseitige Stetigkeit:
%Bei histc() werden jeweils Werte gleich der unteren Intervallgrenze mitgezählt.
%Unterhalb des ersten Abfragewertes werden keine Werte gezählt.
%Deshalb wird ein Intervall von -∞ bis zur ersten Intervallgrenze ergänzt.
F=cumsum(histc(x,[-inf,xi]))/length(x);
F=F(1:end-1);
Fp=zeros([1,length(xi)]);
for i=1:length(xi)
Fp(i)=sum(exp(-lambda)*lambda.^(0:ceil(xi(i))-1)./factorial(0:ceil(xi(i))-1));
end
plot(xi,Fp,'o',xi,F,'.')
rechtsseitige Stetigkeit:
%Für rechtsseitige Stetigkeit werden alle Werte und Intervallgrenzen invertiert und ihre Reihenfolge umgekehrt.
%Auch hier muss ein Intervall ergänzt werden, um alle Werte vor der ersten Intervallgrenze mitzuzählen.
F=cumsum(flip(histc(-x,flip(-[-inf,xi]))))/length(x);
F=F(2:end);
Fp=zeros([1,length(xi)]);
for i=1:length(xi)
Fp(i)=sum(exp(-lambda)*lambda.^(0:floor(xi(i)))./factorial(0:floor(xi(i))));
end
plot(xi,Fp,'o',xi,F,'.')
zum Test nicht alle Möglichkeiten berücksichtigt linksseitige Stetigkeit:
K=10
xi=arange(2,2+K)
plot(xi,'o')
show()
#Bei histogram() werden jeweils Werte gleich der unteren Intervallgrenze mitgezählt.
#Im letzten Intervall werden dann aber Werte beider Intervallgrenzen gezählt.
#Unterhalb des ersten Abfragewertes werden keine Werte gezählt.
#Deshalb werden ein Intervall von -∞ bis zur ersten Intervallgrenze und ein Intervall ab der letzten Intervallgrenze bis +∞ ergänzt.
F=cumsum(histogram(x,bins=append(append(-inf,xi),inf))[0][0:len(xi)])/float(len(x))
Fp=zeros(len(xi))
for i in range(0,len(xi)):
  Fp[i]=sum(exp(-Lambda)*Lambda**(arange(0,int(ceil(xi[i]))))/factorial(arange(0,int(ceil(xi[i])))))

plot(xi,Fp,'o',xi,F,'.')
show()
rechtsseitige Stetigkeit:
K=10
xi=arange(2,2+K)
plot(xi,'o')
show()
#Für rechtsseitige Stetigkeit werden alle Werte und Intervallgrenzen invertiert und ihre Reihenfolge umgekehrt.
#Hier hilft dann, dass von histogram() Werte gleich der unteren Intervallgrenzen mitgezählt werden.
#Es muss dann nur noch ein Intervall zur Unterscheidung der Werte bis zur letzten Intervallgrenze und oberhalb der letzten Intervallgrenze ergänzt werden.
F=cumsum(flip(histogram(-x,bins=flip(-append(-inf,xi)))[0]))/float(len(x))
Fp=zeros(len(xi))
for i in range(0,len(xi)):
  Fp[i]=sum(exp(-Lambda)*Lambda**(arange(0,int(floor(xi[i]))+1))/factorial(arange(0,int(floor(xi[i]))+1)))

plot(xi,Fp,'o',xi,F,'.')
show()
linksseitige Stetigkeit:
K=10;
xi=(2:2+K-1);
plot(xi,'o')
waitforbuttonpress
%Bei histc() werden jeweils Werte gleich der unteren Intervallgrenze mitgezählt.
%Unterhalb des ersten Abfragewertes werden keine Werte gezählt.
%Deshalb wird ein Intervall von -∞ bis zur ersten Intervallgrenze ergänzt.
F=cumsum(histc(x,[-inf,xi]))/length(x);
F=F(1:end-1);
Fp=zeros([1,length(xi)]);
for i=1:length(xi)
Fp(i)=sum(exp(-lambda)*lambda.^(0:ceil(xi(i))-1)./factorial(0:ceil(xi(i))-1));
end
plot(xi,Fp,'o',xi,F,'.')
rechtsseitige Stetigkeit:
K=10;
xi=(2:2+K-1);
plot(xi,'o')
waitforbuttonpress
%Für rechtsseitige Stetigkeit werden alle Werte und Intervallgrenzen invertiert und ihre Reihenfolge umgekehrt.
%Auch hier muss ein Intervall ergänzt werden, um alle Werte vor der ersten Intervallgrenze mitzuzählen.
F=cumsum(flip(histc(-x,flip(-[-inf,xi]))))/length(x);
F=F(2:end);
Fp=zeros([1,length(xi)]);
for i=1:length(xi)
Fp(i)=sum(exp(-lambda)*lambda.^(0:floor(xi(i)))./factorial(0:floor(xi(i))));
end
plot(xi,Fp,'o',xi,F,'.')
Wahrscheinlichkeitsfunktion P(ξk)=λkk!eλξk=k
xi=arange(0,11)
#Bei histogram() werden im letzten Intervall Werte beider Intervallgrenzen gezählt.
#Deshalb wird der letzte Abfragewert noch ein zweites Mal eingetragen, um die Werte am letzten und am vorlestzten Abfragewert zu trennen.
P=histogram(x,bins=append(xi,xi[-1]))[0]/float(len(x))
Pp=zeros(len(xi))
for i in range(0,len(xi)):
  Pp[i]=(xi[i]>=0)*(ceil(xi[i])==floor(xi[i]))*exp(-Lambda)*Lambda**xi[i]/factorial(maximum(0,xi[i]))

plot(xi,Pp,'o',xi,P,'.')
show()
xi=(0:10);
P=histc(x,xi)/length(x);
Pp=zeros([1,length(xi)]);
for i=1:length(xi)
Pp(i)=(xi(i)>=0)*(ceil(xi(i))==floor(xi(i)))*exp(-lambda)*lambda^xi(i)./factorial(max(0,floor(xi(i))));
end
plot(xi,Pp,'o',xi,P,'.')
Beispiel 1: kumulative Verteilungsfunktion und Wahrscheinlichkeitsdichtefunktion einer gleichverteilten Zufallsgröße im Intervall [1,1)
N=10000
x=uniform(-1,1,N)
plot(x,'.')
show()
K=40
xi=0.1*arange(0,K)-2
F=cumsum(histogram(x,bins=append(append(-inf,xi),inf))[0][0:len(xi)])/float(len(x))
plot([-2,-1,1,2],[0,0,1,1],xi,F,'o')
show()
f=histogram(x,bins=xi)[0][0:K-1]/float(len(x))/(xi[1:K]-xi[0:K-1])
plot([-2,-1,-1,1,1,2],[0,0,0.5,0.5,0,0])
step(xi,append(0,f))
show()
N=10000;
x=unifrnd(-1,1,[1,N]);
plot(x,'.')
waitforbuttonpress
K=40;
xi=0.1*(0:K-1)-2;
F=cumsum(histc(x,[-inf,xi]))/length(x);
F=F(1:end-1);
plot([-2,-1,1,2],[0,0,1,1],xi,F,'o')
waitforbuttonpress
f=histc(x,xi)/length(x);
f=f(1:length(xi)-1);
f=f./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(f));
plot([-2,-1,-1,1,1,2],[0,0,0.5,0.5,0,0])
hold on
stairs(xi,[f,0])
hold off
Beispiel 2: kumulative Verteilungsfunktion und Wahrscheinlichkeitsfunktion eines Münzwurfs, wobei nur "Kopf" (0) und "Zahl" (1) als mögliche Ereignisse angenommen werden.
N=10000
x=randint(0,2,N)
plot(x,'.')
show()
K=30
xi=0.1*arange(0,K)-1
F=cumsum(histogram(x,bins=append(append(-inf,xi),inf))[0][0:len(xi)])/float(len(x))
plot([-1,0,0,1,1,2],[0,0,0.5,0.5,1,1],xi,F,'o')
show()
xi=[0,1]
P=histogram(x,bins=append(xi,xi[-1]))[0]/float(len(x))
axes=gca()
axes.set_xlim([-1,2])
axes.set_ylim([0,1])
plot([0,1],[0.5,0.5],'o',xi,P,'.')
show()
N=10000;
x=randi([0,1],[1,N]);
plot(x,'.')
waitforbuttonpress
K=30;
xi=0.1*(0:K-1)-1;
F=cumsum(histc(x,[-inf,xi]))/length(x);
F=F(1:end-1);
plot([-1,0,0,1,1,2],[0,0,0.5,0.5,1,1],xi,F,'o')
waitforbuttonpress
xi=[0,1];
P=histc(x,xi)/length(x);
plot([0,1],[0.5,0.5],'o',xi,P,'.')
axis([-1,2,0,1])
Wahrscheinlichkeitsdichtefunktion von unabhängigen standardnormalverteilten Zufallszahlen, …
N=10000
x=standard_normal(N)
K=100
xi=0.1*arange(0,K)-5
K=len(xi)
f=histogram(x,bins=xi)[0][0:K-1]/float(len(x))/(xi[1:K]-xi[0:K-1])
plot(0.05*arange(-100,100),exp(-(0.05*arange(-100,100))**2/2.0)/(sqrt(2*pi)))
step(xi,append(0,f))
show()
N=10000;
x=randn(1,N);
K=100;
xi=0.1*(0:K-1)-5;
f=histc(x,xi)/length(x);
f=f(1:length(xi)-1);
f=f./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(f));
plot((-5:0.05:5),exp(-(-5:0.05:5).^2/2)/sqrt(2*pi))
hold on
stairs(xi,[f,0])
hold off
ihren Quadraten und …
y=x**2
K=100
xi=0.1*arange(0,K)-5
K=len(xi)
f=histogram(y,bins=xi)[0][0:K-1]/float(len(y))/(xi[1:K]-xi[0:K-1])
step(xi,append(0,f))
show()
y=x.^2;
K=100;
xi=0.1*(0:K-1)-5;
f=histc(y,xi)/length(y);
f=f(1:length(xi)-1);
f=f./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(f));
stairs(xi,[f,0])
ihren Inversen
y=1.0/x
K=100
xi=0.1*arange(0,K)-5
K=len(xi)
f=histogram(y,bins=xi)[0][0:K-1]/float(len(y))/(xi[1:K]-xi[0:K-1])
plot(0.05*arange(-100,100),exp(-1.0/(2*(0.05*arange(-100,100))**2))/((0.05*arange(-100,100))**2*sqrt(2*pi)))
step(xi,append(0,f))
show()
y=1./x;
K=100;
xi=0.1*(0:K-1)-5;
f=histc(y,xi)/length(y);
f=f(1:length(xi)-1);
f=f./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(f));
plot(0.05*(-100:100),exp(-1.0./(2*(0.05*(-100:100)).^2))./((0.05*(-100:100)).^2*sqrt(2*pi)))
hold on
stairs(xi,[f,0])
hold off
Wahrscheinlichkeitsdichtefunktion von zwei unabhängigen standardnormalverteilten Zufallszahlen, …
N=10000
x=standard_normal(N)
y=standard_normal(N)
K=100
xi=0.1*arange(0,K)-5
K=len(xi)
fx=histogram(x,bins=xi)[0][0:K-1]/float(len(x))/(xi[1:K]-xi[0:K-1])
fy=histogram(y,bins=xi)[0][0:K-1]/float(len(y))/(xi[1:K]-xi[0:K-1])
plot(0.05*arange(-100,100),exp(-(0.05*arange(-100,100))**2/2.0)/(sqrt(2*pi)))
step(xi,append(0,fx))
step(xi,append(0,fy))
show()
N=10000;
x=randn(1,N);
y=randn(1,N);
K=100;
xi=0.1*(0:K-1)-5;
fx=histc(x,xi)/length(x);
fx=fx(1:length(xi)-1);
fx=fx./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(fx));
fy=histc(y,xi)/length(y);
fy=fy(1:length(xi)-1);
fy=fy./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(fy));
plot((-5:0.05:5),exp(-(-5:0.05:5).^2/2)/sqrt(2*pi))
hold on
stairs(xi,[fx,0])
stairs(xi,[fy,0])
hold off
ihrer Summen (im Vergleich zu der eines der beiden Summanden), …
z=x+y
K=100
xi=0.1*arange(0,K)-5
K=len(xi)
fz=histogram(z,bins=xi)[0][0:K-1]/float(len(z))/(xi[1:K]-xi[0:K-1])
step(xi,append(0,fz))
step(xi,append(0,fx))
show()
z=x+y;
K=100;
xi=0.1*(0:K-1)-5;
fz=histc(z,xi)/length(z);
fz=fz(1:length(xi)-1);
fz=fz./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(fz));
stairs(xi,[fz,0])
hold on
stairs(xi,[fx,0])
hold off
ihren Produkten (im Vergleich zum Quadrat einer der beiden) …
z=x*y
K=100
xi=0.1*arange(0,K)-5
K=len(xi)
f=histogram(z,bins=xi)[0][0:K-1]/float(len(z))/(xi[1:K]-xi[0:K-1])
f2=histogram(x**2,bins=xi)[0][0:K-1]/float(len(x))/(xi[1:K]-xi[0:K-1])
step(xi,append(0,f))
step(xi,append(0,f2))
show()
z=x.*y;
K=100;
xi=0.1*(0:K-1)-5;
f=histc(z,xi)/length(z);
f=f(1:length(xi)-1);
f=f./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(f));
f2=histc(x.^2,xi)/length(x);
f2=f2(1:length(xi)-1);
f2=f2./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(f2));
stairs(xi,[f,0])
hold on
stairs(xi,[f2,0])
hold off
z. B. Produkt korrelierter Zufallszahlen
y1=sqrt(0.5)*(x+y)
K=100
xi=0.1*arange(0,K)-5
K=len(xi)
fx=histogram(x,bins=xi)[0][0:K-1]/float(len(x))/(xi[1:K]-xi[0:K-1])
fy=histogram(y,bins=xi)[0][0:K-1]/float(len(y))/(xi[1:K]-xi[0:K-1])
fy1=histogram(y1,bins=xi)[0][0:K-1]/float(len(y1))/(xi[1:K]-xi[0:K-1])
step(xi,append(0,fx))
step(xi,append(0,fy))
step(xi,append(0,fy1))
show()
z=x*y
z1=x*y1
fz=histogram(z,bins=xi)[0][0:K-1]/float(len(z))/(xi[1:K]-xi[0:K-1])
fz1=histogram(z1,bins=xi)[0][0:K-1]/float(len(z1))/(xi[1:K]-xi[0:K-1])
step(xi,append(0,fz1))
step(xi,append(0,fz))
show()
y1=sqrt(0.5)*(x+y);
K=100;
xi=0.1*(0:K-1)-5;
fx=histc(x,xi)/length(x);
fx=fx(1:length(xi)-1);
fx=fx./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(fx));
fy=histc(y,xi)/length(y);
fy=fy(1:length(xi)-1);
fy=fy./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(fy));
fy1=histc(y1,xi)/length(y1);
fy1=fy1(1:length(xi)-1);
fy1=fy1./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(fy1));
stairs(xi,[fx,0])
hold on
stairs(xi,[fy,0])
stairs(xi,[fy1,0])
hold off
waitforbuttonpress
z=x.*y;
z1=x.*y1;
fz=histc(z,xi)/length(z);
fz=fz(1:length(xi)-1);
fz=fz./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(fz));
fz1=histc(z1,xi)/length(z1);
fz1=fz1(1:length(xi)-1);
fz1=fz1./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(fz1));
stairs(xi,[fz1,0])
hold on
stairs(xi,[fz,0])
hold off
… und ihren Quotienten (haben relativ oft extreme Werte, im Vergleich zur Normalverteilung)
z=x/y
K=100
xi=0.2*arange(0,K)-10
K=len(xi)
fz=histogram(z,bins=xi)[0][0:K-1]/float(len(z))/(xi[1:K]-xi[0:K-1])
fx=histogram(x,bins=xi)[0][0:K-1]/float(len(x))/(xi[1:K]-xi[0:K-1])
step(xi,append(0,fz))
step(xi,append(0,fx))
yscale('log')
show()
z=x./y;
K=100;
xi=0.2*(0:K-1)-10;
fz=histc(z,xi)/length(z);
fz=fz(1:length(xi)-1);
fz=fz./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(fz));
fx=histc(x,xi)/length(x);
fx=fx(1:length(xi)-1);
fx=fx./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(fx));
stairs(xi,[fz,0])
hold on
stairs(xi,[fx,0])
hold off
set(gca,'YScale','log')
Generieren von Zufallszahlen einer beliebigen Verteilung (hier quadratische Funktion im Intervall [1,1))
N=10000
x=uniform(-1,1,N)
for i in range(0,N):
  while (random()>=1-x[i]**2):
    x[i]=uniform(-1,1)
  

plot(x,'.')
show()
K=21
xi=0.1*arange(0,K)-1
K=len(xi)
f=histogram(x,bins=xi)[0][0:K-1]/float(len(x))/(xi[1:K]-xi[0:K-1])
plot(0.01*arange(-100,101),0.75*(1-(0.01*arange(-100,101))**2))
step(xi,append(0,f))
show()
N=10000;
x=unifrnd(-1,1,[1,N]);
for i=1:N
while (rand()>=1-x(i)^2)
x(i)=unifrnd(-1,1);
end
end
plot(x,'.')
waitforbuttonpress
K=21;
xi=0.1*(0:K-1)-1;
f=histc(x,xi)/length(x);
f=f(1:length(xi)-1);
f=f./reshape(xi(2:length(xi))-xi(1:length(xi)-1),size(f));
plot((-1:0.01:1),0.75*(1-(-1:0.01:1).^2))
hold on
stairs(xi,[f,0])
hold off