LABORATORIO DE PROBABILIDAD I, 2005/06                 	 18/11/05

%--------- PROPUESTA PARA HOY:

Generar permutaciones de los índices 1,...,N con igual probabilidad para las  N!
permutaciones posibles.
El código debe permitir esta variante: para un  n  fijado < N, extraer al azar
n  de esos índices sin repetirlos; lo anterior sería el caso n = N .

Una vez conseguido ese código, se pueden observar muchas variables aleatorias,
por ejemplo: signo de la permutación, número de inversiones, número de índices
que deja fijos... , y compararlos con sus valores esperados, que no es difícil
hallar con argumentos teóricos.


%--------- UN POSIBLE CODIGO PARA EL ALBUM, COMENTADO:

Se parte de la idea sugerida al final del Guión anterior: cada "subcolección" es
uniforme, luego el experimento equivale a 
  1) escoger subcolección, con probabilidades que NO sean proporcionales a su
tamaño;
  2) en la escogida, sacar al azar un cromo (todos equi-IP).

El código que sigue usa valores fijados para los tamaños cN + dN = N de ambas
partes, y repite el experimento R=100 veces para cada uno de varios valores de
la
     p = IP(escoger la parte pequeña)
que son unos mayores, otros menores que la proporción = 1/10 del tamaño de esa
parte.
Produce también:
   * valores teóricos: la IE(X) = N(1 + 1/2 + ... + 1/N)  para el álbum
uniforme, y los dos valores esperados  IE(Xc) , IE(Xd)  de los números de cromos 
Xc, Xd  que habrá que comprar hasta completar cada parte; nótese que el cálculo
de estos valores esperados es como el anterior: sólo hay que incluir el factor 
p  ó  1-p  en la "IP(éxito)" ;
   * estadísticas: la media  m  de la variable  X  en cada serie de R réplicas; 
   * gráficos: los valores de  X  en esas R réplicas, comparados con los valores 
IE(Xc) , IE(Xd)  para ese valor de  p ,  y con el valor observado de  m .

%------------ PARAMETROS:
R= 100; 						% R réplicas de cada experimento
N= 50; c= 1/10; cN= ceil(c*N); dN= N-cN;		% tamaño de la colección y de
								%sus dos partes
Valor_esperado_de_X_con_IP_uniforme = sum(N./(1:N))

%------------ Este bucle externo es para variar el valor de	p = IP(escoger
la parte pequeña)
for k=1:6
	p= c*2.4/k; 	% Tras algunos tanteos (ejecutando el código) se han
				% escogido estos valores;
				% son pocos para permitir mostrar los resultados en un
				% sólo gráfico.
	EXc = sum(cN./(1:cN)) /p; EXd = sum(dN./(1:dN)) /(1-p);	% valores
									%esperados de  Xc , Xd

	f=ones(N,R); s=sum(f); X=0*s; F=s;
				% Valores iniciales: f==1 son los s cromos que faltan; 
				%   X, el número comprado; F es el historial de s.
	%------------- Este bucle produce las R réplicas de la simulación:
	while sum(s)>0
		e= (rand(1,R)<p); 	% Escoge una parte u otra para generar el
						% índice de cromo:
		r= rand(1,R); i= e.*ceil(r*cN) + (1-e).*ceil(cN+r*dN);
		f(i+N*(0:R-1))= 0;   		% Lo anota como ya adquirido, haciendo
							% =0 cada entrada  i(k),k  de f ; 
							% luego actualiza  los valores de las
							%otras variables:
		X=X+(s>0); s=sum(f); F=[F;s];
	end
 	 %------------- 
	m= mean(X);				% valor medio observado; 
	subplot(3,2,k); 				% 		OUTPUT GRAFICO:
	plot(sort(X),1:R,0,0,[EXc,EXd],R/2,'+r',m,R/2,'ob'); 
	gtext( ['p/c = ',num2str(p/c),'    media = ',num2str(m)] , 'Fontsize',14)

	Ec(k)= EXc; Ed(k)= EXd; mX(k)= m; pc1(k)= 100*p;
				% Se guardan valores para ...
