home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Programmierung / SOURCE.mdf / programm / msdos / c / yamp / testgraf.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1993-01-11  |  5.5 KB  |  208 lines
  1.  
  2. /*
  3. YAMP - Yet Another Matrix Program
  4. Version: 1.6
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/11/93  
  7.  
  8. Copyright(c) Mark Von Tress 1993
  9. Portions of this code are (c) 1991 by Allen I. Holub and are used by
  10. permission
  11.  
  12. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  13. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  14. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  15. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  16. FROM USE OF THIS PROGRAM.
  17.  
  18. */
  19.  
  20. #include "virt.h"
  21.  
  22. //required global declaration for the
  23. //  matrix stack object
  24.  
  25. // unsigned int _stklen = STACKLENGTH;
  26.  
  27. MStack *Dispatch = new MStack;
  28.  
  29.  
  30. VMatrix &getx(int N)
  31. // create an x matrix
  32. {
  33.     Dispatch->Inclevel();
  34.     VMatrix x, c1, x2;
  35.     
  36.     c1 = Fill(N, 1, 1.0);
  37.     x = Index(1, N) - ((double) N) *0.5;
  38.     x = Ch(c1, x);
  39.     
  40.     // push x onto the stack
  41.     Dispatch->Push(x);
  42.     // decrement the subroutine nesting level
  43.     // and return the stack top
  44.     return Dispatch->DecReturn();
  45. }
  46.  
  47. VMatrix &gety(VMatrix &x, VMatrix &beta)
  48. // create a y matrix
  49. {
  50.     Dispatch->Inclevel();
  51.     VMatrix y;
  52.     
  53.     y = x*beta;
  54.     srand(123);
  55.     for (int i = 1; i <= y.r; i++) {
  56.         // use sum of 3 uniforms for an approximate
  57.         // normal random variable
  58.         int u = random(100) + random(100) + random(100) + 3;
  59.         y.M(i, 1) = y.m(i, 1) + 5.0*sin(3.14*((double) (i % 8)) / 7.0)
  60.         + ((double) (u - 150)) / 300.0;
  61.     }
  62.     Dispatch->Push(y);
  63.     return Dispatch->DecReturn();
  64. }
  65.  
  66. VMatrix ®ression(VMatrix& x, VMatrix& y)
  67. // do a multiple linear regression
  68. {
  69.     Dispatch->Inclevel();
  70.     VMatrix yx, reg, betahat;
  71.     int N = x.r, p = x.c;
  72.     
  73.     // solve for regression parameters using sweep
  74.     yx = Ch(y, x);
  75.     reg = Sweep(2, p + 1, Tran(yx) *yx);
  76.     // calculate mean square error
  77.     reg.M(1, 1) = reg.m(1, 1) / ((double) (N - p));
  78.     reg.DisplayMat();
  79.     
  80.     
  81.     // solve regression using normal equations
  82.     betahat = Inv(Tran(x) *x) *Tran(x) *y;
  83.     
  84.     Dispatch->Push(betahat);
  85.     return Dispatch->DecReturn();
  86. }
  87. void plotResiduals(VMatrix &resids)
  88. {
  89.     Dispatch->Inclevel();
  90.     VMatrix grf = Ch(Index(1, resids.r), resids);
  91.     GMatrix Agraph(grf, - '%');
  92.     *Agraph.PathToDriver = "C:\\tc\\bgi";
  93.     *Agraph.title = "Residuals for data";
  94.     *Agraph.title2 = "with serial correlations with frequency 0.25";
  95.     *Agraph.yname = "Residuals";
  96.     *Agraph.xname = "Index";
  97.     Agraph.Href(0.0);
  98.     Agraph.Show();
  99.     Dispatch->Cleanstack();
  100.     Dispatch->Declevel();
  101. }
  102. VMatrix &GetSerialCovar(VMatrix &R, VMatrix &spectdens)
  103. {
  104.    // Parameters to CORR in Timeslab
  105.    // correlogram = CORR(x=R,n=R.r,M=R.r-1,Q=2*R.r,
  106.    //                    iopt=1,r0=r0, per = spectdens)
  107.    // covar = correlogram*r0
  108.    Dispatch->Inclevel();
  109.    VMatrix centered, z, covar;
  110.    double n = (double) R.r;
  111.    // center a column vector
  112.    centered = R - Sum(R).m(1, 1) / n;
  113.    // zero pad to length 2n: 2n periodic for full
  114.    // sample spectral density
  115.    centered = Cv(centered, Fill(R.r, R.c, 0));
  116.    // take fft
  117.    z = Fft(centered);
  118.    // take convolution : gives sample spectral density
  119.    spectdens = Sum(z % z / n, COLUMNS);
  120.    // inverse fft for serial correlation (autocorrelation)
  121.    covar = Fft(spectdens, FFALSE);
  122.    // throw away last half.
  123.    covar = Submat(covar, R.r, 2);
  124.    Dispatch->Push(covar);
  125.    return Dispatch->DecReturn();
  126. }
  127.  
  128. void plotCorrelogram(VMatrix &serial)
  129. {
  130.    Dispatch->Inclevel();
  131.    
  132.    int n = serial.r;
  133.    double sigma = serial.m(1, 1);
  134.    VMatrix Correlogram = serial / sigma;
  135.    Correlogram = Submat(Correlogram, n, 1, 2, 1);
  136.    VMatrix graf = Ch(Index(1, Correlogram.r), Correlogram);
  137.    
  138.    GMatrix Agraph(graf);
  139.    *Agraph.PathToDriver = "C:\\tc\\bgi";
  140.    *Agraph.title = "Serial Correlations";
  141.    *Agraph.title2 = "for sample residuals";
  142.    *Agraph.yname = "Serial correlations";
  143.    *Agraph.xname = "Lags";
  144.    Agraph.Href(0.0);
  145.    Agraph.Show();
  146.    
  147.    Dispatch->Cleanstack();
  148.    Dispatch->Declevel();
  149. }
  150.  
  151. void plotPeriodogram(VMatrix &spectdens)
  152. {
  153.    // calculate a standardized periodogram on the log scale
  154.    Dispatch->Inclevel();
  155.    int n = spectdens.r;
  156.    // this works because data is already centered, which
  157.    // forces spectdens.m(1,1) = 0.0;
  158.    double sigma = Sum(spectdens).m(1, 1) / ((double) n);
  159.    VMatrix Periodogram = spectdens / sigma;
  160.    Periodogram = Mlog(Submat(Periodogram, n / 2 + 1, 1, 2, 1));
  161.    
  162.    // frequencies
  163.    double dn =(double) (n / 2 + 1);
  164.    VMatrix graf = Ch(Index(1, Periodogram.r) / dn, Periodogram);
  165.    
  166.    GMatrix Agraph(graf);
  167.    *Agraph.PathToDriver = "C:\\tc\\bgi";
  168.    *Agraph.title = "Periodogram";
  169.    *Agraph.title2 = "for sample residuals";
  170.    *Agraph.yname = "Log periodogram";
  171.    *Agraph.xname = "Frequencies";
  172.    Agraph.Vref(0.25);
  173.    Agraph.Show();
  174.    
  175.    Dispatch->Cleanstack();
  176.    Dispatch->Declevel();
  177. }
  178.  
  179. main()
  180. {
  181.     Dispatch->Inclevel();
  182.     VMatrix x, beta("beta", 2, 1), y, betahat, resids, serial;
  183.     Setwid(15);
  184.     Setdec(10);
  185.     
  186.     beta.M(1, 1) = 1;
  187.     beta.M(2, 1) = 0.5;
  188.     
  189.     x = getx(128);
  190.     y = gety(x, beta);
  191.     
  192.     betahat = regression(x, y);
  193.     betahat.Nameit("Text book betahat");
  194.     (Tran(beta)).DisplayMat();
  195.     (Tran(betahat)).DisplayMat();
  196.     
  197.     resids = y - x*betahat;
  198.     plotResiduals(resids);
  199.     
  200.     VMatrix spectdens;
  201.     serial = GetSerialCovar(resids, spectdens);
  202.     
  203.     plotCorrelogram(serial);
  204.     plotPeriodogram(spectdens);
  205.     
  206.     vclose();
  207. }
  208.