(*  WQ  FS  2013  

 1, 2, oder 3 Wuerfel: Jeweils Augenzahl addieren.
Zufallsvariable addieren bedeutet Verteilungsdichten falten. *)

p[n_Integer]= If[(n>0) && (n<7), 1./6., 0.0];
w1=Table[6*p[k],{k,1,7}] 

p2[k_Integer]= Sum[p[j] * p[k-j],{j,1,6}];
w2=Table[36*p2[kk],{kk,1,13}]

p3[l_Integer]= Sum[p2[k]*p[l-k],{k,1,12}];        
w3=Table[216*p3[k],{k,19}]

grau=GrayLevel[0.5];
hell=GrayLevel[0.8];
weiss=GrayLevel[1.0];
b1=BarChart[{w1},ChartStyle->{grau},PlotRange->All];
b2=BarChart[{w2},ChartStyle->{hell},PlotRange->All];
b3=BarChart[{w3},ChartStyle->{weiss},PlotRange->All,ChartLabels->{Range[19]},AxesLabel->{"Augenzahl","Anzahl"} ];

Show[b3,b2,b1]

(* Analog: The convolution of UnitBox with itself is UnitTriangle *)
Convolve[UnitBox[x],UnitBox[x],x,t]

Plot[%, {t,-1.2,1.2}]

Table[Plot[PDF[UniformSumDistribution[k],x],{x,0,4},Filling->Axis,Exclusions->None,PlotRange->{0,1},PlotLabel->k],{k,1,4}]

(*  Gleichverteilung gegen Binomialverteilung
6 x Wuerfel  <--> 6 x Muenzwurf
Wuerfel habe Punkte von 0 bis 5     *)

w1=Table[p[k],{k,1,6}] 
b6=Table[Binomial[5,k]/2^5,{k,0,5}]

v1=BarChart[{w1},ChartStyle->Opacity[0.5],PlotRange->All,ChartLabels->{Range[6]-1}]

v2=BarChart[{b6},ChartStyle->{grau},PlotRange->All,ChartLabels->{Range[6]-1},
AxesLabel->{"Augenzahl/AnzahlKopf","Anzahl"} ]

Show[v2,v1]

(*  mit Mma wirklich wuerfeln *)

w1=Table[RandomInteger[5]+1,{k,1,1000}];
w2=Table[RandomInteger[5]+1,{k,1,1000}];
w3=Table[RandomInteger[5]+1,{k,1,1000}];

Short[w1]

Histogram[w1]

Histogram[w1+w2]

Histogram[w1+w2+w3]

(* Mma kennt alle Verteilungen, hier die Normalverteilung *)

da1=Histogram[RandomVariate[NormalDistribution[0,1],1000]]

da2=Histogram[RandomVariate[NormalDistribution[3,1/2],1000]]

Show[da2,da1,PlotRange->All]

(* Overlay several probability density function ( PDF ) plots for Poisson distributions: *)

PDF[PoissonDistribution[λ],k]

Show[Histogram[Table[RandomInteger[PoissonDistribution[10 i],1000],{i,5}],{1},"PDF",PlotRange->All],Plot[Evaluate[Table[PDF[PoissonDistribution[10 μ],x],{μ,5}]],{x,0,80},PlotStyle->Thick,PlotRange->All]]

(* PDF direkt *)

PDF[NormalDistribution[0,1],x]  (*  my=0, sigma=1 *)

PDF[NormalDistribution[3,1/2],x]  (*  my=3, sigma=1/2 *)

Plot[PDF[NormalDistribution[0,1],x],{x,-5,5},Filling->Axis]

punkte=Graphics[{PointSize[0.02],Table[Point[{x,N[PDF[NormalDistribution[0,1],x]]}],{x,1,3}]}];
bi=Plot[{PDF[NormalDistribution[0,1],x] ,PDF[NormalDistribution[0,2],x],PDF[NormalDistribution[0,3],x]},{x,-3,7}];
Show[bi,punkte]

