--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/draft/prob/prob.m	Mon Dec 26 19:21:22 2005 -0500
@@ -0,0 +1,135 @@
+# Nombre de résultats calculés
+# de -N à +N
+N = 20;
+# Taille du dé
+M = 10;
+# Nombre de résultats affichés
+# de -R à +R
+R = 20;
+# Nombre de dés
+J = 2;
+# Nombre de résultats affichés
+# de -R à +R
+R = 20;
+# Prob de 1 à N
+P1N = 1:N;
+# Prob de -N à N
+PDN = 0:2*N;
+
+NSG = 0:2*N;
+
+function y = P1_ferme(x,M,J)
+  if ((x < 1) || (x > M))
+    y = 0;
+  else
+    y = 1/M;
+  end;
+end;
+
+function y = P1_limite(x,M,J)
+  if (x < 1)
+    y = 0;
+  elseif (x < M)
+    y = 1/M;
+  elseif (x < J*M)
+    y = P1_limite(x-M,M,J)/M;
+  elseif (x == J*M)
+    y = (1/M)^J;
+  else
+    y = 0;	
+  end;
+end;
+
+function y = P1_ouvert(x,M,J)
+  if (x < 1)
+    y = 0;
+  elseif (x < M)
+    y = 1/M;
+  else
+    y = P1_ouvert(x-M,M)/M;
+  end;
+end;
+
+function y = P1(x,M,J)
+  y = P1_limite(x,M,J);
+end;
+
+function y = diff(x,P1N,N)
+  if (x < 0)
+    c = -x;
+  else
+    c = x;
+  end
+
+  if (c >= N)
+    y = 0;
+    return;
+  end
+
+  y = 0;
+  for i = (c+1):N
+    y = y + P1N(i) * P1N(i-c);
+  end;
+end;
+
+function y = prob_real(N, M, J)
+# usage: y = prob_real(N, M, J)
+# Return a -N:N vector containing
+# the compute value of probability
+# for J M-side dices.
+  for i = 1:N
+    p1(i) = P1(i,M,J);
+  end;
+  for i = -N:N
+    y(i+N+1) = diff(i,p1,N);
+  end;
+end;
+
+function y = prob_theo(N)
+# usage: y = prob_theo(N)
+# Return a -N:N vector containing
+# the theorical value of probability
+  for i = -N:N
+    if (i < 0)
+      y(i+N+1) = 0.5 * 2^(i/3);
+    else
+      y(i+N+1) = 1 - 0.5 * 2^(-i/3);
+    end;
+  end;
+end;
+
+PDN = prob_real(N,M,J);
+IPDN = cumsum(PDN);
+NSG = prob_theo(N);
+
+RANGE=-R:R;
+IRANGE=N-R+1:N+R+1;
+
+# Affichage des résultats
+IPDN(IRANGE)
+
+title("courbe de probabilité");
+plot(RANGE, IPDN(IRANGE), "-.", RANGE, NSG(IRANGE), "-");
+r = scanf("%s","C");
+title("erreur absolue");
+plot(RANGE, abs(IPDN(IRANGE)-NSG(IRANGE))*100);
+r = scanf("%s","C");
+title("erreur relative");
+plot(RANGE, abs(IPDN(IRANGE)-NSG(IRANGE))./NSG(IRANGE)*100);
+r = scanf("%s","C");
+
+function y = diff_sum(P,R)
+  Psize = size(P,1)
+  Rsize = size(R,1)
+  for i = 1:Rsize-1
+    off = R(i+1)-1
+    y(i) = sum(P(R(i):off))
+  end;
+  y(Rsize) = sum(P(R(Rsize):Psize))
+end;
+
+X = PDN(N:-1:1);
+X = X / sum(X);
+S = [1;5;10;15]
+diff_sum(X,S)
+