function glmshow(option)
% glmshow:  A demonstration for generalized linear models

% Copyright 2001
% Peter Dunn

xvalmin = 5;
axismin = 0;
plotmin = 4.5;
xvalmax = 20;
plotmax = 20;
scalefactor = 5;
fontscale = 0.6;


if nargin == 0,

   option = 0;


   % SET UP CANVAS
   figure('Name','Generalized Linear Models Demonstration')

   % Frames
   % * First, the frame
   uicontrol( 'Style','Frame', ...
      'Units','normalized', ...
      'Position',[0.01, 0.02, 0.84, 0.15])


   % Slider
   uicontrol('Style','slider', ...
      'Units','normalized', ...
      'Position',[0.03 0.10, 0.8, 0.05], ...
      'Max',xvalmax,  'Min',xvalmin,  ...
      'Value',5, ...
      'tag','glmshowXValueSlider', ...
      'Callback','glmshow(1)');


   % Distribution
   uicontrol('Style','Text',...
      'String','Distribution:',...
      'HorizontalAlignment','Right', ...
      'FontUnits','normalized',...
      'FontSize',fontscale,...
      'Units','normalized',...
      'Position',[0.02 0.03, 0.13, 0.05]);

   uicontrol('Style','PopupMenu', ...
      'String','normal|gamma|poisson', ...
      'Units','normalized', ...
      'Position',[0.16 0.03, 0.15, 0.05],...
      'Value',1, ...
      'FontUnits','normalized',...
      'FontSize',fontscale,...
      'Callback','glmshow(20)', ...
      'tag','glmshowDistInfo' )


   % Link function
   uicontrol('Style','Text', ...
      'String','Link:', ...
      'FontUnits','normalized',...
      'FontSize',fontscale,...
      'Units','normalized', ...
      'HorizontalAlignment','Right', ...
      'Position',[0.32 0.03, 0.07, 0.05])

   uicontrol('Style','PopUpMenu', ...
      'String','Id|Recip|Log', ...
      'Value',1,...
      'FontUnits','normalized',...
      'FontSize',fontscale,...
      'Units','normalized', ...
      'HorizontalAlignment','Right', ...
      'Position',[0.39 0.03, 0.13, 0.05], ...
      'Callback','glmshow(5)', ...
      'tag','glmshowLink');

   % Dispersion Parameter
   uicontrol('Style','Text', ...
      'String','Dispersion Parameter:', ...
      'FontUnits','normalized',...
      'FontSize',fontscale,...
      'Units','normalized', ...
      'HorizontalAlignment','Right', ...
      'Position',[0.53 0.03, 0.22, 0.05])
   uicontrol('Style','Edit', ...
      'String','1', ...
      'Units','normalized', ...
      'FontUnits','normalized',...
      'FontSize',fontscale,...
      'Position',[0.75 0.03, 0.08, 0.05], ...
      'HorizontalAlignment','Left',... 
      'Callback','glmshow(4)',...
      'tag','glmshowDisPar');
      

   % Close Button
   uicontrol('Style','PushButton', ...
      'Units','normalized', ...
      'Position',[0.87 0.10, 0.1, 0.05], ...
      'String','Close', ...
      'FontUnits','normalized',...
      'FontSize',fontscale,...
      'Callback','close(gcf);');
	   
 

   % Plot the yval-xval curve
   axes('Position',[0.09 0.28 0.89 0.70],'tag','glmshowCurve');
   curveplot = linspace( plotmin, plotmax, 100 );
   
   plot( curveplot, xval2yval( curveplot ), ...
         'LineWidth',2, ...
         'tag','glmshowCurve');
   hold on;

   % Plot lines
   xvalue = get( findobj('tag','glmshowXValueSlider'),'Value');

   % Plot to curve
   yvalue = xval2yval(xvalue);
   plot( xvalue, yvalue, 'ro', ...
         'MarkerFaceColor','r', ...
         'tag','glmshowMean');

   % Curve to distn
   plot( [xvalue 0], [yvalue, yvalue], ...
         'LineStyle',  '--', ...
         'tag','glmshowCurve2Distn' );


   % Plot the distribution

   % Get the information
   xvalue = get( findobj('tag','glmshowXValueSlider'),'Value');
   yvalue = xval2yval( xvalue );
   dispar = get( findobj('tag','glmshowDisPar'), 'String' );
   dispar = str2num( dispar );

   % Get yval distn values
   plotyval = linspace( 0.1, 15, 500 );
   yval = plotdist( plotyval, yvalue, dispar );

   xvalscale = scalefactor / max( yval );

   plot( xvalue - xvalscale*yval, plotyval, 'r-',...
         'lineWidth',2, ...
         'tag', 'glmshowDistn' );
   plot( [xvalue xvalue], [0 15], 'r-', ...
         'tag','glmshowBaseLine' );

   xlabel('x-Value','FontUnits','normalized','FontSize',0.05);
   ylabel('y-Value','FontUnits','normalized','FontSize',0.05);
   set(gca,'FontUnits','normalized','FontSize',0.05);
 

