home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 3.ddi / DEMOS.DI$ / CENSUSEX.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  2.8 KB  |  118 lines
  1. %CENSUSEX Callback script involked by the buttons in CENSUS demo.
  2.  
  3. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  4.  
  5. gold = [.9 .6 .3];
  6. tinv = [12.706 4.303 3.182 2.776 2.571 2.447 2.365 2.306 2.262];
  7.  
  8. if offscale, close(offscale), offscale = 0; end
  9.  
  10. if button == 13
  11.    % Toggle confidence intervals.
  12.    bounds = 1 - bounds;
  13. else
  14.    % Reset the color of the previous button.
  15.    set(buttons(d+1),'back','default')
  16.    % Get the new degree from the button.
  17.    d = button-1;
  18.    % Set the new degree to gold.
  19.    set(buttons(d+1),'back',gold);
  20. end
  21.  
  22. % Set the cursor to a watch to show computation in progress.
  23. point = get(fig,'pointer');
  24. set(fig,'pointer','watch');
  25. drawnow
  26.  
  27. if d <= 9 
  28.    % Polynomial fit or interpolation.  Solve the linear equation
  29.    % involving the last d+1 columns of the Vandermonde matrix.
  30.    [Q,R] = qr(A(:,n-d:n));
  31.    R = R(1:d+1,:);
  32.    Q = Q(:,1:d+1);
  33.    c = R\(Q'*p);    % Same as c = A(:,n-d:n)\p;
  34.    y = polyval(c,x);
  35.    z = polyval(c,10);
  36.    if d < 9
  37.       title(['Polynomial fit, degree = ' int2str(d)]);
  38.    else
  39.       title('Polynomial interpolation, degree = 9')
  40.    end
  41.    if bounds
  42.       RI = inv(R);
  43.       E = zeros(length(x),d+1);
  44.       for j = 1:d+1
  45.          E(:,j) = polyval(RI(:,j),x);
  46.       end
  47.       e = sqrt(1+sum((E.*E)')');
  48.       r = p - A(:,n-d:n)*c;
  49.       if d < 9
  50.          sig = norm(r)/sqrt(9-d) * tinv(9-d);
  51.       else
  52.          sig = 9;   % Just a guess.
  53.       end
  54.       if all(y < sig*e)
  55.          b = boff;
  56.          offscale = text(1950,50,'Bounds off scale.');
  57.       else
  58.          b = [y-sig*e; NaN; y+sig*e];
  59.       end
  60.    end
  61.  
  62. elseif d == 10
  63.    % Cubic spline interpolation
  64.    y = spline(s,p,x);
  65.    z = spline(s,p,10);
  66.    title('Cubic spline interpolation')
  67.    if bounds
  68.       I = eye(n,n);
  69.       E = zeros(length(x),n);
  70.       for j = 1:n
  71.          E(:,j) = spline(s,I(:,j),x);
  72.       end
  73.       e = sqrt(1+sum((E.*E)')');
  74.       sig = 9;
  75.       b = [y-sig*e; NaN; y+sig*e];
  76.    end
  77.  
  78. elseif d == 11
  79.    % Exponential least squares fit
  80.    E = [ones(size(t)) log(t)];
  81.    [Q,R] = qr(E);
  82.    c = R\(Q'*log(p));
  83.    r = log(p) - E*c;
  84.    E = [ones(size(v)) log(v)];
  85.    y = exp(E*c);
  86.    z = exp([1 log(2000)]*c);
  87.    title('Exponential fit')
  88.    if bounds
  89.       E = E/R(1:2,1:2);
  90.       e = sqrt(1+sum((E.*E)')');
  91.       sig = norm(r)/sqrt(n-2) * tinv(n-2);
  92.       b = exp(sig*e);
  93.       b = [y./b; NaN; y.*b];
  94.    end
  95.  
  96. else
  97.    % Done
  98.    clf
  99.    set(fig,'pointer',point)
  100.    disp('Done')
  101.    return
  102. end
  103.  
  104. % Update the plot
  105.  
  106. set(yhandle,'ydata',y);
  107. set(zhandle,'ydata',z);
  108. zz = max(10,min(390,z-10));  % Avoid the edge of the graph.
  109. set(ztext,'pos',[2000 zz],'string',num2str(z))
  110. if bounds
  111.    set(bhandle,'ydata',b);
  112.    set(buttons(13),'back',gold)
  113. else
  114.    set(bhandle,'ydata',boff);
  115.    set(buttons(13),'back','default')
  116. end
  117. set(fig,'pointer',point)
  118.