/**************************************/
/* Apartado (a) :  p=3  a desconocido */
/**************************************/

X: [1.02, 0.87, 0.89, 0.56, 0.39, 1.05, 0.66, 0.64, 0.61, 0.78, 0.50,
0.58, 0.80, 1.45, 1.05, 0.60, 0.54, 0.25, 0.72, 1.23, 0.35, 0.53,
0.61, 0.78, 1.04, 1.23, 0.47, 0.51, 1.45, 1.05, 0.66, 0.85, 0.86,
0.41, 1.38, 0.50, 0.17, 0.97, 0.81, 1.12, 1.27, 1.23, 0.78, 0.64,
0.71, 0.47, 0.77, 0.90, 1.21, 0.40, 0.86, 0.82, 0.67, 0.55, 1.04,
0.98, 0.57, 0.84, 1.09, 0.81, 1.20, 1.72, 0.53, 0.54, 1.18, 0.69,
1.79, 1.09, 0.38, 0.50, 1.06, 1.35, 0.69, 0.76, 1.46] ;

/* subapartado 1 : dibuja la verosimilitud y estima a*/

load (distrib);
f_gama (x, p, a) := pdf_gamma (x, p, 1/a);/* densidad de una gama */
f_gama (x, p, a);/* unit_step = función indicatriz de la semirrecta positiva */
assume (p > 0, a > 0);
declare (p, integer);/* para que no pregunte el siguiente integrate */
integrate (f_gama(x,p,a), x, 0, inf);/* por ser densidad, debe dar 1 */
logL(a) := apply ("+", map (lambda ([Xi], log(f_gama(Xi,3,a))), X));/* veros. */
plot2d (logL(a), [a, 0.1, 10]);/* el máximo está cerca de a=4 */
load (lbfgs);
lbfgs (-logL(a), [a], [4], 0.001, [-1,0]);/* -log porque lbfgs minimiza */

/* subapartado 2 : calcula EMV y analiza sus propiedades */

logL(a,x,n) := sum (log(f_gama(x[i],3,a)), i, 1, n);
logexpand: super;/* por omisión, true; sirve para log(prod) -> sum(log) */
simpsum: true;/* por omisión, false; sirve para simplificar un sum */
solve (diff(logL(a,x,n),a) = 0, a);/* a = 3n/sum(x) */
diff(logL(a,x,n),a,2);/* segunda derivada negativa => máximo */

/* suficiencia */
logL(a,x,n);
/* esta parte depende sólo de a y sum(x)
                                 n
                                ====
                                \
        3 log(a) n - log(2) n +  >    - a x 
                                /          i
                                ====
                                i = 1
y ésta sólo de x:                                
                                 n
                                ====
                                \
                                 >    (log(unit_step(x )) + 2 log(x )
                                /                     i            i 
                                ====
                                i = 1
suficiencia mínima: */
assume(x[i]>0, y[i]>0);
diff (logL(a,x,n) - logL(a,y,n), a);/* no depende de a sii sum(x)=sum(y) */

/* consistencia */

/* asintóticamente centrado */
f_sumaX (x, p, a, n) := f_gama (x, n*p, a);/* por reproductividad */
assume (n > 1);/* para que no pregunten los integrate */
declare (n, integer);
integrate (f_sumaX(x,p,a,n), x, 0, inf);/* debe dar: 1*/
E_aMV (p, a, n) := integrate (n*p/x * f_sumaX(x,p,a,n), x, 0, inf);
gamma_expand: true $ /* para que simplique las gamas */
E_aMV (3, a, n);
limit (E_aMV(3,a,n), n, inf);/* debe dar: a */
/* varianza asintóticamente nula */
E_aMV2 (p, a, n) := integrate ((n*p/x)^2 * f_sumaX(x,p,a,n), x, 0, inf);
V_aMV  (p, a, n) := E_aMV2(p,a,n) - E_aMV(p,a,n)^2;
factor (ratsimp (V_aMV (3, a, n)));/* simplificar como razón de polinomios */
limit (V_aMV(3,a,n), n, inf);/* debe dar: 0 */
/* error cuadrático medio = sesgo² + varianza */
sesgo : ratsimp (E_aMV (3, a, n) - a);
ECM : factor (ratsimp (sesgo^2 + V_aMV(3,a,n)));