end				%   ... mostarlos en la CW al terminar (redondeados):
Percent_1___EXc___EXd___media_observada = round([pc1;Ec;Ed;mX])'


%------------------------------------- 	FIN DEL CODIGO.	Lo que sigue es su output en la CW :

Valor_esperado_de_X_con_IP_uniforme =

  224.9603

Percent_1___EXc___EXd___media_observada =

    24    48   260   257
    12    95   225   232
     8   143   215   230
     6   190   210   248
     5   238   208   279
     4   285   206   298


%--------- COMENTARIOS: 

Puede observarse cómo la media es mínima al pasar por  p=c  , y se hace mayor
con ambos tipos de desequilibrio, pero "más rápidamente" cuando  p<c  (por
ejemplo, es mayor para  p=c/2  que para  p=2c ).

También puede verse que ese valor es algo mayor que los otros dos (es fácil
probar que su valor esperado es mayor que ambos), pero que casi coincide con el
mayor de ellos cuando son muy distintos. 
Esto último es bastante lógico: si una parte es "mucho más difícil" que la otra,
casi siempre se acaba el total al acabar la parte difícil.

Los gráficos muestran todo eso, y algo más: la distribución empírica de  X  en
cada caso.



----------------------- APENDICE -----------------------------
UN POSIBLE CODIGO PARA PERMUTACIONES Y MUESTRAS SIN REPETICION

%--------- Saca una muestra de tamaño n del conjunto {1,...,N}
N= 50; n= 5;  		% parámetros; si n = N, sale una permutación
I=1:n; s=1:N; 		% s  son los valores que se van a permutar
			% al final la muestra será  s(I)
for i=I, 	
	k= fix(rand*(N+1-i)); s([i,i+k])= s([i+k,i]); 
	% intercambia s(i) con s(i+k), i+k escogido al azar en [i,N]
end
disp(s(I)')		% para pruebas: sólo si n es pequeña !!!!


%--------- Saca R muestras de tamaño n del conjunto {1,...,N}
N= 50; R= 400; n= 50; 
I=1:n; u=ones(1,R); s=(1:N)'*u;
for i=I
	k= fix(rand(1,R)*(N+1-i)); 
	for j=1:R
	   s([i,i+k(j)],j)= s([i+k(j),i],j);	% cada  i+k(j)  está en [i,N]
	end
end
%--------- Cuenta inversiones e índices fijos en esas R muestras
inv=zeros(1,R); cf=inv; 
for i=1:n
	cf= cf+( s(i,:)==i*u ); 	% cf cuenta índices fijos
	i1=1:i-1; inv= inv + sum( ones(i-1,1)*s(i,:) < s(i1,:) ,1 );
end

%--------------- ESTADISTICAS:

	%---------- del no. de inversiones
disp('IE y media observada del número de inversiones')
disp([N*(N-1)/4, mean(inv)]) 
plot( sort(inv),1:R, mean(inv),R/2,'+r' ); 
title(['N = ',num2str(N),'  , n = ',num2str(n)],'Fontsize',14)
xlabel('No. de inversiones','Fontsize',14)
ylabel('No. de repeticiones','Fontsize',14)

	%---------- del signo
disp('Media observada del signo, y su SD teórica:')
disp([mean((-1).^inv), 1/sqrt(R)])

	%---------- del no. de puntos fijos; sólo interesa si n=N !!!!!
disp('Media del número de fijos:')
disp(mean(cf))
v=0:max(cf); f=exp(1);
for k=v, 	% divide las frecuencias de cf por las Poisson_1,
			%   que es su aproximación para N grande
	q(k+1)= mean(cf==k)*f; f=f*(k+1); 
end
disp('frecuencia(no. de puntos fijos = k) / poisson_1(k)')	
disp('    k       cociente')
disp([v;q(v+1)]')