diff --git a/dist/prior_corrunif.m b/dist/prior_corrunif.m index 0fabe324..d5085eb7 100644 --- a/dist/prior_corrunif.m +++ b/dist/prior_corrunif.m @@ -35,7 +35,7 @@ % Copyright (c) 2000-2001,2010 Aki Vehtari % Copyright (c) 2009,2015 Jarno Vanhatalo % Copyright (c) 2010 Jaakko Riihimäki -% ------------ 2015 Marcelo Hartmann +% ───────────── 2015, 2018 Marcelo Hartmann % This software is distributed under the GNU General Public % License (version 3 or later); please refer to the file @@ -93,7 +93,7 @@ % set functions p.fh.pak = @prior_corrunif_pak; p.fh.unpak = @prior_corrunif_unpak; - % p.fh.RealToRho = @prior_corrunif_RealToRho; + % p.fh.RealToRho = @prior_corrunif_RealToRho; p.fh.lp = @prior_corrunif_lp; p.fh.lpg = @prior_corrunif_lpg; p.fh.recappend = @prior_corrunif_recappend; diff --git a/gp/demo_epinf.m b/gp/demo_epinf.m index d49e765a..d1c997f4 100644 --- a/gp/demo_epinf.m +++ b/gp/demo_epinf.m @@ -1,12 +1,21 @@ -%DEMO_EPINF Demonstration of input dependent noise and magnitude model +%DEMO_EPINF Demonstration of input dependent noise and magnitude model with EP % -% Description -% Uses toy data sets to demonstrate heteroscedastic noise with input -% dependent noise and magnitude model. Compare with results from -% DEMO_INPUTDEPENDENTNOISE. +% Description +% Uses toy data sets to demonstrate heteroscedastic noise with +% input dependent noise and magnitude model using EP +% inference. Compare with results from +% DEMO_INPUTDEPENDENTNOISE. % -% See also -% DEMO_INPUTDEPENDENTNOISE, LIK_EPGAUSSIAN +% Reference: +% Ville Tolvanen, Pasi Jylänki and Aki Vehtari (2014). Expectation +% propagation for nonstationary heteroscedastic Gaussian process +% regression. In Machine Learning for Signal Processing (MLSP), +% 2014 IEEE International Workshop on, +% DOI:10.1109/MLSP.2014.6958906. +% +% +% See also +% DEMO_INPUTDEPENDENTNOISE, DEMO_STOCHASTICVOLATILITY, LIK_EPGAUSSIAN % % % Copyright (c) Ville Tolvanen 2011-2014 diff --git a/gp/demo_inputdependentnoise.m b/gp/demo_inputdependentnoise.m index 71baf42e..9b291f37 100644 --- a/gp/demo_inputdependentnoise.m +++ b/gp/demo_inputdependentnoise.m @@ -1,11 +1,20 @@ %DEMO_INPUTDEPENDENTNOISE Demonstration of input dependent-noise -% model using Gaussian process prior +% model using Gaussian process prior and Laplace approx % % Description % Uses toy data sets to demonstrate how inferring % heteroscedastic noise with input dependent noise model % differs from standard noise models (Gaussian, Student-t). % +% Reference: +% Ville Tolvanen, Pasi Jylänki and Aki Vehtari (2014). Expectation +% propagation for nonstationary heteroscedastic Gaussian process +% regression. In Machine Learning for Signal Processing (MLSP), +% 2014 IEEE International Workshop on, +% DOI:10.1109/MLSP.2014.6958906. +% +% +% See also DEMO_EPINF, DEMO_STOCHASTICVOLATILITY % % Copyright (c) Ville Tolvanen 2011-2012 diff --git a/gp/demo_multivariategp.m b/gp/demo_multivariategp.m index d141dbda..ce9a8ba3 100644 --- a/gp/demo_multivariategp.m +++ b/gp/demo_multivariategp.m @@ -42,7 +42,7 @@ % and that N1(t)|f1(t) is independent of N2(t)|f2(t). The expression exp(f1(t)) % is the expected value of N1 given f1 and exp(f2(t)) is the expected value % of N2 given f2. The correlation is now introduced through the multivariate -% Gaussian process model assuming the linear model of coregionalization, i.e., +% Gaussian process prior assuming the linear model of coregionalization, i.e., % % [f1(t), f2(t')] ~ MGP(0, K) % @@ -132,9 +132,12 @@ 'R_prior', prior_corrunif('nu', 3), 'corrFun', {k1 k2}); % set gp structure +% for the EP approximation gp = gp_set('lik', lik, 'cf', k, 'latent_method', 'EP'); -% gp = gp_set('lik', lik, 'cf', k, 'latent_method', 'Laplace'); +% For the Laplace approximation +%gp = gp_set('lik', lik, 'cf', k, 'latent_method', 'Laplace'); + % optimzation options opt = optimset('TolFun', 1e-3, 'TolX', 1e-5, 'Display', 'iter'); @@ -150,8 +153,10 @@ % take the parameters, hyperparameters estimates [th ss] = gp_pak(gp); -% correlation estimate -rho = k.fh.RealToRho(th(end - 2 - 2), 1, []) +% correlation matrix estimate +%rho = k.fh.RealToRho(th(end - 2 - 2), 1, []) +Corr = gp.cf{1}.fh.sigma(gp.cf{1}, 'corr'); +Corr{1} % prediction np = 200; diff --git a/gp/demo_regression_additive1.m b/gp/demo_regression_additive1.m index 80349d59..0d5d43cb 100644 --- a/gp/demo_regression_additive1.m +++ b/gp/demo_regression_additive1.m @@ -119,7 +119,7 @@ warning('GPstuff:SuiteSparseMissing',... ['SuiteSparse is not properly installed. (in BECS try ''use suitesparse'')\n' ... 'Using gpcf_sexp (non-compact support) instead of gpcf_ppcs2 (compact support)']); - gpcf2 = gpcf_sexp('lengthScale', 5, 'magnSigma2', 1, 'lengthScale_prior', pl2, 'magnSigma2_prior', pm); + gpcf2 = gpcf_sexp('lengthScale', 5, 'magnSigma2', 1, 'lengthScale_prior', pl2, 'magnSigma2_prior', pm2); end lik = lik_gaussian('sigma2', 0.1, 'sigma2_prior', pn); diff --git a/gp/demo_stochasticvolatility.m b/gp/demo_stochasticvolatility.m new file mode 100644 index 00000000..6cfedf16 --- /dev/null +++ b/gp/demo_stochasticvolatility.m @@ -0,0 +1,125 @@ +%DEMO_STOCHASTICVOLATILITY Demonstration of input dependent-noise +% model using Gaussian process prior, which corresponds to a +% stochastic volatility model +% +% Description +% Demonstrates stochastic volatility model for S&P 500 +% weekly closing value. +% +% Reference: +% Ville Tolvanen, Pasi Jylänki and Aki Vehtari (2014). Expectation +% propagation for nonstationary heteroscedastic Gaussian process +% regression. In Machine Learning for Signal Processing (MLSP), +% 2014 IEEE International Workshop on, +% DOI:10.1109/MLSP.2014.6958906. +% +% +% See also DEMO_INPUTDEPENDENTNOISE, DEMO_EPINF +% +% Copyright (c) Aki Vehtari 2015 + +% This software is distributed under the GNU General Public +% License (version 3 or later); please refer to the file +% License.txt, included with the software, for details. + + +%% Data +% Stochastic volatility model for change in +% S&P 500 stock index weekly closing value from 2001-01-02 to 2014-12-29 +S = which('demo_stochasticvolatility'); +L = strrep(S,'demo_stochasticvolatility.m','demodata/sp500_weekly.csv'); +d=dataset('File',L,'Delimiter',','); +d.Datenum=datenum(d.Date); +% order from oldest to newest +x=flipud(d.Datenum); +x=x(2:end)-x(1); +% normalise +[xn,nd.xmean,nd.xstd]=normdata(x); +% order from oldest to newest +y=flipud(d.Close); +% change in closing price +y=diff(y); +% normalise +[yn,nd.ymean,nd.ystd]=normdata(y); +% +n=numel(y); + +%% Covariance functions +pl = prior_t('s2',1); +pm = prior_t('s2',1); +gpcf1 = gpcf_sexp('lengthScale', 1, 'magnSigma2', 0.5, ... + 'lengthScale_prior', pl, 'magnSigma2_prior', pm); +gpcf2 = gpcf_exp('lengthScale', 1, 'magnSigma2', 0.1, ... + 'lengthScale_prior', pl, 'magnSigma2_prior', pm); + +%% Inference using Laplace approximate + +% Create the likelihood structure. Don't set prior for sigma2 if covariance +% function magnitude for noise process has a prior. +lik = lik_inputdependentnoise('sigma2', 0.1, 'sigma2_prior', prior_fixed()); + +% NOTE! if multiple covariance functions per latent is used, define +% gp.comp_cf as follows: +% gp = gp_set(..., 'comp_cf' {[1 2] [5 6]}; +gp = gp_set('lik', lik, 'cf', {gpcf1 gpcf2}, 'jitterSigma2', 1e-9, 'comp_cf', {[1] [2]}); + +% Set the approximate inference method to Laplace +gp = gp_set(gp, 'latent_method', 'Laplace'); +% For more complex problems, maxiter in latent_opt should be increased. +% gp.latent_opt.maxiter=1e6; + +% Set the options for the optimization +opt=optimset('TolFun',1e-3,'TolX',1e-3,'Derivativecheck','off','Display','iter'); +% Optimize with the scaled conjugate gradient method +gp=gp_optim(gp,xn,yn,'opt',opt); + +% make prediction to the data points +[Ef, Varf,lpyt, Ey, Vary] = gp_pred(gp, xn, yn); +Ef11=Ef(1:n);Ef12=Ef(n+1:end); +Varf11=diag(Varf(1:n,1:n)); +prctmus=[Ef11-1.645*sqrt(Varf11) Ef11 Ef11+1.645*sqrt(Varf11)].*nd.ystd+nd.ymean; +prctys=[Ey-1.645*sqrt(Vary) Ey Ey+1.645*sqrt(Vary)].*nd.ystd+nd.ymean; + +figure; +a4 +% plot mean and 5% and 95% quantiles +subplot(2,1,1) +plot(x,prctmus(:,2),'b',x,prctmus(:,1),'r',x,prctmus(:,3),'r',x,prctys(:,1),'r--',x,prctys(:,3),'r--',x,y,'.') +datetick('x',11) +title('Input dependent noise model -- Laplace approximation'); +legend('Mean','5% f','95% f','5% y','95% y','location','southwest') + + +%% Inference using EP + +% note that lik_epgaussian could be used to include input dependent +% magnitude, too, see demo_epinf +sigma2=1; +lik = lik_epgaussian('sigma2', sigma2, 'sigma2_prior', prior_fixed(), ... + 'int_likparam', true, 'inputparam', true); + +% Set latent options +latent_opt = struct('maxiter',1000, 'df',0.8, 'df2',0.6, 'tol',1e-6, ... + 'parallel', 'on', 'init_prev','off', 'display','off'); + +% NOTE! if multiple covariance functions per latent is used, define +% gp.comp_cf as follows: +% gp = gp_set(..., 'comp_cf' {[1 2] [5 6]}; +gp2 = gp_set('lik', lik, 'cf', {gpcf1 gpcf2}, ... + 'jitterSigma2', 1e-9, 'comp_cf', {1 2}, ... + 'latent_method', 'EP', 'latent_opt', latent_opt); + +gp2 = gp_optim(gp2,xn,yn,'opt',opt, 'optimf', @fminscg); +[Ef, Varf, lpyt, Ey, Vary] = gp_pred(gp2, xn, yn); +prctmus=[Ef(:,1)-1.645*sqrt(Varf(:,1)) ... + Ef(:,1) ... + Ef(:,1)+1.645*sqrt(Varf(:,1))].*nd.ystd+nd.ymean; +prctys=[Ey(:,1)-1.645*sqrt(Vary(:,1)) ... + Ey(:,1) ... + Ey(:,1)+1.645*sqrt(Vary(:,1))].*nd.ystd+nd.ymean; + +subplot(2,1,2) +plot(x,prctmus(:,2),'b',x,prctmus(:,1),'r',x,prctmus(:,3),'r',x,prctys(:,1),'r--',x,prctys(:,3),'r--',x,y,'.') +datetick('x',11) +title('Input dependent noise model -- EP approximation'); +legend('Mean','5% f','95% f','5% y','95% y','location','southwest') diff --git a/gp/demodata/sp500_weekly.csv b/gp/demodata/sp500_weekly.csv new file mode 100644 index 00000000..3b234d49 --- /dev/null +++ b/gp/demodata/sp500_weekly.csv @@ -0,0 +1,731 @@ +Date,Open,High,Low,Close,Volume,AdjClose +2014-12-29,2087.63,2093.55,2057.94,2058.90,2499570000,2058.90 +2014-12-22,2069.28,2092.70,2069.28,2088.77,2391420000,2088.77 +2014-12-15,2005.03,2077.85,1972.56,2070.65,5086390000,2070.65 +2014-12-08,2074.84,2075.78,2002.33,2002.33,3992236000,2002.33 +2014-12-01,2065.78,2079.47,2049.57,2075.37,3657260000,2075.37 +2014-11-24,2065.07,2075.76,2064.75,2067.56,2942725000,2067.56 +2014-11-17,2038.29,2071.46,2034.46,2063.50,3400928000,2063.50 +2014-11-10,2032.01,2046.18,2030.17,2039.82,3234462000,2039.82 +2014-11-03,2018.21,2034.26,2001.01,2031.92,3730468000,2031.92 +2014-10-27,1962.97,2018.19,1951.37,2018.05,3762182000,2018.05 +2014-10-20,1885.62,1965.27,1882.30,1964.58,3589572000,1964.58 +2014-10-13,1905.65,1912.09,1820.66,1886.76,4962132000,1886.76 +2014-10-06,1970.01,1977.84,1906.05,1906.13,4072602000,1906.13 +2014-09-29,1978.96,1985.17,1926.03,1967.90,3761592000,1967.90 +2014-09-22,2009.08,2009.08,1965.99,1982.85,3229072000,1982.85 +2014-09-15,1986.04,2019.26,1978.48,2010.40,3452364000,2010.40 +2014-09-08,2007.17,2007.17,1980.26,1985.54,2946522000,1985.54 +2014-09-02,2004.07,2011.17,1990.10,2007.71,2880167500,2007.71 +2014-08-25,1991.74,2005.04,1990.52,2003.37,2314342000,2003.37 +2014-08-18,1958.36,1994.76,1958.36,1988.40,2562950000,1988.40 +2014-08-11,1933.43,1964.04,1928.29,1955.06,2749490000,1955.06 +2014-08-04,1926.62,1942.92,1904.78,1931.59,3241478000,1931.59 +2014-07-28,1978.25,1984.85,1916.37,1925.15,3483506000,1925.15 +2014-07-21,1976.93,1991.39,1965.77,1978.34,2842770000,1978.34 +2014-07-14,1969.86,1983.94,1955.59,1978.22,3190470000,1978.22 +2014-07-07,1984.22,1984.22,1952.86,1967.57,2938562000,1967.57 +2014-06-30,1960.79,1985.59,1958.22,1985.44,2768790000,1985.44 +2014-06-23,1962.92,1968.17,1944.69,1960.96,3196694000,1960.96 +2014-06-16,1934.84,1963.91,1930.91,1962.87,3250200000,1962.87 +2014-06-09,1948.97,1955.55,1925.78,1936.16,2772774000,1936.16 +2014-06-02,1923.87,1949.44,1915.98,1949.44,2829538000,1949.44 +2014-05-27,1902.01,1924.03,1902.01,1923.57,2965002500,1923.57 +2014-05-19,1876.66,1901.26,1868.14,1900.53,2721034000,1900.53 +2014-05-12,1880.03,1902.17,1862.36,1877.86,3093954000,1877.86 +2014-05-05,1879.45,1889.07,1859.79,1878.48,3222476000,1878.48 +2014-04-28,1865.00,1891.33,1850.61,1881.14,3607606000,1881.14 +2014-04-21,1865.79,1884.89,1859.70,1863.40,3069702000,1863.40 +2014-04-14,1818.18,1869.63,1815.80,1864.85,3336122500,1864.85 +2014-04-07,1863.92,1872.53,1814.36,1815.69,3666776000,1815.69 +2014-03-31,1859.16,1897.28,1859.16,1865.09,3276300000,1865.09 +2014-03-24,1867.67,1875.92,1842.11,1857.62,3355872000,1857.62 +2014-03-17,1842.81,1883.97,1842.81,1866.52,3535628000,1866.52 +2014-03-10,1877.86,1882.35,1839.57,1841.13,3328212000,1841.13 +2014-03-03,1857.68,1883.57,1834.44,1878.04,3502434000,1878.04 +2014-02-24,1836.78,1867.92,1836.78,1859.45,3742346000,1859.45 +2014-02-18,1839.03,1847.50,1824.58,1836.25,3472885000,1836.25 +2014-02-10,1796.20,1841.65,1791.83,1838.63,3348436000,1838.63 +2014-02-03,1782.68,1798.03,1737.92,1797.02,4076028000,1797.02 +2014-01-27,1791.03,1798.77,1770.45,1782.59,3810850000,1782.59 +2014-01-21,1841.05,1849.31,1790.29,1790.29,3936835000,1790.29 +2014-01-13,1841.26,1850.84,1815.52,1838.70,3567970000,1838.70 +2014-01-06,1832.31,1843.23,1823.73,1842.37,3475120000,1842.37 +2013-12-30,1841.47,1849.44,1827.74,1831.37,2615392500,1831.37 +2013-12-23,1822.92,1844.89,1822.92,1841.40,2048590000,1841.40 +2013-12-16,1777.48,1823.75,1767.99,1818.32,3880520000,1818.32 +2013-12-09,1806.21,1811.52,1772.28,1775.32,3217320000,1775.32 +2013-12-02,1806.55,1810.02,1779.09,1805.09,3333712000,1805.09 +2013-11-25,1806.33,1813.55,1800.58,1805.81,2659387500,1805.81 +2013-11-18,1798.82,1804.84,1777.23,1804.76,3162776000,1804.76 +2013-11-11,1769.96,1798.22,1760.64,1798.18,3095290000,1798.18 +2013-11-04,1763.40,1774.54,1746.20,1770.61,3602804000,1770.61 +2013-10-28,1759.42,1775.22,1752.70,1761.64,3535324000,1761.64 +2013-10-21,1745.20,1759.82,1740.50,1759.77,3492870000,1759.77 +2013-10-14,1699.86,1745.31,1692.13,1744.50,3302596000,1744.50 +2013-10-07,1687.15,1703.44,1646.47,1703.20,3226506000,1703.20 +2013-09-30,1687.26,1696.55,1670.36,1690.50,3171168000,1690.50 +2013-09-23,1711.44,1711.44,1687.11,1691.75,3062048000,1691.75 +2013-09-16,1691.70,1729.86,1691.70,1709.91,3731592000,1709.91 +2013-09-09,1656.85,1689.97,1656.85,1687.99,3154566000,1687.99 +2013-09-03,1635.95,1664.83,1633.41,1655.17,3281187500,1655.17 +2013-08-26,1664.29,1669.51,1627.47,1632.97,2739144000,1632.97 +2013-08-19,1655.25,1664.85,1639.43,1663.50,2790186000,1663.50 +2013-08-12,1688.37,1696.81,1652.61,1655.83,3066858000,1655.83 +2013-08-05,1708.01,1709.24,1684.91,1691.42,2982014000,1691.42 +2013-07-29,1690.32,1709.67,1681.86,1709.67,3384048000,1709.67 +2013-07-22,1694.41,1698.78,1676.03,1691.65,3059340000,1691.65 +2013-07-15,1679.59,1693.12,1671.84,1692.09,3122660000,1692.09 +2013-07-08,1634.20,1680.19,1634.20,1680.19,3233274000,1680.19 +2013-07-01,1609.78,1632.07,1604.57,1631.89,2755502500,1631.89 +2013-06-24,1588.77,1620.07,1560.33,1606.28,4078980000,1606.28 +2013-06-17,1630.64,1654.19,1577.70,1592.43,4091850000,1592.43 +2013-06-10,1644.67,1648.69,1608.07,1626.73,3187002000,1626.73 +2013-06-03,1631.71,1646.53,1598.23,1643.38,3631526000,1643.38 +2013-05-28,1652.63,1674.21,1630.74,1630.74,3660690000,1630.74 +2013-05-20,1665.71,1687.18,1635.53,1649.60,3570650000,1649.60 +2013-05-13,1632.10,1667.47,1626.74,1667.47,3395934000,1667.47 +2013-05-06,1614.40,1635.01,1614.21,1633.70,3294078000,1633.70 +2013-04-29,1582.34,1618.46,1581.28,1614.42,3427490000,1614.42 +2013-04-22,1555.25,1592.64,1548.19,1582.24,3450094000,1582.24 +2013-04-15,1588.84,1588.84,1536.03,1555.25,4005162000,1555.25 +2013-04-08,1553.26,1597.35,1548.63,1588.85,3238698000,1588.85 +2013-04-01,1569.18,1573.66,1539.50,1553.28,3398392000,1553.28 +2013-03-25,1556.89,1570.28,1546.22,1569.19,3066520000,1569.19 +2013-03-18,1560.70,1561.56,1538.57,1556.89,3300302000,1556.89 +2013-03-11,1551.15,1563.62,1547.36,1560.70,3614986000,1560.70 +2013-03-04,1518.20,1552.48,1512.29,1551.18,3597796000,1551.18 +2013-02-25,1515.60,1525.84,1485.01,1518.20,3829222000,1518.20 +2013-02-19,1519.79,1530.94,1497.29,1515.60,3920850000,1515.60 +2013-02-11,1517.93,1524.69,1513.61,1519.79,3416520000,1519.79 +2013-02-04,1513.17,1518.31,1495.02,1517.93,3444132000,1517.93 +2013-01-28,1502.96,1514.41,1496.33,1513.17,3780238000,1513.17 +2013-01-22,1485.98,1503.26,1481.16,1502.96,3574670000,1502.96 +2013-01-14,1472.05,1485.98,1463.76,1485.98,3404978000,1485.98 +2013-01-07,1466.47,1472.75,1451.64,1472.05,3600690000,1472.05 +2012-12-31,1402.43,1467.94,1398.11,1466.47,3665237500,1466.47 +2012-12-24,1430.15,1430.15,1401.58,1402.43,2197712500,1402.43 +2012-12-17,1413.54,1448.00,1413.54,1430.15,4108678000,1430.15 +2012-12-10,1418.07,1438.59,1411.88,1413.58,3383768000,1413.58 +2012-12-03,1416.34,1423.73,1398.23,1418.07,3386154000,1418.07 +2012-11-26,1409.15,1419.70,1385.43,1416.18,3390836000,1416.18 +2012-11-19,1359.88,1409.16,1359.88,1409.15,2688502500,1409.15 +2012-11-12,1379.86,1388.81,1343.35,1359.88,3621476000,1359.88 +2012-11-05,1414.02,1433.38,1373.03,1379.85,3602274000,1379.85 +2012-10-31,1410.99,1434.27,1405.95,1414.20,3746493300,1414.20 +2012-10-22,1433.21,1435.46,1403.28,1411.94,3397482000,1411.94 +2012-10-15,1428.75,1464.02,1427.24,1433.19,3692620000,1433.19 +2012-10-08,1460.93,1460.93,1425.53,1428.59,3115478000,1428.59 +2012-10-01,1440.90,1470.96,1439.01,1460.93,3429462000,1460.93 +2012-09-24,1459.76,1463.24,1430.53,1440.67,3394752000,1440.67 +2012-09-17,1465.42,1467.07,1449.98,1460.15,3705514000,1460.15 +2012-09-10,1437.92,1474.51,1428.98,1465.77,4004608000,1465.77 +2012-09-04,1406.54,1437.92,1396.56,1437.92,3564977500,1437.92 +2012-08-27,1411.13,1416.17,1397.01,1406.58,2628268000,1406.58 +2012-08-20,1417.85,1426.68,1398.04,1411.13,2943798000,1411.13 +2012-08-13,1405.87,1418.71,1397.32,1418.16,2824746000,1418.16 +2012-08-06,1391.04,1407.14,1391.04,1405.87,3182784000,1405.87 +2012-07-30,1385.94,1394.16,1354.65,1390.99,3883892000,1390.99 +2012-07-23,1362.34,1389.19,1329.24,1385.97,4031190000,1385.97 +2012-07-16,1356.50,1380.39,1345.07,1362.66,3608082000,1362.66 +2012-07-09,1354.66,1361.54,1325.41,1356.78,3333824000,1356.78 +2012-07-02,1362.33,1374.81,1348.03,1354.68,2801175000,1354.68 +2012-06-25,1334.90,1362.17,1309.27,1362.16,3752304000,1362.16 +2012-06-18,1342.42,1363.46,1324.41,1335.02,4027288000,1335.02 +2012-06-11,1325.72,1343.32,1306.62,1342.84,3715250000,1342.84 +2012-06-04,1278.29,1329.05,1266.74,1325.66,3887776000,1325.66 +2012-05-29,1318.90,1334.93,1277.25,1278.04,4050725000,1278.04 +2012-05-21,1295.73,1328.49,1295.73,1317.82,3765818000,1317.82 +2012-05-14,1351.93,1351.93,1291.98,1295.22,4251866000,1295.22 +2012-05-07,1368.79,1373.91,1343.13,1353.39,3941332000,1353.39 +2012-04-30,1403.26,1415.32,1367.96,1369.10,3833174000,1369.10 +2012-04-23,1378.53,1406.64,1358.79,1403.36,3790184000,1403.36 +2012-04-16,1370.27,1392.76,1365.38,1378.53,3701492000,1378.53 +2012-04-09,1397.45,1397.45,1357.38,1370.26,3818638000,1370.26 +2012-04-02,1408.47,1422.38,1392.92,1398.08,3659032500,1398.08 +2012-03-26,1397.11,1419.15,1391.56,1408.47,3698456000,1408.47 +2012-03-19,1404.17,1414.00,1386.87,1397.11,3682996000,1397.11 +2012-03-12,1370.78,1405.88,1366.69,1404.17,4281244000,1404.17 +2012-03-05,1369.59,1374.76,1340.03,1370.87,3676690000,1370.87 +2012-02-27,1365.20,1378.04,1354.92,1369.63,3782622000,1369.63 +2012-02-21,1361.22,1368.92,1352.28,1365.74,3680180000,1365.74 +2012-02-13,1343.06,1363.40,1340.80,1361.23,3882884000,1361.23 +2012-02-06,1344.32,1354.32,1335.92,1342.64,3861272000,1342.64 +2012-01-30,1316.16,1345.34,1300.49,1344.90,4225678000,1344.90 +2012-01-23,1315.29,1333.47,1306.06,1316.33,4080966000,1316.33 +2012-01-17,1290.22,1315.49,1290.22,1315.38,4121290000,1315.38 +2012-01-09,1277.83,1296.82,1274.55,1289.09,3854788000,1289.09 +2012-01-03,1258.86,1284.62,1258.86,1277.81,3877267500,1277.81 +2011-12-27,1265.02,1269.37,1248.64,1257.60,2257637500,1257.60 +2011-12-19,1219.74,1265.42,1202.37,1265.33,3280102000,1265.33 +2011-12-12,1255.05,1255.05,1209.47,1219.66,4235314000,1219.66 +2011-12-05,1244.33,1267.06,1231.47,1255.19,4034362000,1255.19 +2011-11-28,1158.67,1260.08,1158.67,1244.28,4335660000,1244.28 +2011-11-21,1215.62,1215.62,1158.66,1158.67,3356230000,1158.67 +2011-11-14,1263.85,1264.25,1209.43,1215.65,3865610000,1215.65 +2011-11-07,1253.21,1277.55,1226.64,1263.85,3874182000,1263.85 +2011-10-31,1284.96,1284.96,1215.42,1253.23,4549214000,1253.23 +2011-10-24,1238.72,1292.66,1221.06,1285.09,4912236000,1285.09 +2011-10-17,1224.47,1239.03,1191.48,1238.25,4767664000,1238.25 +2011-10-10,1158.15,1224.61,1158.15,1224.58,4555924000,1224.58 +2011-10-03,1131.21,1171.40,1074.77,1155.46,4514868000,1155.46 +2011-09-26,1136.91,1195.86,1131.07,1131.42,4960282000,1131.42 +2011-09-19,1214.99,1220.39,1114.22,1136.43,5128284000,1136.43 +2011-09-12,1153.50,1220.06,1136.07,1216.01,4913056000,1216.01 +2011-09-06,1173.97,1204.40,1140.13,1154.23,4649140000,1154.23 +2011-08-29,1177.91,1230.71,1170.56,1173.97,4650126000,1173.97 +2011-08-22,1123.55,1190.68,1121.09,1176.80,5309696000,1176.80 +2011-08-15,1178.86,1208.47,1122.05,1123.53,4427032000,1123.53 +2011-08-08,1198.48,1198.48,1101.54,1178.81,3865062000,1178.81 +2011-08-01,1292.59,1307.38,1168.09,1199.38,5268348000,1199.38 +2011-07-25,1344.32,1344.32,1282.86,1292.28,4207194000,1292.28 +2011-07-18,1315.94,1347.00,1295.92,1345.02,4110088000,1345.02 +2011-07-11,1343.31,1343.31,1306.51,1316.14,4153686000,1316.14 +2011-07-05,1339.59,1356.48,1330.92,1343.80,3737600000,1343.80 +2011-06-27,1268.44,1341.01,1267.53,1339.67,3901108000,1339.67 +2011-06-20,1271.50,1298.61,1262.87,1268.45,3977604000,1268.45 +2011-06-13,1271.31,1292.50,1258.07,1271.50,4093202000,1271.50 +2011-06-06,1300.26,1300.26,1268.28,1270.98,3710360000,1270.98 +2011-05-31,1331.10,1345.20,1297.90,1300.16,4051132500,1300.16 +2011-05-23,1333.07,1334.62,1311.80,1331.10,3519106000,1331.10 +2011-05-16,1334.77,1346.82,1318.51,1333.27,3902876000,1333.27 +2011-05-09,1340.20,1359.44,1332.03,1337.77,3907822000,1337.77 +2011-05-02,1365.21,1370.58,1329.17,1340.20,4072744000,1340.20 +2011-04-25,1337.14,1364.56,1331.47,1363.61,3523530000,1363.61 +2011-04-18,1313.35,1337.49,1294.70,1337.38,3983390000,1337.38 +2011-04-11,1329.01,1333.77,1302.42,1319.68,3940338000,1319.68 +2011-04-04,1333.56,1339.46,1322.94,1328.17,3977634000,1328.17 +2011-03-28,1315.45,1337.85,1305.26,1332.41,3659466000,1332.41 +2011-03-21,1281.65,1319.18,1281.65,1313.80,4018022000,1313.80 +2011-03-14,1301.19,1301.19,1249.05,1279.21,4781044000,1279.21 +2011-03-07,1322.72,1327.68,1291.99,1304.28,4133818000,1304.28 +2011-02-28,1321.61,1332.28,1302.58,1321.15,2404496000,1321.15 +2011-02-22,1338.91,1338.91,1294.26,1319.88,1928012500,1319.88 +2011-02-14,1328.73,1344.07,1324.61,1343.01,2517822000,1343.01 +2011-02-07,1311.85,1330.79,1311.74,1329.15,4021990000,1329.15 +2011-01-31,1276.50,1311.00,1276.50,1310.87,4345372000,1310.87 +2011-01-24,1283.29,1302.67,1275.10,1276.34,4631330000,1276.34 +2011-01-18,1293.22,1296.06,1271.26,1283.35,4974835000,1283.35 +2011-01-10,1270.84,1293.24,1262.18,1293.24,4257314000,1293.24 +2011-01-03,1257.62,1278.17,1257.62,1271.50,4731044000,1271.50 +2010-12-27,1254.66,1262.60,1251.48,1257.64,2091158000,1257.64 +2010-12-20,1245.76,1259.39,1241.51,1256.77,2707105000,1256.77 +2010-12-13,1242.52,1246.73,1232.85,1243.91,4454044000,1243.91 +2010-12-06,1223.87,1240.40,1219.50,1240.40,4835082000,1240.40 +2010-11-29,1189.08,1225.57,1173.64,1224.71,4242568000,1224.71 +2010-11-22,1198.07,1198.94,1176.91,1189.40,3205160000,1189.40 +2010-11-15,1200.44,1207.43,1173.00,1199.73,4177436000,1199.73 +2010-11-08,1223.24,1226.84,1194.08,1199.21,4298262000,1199.21 +2010-11-01,1185.71,1227.08,1177.65,1225.85,4798758000,1225.85 +2010-10-25,1184.74,1196.14,1171.70,1183.26,4116414000,1183.26 +2010-10-18,1176.83,1189.43,1159.71,1183.08,4576282000,1183.08 +2010-10-11,1165.32,1184.38,1155.71,1176.19,4449160000,1176.19 +2010-10-04,1144.96,1167.73,1131.87,1165.15,3905616000,1165.15 +2010-09-27,1148.64,1157.16,1132.09,1146.24,4037410000,1146.24 +2010-09-20,1126.57,1148.90,1122.79,1148.67,3884522000,1148.67 +2010-09-13,1113.38,1131.47,1113.38,1125.59,3972432000,1125.59 +2010-09-07,1102.60,1110.88,1091.15,1109.55,3195237500,1109.55 +2010-08-30,1062.90,1105.10,1040.88,1104.51,3718470000,1104.51 +2010-08-23,1073.36,1081.58,1039.70,1064.59,3951328000,1064.59 +2010-08-16,1077.49,1100.14,1063.91,1071.69,3777406000,1071.69 +2010-08-09,1122.80,1129.24,1076.69,1079.25,4064104000,1079.25 +2010-08-02,1107.53,1128.75,1107.17,1121.64,3963460000,1121.64 +2010-07-26,1102.89,1120.95,1088.01,1101.60,4271320000,1101.60 +2010-07-19,1066.85,1103.73,1056.88,1102.66,4580286000,1102.66 +2010-07-12,1077.23,1099.46,1063.32,1064.88,4487664000,1064.88 +2010-07-06,1028.09,1078.16,1018.35,1077.96,4419372500,1077.96 +2010-06-28,1077.50,1082.60,1010.91,1022.58,5100892000,1022.58 +2010-06-21,1122.79,1131.23,1067.89,1076.76,4699712000,1076.76 +2010-06-14,1095.00,1121.01,1089.03,1117.51,4637208000,1117.51 +2010-06-07,1065.84,1092.25,1042.17,1091.60,5369514000,1091.60 +2010-06-01,1087.30,1105.67,1060.50,1064.88,5368597500,1064.88 +2010-05-24,1084.78,1103.52,1040.78,1089.41,5528868000,1089.41 +2010-05-17,1136.52,1148.66,1055.90,1087.69,6528051900,1087.69 +2010-05-10,1122.27,1173.57,1122.27,1135.68,5791750000,1135.68 +2010-05-03,1188.58,1205.13,1065.79,1110.88,7683886000,1110.88 +2010-04-26,1217.07,1219.80,1181.62,1186.69,6310456000,1186.69 +2010-04-19,1192.06,1217.28,1183.68,1217.28,5800096000,1217.28 +2010-04-12,1194.94,1213.92,1186.77,1192.13,5974902000,1192.13 +2010-04-05,1178.71,1194.66,1175.12,1194.37,4461554000,1194.37 +2010-03-29,1167.71,1181.43,1165.77,1178.10,4237947500,1178.10 +2010-03-22,1157.25,1180.69,1152.88,1166.59,4751278000,1166.59 +2010-03-15,1148.53,1169.84,1141.45,1159.90,4588800000,1159.90 +2010-03-08,1138.40,1153.41,1134.90,1149.99,4805318000,1149.99 +2010-03-01,1105.36,1139.38,1105.36,1138.70,4002330000,1138.70 +2010-02-22,1110.00,1112.29,1086.02,1104.49,4194034000,1104.49 +2010-02-16,1079.13,1112.42,1079.13,1109.17,4040725000,1109.17 +2010-02-08,1065.51,1080.04,1056.51,1075.51,4403416000,1075.51 +2010-02-01,1073.89,1104.73,1044.50,1066.19,5082238000,1066.19 +2010-01-25,1092.40,1103.69,1071.59,1073.87,5079534000,1073.87 +2010-01-19,1136.03,1150.45,1090.18,1091.76,5654582400,1091.76 +2010-01-11,1145.96,1150.41,1131.39,1136.03,4363246000,1136.03 +2010-01-04,1116.56,1145.39,1116.56,1144.98,4223070000,1144.98 +2009-12-28,1127.53,1130.38,1114.81,1115.10,2390427500,1115.10 +2009-12-21,1105.31,1126.48,1105.31,1126.48,3013262500,1126.48 +2009-12-14,1107.84,1116.21,1093.88,1102.47,5672874000,1102.47 +2009-12-07,1105.52,1110.72,1085.89,1106.41,4150876000,1106.41 +2009-11-30,1091.07,1119.13,1086.25,1105.98,4535468000,1105.98 +2009-11-23,1094.86,1112.38,1083.74,1091.49,3232000000,1091.49 +2009-11-16,1094.13,1113.69,1086.81,1091.38,4122504000,1091.38 +2009-11-09,1072.31,1105.37,1072.31,1093.48,4218872000,1093.48 +2009-11-02,1036.18,1071.48,1029.38,1069.30,5290226000,1069.30 +2009-10-26,1080.36,1091.75,1033.38,1036.19,6081714000,1036.19 +2009-10-19,1088.22,1101.36,1074.31,1079.60,5118466000,1079.60 +2009-10-12,1071.63,1096.56,1066.71,1087.68,4740370000,1087.68 +2009-10-05,1026.87,1071.51,1025.92,1071.49,4466710000,1071.49 +2009-09-28,1045.38,1069.62,1019.95,1025.21,5210080000,1025.21 +2009-09-21,1067.14,1080.15,1041.17,1044.38,5081302000,1044.38 +2009-09-14,1040.15,1074.77,1035.00,1068.30,6046967900,1068.30 +2009-09-08,1018.67,1048.18,1018.67,1042.73,5137922500,1042.73 +2009-08-31,1025.21,1028.45,991.97,1016.40,5286260000,1016.40 +2009-08-24,1026.59,1039.47,1016.20,1028.93,5744582000,1028.93 +2009-08-17,998.18,1027.59,978.51,1026.13,4664650000,1026.13 +2009-08-10,1008.89,1013.14,992.40,1004.09,5373764000,1004.09 +2009-08-03,990.22,1018.00,990.22,1010.48,6427945900,1010.48 +2009-07-27,978.63,996.68,968.65,987.48,5294932000,987.48 +2009-07-20,942.07,979.79,940.99,979.26,5003300000,979.26 +2009-07-13,879.57,943.96,875.32,940.38,4785464000,940.38 +2009-07-06,894.27,898.72,869.32,879.13,4673382000,879.13 +2009-06-29,919.86,931.92,896.42,896.42,4172432500,896.42 +2009-06-22,918.13,922.00,888.86,918.90,5119916000,918.90 +2009-06-15,942.45,942.45,903.78,921.23,5114026000,921.23 +2009-06-08,938.12,956.23,926.44,946.21,4866352000,946.21 +2009-06-01,923.26,951.69,923.26,940.09,5662470000,940.09 +2009-05-26,887.00,920.02,881.46,919.14,5788812500,919.14 +2009-05-18,886.07,924.60,879.61,887.00,6339728000,887.00 +2009-05-11,922.99,922.99,878.94,882.88,6337752000,882.88 +2009-05-04,879.21,930.17,879.21,929.23,7952024000,929.23 +2009-04-27,862.82,888.70,847.12,877.52,6043558000,877.52 +2009-04-20,868.27,871.80,826.83,866.23,7083169900,866.23 +2009-04-13,855.33,875.63,835.58,869.60,6839301900,869.60 +2009-04-06,839.75,856.91,814.53,856.56,6226187600,856.56 +2009-03-30,809.07,845.61,779.81,842.50,6286869900,842.50 +2009-03-23,772.31,832.98,772.31,815.94,6952819900,815.94 +2009-03-16,758.84,803.24,749.93,768.54,7963276000,768.54 +2009-03-09,680.76,758.29,672.88,756.55,7459435800,756.55 +2009-03-02,729.57,729.57,666.79,683.38,7592844000,683.38 +2009-02-23,773.25,780.12,734.52,735.09,7550775800,735.09 +2009-02-17,818.61,818.61,754.25,770.05,6401515100,770.05 +2009-02-09,868.24,875.01,808.06,826.84,6008821900,826.84 +2009-02-02,823.09,870.75,812.87,868.60,6217632000,868.60 +2009-01-26,832.50,877.86,821.67,825.88,5602004000,825.88 +2009-01-20,849.64,849.64,804.30,831.95,6129762500,831.95 +2009-01-12,890.40,890.40,817.04,850.12,6058756000,850.12 +2009-01-05,929.17,943.85,888.31,890.35,5043904000,890.35 +2008-12-29,872.37,934.73,857.07,931.80,3793110000,931.80 +2008-12-22,887.20,887.37,857.09,872.80,3087105000,872.80 +2008-12-15,881.07,918.85,857.72,887.88,5855972000,887.88 +2008-12-08,882.71,918.57,851.35,879.73,5932454000,879.73 +2008-12-01,888.61,888.61,815.69,876.07,6093950000,876.07 +2008-11-24,801.20,896.25,801.20,896.24,5841565000,896.24 +2008-11-17,873.23,882.29,741.02,800.03,7349040000,800.03 +2008-11-10,936.75,951.95,818.69,873.29,5812934000,873.29 +2008-11-03,968.67,1007.51,899.73,930.99,5296816000,930.99 +2008-10-27,874.28,984.38,845.27,968.75,6460596000,968.75 +2008-10-20,943.51,985.44,852.85,876.77,6037080000,876.77 +2008-10-13,912.75,1044.31,865.83,940.55,7306794000,940.55 +2008-10-06,1097.56,1097.56,839.80,899.22,8403357900,899.22 +2008-09-29,1209.07,1209.07,1098.14,1099.23,6205326000,1099.23 +2008-09-22,1255.37,1255.37,1179.79,1213.27,5327094000,1213.27 +2008-09-15,1250.92,1265.12,1133.50,1255.08,9328214000,1255.08 +2008-09-08,1249.50,1274.42,1211.54,1251.70,6883584000,1251.70 +2008-09-02,1287.83,1303.04,1217.23,1242.31,5017530000,1242.31 +2008-08-25,1290.47,1300.68,1263.21,1282.83,3530036000,1282.83 +2008-08-18,1298.14,1300.22,1261.16,1292.20,4063548000,1292.20 +2008-08-11,1294.42,1313.15,1274.86,1298.20,4534404000,1298.20 +2008-08-04,1253.27,1297.85,1247.45,1296.32,4188240000,1296.32 +2008-07-28,1257.76,1284.93,1234.37,1260.31,5071890000,1260.31 +2008-07-21,1261.82,1291.17,1248.83,1257.76,5663448000,1257.76 +2008-07-14,1241.61,1262.31,1200.44,1260.68,6511124000,1260.68 +2008-07-07,1262.90,1277.36,1225.35,1239.49,5812632000,1239.49 +2008-06-30,1278.06,1292.17,1252.01,1262.90,4850575000,1262.90 +2008-06-23,1319.77,1335.63,1272.00,1278.38,5031320000,1278.38 +2008-06-16,1358.85,1366.59,1314.46,1317.93,4443808000,1317.93 +2008-06-09,1360.83,1370.63,1331.29,1360.03,4526856000,1360.03 +2008-06-02,1399.62,1404.05,1359.90,1360.68,4314358000,1360.68 +2008-05-27,1375.97,1406.32,1373.07,1400.38,3814042500,1400.38 +2008-05-19,1425.28,1440.24,1373.72,1375.93,3905724000,1375.93 +2008-05-12,1389.40,1425.82,1386.20,1425.35,3809532000,1425.35 +2008-05-05,1415.34,1421.57,1384.11,1388.28,3751244000,1388.28 +2008-04-28,1397.96,1422.72,1383.07,1413.90,4066604000,1413.90 +2008-04-21,1387.72,1399.11,1369.84,1397.84,3939778000,1397.84 +2008-04-14,1332.20,1395.90,1324.35,1390.33,3868576000,1390.33 +2008-04-07,1373.69,1386.74,1331.21,1332.83,3663378000,1332.83 +2008-03-31,1315.92,1380.91,1312.81,1370.40,4175550000,1370.40 +2008-03-24,1330.29,1359.68,1312.95,1315.22,4084940000,1315.22 +2008-03-17,1283.21,1341.51,1256.98,1329.51,5630602500,1329.51 +2008-03-10,1293.16,1333.26,1272.66,1288.14,4802348000,1288.14 +2008-03-03,1330.45,1344.19,1282.43,1293.37,4408266000,1293.37 +2008-02-25,1352.75,1388.34,1325.42,1330.63,4046484000,1330.63 +2008-02-19,1355.86,1367.94,1327.04,1353.11,3688347500,1353.11 +2008-02-11,1331.92,1369.23,1320.32,1349.99,3744452000,1349.99 +2008-02-04,1395.38,1395.38,1316.75,1331.29,4035458000,1331.29 +2008-01-28,1330.70,1396.02,1322.26,1395.42,4539542000,1395.42 +2008-01-22,1312.94,1368.56,1270.05,1330.61,5100980000,1330.61 +2008-01-14,1402.91,1417.89,1312.51,1325.19,5006464000,1325.19 +2008-01-07,1414.07,1430.28,1378.70,1401.02,4788802000,1401.02 +2007-12-31,1475.25,1475.83,1411.19,1411.63,3372257500,1411.63 +2007-12-24,1484.55,1498.85,1471.70,1478.49,2016050000,1478.49 +2007-12-17,1465.05,1485.40,1435.65,1484.46,3745900000,1484.46 +2007-12-10,1505.11,1523.57,1467.78,1467.95,3702056000,1467.95 +2007-12-03,1479.63,1510.63,1460.66,1504.66,3415362000,1504.66 +2007-11-26,1440.74,1488.94,1406.10,1481.14,4096428000,1481.14 +2007-11-19,1456.70,1456.70,1415.64,1440.70,3670937500,1440.70 +2007-11-12,1453.66,1492.14,1438.53,1458.74,4095036000,1458.74 +2007-11-05,1505.61,1520.77,1448.51,1453.70,4415684000,1453.70 +2007-10-29,1536.92,1552.76,1492.53,1509.65,3763506000,1509.65 +2007-10-22,1497.79,1535.53,1489.56,1535.28,3716066000,1535.28 +2007-10-15,1562.25,1564.74,1500.26,1500.63,3475220000,1500.63 +2007-10-08,1556.51,1576.09,1546.72,1561.80,2943480000,1561.80 +2007-10-01,1527.29,1561.91,1527.25,1557.59,3011736000,1557.59 +2007-09-24,1525.75,1533.74,1507.13,1526.75,3070800000,1526.75 +2007-09-17,1484.24,1538.74,1471.82,1525.75,3358248000,1525.75 +2007-09-10,1453.50,1489.58,1439.29,1484.25,2851118000,1484.25 +2007-09-04,1473.96,1496.40,1449.07,1453.55,2852217500,1453.55 +2007-08-27,1479.36,1481.47,1432.01,1473.99,2724582000,1473.99 +2007-08-20,1445.94,1479.40,1430.54,1479.37,3053680000,1479.37 +2007-08-13,1453.42,1466.29,1370.60,1445.94,4376236000,1445.94 +2007-08-06,1433.04,1503.89,1427.39,1453.64,5342306000,1453.64 +2007-07-30,1458.93,1488.30,1432.80,1433.06,4510208000,1433.06 +2007-07-23,1534.06,1547.23,1458.95,1458.95,4151786000,1458.95 +2007-07-16,1552.50,1555.90,1529.20,1534.10,3263540000,1534.10 +2007-07-09,1530.43,1555.10,1506.10,1552.50,3066650000,1552.50 +2007-07-02,1504.66,1532.40,1504.66,1530.44,2318562500,1530.44 +2007-06-25,1502.56,1517.53,1484.18,1503.35,3251210000,1503.35 +2007-06-18,1532.90,1537.32,1500.74,1502.56,3217232000,1502.56 +2007-06-11,1507.64,1538.71,1492.65,1532.91,2975814000,1532.91 +2007-06-04,1536.28,1540.53,1487.41,1507.67,3034900000,1507.67 +2007-05-29,1515.55,1540.56,1510.06,1536.34,2953637500,1536.34 +2007-05-21,1522.75,1532.43,1505.18,1515.73,3018380000,1515.73 +2007-05-14,1505.76,1522.75,1498.34,1522.75,2918038000,1522.75 +2007-05-07,1505.57,1513.80,1491.42,1505.85,2805676000,1505.85 +2007-04-30,1494.07,1510.34,1476.70,1505.62,3090694000,1505.62 +2007-04-23,1484.33,1498.02,1473.74,1494.07,2978394000,1494.07 +2007-04-16,1452.84,1484.74,1452.84,1484.35,3001118000,1484.35 +2007-04-09,1443.77,1453.11,1433.91,1452.85,2654060000,1452.85 +2007-04-02,1420.83,1444.88,1416.37,1443.76,2692797500,1443.76 +2007-03-26,1436.11,1437.65,1408.90,1420.86,2837362000,1420.86 +2007-03-19,1386.95,1438.89,1386.95,1436.11,2901376000,1436.11 +2007-03-12,1402.80,1409.34,1363.98,1386.95,3224692000,1386.95 +2007-03-05,1387.11,1410.15,1373.97,1402.84,3123586000,1402.84 +2007-02-26,1451.04,1456.95,1380.87,1387.17,3599964000,1387.17 +2007-02-20,1455.53,1461.57,1448.36,1451.19,2368890000,1451.19 +2007-02-12,1438.00,1457.97,1431.44,1455.54,2527498000,1455.54 +2007-02-05,1448.33,1452.99,1433.44,1438.06,2686990000,1438.06 +2007-01-29,1422.03,1449.33,1418.46,1448.39,2779552000,1448.39 +2007-01-22,1430.47,1440.69,1416.96,1422.18,2783864000,1422.18 +2007-01-16,1430.73,1435.27,1424.21,1430.50,2722427500,1430.50 +2007-01-08,1409.26,1431.23,1403.97,1430.73,2822146000,1430.73 +2007-01-03,1418.03,1429.42,1405.75,1409.71,3117673300,1409.71 +2006-12-26,1410.75,1427.72,1410.45,1418.30,1541112500,1418.30 +2006-12-18,1427.08,1431.81,1410.28,1410.76,2328566000,1410.76 +2006-12-11,1409.81,1431.63,1404.75,1427.09,2707922000,1427.09 +2006-12-04,1396.67,1418.27,1396.67,1409.84,2686182000,1409.84 +2006-11-27,1400.95,1406.30,1377.83,1396.71,2989828000,1396.71 +2006-11-20,1401.17,1407.89,1397.85,1400.95,2053727500,1400.95 +2006-11-13,1380.58,1403.76,1378.80,1401.20,2761356000,1401.20 +2006-11-06,1364.27,1388.92,1364.27,1380.90,2657402000,1380.90 +2006-10-30,1377.30,1381.95,1360.98,1364.30,2692108000,1364.30 +2006-10-23,1368.58,1389.45,1363.94,1377.34,2712532000,1377.34 +2006-10-16,1365.61,1372.87,1356.87,1368.60,2526124000,1368.60 +2006-10-09,1349.58,1366.63,1343.57,1365.62,2365916000,1365.62 +2006-10-02,1335.82,1353.79,1327.10,1349.59,2639458000,1349.59 +2006-09-25,1314.78,1340.28,1311.58,1335.85,2560806000,1335.85 +2006-09-18,1319.85,1328.53,1310.94,1314.78,2409864000,1314.78 +2006-09-11,1298.86,1324.65,1290.93,1319.66,2688896000,1319.66 +2006-09-05,1310.94,1314.67,1292.13,1298.92,2225772500,1298.92 +2006-08-28,1295.09,1312.03,1293.97,1311.01,1952878000,1311.01 +2006-08-21,1302.30,1302.49,1289.82,1295.09,1831910000,1295.09 +2006-08-14,1266.67,1302.30,1266.67,1302.30,2299788000,1302.30 +2006-08-07,1279.31,1283.74,1261.30,1266.74,2293082000,1266.74 +2006-07-31,1278.53,1292.92,1265.71,1279.36,2571830000,1279.36 +2006-07-24,1240.25,1280.42,1240.25,1278.55,2560298000,1278.55 +2006-07-17,1236.20,1262.56,1224.54,1240.29,2475962000,1240.29 +2006-07-10,1265.46,1274.06,1228.45,1236.20,2285754000,1236.20 +2006-07-03,1270.06,1280.38,1263.13,1265.48,1819212500,1265.48 +2006-06-26,1244.50,1276.30,1237.59,1270.20,2367602000,1270.20 +2006-06-19,1251.54,1257.96,1237.17,1244.50,2255366000,1244.50 +2006-06-12,1252.27,1258.64,1219.29,1251.54,2737928000,1251.54 +2006-06-05,1288.16,1288.16,1235.18,1252.30,2682616000,1252.30 +2006-05-30,1280.04,1290.68,1259.38,1288.22,2381012500,1288.22 +2006-05-22,1267.03,1280.54,1245.34,1280.16,2512808000,1280.16 +2006-05-15,1291.19,1297.88,1256.28,1267.03,2648372000,1267.03 +2006-05-08,1325.76,1326.70,1290.38,1291.24,2335326000,1291.24 +2006-05-01,1310.61,1326.53,1303.46,1325.76,2392390000,1325.76 +2006-04-24,1311.28,1316.04,1295.57,1310.61,2435666000,1310.61 +2006-04-17,1289.12,1318.16,1280.74,1311.28,2348590000,1311.28 +2006-04-10,1295.51,1300.74,1282.96,1289.12,1990310000,1289.12 +2006-04-03,1302.88,1314.07,1294.18,1295.50,2285182000,1295.50 +2006-03-27,1302.95,1310.15,1291.84,1294.87,2170618000,1294.87 +2006-03-20,1307.25,1310.88,1295.81,1302.95,2094204000,1302.95 +2006-03-13,1281.58,1310.45,1281.58,1307.25,2274080000,1307.25 +2006-03-06,1287.23,1288.23,1268.42,1281.42,2250934000,1281.42 +2006-02-27,1289.43,1297.57,1278.66,1287.23,2260408000,1287.23 +2006-02-21,1287.24,1294.17,1281.33,1289.43,2100980000,1289.43 +2006-02-13,1266.99,1289.47,1258.34,1287.24,2197072000,1287.24 +2006-02-06,1264.03,1274.56,1253.61,1266.99,2337512000,1266.99 +2006-01-30,1283.72,1287.94,1261.02,1264.03,2485592000,1264.03 +2006-01-23,1261.49,1286.38,1259.42,1283.72,2592450000,1283.72 +2006-01-17,1287.61,1287.79,1260.92,1261.49,2425750000,1261.49 +2006-01-09,1285.45,1294.90,1282.78,1287.61,2321112000,1287.61 +2006-01-03,1248.29,1286.09,1245.74,1285.45,2487450000,1285.45 +2005-12-27,1268.66,1271.83,1246.59,1248.29,1447217500,1248.29 +2005-12-19,1267.32,1270.51,1257.21,1268.66,1888996000,1268.66 +2005-12-12,1259.37,1275.80,1255.52,1267.32,2235374000,1267.32 +2005-12-05,1265.08,1272.89,1250.91,1259.37,2121000000,1259.37 +2005-11-28,1268.25,1268.44,1249.39,1265.08,2280068000,1265.08 +2005-11-21,1248.27,1270.64,1246.90,1268.25,1779777500,1268.25 +2005-11-14,1234.72,1249.58,1226.41,1248.27,2226412000,1248.27 +2005-11-07,1220.14,1235.70,1215.05,1234.72,2063738000,1234.72 +2005-10-31,1198.41,1224.70,1198.41,1220.14,2488110000,1220.14 +2005-10-24,1179.59,1204.01,1178.89,1198.41,2350556000,1198.41 +2005-10-17,1186.57,1197.30,1170.55,1179.59,2408668000,1179.59 +2005-10-10,1195.90,1196.52,1168.20,1186.57,2305280000,1186.57 +2005-10-03,1228.81,1233.34,1181.92,1195.90,2380760000,1195.90 +2005-09-26,1215.29,1229.57,1211.11,1228.81,2075822000,1228.81 +2005-09-19,1237.91,1237.91,1205.35,1215.29,2268336000,1215.29 +2005-09-12,1241.48,1242.60,1224.85,1237.91,2247794000,1237.91 +2005-09-06,1218.02,1243.13,1218.02,1241.48,1986932500,1241.48 +2005-08-29,1205.10,1227.29,1201.07,1218.02,1950290000,1218.02 +2005-08-22,1219.71,1228.96,1204.23,1205.10,1668590000,1205.10 +2005-08-15,1230.40,1236.24,1215.93,1219.71,1721880000,1219.71 +2005-08-08,1226.42,1242.69,1222.67,1230.39,1904968000,1230.39 +2005-08-01,1234.18,1245.86,1225.62,1226.42,1934294000,1226.42 +2005-07-25,1233.68,1245.15,1228.15,1234.18,1877768000,1234.18 +2005-07-18,1227.92,1236.56,1221.13,1233.68,1916710000,1233.68 +2005-07-11,1211.86,1233.16,1211.86,1227.92,1871184000,1227.92 +2005-07-05,1194.44,1212.73,1183.55,1211.86,1885635000,1211.86 +2005-06-27,1191.57,1204.07,1188.30,1194.44,1796724000,1194.44 +2005-06-20,1216.96,1219.59,1191.45,1191.57,1941440000,1191.57 +2005-06-13,1198.11,1219.55,1194.51,1216.96,1876670000,1216.96 +2005-06-06,1196.02,1208.85,1191.09,1198.11,1720456000,1198.11 +2005-05-31,1198.78,1205.64,1191.03,1196.02,1773022500,1196.02 +2005-05-23,1189.28,1199.56,1185.96,1198.78,1627978000,1198.78 +2005-05-16,1154.05,1191.22,1153.64,1189.28,1883610000,1189.28 +2005-05-09,1171.35,1178.87,1146.18,1154.05,1953106000,1154.05 +2005-05-02,1156.85,1178.62,1154.71,1171.35,2031568000,1171.35 +2005-04-25,1152.12,1164.80,1139.19,1156.85,2090184000,1156.85 +2005-04-18,1142.62,1159.95,1136.15,1152.12,2178972000,1152.12 +2005-04-11,1181.20,1190.17,1141.92,1142.62,2119976000,1142.62 +2005-04-04,1172.79,1191.88,1167.72,1181.20,1861984000,1181.20 +2005-03-28,1171.42,1189.80,1163.69,1172.92,2089900000,1172.92 +2005-03-21,1189.65,1189.65,1168.70,1171.42,1975625000,1171.42 +2005-03-14,1200.08,1210.54,1182.78,1189.65,1706090000,1189.65 +2005-03-07,1222.12,1229.11,1198.15,1200.08,1554146000,1200.08 +2005-02-28,1211.37,1224.76,1198.13,1222.12,1665028000,1222.12 +2005-02-22,1201.59,1212.15,1184.16,1211.37,1572115000,1211.37 +2005-02-14,1205.30,1212.44,1197.35,1201.59,1487736000,1201.59 +2005-02-07,1203.03,1208.38,1191.54,1205.30,1465690000,1205.30 +2005-01-31,1171.36,1203.47,1171.36,1203.03,1625228000,1203.03 +2005-01-24,1167.87,1177.50,1163.75,1171.36,1596660000,1171.36 +2005-01-18,1184.52,1195.98,1167.82,1167.87,1607750000,1167.87 +2005-01-10,1186.19,1194.78,1175.64,1184.52,1477400000,1184.52 +2005-01-03,1211.92,1217.80,1182.16,1186.19,1603540000,1186.19 +2004-12-27,1210.13,1217.33,1204.92,1211.92,889520000,1211.92 +2004-12-20,1194.20,1213.66,1193.36,1210.13,1313350000,1210.13 +2004-12-13,1188.00,1207.97,1188.00,1194.20,1761040000,1194.20 +2004-12-06,1191.17,1192.41,1173.79,1188.00,1496380000,1188.00 +2004-11-29,1182.65,1197.46,1172.37,1191.17,1609280000,1191.17 +2004-11-22,1170.34,1186.62,1167.89,1182.65,1118795000,1182.65 +2004-11-15,1184.17,1188.46,1169.19,1170.34,1497040000,1170.34 +2004-11-08,1166.17,1184.17,1162.32,1184.17,1447680000,1184.17 +2004-11-01,1130.20,1170.87,1127.60,1166.17,1665900000,1166.17 +2004-10-25,1095.74,1131.40,1090.29,1130.20,1587360000,1130.20 +2004-10-18,1108.20,1117.96,1094.25,1095.74,1587820000,1095.74 +2004-10-11,1122.14,1127.01,1102.06,1108.20,1388940000,1108.20 +2004-10-04,1131.50,1142.05,1120.19,1122.14,1421640000,1122.14 +2004-09-27,1110.11,1131.64,1101.29,1131.50,1478640000,1131.50 +2004-09-20,1128.55,1131.54,1108.05,1110.11,1288840000,1110.11 +2004-09-13,1123.92,1130.14,1119.82,1128.55,1259360000,1128.55 +2004-09-07,1113.63,1125.26,1113.62,1123.92,1273300000,1123.92 +2004-08-30,1107.77,1120.80,1094.72,1113.63,1033194000,1113.63 +2004-08-23,1098.35,1109.68,1092.82,1107.77,1035120000,1107.77 +2004-08-16,1064.80,1100.26,1064.80,1098.35,1241160000,1098.35 +2004-08-09,1063.97,1079.04,1060.72,1064.80,1264440000,1064.80 +2004-08-02,1101.72,1108.60,1062.23,1063.97,1380380000,1063.97 +2004-07-26,1086.20,1103.73,1078.78,1101.72,1481360000,1101.72 +2004-07-19,1101.39,1116.27,1083.56,1086.20,1492700000,1086.20 +2004-07-12,1112.81,1119.60,1101.07,1101.39,1327060000,1101.39 +2004-07-06,1125.38,1125.38,1108.72,1112.81,1299825000,1112.81 +2004-06-28,1134.43,1144.20,1123.06,1125.38,1356820000,1125.38 +2004-06-21,1135.02,1146.34,1124.37,1134.43,1431640000,1134.43 +2004-06-14,1136.47,1138.96,1122.16,1135.02,1298200000,1135.02 +2004-06-07,1122.50,1142.18,1122.50,1136.47,1209875000,1136.47 +2004-06-01,1120.68,1129.17,1113.32,1122.50,1209350000,1122.50 +2004-05-24,1093.56,1123.95,1090.74,1120.68,1352540000,1120.68 +2004-05-17,1095.70,1105.93,1079.36,1093.56,1360260000,1093.56 +2004-05-10,1098.70,1102.77,1076.32,1095.70,1579360000,1095.70 +2004-05-03,1107.30,1127.74,1098.63,1098.70,1573120000,1098.70 +2004-04-26,1140.60,1146.56,1107.23,1107.30,1631580000,1107.30 +2004-04-19,1134.56,1142.77,1116.03,1140.60,1532860000,1140.60 +2004-04-12,1139.32,1147.78,1120.75,1134.61,1425960000,1134.61 +2004-04-05,1141.81,1150.57,1134.52,1139.32,1367500000,1139.32 +2004-03-29,1108.06,1144.81,1108.06,1141.81,1497700000,1141.81 +2004-03-22,1109.78,1115.27,1087.16,1108.06,1445820000,1108.06 +2004-03-15,1120.57,1125.76,1102.61,1109.78,1483600000,1109.78 +2004-03-08,1156.86,1159.94,1105.87,1120.57,1536120000,1120.57 +2004-03-01,1144.94,1163.23,1143.78,1156.86,1394320000,1156.86 +2004-02-23,1144.11,1151.68,1134.43,1144.94,1441800000,1144.94 +2004-02-17,1145.81,1158.98,1139.00,1144.11,1455325000,1144.11 +2004-02-09,1142.76,1158.89,1138.70,1145.81,1440040000,1145.81 +2004-02-02,1131.13,1142.79,1124.44,1142.76,1551020000,1142.76 +2004-01-26,1141.55,1155.38,1122.38,1131.13,1710520000,1131.13 +2004-01-20,1139.83,1150.51,1134.62,1141.55,1677675000,1141.55 +2004-01-12,1121.86,1139.83,1115.19,1139.83,1607360000,1139.83 +2004-01-05,1108.48,1131.92,1108.48,1121.86,1673340000,1121.86 +2003-12-29,1095.89,1118.85,1095.89,1108.48,1063025000,1108.48 +2003-12-22,1088.66,1098.47,1086.14,1095.89,817782500,1095.89 +2003-12-15,1074.14,1091.06,1068.00,1088.66,1549520000,1088.66 +2003-12-08,1061.50,1074.76,1053.41,1074.14,1358520000,1074.14 +2003-12-01,1058.20,1074.30,1058.20,1061.50,1385780000,1061.50 +2003-11-24,1035.28,1060.63,1035.28,1058.20,1055355000,1058.20 +2003-11-17,1050.35,1050.35,1031.20,1035.28,1331060000,1035.28 +2003-11-10,1053.21,1063.65,1043.46,1050.35,1298900000,1050.35 +2003-11-03,1050.71,1062.39,1044.88,1053.21,1418400000,1053.21 +2003-10-27,1028.91,1053.09,1028.91,1050.71,1538440000,1050.71 +2003-10-20,1039.32,1048.57,1018.32,1028.91,1468480000,1028.91 +2003-10-13,1038.06,1053.79,1036.57,1039.32,1320640000,1039.32 +2003-10-06,1029.85,1048.28,1026.27,1038.06,1250920000,1038.06 +2003-09-29,996.85,1039.31,990.36,1029.85,1472620000,1029.85 +2003-09-22,1036.30,1036.30,996.08,996.85,1427800000,996.85 +2003-09-15,1018.63,1040.29,1013.59,1036.30,1382022000,1036.30 +2003-09-08,1021.39,1032.41,1007.71,1018.63,1373760000,1018.63 +2003-09-02,1008.01,1029.34,1005.73,1021.39,1516300000,1021.39 +2003-08-25,993.06,1008.85,983.57,1008.01,1062420000,1008.01 +2003-08-18,990.67,1011.01,990.67,993.06,1271000000,993.06 +2003-08-11,977.59,992.50,974.21,990.67,1037294000,990.67 +2003-08-04,980.15,985.75,960.84,977.59,1327460000,977.59 +2003-07-28,998.68,1004.59,978.86,980.15,1445600000,980.15 +2003-07-21,993.32,998.89,975.63,998.68,1402620000,998.68 +2003-07-14,998.14,1015.41,978.60,993.32,1531220000,993.32 +2003-07-07,985.70,1010.43,983.63,998.14,1458240000,998.14 +2003-06-30,976.22,995.00,962.10,985.70,1335650000,985.70 +2003-06-23,995.69,995.69,973.80,976.22,1380160000,976.22 +2003-06-16,988.61,1015.33,988.61,995.69,1508520000,995.69 +2003-06-09,987.76,1002.74,972.59,988.61,1385420000,988.61 +2003-06-02,963.59,1007.69,963.59,987.76,1652340000,987.76 +2003-05-27,933.22,965.38,927.33,963.59,1616400000,963.59 +2003-05-19,944.30,944.30,912.05,933.22,1397660000,933.22 +2003-05-12,933.41,948.65,929.30,944.30,1442580000,944.30 +2003-05-05,930.08,939.61,919.72,933.41,1466700000,933.41 +2003-04-28,898.81,930.56,898.81,930.08,1507782000,930.08 +2003-04-21,893.58,919.74,886.70,898.81,1480200000,898.81 +2003-04-14,868.30,896.77,868.30,893.58,1402350000,893.58 +2003-04-07,878.85,904.89,862.76,868.30,1288000000,868.30 +2003-03-31,863.50,885.89,843.68,878.85,1425520000,878.85 +2003-03-24,895.79,895.79,858.09,863.50,1281200000,863.50 +2003-03-17,833.27,895.90,827.17,895.79,1610346000,895.79 +2003-03-10,828.89,841.39,788.90,833.27,1532180000,833.27 +2003-03-03,841.15,852.34,811.23,828.89,1293180000,828.89 +2003-02-24,848.17,848.17,818.54,841.15,1349680000,841.15 +2003-02-18,834.89,852.87,831.48,848.17,1229675000,848.17 +2003-02-10,829.69,843.02,806.29,834.89,1339920000,834.89 +2003-02-03,855.70,864.64,826.70,829.69,1373720000,829.69 +2003-01-27,861.40,868.72,840.34,855.70,1515846000,855.70 +2003-01-21,901.78,906.00,859.71,861.40,1553837500,861.40 +2003-01-13,927.57,935.05,899.02,901.78,1420120000,901.78 +2003-01-06,908.59,932.89,908.32,927.57,1498880000,927.57 +2002-12-30,875.40,911.25,869.45,908.59,1126575000,908.59 +2002-12-23,895.74,903.89,873.62,875.40,762477500,875.40 +2002-12-16,889.48,911.22,880.32,895.76,1427646000,895.76 +2002-12-09,912.23,912.23,888.48,889.48,1295720000,889.48 +2002-12-02,936.31,954.28,895.96,912.23,1436120000,912.23 +2002-11-25,930.55,941.82,912.10,936.31,1277840000,936.31 +2002-11-18,909.83,937.28,893.09,930.55,1635840000,930.55 +2002-11-11,894.74,910.21,872.05,909.83,1374520000,909.83 +2002-11-04,900.96,925.66,891.62,894.74,1517480000,894.74 +2002-10-28,897.65,907.44,867.91,900.96,1485260000,900.96 +2002-10-21,884.39,902.94,873.06,897.65,1526214000,897.65 +2002-10-14,835.32,886.68,828.37,884.39,1588958000,884.39 +2002-10-07,800.58,843.27,768.63,835.32,1868864000,835.32 +2002-09-30,827.37,851.93,794.10,800.58,1736420000,800.58 +2002-09-23,845.39,856.60,817.38,827.37,1572028000,827.37 +2002-09-16,889.81,902.68,839.09,845.39,1453560000,845.39 +2002-09-09,893.92,924.02,877.05,889.81,1125240000,889.81 +2002-09-03,916.07,916.07,870.50,893.92,1311925000,893.92 +2002-08-26,940.86,955.82,903.33,916.07,1134440000,916.07 +2002-08-19,928.77,965.00,927.21,940.86,1281180000,940.86 +2002-08-12,908.64,935.38,876.20,928.77,1327680000,928.77 +2002-08-05,864.24,913.95,833.44,908.64,1474320000,908.64 +2002-07-29,852.84,911.64,852.84,864.24,1772880000,864.24 +2002-07-22,847.76,854.13,775.68,852.84,2337088000,852.84 +2002-07-15,921.39,926.52,842.07,847.75,2275080000,847.75 +2002-07-08,989.03,993.56,900.94,921.39,1607616000,921.39 +2002-07-01,989.82,994.46,934.87,989.03,1368925000,989.03 +2002-06-24,989.14,1005.88,952.92,989.82,1821238000,989.82 +2002-06-17,1007.27,1040.83,985.65,989.14,1330540000,989.14 +2002-06-10,1027.53,1039.04,981.63,1007.27,1437764000,1007.27 +2002-06-03,1067.14,1070.74,1012.49,1027.53,1406760000,1027.53 +2002-05-28,1083.82,1085.98,1054.26,1067.14,1160550000,1067.14 +2002-05-20,1106.59,1106.59,1075.64,1083.82,1080980000,1083.82 +2002-05-13,1054.99,1106.59,1053.90,1106.59,1290860000,1106.59 +2002-05-06,1073.43,1088.92,1048.96,1054.99,1260840000,1054.99 +2002-04-29,1076.32,1091.42,1063.46,1073.43,1408640000,1073.43 +2002-04-22,1125.17,1125.17,1076.31,1076.32,1367020000,1076.32 +2002-04-15,1111.01,1133.00,1099.41,1125.17,1276580000,1125.17 +2002-04-08,1122.73,1131.76,1102.42,1111.01,1313260000,1111.01 +2002-04-01,1147.39,1147.84,1119.49,1122.73,1168260000,1122.73 +2002-03-25,1148.70,1154.45,1131.61,1147.39,1152300000,1147.39 +2002-03-18,1166.16,1173.94,1139.48,1148.70,1262380000,1148.70 +2002-03-11,1164.31,1173.03,1151.01,1166.16,1314260000,1166.16 +2002-03-04,1131.78,1172.76,1130.93,1164.31,1522860000,1164.31 +2002-02-25,1089.84,1131.79,1089.84,1131.78,1383820000,1131.78 +2002-02-19,1104.18,1104.18,1074.36,1089.84,1355350000,1089.84 +2002-02-11,1096.22,1124.72,1094.68,1104.18,1220240000,1104.18 +2002-02-04,1122.20,1122.20,1077.78,1096.22,1539040000,1096.22 +2002-01-28,1133.28,1138.63,1081.66,1122.20,1588520000,1122.20 +2002-01-22,1127.58,1139.50,1117.43,1133.28,1422175000,1133.28 +2002-01-14,1145.60,1148.81,1124.45,1127.58,1373840000,1127.58 +2002-01-07,1172.51,1176.97,1145.45,1145.60,1306000000,1145.60 +2001-12-31,1161.02,1176.55,1136.23,1172.51,1256625000,1172.51 +2001-12-24,1144.89,1164.64,1144.62,1161.02,756117500,1161.02 +2001-12-17,1123.09,1152.44,1122.66,1144.89,1456760000,1144.89 +2001-12-10,1158.31,1158.31,1114.53,1123.09,1370780000,1123.09 +2001-12-03,1139.45,1173.62,1125.78,1158.31,1404560000,1158.31 +2001-11-26,1150.34,1163.38,1125.51,1139.45,1312160000,1139.45 +2001-11-19,1138.65,1152.45,1129.78,1150.34,1021650000,1150.34 +2001-11-12,1120.31,1148.28,1098.32,1138.65,1319400000,1138.65 +2001-11-05,1087.20,1135.75,1087.20,1120.31,1329260000,1120.31 +2001-10-29,1104.61,1104.61,1053.61,1087.20,1239060000,1087.20 +2001-10-22,1073.48,1110.61,1065.64,1104.61,1273620000,1104.61 +2001-10-15,1091.65,1107.12,1057.24,1073.48,1249040000,1073.48 +2001-10-08,1071.37,1099.16,1052.76,1091.65,1311036000,1091.65 +2001-10-01,1040.94,1084.12,1026.76,1071.38,1405360000,1071.38 +2001-09-24,965.80,1040.94,965.80,1040.94,1595600000,1040.94 +2001-09-10,1085.78,1096.94,944.75,965.80,1950081600,965.80 +2001-09-04,1133.58,1155.40,1082.12,1085.78,1336700000,1085.78 +2001-08-27,1184.93,1186.85,1124.87,1133.58,974100000,1133.58 +2001-08-20,1161.97,1185.15,1153.34,1184.93,1015860000,1184.93 +2001-08-13,1190.16,1198.79,1156.07,1161.97,979500000,1161.97 +2001-08-06,1214.35,1214.35,1169.55,1190.16,1002680000,1190.16 +2001-07-30,1205.82,1226.27,1200.41,1214.35,1107360000,1214.35 +2001-07-23,1210.85,1215.22,1165.54,1205.82,1139100000,1205.82 +2001-07-16,1215.68,1225.04,1196.14,1210.85,1221720000,1210.85 +2001-07-09,1190.59,1218.54,1168.46,1215.68,1241860000,1215.68 +2001-07-02,1224.42,1239.78,1188.74,1190.59,935502500,1190.59 +2001-06-25,1225.35,1237.29,1204.64,1224.38,1314152000,1224.38 +2001-06-18,1214.36,1240.24,1207.71,1225.35,1276524000,1225.35 +2001-06-11,1264.96,1264.96,1203.03,1214.36,1189730000,1214.36 +2001-06-04,1260.67,1286.62,1256.36,1264.96,966200000,1264.96 +2001-05-29,1277.89,1278.42,1245.96,1260.67,1106550000,1260.67 +2001-05-21,1291.96,1315.93,1276.42,1277.89,1099780000,1277.89 +2001-05-14,1245.67,1296.48,1241.02,1291.96,1164340000,1291.96 +2001-05-07,1266.61,1270.00,1240.79,1245.67,1010120000,1245.67 +2001-04-30,1253.05,1272.93,1232.00,1266.61,1202060000,1266.61 +2001-04-23,1242.98,1253.07,1207.38,1253.05,1173840000,1253.05 +2001-04-16,1183.50,1253.71,1167.38,1242.98,1353580000,1242.98 +2001-04-09,1128.43,1183.51,1126.38,1183.50,1201175000,1183.50 +2001-04-02,1160.33,1169.51,1091.99,1128.43,1340278000,1128.43 +2001-03-26,1139.83,1183.35,1136.26,1160.33,1255380000,1160.33 +2001-03-19,1150.53,1180.56,1081.19,1139.83,1359450000,1139.83 +2001-03-12,1233.42,1233.42,1148.64,1150.53,1358072000,1150.53 +2001-03-05,1234.18,1267.42,1228.42,1233.42,1070640000,1233.42 +2001-02-26,1245.86,1272.76,1214.50,1234.18,1211820000,1234.18 +2001-02-20,1301.53,1307.16,1215.44,1245.86,1229475000,1245.86 +2001-02-12,1314.76,1336.62,1293.18,1301.53,1135100000,1301.53 +2001-02-05,1349.47,1363.55,1309.98,1314.76,1082720000,1314.76 +2001-01-29,1354.92,1383.37,1348.72,1349.47,1133080000,1349.47 +2001-01-22,1342.54,1369.75,1333.84,1354.95,1212320000,1354.95 +2001-01-16,1318.32,1354.55,1313.33,1342.54,1351900000,1342.54 +2001-01-08,1298.35,1333.21,1276.29,1318.55,1258100000,1318.55 +2001-01-02,1320.28,1350.24,1274.62,1298.35,1642975000,1298.35 diff --git a/gp/gpcf_covar.m b/gp/gpcf_covar.m index 38c3fa0a..7dc2e8c6 100644 --- a/gp/gpcf_covar.m +++ b/gp/gpcf_covar.m @@ -1,23 +1,29 @@ function gpcf = gpcf_covar(varargin) -% GPCF_COVAR creates a covariance matrix in the parameterization of -% correlation and variances for the multivariates Gaussian process model -% based on coregionalization models (LMC). See Gelfand et al. (2004) in -% nonstationary multivariate process modelling through spatially varying -% coregionalization. +% GPCF_COVAR creates a gpstuff covariance structure for the multivariate +% Gaussian Process model based on the coregionalization models (LMC). +% See Gelfand et al. (2004) in nonstationary multivariate process modelling +% through spatially varying coregionalization. Sociedad de Estadistica e +% Investigacion Operativa, Test. +% +% The coregionalization matrix is a covariance matrix and its +% parametrization is given via variances and the correlation matrix. The +% correlation matrix is parametrizated with the transformation proposed by +% Lewandowski et al. (2009) in generating correlation matrices based on +% vines and extended onion method. Journal of multivariate analysis. % % Description : % -% - GPCF = GPCF_COVAR('PARAM1',VALUE1,'PARAM2,VALUE2,...) +% ─ GPCF = GPCF_COVAR('PARAM1',VALUE1,'PARAM2,VALUE2,...) % creates correlation matrix structure in which % the named parameters have the specified values. Any % unspecified parameters are set to default values. % -% - GPCF = GPCF_COVAR(GPCF,'PARAM1',VALUE1,'PARAM2,VALUE2,...) +% ─ GPCF = GPCF_COVAR(GPCF,'PARAM1',VALUE1,'PARAM2,VALUE2,...) % modify a covariance function structure with the named % parameters altered with the specified values. % -% - Parameters for covariance matrix [default] -% R_prior - prior for correlation matrix [prior_R] +% ─ Parameters (of the gpstuff structure) for the covariance structure [default] +% R_prior - prior for correlation matrix [prior_corrunif()] % V_prior - prior for the variances. Must be a structure % with each component being also a structure (the % priors). Otherwise, just use prior_fixed. @@ -31,11 +37,11 @@ % using that specific type of covariance function % in the gpstructure (gp_set). % -% numberClass - number of classes (categories, species, ...) +% numberClass - number of classes (GPs, categories, species, ...) % being modelled (this is required to initialize % other parameters) % -% degreeFreedom_prior - prior for degree of freedoms 'nu' [prior_corrunif] +% degreeFreedom_prior - prior for degree of freedoms 'nu' [prior_corrunif()] % classVariables - value defining which column of x is used to % identify the class variables. If this is not % given returns an error. They have to be given @@ -50,10 +56,10 @@ % * See also % GP_SET, GPCF_*, PRIOR_* % -% Copyright (c) 2007-2010, 2015 Jarno Vanhatalo -% Copyright (c) 2010 Aki Vehtari -% Copyright (c) 2014 Arno Solin and Jukka Koskenranta -% ̣̣̣------------ 2015 Marcelo Hartmann +% Copyright (c) 2018 Marcelo Hartmann. + +% If you use this file. Learn to recognize its author. Copyright is good but +% misleading. % This software is distributed under the GNU General Public % License (version 3 or later); please refer to the file @@ -70,7 +76,7 @@ ip.addParamValue('R', 0.01, @(x) isvector(x) && ~any(abs(x) > 1)); ip.addParamValue('V', 1, @(x) isvector(x) && ~any(x < 0)); ip.addParamValue('R_prior', prior_corrunif(), @(x) isstruct(x) || isempty(x)); - ip.addParamValue('V_prior', prior_t(), @(x) isstruct(x) || isempty(x)); + ip.addParamValue('V_prior', prior_t('s2', 4), @(x) isstruct(x) || isempty(x)); ip.addParamValue('corrFun', {}, @(x) ~isempty(x) && iscell(x)); ip.addParamValue('classVariables', [], @(x) ~isempty(x) && mod(x, 1) == 0 && x > 0); ip.addParamValue('numberClass', [], @(x) mod(x, 1) == 0 && x > 1); @@ -92,8 +98,10 @@ if init % Set the function handles to the subfunctions - gpcf.fh.RealToRho = @gpcf_covar_RealToRho; - gpcf.fh.RhoToReal = @gpcf_covar_RhoToReal; + gpcf.fh.realtoz = @gpcf_covar_realtoz; + gpcf.fh.ztoreal = @gpcf_covar_ztoreal; + gpcf.fh.rhotoreal = @gpcf_covar_rhotoreal; + gpcf.fh.realtorho = @gpcf_covar_realtorho; gpcf.fh.pak = @gpcf_covar_pak; gpcf.fh.unpak = @gpcf_covar_unpak; gpcf.fh.lp = @gpcf_covar_lp; @@ -115,6 +123,14 @@ gpcf.vectSize = (gpcf.numberClass^2 - gpcf.numberClass)/2; gpcf.aValue = ip.Results.aValue; + % create indexes which will easy our lives from now on + seq = 1:gpcf.vectSize; + i = ceil(0.5 + 0.5 * sqrt(1 + 8 * seq)); % ith column + j = seq - (i - 2).*(i - 1)/2; % jth column + ind1 = (j - 1) * gpcf.numberClass + i; % linear indexes for lower triangular + ind2 = (i - 1) * gpcf.numberClass + j; % linear indexes for upper triangular + gpcf.index = [i; j; ind1; ind2; seq]'; + % Initialize correlations functions if init || ~ismember('corrFun', ip.UsingDefaults) gpcf.corrFun = ip.Results.corrFun; @@ -176,7 +192,10 @@ else var = [repmat('var', gpcf.numberClass , 1), num2str((1:gpcf.numberClass)')]; for ind = 1:size(var, 1) - Vpriors.(var(ind, :)).p = ip.Results.V_prior; + var_aux = var(ind, :); + var_aux = var_aux(~isspace(var_aux)); + Vpriors.(var_aux).p = ip.Results.V_prior; + %Vpriors.(var(ind, :)).p = ip.Results.V_prior; end gpcf.p.V = Vpriors; gpcf.varfields = (fields(gpcf.p.V)'); @@ -185,19 +204,19 @@ end -function y = gpcf_covar_RealToRho(x, a, d) +function y = gpcf_covar_realtoz(x, a, d) % Description : -% Transforms the real line to the interval (-1, 1) using modified -% logistic distribution function +% · Transforms the real line to the interval (-1, 1) using the shifted +% logistic function % -% - rho = 2/(1+exp(-ax)) - 1 -% - x is the vector on the real line -% - d equals to [], 1 or 2, and indicates whether the output corresponds +% · rho = 2 ./ (1 + exp(-a*x)) - 1; +% · x is the vector on the real line +% · d equals to [], 1 or 2, and indicates whether the output corresponds % no derivatives, first or second derivative. if isempty(d) % modified logistic distribution function - y = 2 ./ (1 + exp(-a*x)) - 1; + y = 2 ./ (1 + exp(-a*x)) - 1; elseif d == 1 % repeated terms eminax = exp(-a*x); @@ -215,20 +234,148 @@ end -function y = gpcf_covar_RhoToReal(x, a) +function y = gpcf_covar_ztoreal(x, a) % Description : -% Transforms the interval (-1, 1) to the real line, using inverse -% modified logistic distribution function +% · Transforms the interval (-1, 1) to the real line, using the inverse +% hyperbolic tangent function % -% - y = - 1/a * log(2/(x+1) - 1) -% - x in (-1, 1)^(dim(x)) +% · y = - 1/a * log(2/(x+1) - 1) +% · x in (-1, 1)^(dim(x)) -if any(abs(x) > 1) +if any(abs(x) >= 1) error('domain error'); else y = -1/a * log(2./(x+1) - 1); end +end + +function y = gpcf_covar_rhotoreal(gpcf) +% Description : +% · Transforms the matrix R to the matrix Z (composed by the elements +% gamma on the real line) + +ind1 = gpcf.index(:, 3); +ind2 = gpcf.index(:, 4); + +Z = zeros(gpcf.numberClass, gpcf.numberClass); +R = eye(gpcf.numberClass); +R([ind1, ind2]) = [gpcf.R'; gpcf.R']; + +W = chol(R); +Z(1, 2:gpcf.numberClass) = W(1, 2:gpcf.numberClass); +for i1 = 2:(gpcf.numberClass-1) + for j1 = (i1+1):gpcf.numberClass + Z(i1, j1) = W(i1, j1) .* exp(-0.5*sum(log(1-Z(1:(i1-1), j1).^2))); + % Z(i1, j1) = W(i1, j1) .* 1/sqrt(prod(1-Z(1:(i1-1), j1).^2)); + end +end + +y = gpcf.fh.ztoreal(Z(ind2), gpcf.aValue); +end + + +function y = gpcf_covar_realtorho(gpcf, w, d) +% Description : +% · Transforms the vector gamma on the real line to the matrix R (in a vector); + +% · z = 2 ./ (1 + exp(-a*x)) - 1; +% · x is the vector on the real line +% · d equals to [], 1 or 2, and indicates whether the output corresponds +% no derivatives, first or second derivative. + +i = gpcf.index(:, 1); +j = gpcf.index(:, 2); +ind2 = gpcf.index(:, 4); +seq = gpcf.index(:, end); +i2 = gpcf.vectSize; + +Z = zeros(gpcf.numberClass, gpcf.numberClass); +W = eye(gpcf.numberClass, gpcf.numberClass); + +Z(ind2) = gpcf.fh.realtoz(w(1:i2), gpcf.aValue, []); +W(1, 2:gpcf.numberClass) = Z(1, 2:gpcf.numberClass); + +for i1 = 2:gpcf.numberClass + for j1 = i1:gpcf.numberClass + ztmp = 0.5.*sum(log(1 - Z(1:(i1-1), j1).^2)); + W(i1, j1) = ((i1 == j1) + (i1 ~= j1)*Z(i1, j1)) * exp(ztmp); + % ztmp = sqrt(1-Z(1:(i1-1), j1).^2); + % W(i1, j1) = ((i1 == j1) + ~(i1 == j1)*Z(i1, j1)) * prod(ztmp); + end +end + +if isempty(d) + % correlation matrix + R = W'*W; + y = R(ind2); + +elseif d == 1 + % indexes + ind3 = [j i seq]; + + % derivatives of the cholesk decomposition w.r.t values on the real line + dW = zeros(gpcf.numberClass, gpcf.numberClass, i2); + dL = zeros(gpcf.numberClass, gpcf.numberClass, i2); + z = Z(ind2); + dz = gpcf.fh.realtoz(w, gpcf.aValue, 1); + + fz = sqrt(1 - z.^2); + dfz = - 1./sqrt(1 - z.^2) .* z .* dz; + + for k = 1:i2 + % take the column + j1 = ind3(ind3(:, 3) == k, 2); + + % take the indexes for that column + ind4 = ind3(ind3(:, 2) == j1, end); + + fzaux = fz(ind4); + dfzaux = dfz(ind4); + + for i1 = 1:j1 + if i1 == 1 + dW(i1, j1, k) = (k == ind4(1))*dz(k); + + elseif i1 ~= j1 + kaux = ind4(i1); + + if kaux < k + continue + else + ind5 = ind4(1:(i1-1)); + fzaux2 = fz(ind5); + dfzaux2 = dfz(ind5); + e = (k == ind5); + + if kaux == k + %dW(i1, j1, k) = dz(k).*prod(fzaux2(~e)); + dW(i1, j1, k) = dz(k).*exp(sum(log(fzaux2(~e)))); + end + + if kaux > k + %dW(i1, j1, k) = z(kaux).*dfzaux2(e)*prod(fzaux2(~e)); + dW(i1, j1, k) = z(kaux).*dfzaux2(e)*exp(sum(log(fzaux2(~e)))); + end + end + + elseif i1 == j1 + e = (k == ind4); + %dW(i1, j1, k) = dfzaux(e)*prod(fzaux(~e)); + dW(i1, j1, k) = dfzaux(e)*exp(sum(log(fzaux(~e)))); + + end + end + dL(:, :, k) = dW(:, :, k)'; + + end + y = dL; + +elseif d == 2 + % lower cholesk decompostion of R + y = W'; + +end end @@ -237,12 +384,12 @@ % keeping the order of input below. % % Description : -% W = GPCF_COVAR_PAK(GPCF) takes the non-diagonal elements of +% · W = GPCF_COVAR_PAK(GPCF) takes the non-diagonal elements of % the correlation matrix and put them into a single row vector w % This is a mandatory subfunction used for example % in energy and gradient computations. % -% w = [rho_(1, 2), ..., rho_(1, J), rho_(2, 3), ..., rho_(2, J), ..., rho_((J-1), J)] +% · w = [rho_(1, 2), ..., rho_(1, J), rho_(2, 3), ..., rho_(2, J), ..., rho_((J-1), J)] % % * See also % GPCF_covar_UNPAK @@ -259,12 +406,11 @@ % pak correlations if ~isempty(gpcf.p.R) - seq = 1:gpcf.vectSize; - i = ceil(0.5 + 0.5 * sqrt(1 + 8 * seq)); - j = seq - (i - 2).*(i - 1)/2; - - w = [w gpcf.fh.RhoToReal(gpcf.R, gpcf.aValue)]; - S = [repmat('rhoreal.', gpcf.vectSize, 1), num2str(j'), num2str(i')]; + i = gpcf.index(:, 1); + j = gpcf.index(:, 2); + waux = gpcf.fh.rhotoreal(gpcf); + w = [w waux']; + S = [repmat('realcorr.', gpcf.vectSize, 1), num2str(j), num2str(i)]; s = [s; cellstr(S)]; h = [h 1]; @@ -302,13 +448,13 @@ % structure % % Description : -% [GPCF, W] = GPCF_COVAR_UNPAK(GPCF, W) takes a covariance -% function structure GPCF and a hyperparameter vector W, -% and returns a covariance function structure identical to -% the input, except that the covariance hyperparameters have -% been set to the values in W. Deletes the values set to GPCF -% from W and returns the modified W. This is a mandatory -% subfunction used for example in energy and gradient computations. +% · [GPCF, W] = GPCF_COVAR_UNPAK(GPCF, W) takes a covariance +% function structure GPCF and a hyperparameter vector W, +% and returns a covariance function structure identical to +% the input, except that the covariance hyperparameters have +% been set to the values in W. Deletes the values set to GPCF +% from W and returns the modified W. This is a mandatory +% subfunction used for example in energy and gradient computations. % % * See also % GPCF_covar_PAK @@ -316,16 +462,15 @@ % unpak parameters of correlations functions for i = 1:gpcf.numberClass cf = gpcf.corrFun{i}; - [gpcf.corrFun{i} w] = cf.fh.unpak(cf, w); + [gpcf.corrFun{i}, w] = cf.fh.unpak(cf, w); end gpp = gpcf.p; % unpak for Σ if ~ isempty(gpp.R) - i2 = gpcf.vectSize; - gpcf.R = gpcf.fh.RealToRho(w(1:i2), gpcf.aValue, []); - w = w((i2 + 1):end); + gpcf.R = gpcf.fh.realtorho(gpcf, w(1:gpcf.vectSize), [])'; + w = w((gpcf.vectSize + 1):end); % unpak hyperparameters of R [p, w] = gpcf.p.R.fh.unpak(gpcf.p.R, w); @@ -352,8 +497,8 @@ % GPCF_COVAR_LP Evaluate the log prior of covariance function parameters % % Description : -% LP = GPCF_COVAR_LP(GPCF, X, T) takes a correlation function -% structure GPCF and evaluates log-prior(R) +% · LP = GPCF_COVAR_LP(GPCF, X, T) takes a correlation function +% structure GPCF and evaluates log-prior(R) % % * See also % GPCF_COVAR_PAK, GPCF_covar_UNPAK, GPCF_covar_LPG, GP_E @@ -370,14 +515,18 @@ % log-prior for correlations if ~isempty(gpp.R) % mapping to the real line - x = gpcf.fh.RhoToReal(gpcf.R, gpcf.aValue); - drho = gpcf.fh.RealToRho(x, gpcf.aValue, 1); + x = gpcf.fh.rhotoreal(gpcf); + z = gpcf.fh.realtoz(x, gpcf.aValue, []); + dz = gpcf.fh.realtoz(x, gpcf.aValue, 1); % add log-prior on rho parameterization lp = lp + gpp.R.fh.lp(gpcf.R, gpp.R); - % add log(|det J|) on the real line parameterization - lp = lp + sum(log(drho)); + % add log(|det J|) on the real line parameterization. Equation 11 in + % Lewandowski et al (2009). + pwr = gpcf.numberClass - gpcf.index(:, 2) - 1; + lp = lp + 0.5 .* sum(pwr .* log(1 - z.^2)) + sum(log(dz)); + end % log-prior for variances @@ -397,11 +546,11 @@ % to the parameters. % % Description : -% LPG = GPCF_COVAR_LPG(GPCF) takes a covariance function -% structure GPCF and returns LPG = d log (p(th))/dth, where th -% is the parametric vector of correlations and variances. -% This is a mandatory subfunction used for example in gradient -% computations. +% · LPG = GPCF_COVAR_LPG(GPCF) takes a covariance function +% structure GPCF and returns LPG = d log (p(th))/dth, where th +% is the parametric vector of correlations and variances. +% This is a mandatory subfunction used for example in gradient +% computations. % % * See also % GPCF_covar_PAK, GPCF_covar_UNPAK, GPCF_covar_LP, GP_G @@ -418,16 +567,28 @@ % for correlations if ~isempty(gpcf.p.R) % mapping to the real line - x = gpcf.fh.RhoToReal(gpcf.R, gpcf.aValue); + real = gpcf.fh.rhotoreal(gpcf); + z = gpcf.fh.realtoz(real, gpcf.aValue, [])'; + dz = gpcf.fh.realtoz(real, gpcf.aValue, 1)'; + d2z = gpcf.fh.realtoz(real, gpcf.aValue, 2)'; + pwr = (gpcf.numberClass - gpcf.index(:, 2) - 1)'; + + % cholesk decomposition + L = gpcf.fh.realtorho(gpcf, real, 2); + + % transforming to first derivative of L w.r.t gammas + dL = gpcf.fh.realtorho(gpcf, real, 1); - % mapping to the first derivative of rho - drho = gpcf.fh.RealToRho(x, gpcf.aValue, 1); + b = zeros(gpcf.vectSize, gpcf.vectSize); + for k = 1:size(dL, 3) + dRaux = dL(:, :, k)*L' + L*dL(:, :, k)'; + b(:, k) = dRaux(gpcf.index(:, 4)); + end - % mapping to the second derivative of rho - d2rho = gpcf.fh.RealToRho(x, gpcf.aValue, 2); + % dlogp(R(gamma))/dgamma + dlpdgm = gpp.R.fh.lpg(gpcf.R, gpp.R) * b; - % building the grad vector on the real line parameterization - lpgs = gpp.R.fh.lpg(gpcf.R, gpp.R) .* drho + d2rho./drho; + lpgs = dlpdgm - pwr.*(z.*dz./(1 - z.^2)) + d2z./dz; lpg = [lpg lpgs]; end @@ -448,51 +609,54 @@ % coregionalization and its Cholesky decomposition. % % Description : -% C = GP_COVAR_COV(GP, i1) takes the gpcf structure and -% returns the covariance matrix Σ and its Cholesky decomposition. -% Every element (i, j) of Σ contains the covariance (parameterised by -% correlation and variances) between class i and class j. This is a -% mandatory subfunction. +% · C = GP_COVAR_COV(GP, i1) takes the gpcf structure and +% returns the covariance matrix Σ and its Cholesky decomposition. +% Every element (i, j) of Σ contains the covariance (parameterised by +% correlation and variances) between class i and class j. This is a +% mandatory subfunction. % % * See also % GPCF_COVAR_TRCOV, GPCF_COVAR_TRVAR, GP_COV, GP_TRCOV -% building covariance matrix in correlation parameterization -seq = 1 : gpcf.vectSize; - % linear indexes -i = ceil(0.5 + 0.5 * sqrt(1 + 8 * seq)); -j = seq - (i - 2) .* (i - 1) / 2; -ind1 = (j - 1) * gpcf.numberClass + i; -ind2 = (i - 1) * gpcf.numberClass + j; +ind1 = gpcf.index(:, 3); +ind2 = gpcf.index(:, 4); % corr matrix R = eye(gpcf.numberClass); % filling elements in lower and upper part -R([ind1, ind2]) = [gpcf.R'; gpcf.R']; -S = diag(sqrt(gpcf.V)); +R([ind1, ind2]) = [gpcf.R, gpcf.R]; % matrix Σ -Sig = S * R * S; +Sig = bsxfun(@times, bsxfun(@times, sqrt(gpcf.V)', R), sqrt(gpcf.V)); +% Sig = S * R * S; % other options if nargin == 2 - % empty matrix - C{1} = Sig; - - % calcule cholesk decomposition - [L, notpositivedefinite] = chol(Sig, 'lower'); - - % test if it is whether is positive-definite or not - if notpositivedefinite - L = NaN(gpcf.numberClass, gpcf.numberClass); - L(ind2) = 0; + switch type + case 'cov' + % empty matrix + C{1} = Sig; + + % calcule cholesk decomposition of Sigma + [L, notpositivedefinite] = chol(Sig, 'lower'); + % test whether is positive-definite or not + if notpositivedefinite + L = NaN(gpcf.numberClass, gpcf.numberClass); + L(ind2) = 0; + end + + % take cholesky decomposition + C{2} = L; + + case 'corr' + x = gpcf.fh.rhotoreal(gpcf); + L = gpcf.fh.realtorho(gpcf, x, 2); + + C{1} = L*L'; + C{2} = L; end - - % take cholesky decomposition - C{2} = L; - else C = Sig; @@ -504,12 +668,12 @@ % GP_COVAR_COV Evaluate the covariance matrix between two input vectors % % Description -% C = GP_COVAR_COV(GP, TX, X) takes in correlation structure -% and two matrixes TX and X that contain input vectors to GP. -% Returns covariance matrix C. -% Every element ij of C contains correlation between inputs i -% in TX and j in X. This is a mandatory subfunction used for -% example in prediction and energy computations. +% · C = GP_COVAR_COV(GP, TX, X) takes in correlation structure +% and two matrixes TX and X that contain input vectors to GP. +% Returns covariance matrix C. +% Every element ij of C contains correlation between inputs i +% in TX and j in X. This is a mandatory subfunction used for +% example in prediction and energy computations. % % % * See also @@ -545,14 +709,13 @@ nb1 = find(diff([-inf xClass1' inf])); nb2 = find(diff([-inf xClass2' inf])); -% if ~all([na1 na2] == gpcf.numberClass) -% error('more/less classes than given to the model'); -% end - % chol(Σ) - M = gpcf.fh.sigma(gpcf, 'chol'); + M = gpcf.fh.sigma(gpcf, 'cov'); L = M{2}; + % M = gpcf.fh.sigma(gpcf, 'corr'); + % L = bsxfun(@times, sqrt(gpcf.V)', M{2}); + C = zeros(n1, n2); for k = 1:gpcf.numberClass @@ -579,17 +742,17 @@ % GPCF_COVAR_CFG Evaluate the derivarive of the covar matrix % % Description : -% dRff = GPCF_COVAR_CFG(GPCF, X) takes a -% covariance structure GPCF, a matrix of inputs -% vectors and returns dSRS/drho and/or dSRS/dV, the gradients of -% SRS matrix -% -% dRff = GPCF_COVAR_CFG(GPCF, X, X2) takes a -% covariance structure GPCF, a matrix X of input -% vectors and returns dSRS/drho and/or dSRS/dV, the gradients of -% SRS matrix with respect to rho (cell array with matrix -% elements). This subfunction is needed when using sparse -% approximations (e.g. FIC). +% · dRff = GPCF_COVAR_CFG(GPCF, X) takes a +% covariance structure GPCF, a matrix of inputs +% vectors and returns dSRS/drho and/or dSRS/dV, the gradients of +% SRS matrix +% +% · dRff = GPCF_COVAR_CFG(GPCF, X, X2) takes a +% covariance structure GPCF, a matrix X of input +% vectors and returns dSRS/drho and/or dSRS/dV, the gradients of +% SRS matrix with respect to rho (cell array with matrix +% elements). This subfunction is needed when using sparse +% approximations (e.g. FIC). % % * See also % GPCF_COVAR_PAK, GPCF_COVAR_UNPAK, GPCF_COVAR_LP, GP_G @@ -597,16 +760,14 @@ dCRV = {}; % covariance matrix Σ and chol(Σ) -M = gpcf_covar_sigma(gpcf, 1); -Sig = M{1}; L = M{2}; +M = gpcf.fh.sigma(gpcf, 'corr'); % R = M{1}; +L = M{2}; -% auxiliar indexes -seq = 1 : gpcf.vectSize; -ii = ceil(0.5 + 0.5 * sqrt(1 + 8 * seq)); -jj = seq - (ii - 2) .* (ii - 1) / 2; - -% for correlations and variances; -ind3 = [jj' ii' seq']; +% covariance matrix +%Sig = sVL * sVL'; + +% \sqrt(diag(var)) times L +sVL = bsxfun(@times, sqrt(gpcf.V)', L); % take class variables xC = x(:, gpcf.classVariables); @@ -624,36 +785,19 @@ % number of observations n = length(xC); -if nargin == 5 - % Use memory save option - savememory = 1; -% if i1 == 0 -% i = 0; j = 0; -% % Return number of hyperparameters -% if ~isempty(gpcf.p.R) -% i = length(gpcf.R); -% end -% if ~isempty(gpcf.p.V) -% j = length(gpcf.V); -% end -% dCRV = i + j; -% return -% end -else - savememory = 0; -end - -% Evaluate the values (dSRS/drho_12, dSRS/drho_13, ..., dSRS/drho_1J, dSRS/ -% drho_23, ..., ..., dSRS/drho_2J, ..., dSRS/drho_(J-1)J and/or dSRS/dV_1 -% ... dSRS/dV_J +% use savememory option +savememory = nargin == 5; if nargin == 2 || (isempty(x2) && isempty(mask)) + % Evaluate the values (dSRS/dgamma_12, dSRS/dgamma_13, ..., dSRS/dgamma_(J-1)J + % and/or dSRS/dV_1 ... dSRS/dV_J + % build the big Tj's T = cell(1, gpcf.numberClass); for k = 1:gpcf.numberClass - Taux = L(:, k) * L(:, k)'; + Taux = sVL(:, k) * sVL(:, k)'; for i = 1:gpcf.numberClass inb1 = nb(i) : nb(i + 1) - 1; @@ -668,7 +812,7 @@ end end - % matrices + % correlations matrices dC = {}; for i1 = 1:nc @@ -685,106 +829,52 @@ end end - % correlations and variances in Σ = SRS - dR = {}; - dV = {}; + % correlations and variances in Σ = sqrt(S) L L sqrt(S) + sVdL = []; + dsVL = []; - if ~isempty(gpcf.p.R) || ~isempty(gpcf.p.V); - - dA = zeros(gpcf.numberClass, gpcf.numberClass); - - % for correlations in Σ + if ~isempty(gpcf.p.R) || ~isempty(gpcf.p.V) + + % for gammas in Σ if ~isempty(gpcf.p.R) - % creating cell of matrices - dR = cell(1, gpcf.vectSize); + % transforming to the real line + real = gpcf.fh.rhotoreal(gpcf); - % building the matrices - for k = 1:gpcf.vectSize - % transforming to the real line - real = gpcf.fh.RhoToReal(gpcf.R(k), gpcf.aValue); - - % first derivatives w.r.t correlations x variances - drhoV = gpcf.fh.RealToRho(real, gpcf.aValue, 1) * ... - prod(sqrt(gpcf.V(ind3(k, 1:2)))); - - % indexes - u = ind3(k, 1:2); - - % derivatives - dA(u(1), u(2)) = drhoV; - dA(u(2), u(1)) = drhoV; - - % derivative matrix and put zeros back to dA - dR{k} = dA; - dA(:, :) = 0; - end + % transforming to first derivative of L w.r.t gammas + dL = gpcf.fh.realtorho(gpcf, real, 1); + + % variances times first derivatives w.r.t correlations + % \sqrt(diag(var)) (dLj/dΘ) + sVdL = bsxfun(@times, sqrt(gpcf.V)', dL); end % for variances in Σ - if ~isempty(gpcf.p.V) - % creating cell of matrices - dV = cell(1, nc); - - % building the matrices - for k = 1:gpcf.numberClass - % takes shared correlations - aux = ind3(logical((ind3(:, 1) == k) + (ind3(:, 2) == k)), 1:2); - aux = [k k; aux]; - - % get linear indexes - ind = (aux(:, 1) - 1)*nc + aux(:, 2); - - % take the elements of Sig - u = Sig(ind); - - % derivative w.r.t variances - dA(aux(1, 1), aux(1, 2)) = u(1); - - % loop for building the matrix - for kk = 2:nc - dA(aux(kk, 1), aux(kk, 2)) = 0.5 * u(kk); - dA(aux(kk, 2), aux(kk, 1)) = 0.5 * u(kk); - end - - % derivative matrix and put zeros back to dA - dV{k} = dA; - dA(:, :) = 0; + if ~isempty(gpcf.p.V) + dV = 0.5 .* sqrt(gpcf.V); + dsVLm = bsxfun(@times, dV', L); % take the lines ... + dsVL = zeros(nc, nc, nc); + + for k = 1:nc + dsVL(k, :, k) = dsVLm(k, :); end end - % derivatives w.r.t to correlations or/and variances; - dSig = [dR dV]; + % derivatives w.r.t to correlations or/and variances + dSig = cat(3, sVdL, dsVL); if ~isempty(dSig) - % evaluate derivatives of cholesky decomposition - % see Bayesian filtering and smoothing by Simo Särkkã page 211 - - % indexes for the upper triangular matrix - ind = (ii - 1) * nc + jj; - + nSig = length(dSig); dL = cell(1, nSig); % dCorr = zeros(n, n); for k = 1:nSig - % L^-1 x dΣ/dθ x L^-t - M = (L \ dSig{k}) / L'; - - % apply function ø(·) - M(1 : nc + 1 : nc^2) = M(1 : nc + 1 : nc^2) * 0.5; - M(ind) = 0; - - % compute the derivatives of cholesky decomposition - dA = L * M; - dA(abs(dA) < eps) = 0; - dLs = dA; - % sum_j R_j kron (daj/dΘ * aj' + aj * daj/dΘ') dL{k} = zeros(n, n); for j = 1:nc - aj = repelem(L(:, j).', diffnb); - daj = repelem(dLs(:, j), diffnb); + aj = repelem(sVL(:, j).', diffnb); + daj = repelem(dSig(:, j, k), diffnb); dd = daj * aj; cf = gpcf.corrFun{j}; @@ -806,7 +896,6 @@ error('nargin == 4 || nargin == 5 not implemented yet'); end - end @@ -814,12 +903,12 @@ % GP_COVAR_TRCOV Evaluate training covariance matrix of inputs % % Description : -% C = GP_COVAR_TRCOV(GP, TX) takes in correlation structure -% matrix TX that contains training input vectors. -% Returns covariance matrix C. Every element ij of C contains -% the correlation between inputs i and j in TX. -% This is a mandatory subfunction used for example in -% prediction and energy computations. +% · C = GP_COVAR_TRCOV(GP, TX) takes in correlation structure +% matrix TX that contains training input vectors. +% Returns covariance matrix C. Every element ij of C contains +% the correlation between inputs i and j in TX. +% This is a mandatory subfunction used for example in +% prediction and energy computations. % % * See also % GPCF_COVAR_COV, GPCF_COVAR_TRVAR, GP_COV, GP_TRCOV @@ -831,7 +920,7 @@ % checking if ~issorted(xC) - error('class variable should increasing downwards'); + error('class variable should increase downwards'); end % number of observations @@ -842,6 +931,7 @@ % checking if max(a) > gpcf.numberClass +% if na ~= gpcf.numberClass || max(a) ~= gpcf.numberClass error('more or less classes than given to the model'); end @@ -850,8 +940,6 @@ % number of observations in each class nb = find(diff([-inf xC' inf])); -% nb = diff([0 xC' xC(end) + 1]) .* (1:(n + 1)); -% nb = nb(nb ~= 0); % Try to use the C implementation % C = trcov_corr(gpcf, x1); @@ -860,11 +948,11 @@ if isnan(C) % take chol(Σ) - M = gpcf.fh.sigma(gpcf, 'chol'); + M = gpcf.fh.sigma(gpcf, 'cov'); L = M{2}; % full matrices - C = zeros(n, n); K = zeros(n, n); + C = zeros(n, n); K = zeros(n, n); for k = 1:gpcf.numberClass % T_k = a_k * a_k' matrices, @@ -896,29 +984,64 @@ % GP_COVAR_VECTOR Evaluate training variance vector % % Description: -% C = GP_COVAR_TRVAR(GPCF, TX) takes in correlation structure -% of and matrix TX that contains training inputs. -% Returns correlation vector C, which are ones. -% Every element i of C contains the correlation of the input in TX. -% This is a mandatory subfunction used for example in prediction and -% energy computations. +% · C = GP_COVAR_TRVAR(GPCF, TX) takes in correlation structure +% of and matrix TX that contains training inputs. +% Returns correlation vector C, which are ones. +% Every element i of C contains the correlation of the input in TX. +% This is a mandatory subfunction used for example in prediction and +% energy computations. % % * See also: % GPCF_COVAR_COV, GP_COV, GP_TRCOV - % take indicator class - x = x(:, gpcf.classVariables); - - % information of the classes in the data - a = unique(x); - - % locate each class - nb = find(diff([-inf x' inf])); - - % number of observation in each class - diffnb = diff(nb); - - C = repelem(gpcf.V(a), diffnb)'; +% takes the class variables +xC = x(:, gpcf.classVariables); + +% checking +if ~issorted(xC) + error('class variable should increase downwards'); +end + +% number of observations +n = size(xC, 1); + +% getting the information of the classes in the data +a = unique(xC); % na = size(a, 1); + +% checking +if max(a) > gpcf.numberClass + error('more or less classes than given to the model'); +end + +% number of observations in each class +nb = find(diff([-inf xC' inf])); +diffnb = diff(nb); + +M = gpcf.fh.sigma(gpcf, 'cov'); +L = M{2}; + +% full vector; +C = zeros(n, 1); + +for k = 1:gpcf.numberClass + % T_k = diag(a_k * a_k') matrices, + diagTk = repelem(L(:, k).^2, diffnb); + + cf = gpcf.corrFun{k}; + C = C + cf.fh.trvar(cf, x) .* diagTk; +end + +C(abs(C) < eps) = 0; + +% % take indicator class +% xC = x(:, gpcf.classVariables); +% % information of the classes in the data +% a = unique(xC); +% % locate each class +% nb = find(diff([-inf xC' inf])); +% % number of observation in each class +% diffnb = diff(nb); +% C = repelem(gpcf.V(a), diffnb)'; end @@ -926,12 +1049,12 @@ % RECAPPEND Record append % % Description: -% RECCF = GPCF_COVAR_RECAPPEND(RECCF, RI, GPCF) takes a -% correlation structure record structure RECCF, record index RI -% and correlation structure GPCF with the current MCMC -% samples of the parameters. Returns RECCF which contains -% all the old samples and the current samples from GPCF. -% This subfunction is needed when using MCMC sampling (gp_mc). +% · RECCF = GPCF_COVAR_RECAPPEND(RECCF, RI, GPCF) takes a +% correlation structure record structure RECCF, record index RI +% and correlation structure GPCF with the current MCMC +% samples of the parameters. Returns RECCF which contains +% all the old samples and the current samples from GPCF. +% This subfunction is needed when using MCMC sampling (gp_mc). % % * See also: % GP_MC and GP_MC -> RECAPPEND @@ -950,8 +1073,10 @@ reccf.V = []; % Set the function handles - reccf.fh.RealToRho = @gpcf_covar_RealToRho; - reccf.fh.RhoToReal = @gpcf_covar_RhoToReal; + reccf.fh.realtoz = @gpcf_covar_realtoz; + reccf.fh.ztoreal = @gpcf_covar_ztoreal; + reccf.fh.rhotoreal = @gpcf_covar_rhotoreal; + reccf.fh.realtorho = @gpcf_covar_realtorho; reccf.fh.pak = @gpcf_covar_pak; reccf.fh.unpak = @gpcf_covar_unpak; reccf.fh.lp = @gpcf_covar_lp; @@ -982,6 +1107,7 @@ reccf.vectSize = ri.vectSize; reccf.aValue = ri.aValue; reccf.varfields = ri.varfields; + reccf.index = ri.index; end else diff --git a/gp/lik_binomial.m b/gp/lik_binomial.m index 37fb097c..f36dba02 100644 --- a/gp/lik_binomial.m +++ b/gp/lik_binomial.m @@ -320,22 +320,60 @@ end if nargout > 1 - nt=length(Ef); - Ey=zeros(nt,1); - EVary = zeros(nt,1); - VarEy = zeros(nt,1); - for i1=1:nt - ci = sqrt(Varf(i1)); - F = @(x)zt(i1)./(1+exp(-x)).*norm_pdf(x,Ef(i1),sqrt(Varf(i1))); - Ey(i1) = quadgk(F,Ef(i1)-6*ci,Ef(i1)+6*ci); + p1 = 0.4353; p2 = 0.5647; + c1S = 2.2967^2; c2S = 1.3017^2; - F2 = @(x)zt(i1)./(1+exp(-x)).*(1-1./(1+exp(-x))).*norm_pdf(x,Ef(i1),sqrt(Varf(i1))); - EVary(i1) = quadgk(F2,Ef(i1)-6*ci,Ef(i1)+6*ci); + z1 = Ef ./ sqrt(c1S + Varf); + z2 = Ef ./ sqrt(c2S + Varf); - F3 = @(x)(zt(i1)./(1+exp(-x))).^2.*norm_pdf(x,Ef(i1),sqrt(Varf(i1))); - VarEy(i1) = quadgk(F3,Ef(i1)-6*ci,Ef(i1)+6*ci) - Ey(i1).^2; - end - Vary = EVary+VarEy; + % unconditional expectation + p = (p1 .* normcdf(z1) + p2 .* normcdf(z2)); + Ey = zt .* p; + + % unconditional variance + if any(zt ~= 1) + nt = length(Ef); + zz11 = zeros(nt, 1); + zz12 = zeros(nt, 1); + zz22 = zeros(nt, 1); + + for i = 1:length(Ef) + zz11(i) = mvncdf([Ef(i) Ef(i)], [0 0], ... + [Varf(i) + c1S Varf(i); Varf(i) Varf(i) + c1S]); + + zz12(i) = mvncdf([Ef(i) Ef(i)], [0 0], ... + [Varf(i) + c1S Varf(i); Varf(i) Varf(i) + c2S]); + + zz22(i) = mvncdf([Ef(i) Ef(i)], [0 0], ... + [Varf(i) + c2S Varf(i); Varf(i) Varf(i) + c2S]); + + end + Ey12 = p1^2 * zz11 + 2*p1*p2 * zz12 + p2^2 * zz22; + + else + Ey12 = 0; + + end + + % unconditional variance + Vary = zt .* (p - Ey12) + zt.^2 .* (Ey12 - p.^2); + + % nt=length(Ef); + % Ey=zeros(nt,1); + % EVary = zeros(nt,1); + % VarEy = zeros(nt,1); + % for i1=1:nt + % ci = sqrt(Varf(i1)); + % F = @(x)zt(i1)./(1+exp(-x)).*norm_pdf(x,Ef(i1),sqrt(Varf(i1))); + % Ey(i1) = quadgk(F,Ef(i1)-6*ci,Ef(i1)+6*ci); + % + % F2 = @(x)zt(i1)./(1+exp(-x)).*(1-1./(1+exp(-x))).*norm_pdf(x,Ef(i1),sqrt(Varf(i1))); + % EVary(i1) = quadgk(F2,Ef(i1)-6*ci,Ef(i1)+6*ci); + % + % F3 = @(x)(zt(i1)./(1+exp(-x))).^2.*norm_pdf(x,Ef(i1),sqrt(Varf(i1))); + % VarEy(i1) = quadgk(F3,Ef(i1)-6*ci,Ef(i1)+6*ci) - Ey(i1).^2; + % end + % Vary = EVary+VarEy; end nt=length(yt); diff --git a/gp/lik_epgaussian.m b/gp/lik_epgaussian.m index 71ed8c49..79d20b7e 100644 --- a/gp/lik_epgaussian.m +++ b/gp/lik_epgaussian.m @@ -8,7 +8,7 @@ % See also % GP_SET, LIK_* % - +% % Copyright (c) 2013 Ville Tolvanen % This software is distributed under the GNU General Public diff --git a/gp/lik_liks.m b/gp/lik_liks.m index 2d025070..ee2ca771 100644 --- a/gp/lik_liks.m +++ b/gp/lik_liks.m @@ -44,7 +44,7 @@ % Copyright (c) 2011 Aki Vehtari % Copyright (c) 2012 Ville Tolvanen % Copyright (c) 2015-2017 Jarno Vanhatalo -% ------------- 2015-2017 Marcelo Hartmann +% Copyright (c) 2015-2018 Marcelo Hartmann % Copyright (c) 2017 Ville Tolvanen % This software is distributed under the GNU General Public