/* distribución : método delta                                             */
/* aMV = p/X_                                                              */
/* X_ = gama(p.n,a)/n = gama(p.n,a.n) = (TCL) = N (p/a, raíz(p/n)) =       */
/*                                            = N (p/a, 1/a raíz(p/n))     */
/* X_/p = gama(p.n,a.n.p) = (TCL) = N(1/a, 1/(a*raíz(p.n)))                */
/* aMV = p/X_ = g (X_/p)   con   g(x) = 1/x   =>   |g'(x)| = 1/x**2        */
/* aMV = (delta) = N(g(1/a), |g'(1/a)| . 1/(a*raíz(p.n))) =                */
/*     = N (a, a**2 / (a * raíz(p.n))) = N (a, a / raíz(p.n))              */
g(x) := 1/x;
define (gprima(x), diff(g(x),x));
abs (gprima (1/a)) * 1/(a*sqrt(p*n));
/* otro enfoque: aMV = gama-inversa (1/a, 1/(a*raíz(p.n)))     */
/* véase https://fr.wikipedia.org/wiki/Loi_inverse-gamma       */
/* esperanza = npa / (np - 1) -> a  si n -> inf                */
/* varianza  = (npa)**2 / (np-1)**2 / (np-2) -> 0  si n -> inf */

/* subapartado 3 : método de los momentos                     */
/* X = gama (p, a)   esperanza = p/a   aMM = p/media(X) = aMV */


/**************************************/
/* Apartado (b) :  p desconocido  a=5 */
/**************************************/

kill(all);
X: [1.02, 0.87, 0.89, 0.56, 0.39, 1.05, 0.66, 0.64, 0.61, 0.78, 0.50,
0.58, 0.80, 1.45, 1.05, 0.60, 0.54, 0.25, 0.72, 1.23, 0.35, 0.53,
0.61, 0.78, 1.04, 1.23, 0.47, 0.51, 1.45, 1.05, 0.66, 0.85, 0.86,
0.41, 1.38, 0.50, 0.17, 0.97, 0.81, 1.12, 1.27, 1.23, 0.78, 0.64,
0.71, 0.47, 0.77, 0.90, 1.21, 0.40, 0.86, 0.82, 0.67, 0.55, 1.04,
0.98, 0.57, 0.84, 1.09, 0.81, 1.20, 1.72, 0.53, 0.54, 1.18, 0.69,
1.79, 1.09, 0.38, 0.50, 1.06, 1.35, 0.69, 0.76, 1.46] ;
load(distrib);
f_gama (x, p, a) := pdf_gamma (x, p, 1/a);/* densidad de una gama */
logL(p) := apply ("+", map (lambda ([Xi], log(f_gama(Xi,p,5))), X));/* verosimilitud */
plot2d (logL(p), [p, 0.1, 10]);/* el máximo está cerca de p=4 */
lbfgs (-logL(p), [p], [4], 0.001, [-1,0]);/* -log porque lbfgs minimiza */
/* máximo en p=4,3 */
logL(p,x,n) := sum (log(f_gama(x[i],p,5)), i, 1, n);
logexpand: super;/* por omisión, true; sirve para log(prod) -> sum(log) */
simpsum: true;/* por omisión, false; sirve para simplificar un sum */
solve (diff(logL(p,x,n),p) = 0, a);/* no puede encontrar solución: no hay EMV */
/* método de los momentos */
/* X = gama (p, a)   esperanza = p/a   aMM = a.media(X) */
5 * apply("+",X) / length(X); /* 4,166 */
/* a.media(X) = gama (np, na/a) = gama (np, n) -> esperanza = np/n = p => centrado */
/* var(a.media(X)) = a**2 . var(media(X)) = a**2 . var(gama(np,na)) =
   a**2 . np / (na)**2 = p/n -> 0 si n -> inf => consistente */
/* suficiencia mínima: */
assume(x[i]>0, y[i]>0);
diff (logL(p,x,n) - logL(p,y,n), p);
/* = 0  sii  no depende de p  sii  sum(log x)=sum( logy)  sii  prod(x)=prod(y) */
/* luego prod(x) sería el estadígrafo suficiente mínimo para p */
/* y sum(x), media(x), a.media(x) no inducen una subpartición de la de prod(x) */
/* luego no son suficientes */


