Changeset 61e5bcf in tol
- Timestamp:
- Oct 20, 2008, 8:09:43 AM (18 years ago)
- Branches:
- 1941-docker-image-to-test, devel, master
- Children:
- f22ae37
- Parents:
- 26eb0b1
- File:
-
- 1 edited
-
tol/btol/bmodel/estim.cpp (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
tol/btol/bmodel/estim.cpp
r26eb0b1 r61e5bcf 37 37 #include <tol/tol_btxtgra.h> 38 38 #include <tol/tol_bdtegra.h> 39 //#include <contrib/fftw/bfft.h> 39 40 40 41 //#define _USE_GSL_MARQUARDT_ … … 1502 1503 contens_.SetIndexByName(); 1503 1504 } 1505 /* 1506 //-------------------------------------------------------------------- 1507 void Taper_TukeyHanning(double rho, BArray<double>& h) 1508 //-------------------------------------------------------------------- 1509 { 1510 static double pi = BDat::Pi().Value(); 1511 int n = h.Size(); 1512 int l = (int)rho*n; 1513 int t; 1514 h.AllocBuffer(n); 1515 double* H = h.GetBuffer(); 1516 double HS2 = 0; 1517 for(t=0; t<l; t++, H++) 1518 { 1519 *H = 0.5*(1-cos(pi*(t+0.5)))/l; 1520 HS2 += *H * *H; 1521 } 1522 for(t; t<n-l; t++, H++) 1523 { 1524 *H = 1; 1525 HS2 += *H * *H; 1526 } 1527 for(t; t<n; t++, H++) 1528 { 1529 *H = 0.5*(1-cos(pi*(n-t-0.5)))/l; 1530 HS2 += *H * *H; 1531 } 1532 H = h.GetBuffer(); 1533 double HS = sqrt(HS2); 1534 for(t=0; t<n; t++, H++) 1535 { 1536 *H/=HS; 1537 } 1538 } 1539 1540 /* 1541 //-------------------------------------------------------------------- 1542 struct harmonic_power 1543 //-------------------------------------------------------------------- 1544 { 1545 bool calculated_; 1546 double exp_2_pi_i_ 1547 } 1548 1549 1550 //-------------------------------------------------------------------- 1551 complex eval_hasrmonic_pol(const BPol& pol, int T, int k) 1552 //-------------------------------------------------------------------- 1553 { 1554 register int n = pol->Size(); 1555 register int i, deg, nextDeg; 1556 register double y = 0; 1557 register const BMonome<BDat>* aux = pol.Buffer()+(n-1); 1558 deg = aux.Degree(); 1559 z = exp(2*pi*) 1560 for(i=n-1; i>=0; i--) 1561 { 1562 y += aux.Coef(); 1563 aux--; 1564 nextDeg = (i>0)?aux.Degree():0 1565 y *= z^(deg-nextDeg); 1566 deg = nextDeg; 1567 } 1568 return(y); 1569 } 1570 1571 //-------------------------------------------------------------------- 1572 double power_spectrum(const BPol& ar, const BPol& ma, int n, int k) 1573 //-------------------------------------------------------------------- 1574 { 1575 1576 } 1577 1578 //-------------------------------------------------------------------- 1579 DeclareContensClass(BSet, BSetTemporary, BSetARIMATaperedWhittelEval); 1580 DefExtOpr(1, BSetARIMATaperedWhittelEval, "ARIMATaperedWhittelEval", 2, 4, "Ratio arma Matrix Real Matrix", 1581 "(Set arima, Matrix output [, Real sigma=1.0, Matrix taper=TukeyHanningTaper(2/Sqrt(T))])", 1582 I2("Calculates the tapered Whittle aproximation to ARIMA log-likelihood upon on spectrañ density.\n" 1583 "See more on http://epub.wu-wien.ac.at/dyn/virlib/wp/mediate/epub-wu-01_9dc.pdf?ID=epub-wu-01_9dc" 1584 "NOTE: If output matrix is well differenced arima must have difference " 1585 "polynomial 1 in all seasonalities.\n" 1586 "If all works fine, returns following information:\n", 1587 "Calcula el logaritmo de la verosimilitud aproximada de un modelo ARIMA mediante el método " 1588 "de Wittle disminuido (tapered) basado en la densidad espectral.\n" 1589 "Ver más detalles en http://epub.wu-wien.ac.at/dyn/virlib/wp/mediate/epub-wu-01_9dc.pdf?ID=epub-wu-01_9dc" 1590 "NOTA: Si la matriz de output ya está diferenciada hay que introducir" 1591 "en arima el polinomio de diferencias 1 en todas las estacionalidades.\n" 1592 "Si todo funciona correctamente devuelve lo siguiente:\n")+ 1593 "Matrix z:"+ I2("Noise", 1594 "Ruido")+"\n"+ 1595 "Matrix w:"+ I2("Differenced noise", 1596 "Ruido diferenciado")+"\n"+ 1597 "Real logLH:"+ I2("-Sum {Log(S(j)) + (1/(2*PI)*I(j)/S(j)); j=1...T} : Logarithm of spectral density function", 1598 "-Sum {Log(S(j)) + (1/(2*PI)*I(j)/S(j)); j=1...T} : Logaritmo de la función de densidad espectral")+"\n"+ 1599 BOperClassify::Sthocastic_); 1600 //-------------------------------------------------------------------- 1601 void BSetARIMATaperedWhittelEval::CalcContens() 1602 //-------------------------------------------------------------------- 1603 { 1604 BInt i, p, q; 1605 BRat& arma = Rat(Arg(1)); 1606 BMat& z = Mat(Arg(2)); 1607 BDat sigma = (Arg(3))?Dat(Arg(3)):1.0; 1608 BArray<BARIMAFactor> factor; 1609 factor.ReallocBuffer(s.Card()); 1610 1611 for(i=p=q=0; i<s.Card(); i++) 1612 { 1613 BSet& si = Set(s[i+1]); 1614 factor(i).s_ = (BInt)Real(si.Field("Periodicity")); 1615 factor(i).dif_ = Pol(si.Field("DIF")); 1616 factor(i).ar_ = Pol(si.Field("AR")); 1617 factor(i).ma_ = Pol(si.Field("MA")); 1618 p += factor(i).ar_.Degree(); 1619 q += factor(i).ma_.Degree(); 1620 } 1621 if(!(p+q)) 1622 { 1623 } 1624 BARIMA model; 1625 model.PutFactors(factor); 1626 model.PutOutputData(z); 1627 model.CalcLikelihood_Almagro(sigma); 1628 1629 BInt N_ = model.w_.Rows(); 1630 BList* aux = NULL; 1631 BList* lst = NULL; 1632 LstFastAppend(lst,aux,NewFND(BMat, z, model.z_)); 1633 LstFastAppend(lst,aux,NewFND(BMat, w, model.w_)); 1634 LstFastAppend(lst,aux,NewFND(BDat, ACOVDetN, model.aCovDetN_)); 1635 LstFastAppend(lst,aux,NewFND(BDat, wtACOViw, model.wtCoviw_)); 1636 LstFastAppend(lst,aux,NewFND(BDat, aS2, model.ata_)); 1637 LstFastAppend(lst,aux,NewFND(BDat, logLH, model.logLikelihood_)); 1638 contens_.RobElement(lst); 1639 contens_.SetIndexByName(); 1640 } 1641 */ 1642 1504 1643 1505 1644 //--------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
