close Warning: Can't synchronize with repository "(default)" (/var/svn/tolp does not appear to be a Subversion repository.). Look in the Trac log for more information.

Ticket #432: schur.tol

File schur.tol, 6.4 KB (added by Jorge, 19 years ago)

for bug 432

Line 
1//////////////////////////////////////////////////////////////////////////////
2// FILE   : schur.tol
3// PURPOSE: Funciones para la generacion aleatoria de polinomios de Schur
4//          bajo leyes normales multivariantes
5//////////////////////////////////////////////////////////////////////////////
6
7//////////////////////////////////////////////////////////////////////////////
8// INCLUDES
9//////////////////////////////////////////////////////////////////////////////
10//Set Include(PathIniSadd(0));
11
12//////////////////////////////////////////////////////////////////////////////
13// FUNCTIONS
14//////////////////////////////////////////////////////////////////////////////
15
16//////////////////////////////////////////////////////////////////////////////
17// PROCEDURES
18//////////////////////////////////////////////////////////////////////////////
19
20//////////////////////////////////////////////////////////////////////////////
21Real SVDDet(Matrix cov)
22//////////////////////////////////////////////////////////////////////////////
23{
24  Set svd = SVD(cov);
25  Matrix d = svd[2];
26  MatProd(SubDiag(d,0))
27};
28
29
30//////////////////////////////////////////////////////////////////////////////
31Polyn SchurRestrictionPolyn(Polyn p, Real k)
32//////////////////////////////////////////////////////////////////////////////
33{
34  Real n = Degree(p);
35  If(Or(EQ(k,n), GT(k,n), LT(k,0)), p,
36  {
37    Polyn p_ = ChangeBF(p)*(B^Degree(p));
38    Real p_0 = EvalPol(p_, 0);
39    Real p0 = EvalPol(p,0);
40 
41    Polyn sAux = p_0*p-p0*p_;
42    Polyn s = Quotient(sAux/B);
43    SchurRestrictionPolyn(s, k) 
44  })
45};
46//////////////////////////////////////////////////////////////////////////////
47PutDescription("Genera el polinomio de Schur de grado k de un polinomio p.
48El grado k no puede ser mayor que el grado del polinomio p y no inferior a 0",
49SchurRestrictionPolyn);
50//////////////////////////////////////////////////////////////////////////////
51
52//////////////////////////////////////////////////////////////////////////////
53Matrix SchurReduxMatrix(Polyn p)
54//////////////////////////////////////////////////////////////////////////////
55{
56  Real n         = Degree(p);
57  Real a0        = Coef(p,0);
58  Matrix In_1    = Diag(n-1, 1);
59  Matrix In_1Est = Diag(n-1, 1, 0);
60  Matrix redux = RProd(In_1-RProd(In_1Est, a0), (1-a0^2)^(-1));
61  redux
62};
63//////////////////////////////////////////////////////////////////////////////
64PutDescription("Genera la matriz de restricciones k-esima para aplicarlas a
65los coeficientes de un polinomio y obtener el polinomio k-esimo de Schur.",
66SchurReduxMatrix);
67//////////////////////////////////////////////////////////////////////////////
68
69//////////////////////////////////////////////////////////////////////////////
70Polyn SchurReduxPolyn(Polyn p)
71//////////////////////////////////////////////////////////////////////////////
72{
73  Real n          = Degree(p);
74  Matrix m        = SubRow(PolMat(p, n, 1), Range(2, n, 1));
75  Matrix matRedux = SchurReduxMatrix(p);
76  MatPol(Tra(matRedux*m))+B^(n-1)
77};
78
79//////////////////////////////////////////////////////////////////////////////
80Polyn SchurNormPolyn(Polyn p)
81//////////////////////////////////////////////////////////////////////////////
82{
83  p/Coef(p, Degree(p))
84};
85
86//////////////////////////////////////////////////////////////////////////////
87Set SchurDraw(Polyn p, Matrix cov)
88//////////////////////////////////////////////////////////////////////////////
89{
90  Real n = Degree(p);
91  Text WriteLn("Grado "<<n);
92
93  If(LE(n,1), SetOfReal(DrawTruncatedNormal(Coef(p,0), MatDat(cov,1,1), -1, 1)) ,
94  {
95    Matrix nuVector = PolMat(p, n, 1);
96    Real nu         = MatDat(nuVector, 1, 1);
97    Matrix cov22 = Sub(cov, 2, 2, n-1, n-1);
98    Matrix cov12 = Sub(cov, 2, 1, 1, n-1);
99    Matrix cov11 = Sub(cov, 1, 1, 1, 1);
100    Real sigma = SqRt(SVDDet(cov22)/SVDDet(cov));
101    Real a0 = DrawTruncatedNormal(nu, sigma, -1, 1);
102
103    Matrix nu2 = SubRow(nuVector, Range(2, n, 1));
104    Matrix redux = SchurReduxMatrix(p-Coef(p,0)+a0);
105    Matrix newPol =
106     redux*(nu2+Tra(cov12)*InverseDiag(cov11)*(Diag(1,a0)-Diag(1,nu)));
107    Matrix newCov =
108     Tra(redux)*(cov22+Tra(cov12)*InverseDiag(cov11)*cov12)*redux;
109    SetOfReal(a0)<<SchurDraw(MatPol(Tra(newPol))+B^(n-1), newCov)
110  })
111};
112//////////////////////////////////////////////////////////////////////////////
113PutDescription("Muestrea un polinomio de Schur mediante el metodo Almagro
114condicionado a una media dada por p (polinomio normalizado en el grado maximo)
115y su matriz de covarianzas cov. Los parametros que se muestrean seran los
116prefijados como distintos de cero en el polinomio p",
117SchurDraw);
118//////////////////////////////////////////////////////////////////////////////
119
120
121
122/*
123//////////////////////////////////////////////////////////////////////////////
124Polyn DrawSchurPolyn(Matrix beta, Matrix S, Real d, Real p)
125//////////////////////////////////////////////////////////////////////////////
126{
127  Real a0 = 1;
128  Polyn p = 
129  Set coefSet = For(1, d-1, Real(Real k)
130  {
131   
132
133
134  }); 
135
136
137
138};
139//////////////////////////////////////////////////////////////////////////////
140PutDescription("Genera un polinomio de Schur bajo la hipotesis de que sus
141coeficientes siguen una normal multivariante de media beta y covarinza S. El
142polinomio es de la forma 1-a1*B^p-a2*B^(2*p)-...-ad*B^(d*p)",
143DrawSchurPolyn);
144//////////////////////////////////////////////////////////////////////////////
145*/
146
147
148
149Polyn pol = 0.1-2*B^6+B^7;
150Matrix cov = SetDiag(SetOfReal(.5,0.00001,0.00001,0.00001,0.00001,0.00001,.5));
151 Set a = SchurDraw(pol,cov);
152
153
154
155/*
156  Polyn p = pol;
157  Real n = Degree(p);
158
159    Matrix nuVector = PolMat(p, n, 1);
160    Real nu    = MatDat(nuVector, 1, 1);
161    Matrix cov22 = Sub(cov, 2, 2, n-1, n-1);
162    Matrix cov12 = Sub(cov, 2, 1, 1, n-1);
163    Matrix cov11 = Sub(cov, 1, 1, 1, 1);
164    Real sigma = SqRt(SVDDet(cov22)/SVDDet(cov));
165    Real a0 = DrawTruncatedNormal(nu, sigma, -1, 1);
166
167    Matrix nu2 = SubRow(nuVector, Range(2, n, 1));
168    Matrix redux = SchurReduxMatrix(p-Coef(p,0)+a0);
169    Matrix newPol =
170     redux*(nu2+Tra(cov12)*InverseDiag(cov11)*(Diag(1,a0)-Diag(1,nu)));
171    Matrix newCov =
172     Tra(redux)*(cov22+Tra(cov12)*InverseDiag(cov11)*cov12)*redux;
173    Polyn MatPol(newPol)
174//    Set SetOfReal(a0)<<SchurDraw(MatPol(newPol), newCov)
175 
176*/
177
178/*
179  Real n    = Degree(pol);
180  Polyn pol2 = pol-Coef(pol,0);
181  Matrix m = PolMat(pol2, n, 1);
182
183Matrix SchurReduxMatrix(pol);
184Polyn SchurReduxPolyn(pol);
185Polyn SchurNormPolyn(SchurRestrictionPolyn(pol, 6));
186*/
187/*
188Real IsStationary(1-2*B-3*B^2);
189Real IsStationary(1-2*B^6-3*B^7);
190*/