| 1 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 2 | Real SVDDet_(Matrix cov) |
|---|
| 3 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 4 | { |
|---|
| 5 | If(Rows(cov)==2, |
|---|
| 6 | { |
|---|
| 7 | MatDat(cov,1,1)*MatDat(cov,2,2)-MatDat(cov,1,2)*MatDat(cov,2,1) |
|---|
| 8 | }, |
|---|
| 9 | { |
|---|
| 10 | Set svd = SVD(cov); |
|---|
| 11 | Matrix d = svd[2]; |
|---|
| 12 | MatProd(SubDiag(d,0)) |
|---|
| 13 | }) |
|---|
| 14 | }; |
|---|
| 15 | |
|---|
| 16 | |
|---|
| 17 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 18 | Real SVDDet(Matrix cov) |
|---|
| 19 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 20 | { |
|---|
| 21 | Set svd = SVD(cov); |
|---|
| 22 | Matrix d = svd[2]; |
|---|
| 23 | MatProd(SubDiag(d,0)) |
|---|
| 24 | }; |
|---|
| 25 | |
|---|
| 26 | |
|---|
| 27 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 28 | Matrix SchurReduxMatrix(Polyn p) |
|---|
| 29 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 30 | { |
|---|
| 31 | Real n_ = Degree(p); |
|---|
| 32 | Real a0_ = Coef(p,0); |
|---|
| 33 | Matrix In_1 = Diag(n_-1, 1); |
|---|
| 34 | Matrix In_1Est = Diag(n_-1, 1, 0); |
|---|
| 35 | Matrix redux_ = RProd(In_1-RProd(In_1Est, a0_), (1-a0_^2)^(-1)); |
|---|
| 36 | redux_}; |
|---|
| 37 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 38 | PutDescription("Genera la matriz de restricciones k-esima para aplicarlas a |
|---|
| 39 | los coeficientes de un polinomio y obtener el polinomio k-esimo de Schur.", |
|---|
| 40 | SchurReduxMatrix); |
|---|
| 41 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 42 | |
|---|
| 43 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 44 | Polyn SchurReduxPolyn(Polyn p) |
|---|
| 45 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 46 | { |
|---|
| 47 | Real n__ = Degree(p); |
|---|
| 48 | Matrix m = SubRow(PolMat(p, n__, 1), Range(2, n__, 1)); |
|---|
| 49 | Matrix matRedux = SchurReduxMatrix(p); |
|---|
| 50 | MatPol(Tra(matRedux*m))+B^(n__-1) |
|---|
| 51 | }; |
|---|
| 52 | |
|---|
| 53 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 54 | Polyn SchurNormPolyn(Polyn p) |
|---|
| 55 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 56 | { |
|---|
| 57 | p/Coef(p, Degree(p)) |
|---|
| 58 | }; |
|---|
| 59 | |
|---|
| 60 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 61 | Set SchurDraw(Polyn p, Matrix cov) |
|---|
| 62 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 63 | { |
|---|
| 64 | Real n. = Degree(p); |
|---|
| 65 | If(LE(n.,1), SetOfReal(DrawTruncatedNormal(Coef(p,0), MatDat(cov,1,1), -1, 1)) , |
|---|
| 66 | { |
|---|
| 67 | Matrix nuVector = PolMat(p, n., 1); |
|---|
| 68 | Real nu = MatDat(nuVector, 1, 1); |
|---|
| 69 | Matrix cov22 = Sub(cov, 2, 2, n.-1, n.-1); |
|---|
| 70 | Matrix cov12 = Sub(cov, 2, 1, 1, n.-1); |
|---|
| 71 | Matrix cov11 = Sub(cov, 1, 1, 1, 1); |
|---|
| 72 | Real sigma = SqRt(SVDDet(cov22)/SVDDet(cov)); |
|---|
| 73 | Real a0. = DrawTruncatedNormal(nu, sigma, -1, 1); |
|---|
| 74 | Matrix nu2 = SubRow(nuVector, Range(2, n., 1)); |
|---|
| 75 | Real difp0a0 = -Coef(p,0)+a0.; |
|---|
| 76 | Matrix redux = SchurReduxMatrix(p+difp0a0); |
|---|
| 77 | Matrix cov11i = InverseDiag(cov11); |
|---|
| 78 | Matrix cov12t = Tra(cov12); |
|---|
| 79 | Matrix cov12t11i= cov12t*cov11i; |
|---|
| 80 | Matrix newPol = redux*(nu2+cov12t11i*(Diag(1,a0.)-Diag(1,nu))); |
|---|
| 81 | //Matrix difa0nu = Diag(1,a0.)-Diag(1,nu); |
|---|
| 82 | //Matrix newPol = redux*(nu2+cov12t11i*difa0nu); |
|---|
| 83 | Matrix newCov = Tra(redux)*(cov22+cov12t11i*cov12)*redux; |
|---|
| 84 | Matrix aux1 = newCov; |
|---|
| 85 | Matrix aux2 = newPol; |
|---|
| 86 | Polyn aux3 = MatPol(Tra(aux2))+B^(n.-1); |
|---|
| 87 | Set aux4 = SchurDraw(aux3, aux1); |
|---|
| 88 | Set aux5 = SetOfReal(a0.); |
|---|
| 89 | Set aux6 = aux5<<aux4; |
|---|
| 90 | aux6 |
|---|
| 91 | }) |
|---|
| 92 | }; |
|---|
| 93 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 94 | PutDescription("Muestrea un polinomio de Schur mediante el metodo Almagro |
|---|
| 95 | condicionado a una media dada por p (polinomio normalizado en el grado maximo) |
|---|
| 96 | y su matriz de covarianzas cov. Los parametros que se muestrean seran los |
|---|
| 97 | prefijados como distintos de cero en el polinomio p", |
|---|
| 98 | SchurDraw); |
|---|
| 99 | ////////////////////////////////////////////////////////////////////////////// |
|---|
| 100 | |
|---|
| 101 | |
|---|
| 102 | |
|---|
| 103 | |
|---|