/************************************/
/* Apartado (c) :  p,a desconocidos */
/************************************/

kill(all);
X: [1.02, 0.87, 0.89, 0.56, 0.39, 1.05, 0.66, 0.64, 0.61, 0.78, 0.50,
0.58, 0.80, 1.45, 1.05, 0.60, 0.54, 0.25, 0.72, 1.23, 0.35, 0.53,
0.61, 0.78, 1.04, 1.23, 0.47, 0.51, 1.45, 1.05, 0.66, 0.85, 0.86,
0.41, 1.38, 0.50, 0.17, 0.97, 0.81, 1.12, 1.27, 1.23, 0.78, 0.64,
0.71, 0.47, 0.77, 0.90, 1.21, 0.40, 0.86, 0.82, 0.67, 0.55, 1.04,
0.98, 0.57, 0.84, 1.09, 0.81, 1.20, 1.72, 0.53, 0.54, 1.18, 0.69,
1.79, 1.09, 0.38, 0.50, 1.06, 1.35, 0.69, 0.76, 1.46] ;
/* método de los momentos */
/* E(X) = p/a  V(X) = p/a**2  =>  pMM = media(X)**2/var(X)  aMM = media(X)/var(X) */
media(x) := apply("+", x) / length(x);
var(x) := media(x**2) - media(x)**2;
pMM: media(X)**2/var(X);
aMM: media(X)/var(X);   
/* máxima verosimilitud */
load(distrib);
f_gama (x, p, a) := pdf_gamma (x, p, 1/a);/* densidad de una gama */
logL(p,a) := apply ("+", map (lambda ([Xi], log(f_gama(Xi,p,a))), X));/* verosimilitud */
/* maximizando */
load(lbfgs);
lbfgs (-logL(p,a), [p,a], [pMM,aMM], 0.001, [-1,0]);/* -log porque lbfgs minimiza */
/* igualando a cero las derivadas */
load(mnewton);
mnewton ([diff(logL(p,a), a) = 0, 
  diff(logL(p,a), p) = 0],
  [p, a], [pMM, aMM]);
/* superficie */
plot3d(logL(p,a), [a,5,7], [p,6,8]);
/* error cuadrático medio */
EC(p,a,X) := block([pMM, aMM, pMV, aMV, sol],
  pMM: media(X)**2/var(X),
  aMM: media(X)/var(X),
  sol: lbfgs (-logL(p_,a_), [p_,a_], [pMM,aMM], 0.001, [-1,0]),
  pMV: rhs(sol[1]),
  aMV: rhs(sol[2]),
  ['pMV = pMV,
  'aMV = aMV,
  'pMM = pMM,
  'aMM = aMM,
  EC_MV = (pMV-p)**2 + (aMV-a)**2,
  EC_MM = (pMM-p)**2 + (aMM-a)**2]) $
azar_gamma(p,a,n) := random_gamma(p,1/a,n);/* parametrización habitual */
ECazar(p,a,n) := block([X],
  X: azar_gamma(p,a,n),
  EC(p,a,X)) $
ECM(p,a,n,N) := block([aux],
  aux: map (media,
    transpose (apply (matrix,
        map (lambda ([lista], map(rhs, lista)),/* elimina etiquetas ECM_MM y ECM_MV */
          map(lambda([N], ECazar(p,a,n)),/* devuelve [ECM_MM = ..., ECM_MV = ...] */
            makelist(i,i,1,N)))))),
  $/* la i no se usa */
 
map (lambda ([x], round (100*x)), ECM(6,5,75,1000)); /* 627, 524, 630, 526, 200, 237 */
map (lambda ([x], round (100*x)), ECM(6,7,75,1000)); /* 626, 735, 630, 740, 272, 318 */
map (lambda ([x], round (100*x)), ECM(7,5,75,1000)); /* 729, 522, 734, 525, 244, 288 */
map (lambda ([x], round (100*x)), ECM(7,7,75,1000)); /* 726, 728, 729, 732, 344, 398 */