else

   if option==1,
 
      % SLIDER CHANGES

      % Get new value of xval
      xvalue = get( findobj('tag','glmshowXValueSlider'),'Value');

      % Get current dispersion parameter value
      disparvalue = get( findobj('tag','glmshowDisPar'),'String');
      disparvalue = str2num( disparvalue );

      % Determine corresponding value of yval
      yvalue = xval2yval( xvalue );

      % Update the information frame
      set( findobj('tag','glmshowXYValues'), ...
         'String', ['(',num2str(rounddp(xvalue,2)),',',...
                        num2str(rounddp(yvalue,2)),')'] );

      % Update the plots on the graph
      glmshow( 2 );

      % Update distribution
      glmshow( 3 );
      drawnow;

   elseif option==2,

      % PLOT THE DENSITY-YIELD GRAPH LINES

      % Get the information
      xvalue = get( findobj('tag','glmshowXValueSlider'),'Value');


      % Curve to distn
      delete( findobj('tag','glmshowCurve2Distn') )
      yvalue = xval2yval( xvalue );
      plot( [xvalue 0], [yvalue, yvalue], ...
            'LineStyle', '--', ...
            'tag','glmshowCurve2Distn' );       

   elseif option==3,

      % PLOT THE DISTRIBUTION

      % Get the information
      xvalue = get( findobj('tag','glmshowXValueSlider'),'Value');
      yvalue = xval2yval( xvalue );
      dispar = get( findobj('tag','glmshowDisPar'), 'String' );
      dispar = str2num( dispar );

      % Determine distribution
      dist = get( findobj('tag','glmshowDistInfo'), 'Value');

      % Delete old distribution
      delete( findobj('tag','glmshowDistn') )

      if ( dist==1 ) | (dist==2 ) , % Continuous

         % Get yval distn values
         plotyval = linspace( 0.1, 15, 500 );
         yval = plotdist( plotyval, yvalue, dispar );

         xvalscale = scalefactor / max( yval );

         % Plot the distribution
         xvalscale = scalefactor / max( yval );
         plot( xvalue - xvalscale*yval, plotyval, 'r-', ...
               'LineWidth',2, ...
               'tag','glmshowDistn' );

      elseif dist==3, % Poisson
 
         % Get yval distn values
         plotyval = linspace( 0, 14, 15 );
         yval = plotdist( plotyval, yvalue, dispar );

         xvalscale = scalefactor / max( yval(:,2) );

         plot( xvalue - xvalscale*yval(:,2), yval(:,1), 'r-',...
               'LineWidth',2, ...
               'tag','glmshowDistn');
 
      end;

      % Plot the baseline
      delete( findobj('tag','glmshowBaseLine') );
      plot( [xvalue xvalue], [0 15], 'r-', ...
            'tag','glmshowBaseLine' );

         % Plot to curve
         delete( findobj('tag','glmshowMean') )
         plot( xvalue, xval2yval(xvalue), 'ro', ...
            'MarkerFaceColor','r', ...
            'tag', 'glmshowMean');



   elseif option==4,
 
      % THE DISPERSION PARAMETER IS CHANGED

      % Get new value
      disparvalue = get( findobj('tag','glmshowDisPar'),'String');
      disparvalue = str2num( disparvalue );

      % Update the plots on the graph
      glmshow( 2 );

      % Update distribution
      glmshow( 3 );
      drawnow;

   elseif option==5,

      % LINK FUNCTION CHANGES

      % Update the plots on the graph
      glmshow( 2 );

      % Update distribution
      glmshow( 3 );

      % Update curve
      glmshow(10);
      drawnow;

   elseif option==10,

      % PLOT CURVE
      delete( findobj('tag','glmshowCurve') );
      curveplot = linspace( plotmin, plotmax, 100 );
 
      % Plot lines
      xvalue = get( findobj('tag','glmshowXValueSlider'),'Value');

      plot( curveplot, xval2yval( curveplot ), ...
            'LineWidth',2, ...
            'tag','glmshowCurve');

      % Plot to curve
      yvalue = xval2yval(xvalue);
      plot( xvalue, yvalue, 'ro', ...
            'MarkerFaceColor','r', ...
            'tag','glmshowMean');

     % Curve to distn
       plot( [xvalue 0], [yvalue, yvalue], ...
            'LineStyle',  '--', ...
            'tag','glmshowCurve2Distn' );

   elseif option == 20,

      % DISTRIBUTION CHANGES

      % Get new distribution
      dist = get( findobj('tag','glmshowDistInfo'), 'Value');

      if dist==1, % normal

         disparvalue = 1;
         % Also make changeable
         set( findobj('tag','glmshowDisPar'), ...
            'Enable','on');

      elseif dist==2, % gamma

         disparvalue = 0.005;
         % Also make changeable
         set( findobj('tag','glmshowDisPar'), ...
            'Enable','on');

      else, % Poisson

         disparvalue = 1;

         % Also make unchangeable
         set( findobj('tag','glmshowDisPar'), ...
            'Enable','off');

      end;

      % Set new disparvalue
      set( findobj('tag','glmshowDisPar'), ...
         'String',num2str(disparvalue) );

      % Plot curve
      glmshow(10);

      % Plot distribution
      glmshow(3);
   end;