(*  Die Markierungen sind bei 1 sigma, 2 sigma und 3 sigma 
   entsprechend  2/3,  95%,  und  99.5%  der Flaeche *)

(* Noch mal Faltung *)

p[n_Integer]= If[(n>0) && (n<7), 1./6., 0.0];

DiscretePlot[p[n],{n,0,6}]

DiscreteConvolve[p[n],p[n],n,m]

DiscretePlot[%,{m,0,13}]

(* Faltung fuer stetige Dichten  *)

Convolve[PDF[NormalDistribution[3,1/2],x],PDF[NormalDistribution[0,1],x],x,y]

(*  Fuer die Normalverteilung gilt: 
Ihre Faltung mit sich selbt gibt wieder eine Normalverteilung mit
 my= my_1 + my_2 und sigma^2=sigma_1^2 + sigma_2^2 *)

sigma=Sqrt[1/4 +1] ;  PDF[NormalDistribution[3,Sqrt[5]/2],x]

(*  Eine analoge Relation gilt fuer die Binomialverteilung *)

PDF[BinomialDistribution[n,p],k]

Mean[BinomialDistribution[n,p]]
Variance[BinomialDistribution[n,p]]

DiscretePlot[Evaluate@Table[PDF[BinomialDistribution[40,p],k],{p,{0.1,0.5,0.7}}],{k,40},PlotRange->All,PlotMarkers->Automatic]

DiscretePlot[CDF[BinomialDistribution[40,.7],x],{x,10,40},ExtentSize->Right,ExtentMarkers->{"Filled","Empty"}]

CDF[BinomialDistribution[5,.6],x]//PiecewiseExpand

RandomVariate[BernoulliDistribution[0.75],10]

(*  Die Wahrscheinlichkeit eines Tores beim Elfmeter *)

RandomVariate[BernoulliDistribution[0.8],10]/.{0->"Gehalten",1->"Tor"}

(*A packet consisting of a string of n symbols is transmitted over a noisy channel. Each symbol has probability(10^(-4))of incorrect transmission. Find  n such that the probability of incorrect packet transmission is less than (10^(-3)): *)
Probability[x≥1,x\[Distributed]BinomialDistribution[n,(10^(-4))]]

Reduce[%≤(10^(-3)),n,Reals] // N

Compute the same limit using a Poisson approximation:

Probability[x≥1,x\[Distributed]PoissonDistribution[n (10^(-4))]]

Reduce[%≤(10^(-3)),n,Reals] // N

(*  Faltung von Binomial Dichten *)

PDF[BinomialDistribution[n,p],k]

bi1[n_Integer,k_Integer,p_Real]= If[(n>=0)&&(k>=0) && (k<=n),(((1-p))^(-k+n)) (p^k) Binomial[n,k] , 0.0];

(* DiscreteConvolve[bi1[10,k,0.7],bi1[20,k,0.7],k,r] 
  geht nicht ?? Nutze direkte Faltungsformel *)

convBi[r_]:=Sum[bi1[28,k,0.7]*bi1[22,r-k,0.7],{k,0,50}];

DiscretePlot[convBi[r],{r,20,50}]

(*  Vergleich der beiden Faelle: ist natuerlich kein Beweis! *)

DiscretePlot[PDF[BinomialDistribution[50,0.7],k],{k,20,50}]

(* Wenn x≃BinomialDistribution[n,p] und y≃BinomialDistribution[m,p]
verteilt sind, dann ist eine Zufallsgroesse z=x+y nach 
z≃BinomialDistribution[n+m,p] verteilt *)

(*  ErlangDistribution *)

(* Die ErlangDistribution[k,λ] is the convolution of  k ExponentialDistribution[λ]  PDFs: *)

Convolve[PDF[ErlangDistribution[k,λ],x],PDF[ExponentialDistribution[λ],x],x,y]

PDF[ErlangDistribution[k+1,λ],y]

Table[Plot[PDF[ErlangDistribution[k,1],x],{x,0,8},Filling->Axis,Exclusions->None,PlotRange->{0,1},PlotLabel->k],{k,1,4}]