## Copyright (C) 2004  Rafael Laboissiere
##
## This program is free software. You can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## This program is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; see the file COPYING.  If not, write to the
## Free Software Foundation, 59 Temple Place - Suite 330, Boston, MA
## 02111-1307, USA.

## -*- texinfo -*-
## @deftypefn {Function File} {[@var{L}, @var{P}, @var{Q}, @var{y}, @var{threshold}, @var{slope}, @var{chisq}] =} probit_likelihood_function (@var{x}, @var{control_var}, @var{batch_size}, @var{success_count})
## @end deftypefn
##
## @paragraphindent 0
## This function is primarily a helper function for the maximum
## likelihood optimization in the @samp{probit_analysis} function, but
## it can also be used stand alone.
##
## The input argument @var{x} is a vector of length @code{n}, whose
## first @code{n - 1} elements are the thresholds and the last one is
## the slope for the fitted normal cumulative functions.  Input
## arguments @var{control_var}, @var{batch_size}, and
## @var{success_count} have the same meaning as in
## @samp{probit_analysis}.
##
## Output arguments:
##
## @table @code
## @item @var{L}
## Negated log likelihood computed using the fitted probabilities and the
## variables @var{batch_size}, @var{success_count}, which determine
## binomial distributions.  (The likelohood is negated such that this
## function can be used with @samp{fmins}.
## @item @var{P}
## Cell vector containing the fitted probability values for each
## experimental condition.  Varies between 0 and 1.
## @item @var{Q}
## Simply @code{@var{Q}@{i@} = 1 - @var{P}@{i@}}.
## @item @var{y}
## Cell vector containing the values of @code{@var{slope}(i) * (@var{control_var}@{i@} - @var{threshold})}.
## @item @var{threshold}
## The first @code{n - 1} elements of @var{x}.
## @item @var{slope}
## The last element of @var{x}.
## @item @var{chisq}
## The test statistic chi square for the fitted curves.
## @end table
##
## This function may be useful to compute values of @var{P} for plotting
## the fitted curves.  In this case, drop @var{batch_size} and
## @var{success_count} from the input argument list, drop @var{chisq}
## from the output argument list, and use a fine-grained @var{control_var}.
##
## CAVEAT: Since this funcion is primarily intended for use through
## fmins from @samp{probit_analysis}, no input argument checks are
## realized, for the sake of efficiency.
##
## @seealso{probit_analysis}

## Author: Rafael Laboissiere <rafael@laboissiere.net>
## Description: Finney's probit analysis -- likelihood function
## $Id$

function [L, P, Q, y, threshold, slope, chisq] \
      = probit_likelihood_function (x, control_var, batch_size, success_count)

  wdbz = warn_divide_by_zero;
  warn_divide_by_zero = 0;

  n = length (control_var);

  threshold = x ([1 : n]);
  slope = x (n + 1);

  a = - threshold * slope;
  b = slope;

  P = cell (1, n);
  Q = cell (1, n);
  y = cell (1, n);
  L = 0;
  chisq = 0;

  for i = 1 : n

    ## Eq. 3.6, page 25
    y {i} = a (i) + b * control_var {i};

    ## Eq. 3.1, page 22
    P {i} = (erf (y {i} / sqrt (2)) + 1) / 2;
    Q {i} = 1 - P {i};

    if nargin > 2
      ## Eq. 4.1, page 53
      L = L - sum (success_count {i} .* log (P {i}) \
             + (batch_size {i} - success_count {i}) .* log (Q {i}));
    endif

    if nargout > 6
      ## Table. 3.4, page 32
      chisq = chisq \
          + sum (((success_count {i} - batch_size {i} .* P {i}) .^ 2) \
                ./ (batch_size {i} .* P {i} .* Q {i}));
    endif

  endfor

  warn_divide_by_zero = wdbz;

endfunction