end;

%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% SUBFUNCTIONS %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%

%%%   DEN2YIELD   %%%
function y = xval2yval(xval)

% Get link function
link = get( findobj('tag','glmshowLink'), 'Value');

if ( link==1 ), % id

%   y = 150 - 3.25 * xval;
   y = 13 - 0.6 * xval;

elseif ( link==2), % Recip

%   eta = -0.0137 + 0.0047 * xval;
   eta = -0.2 + 0.06*xval;
   y = 1 ./ eta;

elseif ( link==3 ), % log

%    y = exp( 1.1769 + 0.08*xval);
   eta = 3.07 - 0.1535*xval;
   y = exp( eta );

end;



%%% PLOTDIST %%%
function xval = plotdist( y, mu, disparvalue )

% Get distribution specified
dist = get( findobj('tag','glmshowDistInfo'), 'Value');

if dist==1, % normal

   xval = normaldist( y, mu, disparvalue );

elseif  dist==2, % gamma

   xval = gammadist( y, mu, disparvalue );

elseif dist==3, % Poison

   xval = zeros( 5*length(y) + 1, 2 );
   tmpxval = poissondist( y, mu, disparvalue );
   [dump1, dump2, dump3, x, y] = makebars( y,  tmpxval, 'histc' );
   xval(:,1) = x;
   xval(:,2) = y;

end;



%%%   GAMMADIST   %%%

function xval = gammadist( y, mu, disparvalue )

alpha = 1 / disparvalue;

logxval = ( (alpha-1) * log(y) ) + ( alpha * log(alpha) ) - ...
             ( y * alpha / mu ) - ( alpha * log(mu) ) - ...
              gammaln( alpha );
xval = exp( logxval );



%%% NORMALDIST %%%
function xval = normaldist( y, mu, disparvalue)

z = ( y - mu )./disparvalue;
front = sqrt( 2*pi.*disparvalue);
xval = (1./front) .* exp( -0.5.*z.^2) ;



%% POISSONDIST %%%
function xval = poissondist( y, mu, disparvalue)

logxval = -mu + y.*log(mu) - gammaln( y+1);
xval = exp( logxval);



%%% ROUNDDP %%%
function out = rounddp( in, places )
% Rounds to  places  decimal places

out = round( in*10^places) / ( 10^places );