Logo Search packages:      
Sourcecode: r-noncran-lindsey version File versions  Download package

lcr.c

/*
 *  Ordinal: A Library of Ordinal Models
 *  Copyright (C) 1998, 1999, 2000, 2001  J.K. Lindsey and P.J. Lindsey
 *
 *  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 of the License, 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; if not, write to the Free Software
 *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 *  SYNOPSIS
 *
 * void lbr(int y[],int freq[],double cov[],int *ncv,double beta[],
 *          double coef[],int *n,double totncv[],int *iter,int *iterlim,
 *          double *steptol,double sd[],double fd[],int *transform,
 *          int *link,double lambda[],int *rank,double aux[],int pivot[],
 *          double work[],double diff[],int *info,double diag[],
 *          double hessian[],double *loglik,double fit[],int pred[],
 *          double cpred[])
 * void mllbr(int y[],int freq[],double cov[],int *ncv,int *n,double beta[],
 *            int *transform,int *link,double lambda[],double *loglik)
 * void fvbr(int y[],int freq[],double cov[],int *ncv,int *n,double beta[],
 *           int *transform,int *link,double lambda[],double fit[],int pred[],
 *           double cpred[])
 *
 * void lmr(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],
 *          double coef[],int *n,int totnc[],double totncv[],int *iter,
 *          int *iterlim,double *steptol,double sd[],double fd[],double reg[],
 *          int *transform,double p[],int *rank,double aux[],int pivot[],
 *          double work[],double diff[],int *info,double diag[],
 *          double hessian[],double *loglik,double fit[],int pred[],
 *          double cpred[])
 * void mllmr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,
 *            double beta[],int *transform,double *loglik,double top[])
 * void fvmr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,
 *           double beta[],int *transform,double p[],double fit[],int pred[],
 *           double cpred[])
 *
 * void lsmr(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],
 *           double coef[],int *n,int totnc[],double totncv[],int *iter,
 *           int *iterlim,double *steptol,double sd[],double fd[],double reg[],
 *           int *transform,double p[],int *rank,double aux[],int pivot[],
 *           double work[],double diff[],int *info,double diag[],
 *           double hessian[],double *loglik,double fit[],int pred[],
 *           double cpred[])
 * void mllsmr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,
 *             double beta[],int *transform,double *loglik,double top[])
 * void fvsmr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,
 *            double beta[],int *transform,double p[],double fit[],int pred[],
 *            double cpred[])
 *
 * void lucrr(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],
 *            double coef[],int *n,int totnc[],double totncv[],int *iter,
 *            int *iterlim,double *steptol,double sd[],double fd[],double reg[],
 *            int *transform,double p[],double q[],int *rank,double aux[],
 *            int pivot[],double work[],double diff[],int *info,double diag[],
 *            double hessian[],double *loglik,double fit[],int pred[],
 *            double cpred[])
 * void mllucrr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,
 *              double beta[],int *transform,double *loglik,double top[],
 *              double bot[])
 * void fvucrr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,
 *             double beta[],int *transform,double p[],double q[],double fit[],
 *             int pred[],double cpred[])
 *
 * void ldcrr(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],
 *            double coef[],int *n,int totnc[],double totncv[],int *iter,
 *            int *iterlim,double *steptol,double sd[],double fd[],double reg[],
 *            int *transform,double p[],double q[],int *rank,double aux[],
 *            int pivot[],double work[],double diff[],int *info,double diag[],
 *            double hessian[],double *loglik,double fit[],int pred[],
 *            double cpred[])
 * void mlldcrr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,
 *              double beta[],int *transform,double *loglik,double top[],
 *              double bot[])
 * void fvdcrr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,
 *             double beta[],int *transform,double p[],double q[],double fit[],
 *             int pred[],double cpred[])
 *
 * void lpor(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],
 *           double coef[],int *n,int totnc[],double totncv[],int *iter,
 *           int *iterlim,double *steptol,double sd[],double fd[],double reg[],
 *           int *transform,int *link,double p[],double lambda[],double q[],
 *           int *rank,double aux[],int pivot[],double work[],double diff[],
 *           int *info,double diag[],double hessian[],double *loglik,
 *           double fit[],int pred[],double cpred[])
 * void mllpor(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,
 *             double beta[],int *transform,int *link,double lambda[],
 *             double *loglik,double p[],double q[])
 * void fvpor(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,
 *            double beta[],int *transform,int *link,double lambda[],
 *            double p[],double q[],double fit[],int pred[],double cpred[])
 *
 * void lacr(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],
 *           double coef[],int *n,int totnc[],double totncv[],int *iter,
 *           int *iterlim,double *steptol,double sd[],double fd[],double reg[],
 *           int *transform,double p[],int *rank,double aux[],int pivot[],
 *           double work[],double diff[],int *info,double diag[],
 *           double hessian[],double *loglik,double fit[],int pred[],
 *           double cpred[])
 * void mllacr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,
 *             double beta[],int *transform,double *loglik,double top[])
 * void fvacr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,
 *            double beta[],int *transform,double p[],double fit[],
 *            int pred[],double cpred[])
 *
 *  DESCRIPTION
 *
 *    Functions to compute linear regresions, the minus log-likelihood, and the fitted
 * values for the following categorical distributions:
 * Binary, Multinomial, Simplified Multinomial, Upwards Continuation-ratio,
 * Downwards Continuation-ratio, Proportional Odds, and Adjacent Categories.
 *
 */

#include <stdio.h>
#include <math.h>
#include <stddef.h>
#include "R.h"
#include "Rmath.h"

extern int dqrdc2_(double x[],int *ldx,int *n,int *p,double *tol,int *k,double qraux[],int jpvt[],double work[]);
extern int dqrcf_(double x[],int *n,int *k,double qraux[],double y[],int *ny,double b[],int *info);
/* Links:
logit, probit, loglog, cloglog, uniform, log-normal, exponential, Pareto, Cauchy, Laplace, Levy, simplex,
gamma, Weibull, inverse_Gauss, t, chi-square, gen_logistic, gen_extreme_value, Box-Cox, power_exponential,
Burr, Hjorth, beta, stable, gen_gamma, gen_Weibull, gen_inverse_Gauss, F, nc_t, nc_chi-square, Tukey,
nc_beta, nc_F
*/
extern double pnorm(double x,double mean,double sd,int lower_tail,int log_p);
extern double punif(double x,double a,double b,int lower_tail,int log_p);
extern double plnorm(double x,double logmean,double logsd,int lower_tail,int log_p);
extern double pexp(double x,double scale,int lower_tail,int log_p);
extern double ppareto(double q,double m,double s);
extern double pcauchy(double x,double location,double scale,int lower_tail,int log_p);
extern double plaplace(double q,double m,double s);
extern double plevy(double y,double m,double s);
extern double psimplex(double y,double m,double s);
extern double pgamma(double x,double alph,double scale,int lower_tail,int log_p);
extern double pweibull(double x,double shape,double scale,int lower_tail,int log_p);
extern double pinvgauss(double q,double m,double s);
extern double pt(double x,double n,int lower_tail,int log_p);
extern double pchisq(double x,double df,int lower_tail,int log_p);
extern double pglogis(double q,double m,double s,double f);
extern double pgextval(double q,double s,double m,double f);
extern double pboxcox(double q,double m,double s,double f);
extern double ppowexp(double y,double m,double s,double f);
extern double pburr(double q,double m,double s,double f);
extern double phjorth(double q,double m,double s,double f);
extern double pbeta(double x,double pin,double qin,int lower_tail,int log_p);
extern double pstable(double y,double loc,double disp,double beta,double alpha);
extern double pggamma(double q,double s,double m,double f);
extern double pgweibull(double q,double s,double m,double f);
extern double pginvgauss(double y,double m,double s,double f);
extern double pf(double x,double n1,double n2,int lower_tail,int log_p);
extern double pnt(double t,double df,double delta,int lower_tail,int log_p);
extern double pnchisq(double x,double f,double theta,int lower_tail,int log_p);
extern double ptukey(double q,double rr,double cc,double df,int lower_tail,int log_p);
extern double pnbeta(double x,double a,double b,double lambda,int lower_tail,int log_p);
extern double pnf(double x,double n1,double n2,double ncp,int lower_tail,int log_p);

/* fd[] contains score vector (first derivative). */
/* sd[] contains minus Fisher's information matrix (second derivative). */
/* hessian contains the hessian to obtain the standard errors of the converged model. */
/* loglik contains the log-likelihood of the converged model. */

/* Linear binary regresion. */
void lbr(int y[],int freq[],double cov[],int *ncv,double beta[],double coef[],int *n,double totncv[],int *iter,int *iterlim,double *steptol,double sd[],double fd[],
         int *transform,int *link,double lambda[],int *rank,double aux[],int pivot[],double work[],double diff[],int *info,double diag[],double hessian[],
         double *loglik,double fit[],int pred[],double cpred[]);
/* Minus log-likelihood for multinomial regresion. */
void mllbr(int y[],int freq[],double cov[],int *ncv,int *n,double beta[],int *transform,int *link,double lambda[],double *loglik);
/* Fitted values for multinomial regresion. */
void fvbr(int y[],int freq[],double cov[],int *ncv,int *n,double beta[],int *transform,int *link,double lambda[],double fit[],int pred[],double cpred[]);

/* Linear multinomial regression. */
void lmr(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],double coef[],int *n,int totnc[],double totncv[],int *iter,int *iterlim,double *steptol,
         double sd[],double fd[],double reg[],int *transform,double p[],int *rank,double aux[],int pivot[],double work[],double diff[],int *info,double diag[],
         double hessian[],double *loglik,double fit[],int pred[],double cpred[]);
/* Minus log-likelihood for multinomial regresion. */
void mllmr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double *loglik,double top[]);
/* Fitted values for multinomial regresion. */
void fvmr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double p[],double fit[],int pred[],double cpred[]);

/* Linear simplified multinomial regression. */
void lsmr(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],double coef[],int *n,int totnc[],double totncv[],int *iter,int *iterlim,double *steptol,
          double sd[],double fd[],double reg[],int *transform,double p[],int *rank,double aux[],int pivot[],double work[],double diff[],int *info,double diag[],
          double hessian[],double *loglik,double fit[],int pred[],double cpred[]);
/* Minus log-likelihood for simplified multinomial regresion. */
void mllsmr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double *loglik,double top[]);
/* Fitted values for simplified multinomial regresion. */
void fvsmr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double p[],double fit[],int pred[],double cpred[]);

/* Linear upwards continuation-ratio regresion. */
void lucrr(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],double coef[],int *n,int totnc[],double totncv[],int *iter,int *iterlim,double *steptol,
           double sd[],double fd[],double reg[],int *transform,double p[],double q[],int *rank,double aux[],int pivot[],double work[],double diff[],int *info,
           double diag[],double hessian[],double *loglik,double fit[],int pred[],double cpred[]);
/* Minus log-likelihood for upwards continuation-ratio regresion. */
void mllucrr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double *loglik,double top[],double bot[]);
/* Fitted values for upwards continuation-ratio regresion. */
void fvucrr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double p[],double q[],double fit[],int pred[],
            double cpred[]);

/* Linear downwards continuation-ratio regresion. */
void ldcrr(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],double coef[],int *n,int totnc[],double totncv[],int *iter,int *iterlim,double *steptol,
           double sd[],double fd[],double reg[],int *transform,double p[],double q[],int *rank,double aux[],int pivot[],double work[],double diff[],int *info,
           double diag[],double hessian[],double *loglik,double fit[],int pred[],double cpred[]);
/* Minus log-likelihood for downwards continuation-ratio regresion. */
void mlldcrr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double *loglik,double top[],double bot[]);
/* Fitted values for downwards continuation-ratio regresion. */
void fvdcrr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double p[],double q[],double fit[],int pred[],
            double cpred[]);

/* Linear proportional odds regresion. */
void lpor(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],double coef[],int *n,int totnc[],double totncv[],int *iter,int *iterlim,double *steptol,
          double sd[],double fd[],double reg[],int *transform,int *link,double p[],double lambda[],double q[],int *rank,double aux[],int pivot[],double work[],
          double diff[],int *info,double diag[],double hessian[],double *loglik,double fit[],int pred[],double cpred[]);
/* Minus log-likelihood for proportional odds regresion. */
void mllpor(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,int *link,double lambda[],double *loglik,double p[],double q[]);
/* Fitted values for proportional odds regresion. */
void fvpor(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,int *link,double lambda[],double p[],double q[],double fit[],
           int pred[],double cpred[]);

/* Linear adjacent categories regresion. */
void lacr(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],double coef[],int *n,int totnc[],double totncv[],int *iter,int *iterlim,double *steptol,
          double sd[],double fd[],double reg[],int *transform,double p[],int *rank,double aux[],int pivot[],double work[],double diff[],int *info,double diag[],
          double hessian[],double *loglik,double fit[],int pred[],double cpred[]);
/* Minus log-likelihood for adjacent categories regresion. */
void mllacr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double *loglik,double top[]);
/* Fitted values for adjacent categories regresion. */
void fvacr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double p[],double fit[],int pred[],double cpred[]);

/* Linear binary regresion. */
void lbr(int y[],int freq[],double cov[],int *ncv,double beta[],double coef[],int *n,double totncv[],int *iter,int *iterlim,double *steptol,double sd[],double fd[],
         int *transform,int *link,double lambda[],int *rank,double aux[],int pivot[],double work[],double diff[],int *info,double diag[],double hessian[],
         double *loglik,double fit[],int pred[],double cpred[]) {
  int ny,ncncv,i,totnc,j,k;
  double tol,reg,p,q;

  ny=1;
  tol=1.;
  ncncv=*ncv+1;
  for(i=0;i<ncncv;i++)
    beta[i]=coef[i];
  for(i=0,totnc=0;i<*n;i++)
    totnc+=y[i]*freq[i];
  for(i=0;i<*n;i++)
    for(j=0;j<*ncv;j++)
      totncv[j]+=y[i]*freq[i]*cov[*n*j+i];
  while(*iter<*iterlim&&tol>*steptol) {
    for(i=0;i<ncncv;i++) {
      for(j=0;j<ncncv;j++)
        sd[ncncv*i+j]=0.; /* Initialize: I_{aa}, I_{ab_l}, and I_{b_lb_m} */
      fd[i]=0.; /* Initialize: S_{a} and S_{b_l} */
    }
    fd[0]=totnc; /* First term of S_{a} */
    for(i=0;i<*ncv;i++)
      fd[1+i]=totncv[i]; /* First term of S_{b_l} */
    for(i=0;i<*n;i++) {
      for(j=0,reg=beta[0];j<*ncv;j++)
        reg+=beta[1+j]*cov[*n*j+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        reg=reg*reg;
        break;
      case 3: /* log transform */
        reg=exp(reg);
        break;
      case 4: /* logit transform */
        reg=exp(reg);
        reg=reg/(1+reg);
        break;
      case 5: /* cloglog transform */
        reg=1.0-exp(-exp(reg));
        break;
      case 6: /* square transform */
        reg=sqrt(reg);
        break;
      case 7: /* exp transform */
        reg=log(reg);
        break;
      case 8: /* loglog transform */
        reg=exp(-exp(reg));
        break;
      }
      switch(*link) {
      case 1: /* logistic distribution */
        p=exp(reg);
        p=p/(1+p);
        break;
      case 2: /* normal distribution */
        p=pnorm(reg,0.0,1.0,1,0);
        break;
      case 3: /* inverse log-log */
        p=exp(-exp(reg));
        break;
      case 4: /* inverse complementary log-log */
        p=1.0-exp(-exp(reg));
        break;
      case 5: /* uniform distribution */
        p=punif(reg,0.0,1.0,1,0);
        break;
      case 6: /* log-normal distribution */
        p=plnorm(reg,0.0,1.0,1,0);
        break;
      case 7: /* exponential distribution */
        p=pexp(reg,1.0,1,0);
        break;
      case 8: /* Pareto distribution */
        p=ppareto(reg,1.0,2.0);
        break;
      case 9: /* Cauchy distribution */
        p=pcauchy(reg,0.0,1.0,1,0);
        break;
      case 10: /* Laplace distribution */
        p=plaplace(reg,0.0,1.0);
        break;
      case 11: /* Levy distribution */
        p=plevy(reg,0.0,1.0);
        break;
      case 12: /* simplex distribution */
        p=psimplex(reg,0.5,1.0);
        break;
      case 13: /* gamma distribution */
        p=pgamma(reg,lambda[0],1.0/lambda[0],1,0);
        break;
      case 14: /* Weibull distribution */
        p=pweibull(reg,lambda[0],1.0,1,0);
        break;
      case 15: /* inverse Gauss distribution */
        p=pinvgauss(reg,1.0,lambda[0]);
        break;
      case 16: /* t distribution */
        p=pt(reg,lambda[0],1,0);
        break;
      case 17: /* chi-square distribution */
        p=pchisq(reg,lambda[0],1,0);
        break;
      case 18: /* generalized logistic distribution */
        p=pglogis(reg,0.0,1.0,lambda[0]);
        break;
      case 19: /* generalized extreme value distribution */
        p=pgextval(reg,1.0,1.0,lambda[0]);
        break;
      case 20: /* Box-Cox distribution */
        p=pboxcox(reg,1.0,1.0,lambda[0]);
        break;
      case 21: /* power exponential distribution */
        p=ppowexp(reg,0.0,1.0,lambda[0]);
        break;
      case 22: /* Burr distribution */
        p=pburr(reg,1.0,1.0,lambda[0]);
        break;
      case 23: /* Hjorth distribution */
        p=phjorth(reg,1.0,1.0,lambda[0]);
        break;
      case 24: /* beta distribution */
        p=pbeta(reg,lambda[0],lambda[1],1,0);
        break;
      case 25: /* stable distribution */
        p=pstable(reg,0.0,1.0,lambda[0],lambda[1]);
        break;
      case 26: /* generalized gamma distribution */
        p=pggamma(reg,lambda[0],1.0,lambda[1]);
        break;
      case 27: /* generalized Weibull distribution (Mudholkar et al, 1995) */
        p=pgweibull(reg,lambda[0],1.0,lambda[1]);
        break;
      case 28: /* generalized inverse Gauss distribution */
        p=pginvgauss(reg,1.0,lambda[0],lambda[1]);
        break;
      case 29: /* F distribution */
        p=pf(reg,lambda[0],lambda[1],1,0);
        break;
      case 30: /* non-central t distribution */
        p=pnt(reg,lambda[0],lambda[1],1,0);
        break;
      case 31: /* non-central chi-square distribution */
        p=pnchisq(reg,lambda[0],lambda[1],1,0);
        break;
      case 32: /* studentized range distribution */
        p=ptukey(reg,lambda[2],lambda[1],lambda[0],1,0);
        break;
      case 33: /* non-central beta distribution */
        p=pnbeta(reg,lambda[0],lambda[1],lambda[2],1,0);
        break;
      case 34: /* non-central F distribution */
        p=pnf(reg,lambda[0],lambda[1],lambda[2],1,0);
        break;
      }
      q=1.-p;
      fd[0]-=freq[i]*p; /* Second term of S_{a} */
      sd[0]-=freq[i]*p*q; /* I_{aa} */
      for(j=0;j<*ncv;j++) {
        fd[(j+1)]-=freq[i]*cov[*n*j+i]*p; /* Second term of  S_{b_l} */
        sd[j+1]-=freq[i]*cov[*n*j+i]*p*q; /* I_{ab_l} */
        for(k=j;k<*ncv;k++)
          sd[ncncv*(j+1)+(k+1)]-=freq[i]*cov[*n*j+i]*cov[*n*k+i]*p*q; /* I_{b_lb_m} */
      }
    }
    for(i=0;i<(ncncv-1);i++)
      for(j=(i+1);j<ncncv;j++)
        sd[ncncv*j+i]=sd[ncncv*i+j]; /* Complete symmetry of second derivative matrix */
    dqrdc2_(sd,&ncncv,&ncncv,&ncncv,steptol,rank,aux,pivot,work);
    dqrcf_(sd,&ncncv,rank,aux,fd,&ny,diff,info);
    for(i=0,tol=0.;i<ncncv;i++) {
      if(fabs(diff[i])>tol)
        tol=fabs(diff[i]);
      beta[i]-=diff[i];
    }
    (*iter)++;
    if(*info!=0)
      break;
    if(*rank!=ncncv)
      break;
  }
  dqrcf_(sd,&ncncv,rank,aux,diag,&ncncv,hessian,info);
  for(i=0;i<*n;i++) {
    for(j=0,reg=beta[0];j<*ncv;j++)
      reg+=beta[1+j]*cov[*n*j+i];
    switch(*transform) {
    case 2: /* sqrt transform */
      reg=reg*reg;
      break;
    case 3: /* log transform */
      reg=exp(reg);
      break;
    case 4: /* logit transform */
      reg=exp(reg);
      reg=reg/(1+reg);
      break;
    case 5: /* cloglog transform */
      reg=1.0-exp(-exp(reg));
      break;
    case 6: /* square transform */
      reg=sqrt(reg);
      break;
    case 7: /* exp transform */
      reg=log(reg);
      break;
    case 8: /* loglog transform */
      reg=exp(-exp(reg));
      break;
    }
    switch(*link) {
    case 1: /* logistic distribution */
      p=exp(reg);
      p=p/(1+p);
      break;
    case 2: /* normal distribution */
      p=pnorm(reg,0.0,1.0,1,0);
      break;
    case 3: /* inverse log-log */
      p=exp(-exp(reg));
      break;
    case 4: /* inverse complementary log-log */
      p=1.0-exp(-exp(reg));
      break;
    case 5: /* uniform distribution */
      p=punif(reg,0.0,1.0,1,0);
      break;
    case 6: /* log-normal distribution */
      p=plnorm(reg,0.0,1.0,1,0);
      break;
    case 7: /* exponential distribution */
      p=pexp(reg,1.0,1,0);
      break;
    case 8: /* Pareto distribution */
      p=ppareto(reg,1.0,2.0);
      break;
    case 9: /* Cauchy distribution */
      p=pcauchy(reg,0.0,1.0,1,0);
      break;
    case 10: /* Laplace distribution */
      p=plaplace(reg,0.0,1.0);
      break;
    case 11: /* Levy distribution */
      p=plevy(reg,0.0,1.0);
      break;
    case 12: /* simplex distribution */
      p=psimplex(reg,0.5,1.0);
      break;
    case 13: /* gamma distribution */
      p=pgamma(reg,lambda[0],1.0/lambda[0],1,0);
      break;
    case 14: /* Weibull distribution */
      p=pweibull(reg,lambda[0],1.0,1,0);
      break;
    case 15: /* inverse Gauss distribution */
      p=pinvgauss(reg,1.0,lambda[0]);
      break;
    case 16: /* t distribution */
      p=pt(reg,lambda[0],1,0);
      break;
    case 17: /* chi-square distribution */
      p=pchisq(reg,lambda[0],1,0);
      break;
    case 18: /* generalized logistic distribution */
      p=pglogis(reg,0.0,1.0,lambda[0]);
      break;
    case 19: /* generalized extreme value distribution */
      p=pgextval(reg,1.0,1.0,lambda[0]);
      break;
    case 20: /* Box-Cox distribution */
      p=pboxcox(reg,1.0,1.0,lambda[0]);
      break;
    case 21: /* power exponential distribution */
      p=ppowexp(reg,0.0,1.0,lambda[0]);
      break;
    case 22: /* Burr distribution */
      p=pburr(reg,1.0,1.0,lambda[0]);
      break;
    case 23: /* Hjorth distribution */
      p=phjorth(reg,1.0,1.0,lambda[0]);
      break;
    case 24: /* beta distribution */
      p=pbeta(reg,lambda[0],lambda[1],1,0);
      break;
    case 25: /* stable distribution */
      p=pstable(reg,0.0,1.0,lambda[0],lambda[1]);
      break;
    case 26: /* generalized gamma distribution */
      p=pggamma(reg,lambda[0],1.0,lambda[1]);
      break;
    case 27: /* generalized Weibull distribution (Mudholkar et al, 1995) */
      p=pgweibull(reg,lambda[0],1.0,lambda[1]);
      break;
    case 28: /* generalized inverse Gauss distribution */
      p=pginvgauss(reg,1.0,lambda[0],lambda[1]);
      break;
    case 29: /* F distribution */
      p=pf(reg,lambda[0],lambda[1],1,0);
      break;
    case 30: /* non-central t distribution */
      p=pnt(reg,lambda[0],lambda[1],1,0);
      break;
    case 31: /* non-central chi-square distribution */
      p=pnchisq(reg,lambda[0],lambda[1],1,0);
      break;
    case 32: /* studentized range distribution */
      p=ptukey(reg,lambda[2],lambda[1],lambda[0],1,0);
      break;
    case 33: /* non-central beta distribution */
      p=pnbeta(reg,lambda[0],lambda[1],lambda[2],1,0);
      break;
    case 34: /* non-central F distribution */
      p=pnf(reg,lambda[0],lambda[1],lambda[2],1,0);
      break;
    }
    q=1.-p;
    if(y[i]!=0) {
      fit[i]=p;
      *loglik+=freq[i]*log(p);
    }
    else {
      fit[i]=q;
      *loglik+=freq[i]*log(q);
    }
    pred[i]=p<q?0:1;
    cpred[i]=q;
  }
}

/* Minus log-likelihood for binary regresion. */
void mllbr(int y[],int freq[],double cov[],int *ncv,int *n,double beta[],int *transform,int *link,double lambda[],double *loglik) {
  int i,j;
  double p,q;

  for(i=0;i<*n;i++) {
    for(j=0,p=beta[0];j<*ncv;j++)
      p+=beta[1+j]*cov[*n*j+i];
    switch(*transform) {
    case 2: /* sqrt transform */
      p=p*p;
      break;
    case 3: /* log transform */
      p=exp(p);
      break;
    case 4: /* logit transform */
      p=exp(p);
      p=p/(1+p);
      break;
    case 5: /* cloglog transform */
      p=1.0-exp(-exp(p));
      break;
    case 6: /* square transform */
      p=sqrt(p);
      break;
    case 7: /* exp transform */
      p=log(p);
      break;
    case 8: /* loglog transform */
      p=exp(-exp(p));
      break;
    }
    switch(*link) {
    case 1: /* logistic distribution */
      p=exp(p);
      p=p/(1+p);
      break;
    case 2: /* normal distribution */
      p=pnorm(p,0.0,1.0,1,0);
      break;
    case 3: /* inverse log-log */
      p=exp(-exp(p));
      break;
    case 4: /* inverse complementary log-log */
      p=1.0-exp(-exp(p));
      break;
    case 5: /* uniform distribution */
      p=punif(p,0.0,1.0,1,0);
      break;
    case 6: /* log-normal distribution */
      p=plnorm(p,0.0,1.0,1,0);
      break;
    case 7: /* exponential distribution */
      p=pexp(p,1.0,1,0);
      break;
    case 8: /* Pareto distribution */
      p=ppareto(p,1.0,2.0);
      break;
    case 9: /* Cauchy distribution */
      p=pcauchy(p,0.0,1.0,1,0);
      break;
    case 10: /* Laplace distribution */
      p=plaplace(p,0.0,1.0);
      break;
    case 11: /* Levy distribution */
      p=plevy(p,0.0,1.0);
      break;
    case 12: /* simplex distribution */
      p=psimplex(p,0.5,1.0);
      break;
    case 13: /* gamma distribution */
      p=pgamma(p,lambda[0],1.0/lambda[0],1,0);
      break;
    case 14: /* Weibull distribution */
      p=pweibull(p,lambda[0],1.0,1,0);
      break;
    case 15: /* inverse Gauss distribution */
      p=pinvgauss(p,1.0,lambda[0]);
      break;
    case 16: /* t distribution */
      p=pt(p,lambda[0],1,0);
      break;
    case 17: /* chi-square distribution */
      p=pchisq(p,lambda[0],1,0);
      break;
    case 18: /* generalized logistic distribution */
      p=pglogis(p,0.0,1.0,lambda[0]);
      break;
    case 19: /* generalized extreme value distribution */
      p=pgextval(p,1.0,1.0,lambda[0]);
      break;
    case 20: /* Box-Cox distribution */
      p=pboxcox(p,1.0,1.0,lambda[0]);
      break;
    case 21: /* power exponential distribution */
      p=ppowexp(p,0.0,1.0,lambda[0]);
      break;
    case 22: /* Burr distribution */
      p=pburr(p,1.0,1.0,lambda[0]);
      break;
    case 23: /* Hjorth distribution */
      p=phjorth(p,1.0,1.0,lambda[0]);
      break;
    case 24: /* beta distribution */
      p=pbeta(p,lambda[0],lambda[1],1,0);
      break;
    case 25: /* stable distribution */
      p=pstable(p,0.0,1.0,lambda[0],lambda[1]);
      break;
    case 26: /* generalized gamma distribution */
      p=pggamma(p,lambda[0],1.0,lambda[1]);
      break;
    case 27: /* generalized Weibull distribution (Mudholkar et al, 1995) */
      p=pgweibull(p,lambda[0],1.0,lambda[1]);
      break;
    case 28: /* generalized inverse Gauss distribution */
      p=pginvgauss(p,1.0,lambda[0],lambda[1]);
      break;
    case 29: /* F distribution */
      p=pf(p,lambda[0],lambda[1],1,0);
      break;
    case 30: /* non-central t distribution */
      p=pnt(p,lambda[0],lambda[1],1,0);
      break;
    case 31: /* non-central chi-square distribution */
      p=pnchisq(p,lambda[0],lambda[1],1,0);
      break;
    case 32: /* studentized range distribution */
      p=ptukey(p,lambda[2],lambda[1],lambda[0],1,0);
      break;
    case 33: /* non-central beta distribution */
      p=pnbeta(p,lambda[0],lambda[1],lambda[2],1,0);
      break;
    case 34: /* non-central F distribution */
      p=pnf(p,lambda[0],lambda[1],lambda[2],1,0);
      break;
    }
    q=1.-p;
    if(y[i]!=0)
      *loglik+=freq[i]*log(p);
    else
      *loglik+=freq[i]*log(q);
  }
}

/* Fitted values for binary regresion. */
void fvbr(int y[],int freq[],double cov[],int *ncv,int *n,double beta[],int *transform,int *link,double lambda[],double fit[],int pred[],double cpred[]) {
  int i,j;
  double p,q;

  for(i=0;i<*n;i++) {
    for(j=0,p=beta[0];j<*ncv;j++)
      p+=beta[1+j]*cov[*n*j+i];
    switch(*transform) {
    case 2: /* sqrt transform */
      p=p*p;
      break;
    case 3: /* log transform */
      p=exp(p);
      break;
    case 4: /* logit transform */
      p=exp(p);
      p=p/(1+p);
      break;
    case 5: /* cloglog transform */
      p=1.0-exp(-exp(p));
      break;
    case 6: /* square transform */
      p=sqrt(p);
      break;
    case 7: /* exp transform */
      p=log(p);
      break;
    case 8: /* loglog transform */
      p=exp(-exp(p));
      break;
    }
    switch(*link) {
    case 1: /* logistic distribution */
      p=exp(p);
      p=p/(1+p);
      break;
    case 2: /* normal distribution */
      p=pnorm(p,0.0,1.0,1,0);
      break;
    case 3: /* inverse log-log */
      p=exp(-exp(p));
      break;
    case 4: /* inverse complementary log-log */
      p=1.0-exp(-exp(p));
      break;
    case 5: /* uniform distribution */
      p=punif(p,0.0,1.0,1,0);
      break;
    case 6: /* log-normal distribution */
      p=plnorm(p,0.0,1.0,1,0);
      break;
    case 7: /* exponential distribution */
      p=pexp(p,1.0,1,0);
      break;
    case 8: /* Pareto distribution */
      p=ppareto(p,1.0,2.0);
      break;
    case 9: /* Cauchy distribution */
      p=pcauchy(p,0.0,1.0,1,0);
      break;
    case 10: /* Laplace distribution */
      p=plaplace(p,0.0,1.0);
      break;
    case 11: /* Levy distribution */
      p=plevy(p,0.0,1.0);
      break;
    case 12: /* simplex distribution */
      p=psimplex(p,0.5,1.0);
      break;
    case 13: /* gamma distribution */
      p=pgamma(p,lambda[0],1.0/lambda[0],1,0);
      break;
    case 14: /* Weibull distribution */
      p=pweibull(p,lambda[0],1.0,1,0);
      break;
    case 15: /* inverse Gauss distribution */
      p=pinvgauss(p,1.0,lambda[0]);
      break;
    case 16: /* t distribution */
      p=pt(p,lambda[0],1,0);
      break;
    case 17: /* chi-square distribution */
      p=pchisq(p,lambda[0],1,0);
      break;
    case 18: /* generalized logistic distribution */
      p=pglogis(p,0.0,1.0,lambda[0]);
      break;
    case 19: /* generalized extreme value distribution */
      p=pgextval(p,1.0,1.0,lambda[0]);
      break;
    case 20: /* Box-Cox distribution */
      p=pboxcox(p,1.0,1.0,lambda[0]);
      break;
    case 21: /* power exponential distribution */
      p=ppowexp(p,0.0,1.0,lambda[0]);
      break;
    case 22: /* Burr distribution */
      p=pburr(p,1.0,1.0,lambda[0]);
      break;
    case 23: /* Hjorth distribution */
      p=phjorth(p,1.0,1.0,lambda[0]);
      break;
    case 24: /* beta distribution */
      p=pbeta(p,lambda[0],lambda[1],1,0);
      break;
    case 25: /* stable distribution */
      p=pstable(p,0.0,1.0,lambda[0],lambda[1]);
      break;
    case 26: /* generalized gamma distribution */
      p=pggamma(p,lambda[0],1.0,lambda[1]);
      break;
    case 27: /* generalized Weibull distribution (Mudholkar et al, 1995) */
      p=pgweibull(p,lambda[0],1.0,lambda[1]);
      break;
    case 28: /* generalized inverse Gauss distribution */
      p=pginvgauss(p,1.0,lambda[0],lambda[1]);
      break;
    case 29: /* F distribution */
      p=pf(p,lambda[0],lambda[1],1,0);
      break;
    case 30: /* non-central t distribution */
      p=pnt(p,lambda[0],lambda[1],1,0);
      break;
    case 31: /* non-central chi-square distribution */
      p=pnchisq(p,lambda[0],lambda[1],1,0);
      break;
    case 32: /* studentized range distribution */
      p=ptukey(p,lambda[2],lambda[1],lambda[0],1,0);
      break;
    case 33: /* non-central beta distribution */
      p=pnbeta(p,lambda[0],lambda[1],lambda[2],1,0);
      break;
    case 34: /* non-central F distribution */
      p=pnf(p,lambda[0],lambda[1],lambda[2],1,0);
      break;
    }
    q=1.-p;
    if(y[i]!=0)
      fit[i]=p;
    else
      fit[i]=q;
    pred[i]=p<q?0:1;
    cpred[i]=q;
  }
}

/* Linear multinomial regression. */
void lmr(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],double coef[],int *n,int totnc[],double totncv[],int *iter,int *iterlim,double *steptol,
         double sd[],double fd[],double reg[],int *transform,double p[],int *rank,double aux[],int pivot[],double work[],double diff[],int *info,double diag[],
         double hessian[],double *loglik,double fit[],int pred[],double cpred[]) {
  int ny,ncncv,i,j,k,l,m;
  double tol,bot,q;

  ny=1;
  tol=1.;
  ncncv=(*nc-1)*(*ncv+1);
  for(i=0;i<ncncv;i++)
    beta[i]=coef[i];
  for(i=0;i<*n;i++)
    if(y[i]!=0)
      totnc[y[i]-1]+=freq[i];
  for(i=0;i<*n;i++)
    for(j=0;j<*ncv;j++)
      if(y[i]!=0)
        totncv[*ncv*(y[i]-1)+j]+=freq[i]*cov[*n*j+i];
  while(*iter<*iterlim&&tol>*steptol) {
    for(i=0;i<ncncv;i++) {
      for(j=0;j<ncncv;j++)
        sd[ncncv*i+j]=0.; /* Initialize: I_{a_ia_i}, I_{a_ia_l}, I_{a_ib_{il}}, I_{a_ib_{lm}}, I_{b_{il}b_{im}}, and I_{b_{im}b_{ls}} */
      fd[i]=0.; /* Initialize: S_{a_i} and S_{b_{il}} */
    }
    for(i=0;i<(*nc-1);i++) {
      fd[(*ncv+1)*i]=totnc[i]; /* First term of S_{a_i} */
      for(j=0;j<*ncv;j++)
        fd[1+(*ncv+1)*i+j]=totncv[*ncv*i+j]; /* First term of S_{b_{il}} */
    }
    for(i=0;i<*n;i++) {
      for(j=0,bot=1;j<(*nc-1);j++) {
        for(k=0,reg[j]=beta[(*ncv+1)*j];k<*ncv;k++)
          reg[j]+=beta[1+(*ncv+1)*j+k]*cov[*n*k+i];
        switch(*transform) {
        case 2: /* sqrt transform */
          reg[j]=reg[j]*reg[j];
          break;
        case 3: /* log transform */
          reg[j]=exp(reg[j]);
          break;
        case 4: /* logit transform */
          reg[j]=exp(reg[j]);
          reg[j]=reg[j]/(1+reg[j]);
          break;
        case 5: /* cloglog transform */
          reg[j]=1.0-exp(-exp(reg[j]));
          break;
        case 6: /* square transform */
          reg[j]=sqrt(reg[j]);
          break;
        case 7: /* exp transform */
          reg[j]=log(reg[j]);
          break;
        case 8: /* loglog transform */
          reg[j]=exp(-exp(reg[j]));
          break;
        }
        bot+=exp(reg[j]);
      }
      for(j=0;j<(*nc-1);j++)
        p[j]=exp(reg[j])/bot;
      q=1/bot;
      for(j=0;j<(*nc-1);j++) {
        fd[(*ncv+1)*j]-=freq[i]*p[j]; /* Second term of S_{a_i} */
        sd[(1+ncncv)*(*ncv+1)*j]-=freq[i]*p[j]*(1-p[j]); /* I_{a_ia_i} */
        for(k=0;k<*ncv;k++) {
          fd[1+(*ncv+1)*j+k]-=freq[i]*cov[*n*k+i]*p[j]; /* Second term of S_{b_{il}} */
          sd[1+(1+ncncv)*(*ncv+1)*j+k]-=freq[i]*cov[*n*k+i]*p[j]*(1-p[j]); /* I_{a_ib_{il}} */
          sd[(ncncv+1)*(k+1+(*ncv+1)*j)]-=freq[i]*cov[*n*k+i]*cov[*n*k+i]*p[j]*(1-p[j]); /* I_{b_{il}b_{im}} */
        }
        for(k=0;k<(*ncv-1);k++)
          for(l=k+1;l<*ncv;l++)
            sd[(ncncv+1)*(1+j*(*ncv+1))+k*ncncv+l]-=freq[i]*cov[*n*k+i]*cov[*n*l+i]*p[j]*(1-p[j]); /* I_{b_{il}b_{im}} */
      }
      for(j=0;j<(*nc-2);j++) {
        for(k=(j+1);k<(*nc-1);k++) {
          sd[ncncv*(*ncv+1)*j+(*ncv+1)*k]+=freq[i]*p[j]*p[k]; /* I_{a_ia_l} */
          for(l=0;l<*ncv;l++)
            sd[1+(*ncv+1)*(ncncv*j+k)+l]+=freq[i]*cov[*n*l+i]*p[j]*p[k]; /* I_{a_ib_{lm}} part one */
        }
        for(k=0;k<*ncv;k++)
          for(l=(j+1);l<(*nc-1);l++) {
            sd[ncncv*(1+k)+(*ncv+1)*(l+ncncv*j)]+=freq[i]*cov[*n*k+i]*p[l]*p[j]; /* I_{a_ib_{lm}} part two */
            for(m=0;m<*ncv;m++)
              sd[1+ncncv*(1+k)+(*ncv+1)*(l+ncncv*j)+m]+=freq[i]*cov[*n*k+i]*cov[*n*m+i]*p[j]*p[l]; /* I_{b_{im}b_{ls}} */
          }
      }
    }
    for(i=0;i<(ncncv-1);i++)
      for(j=(i+1);j<ncncv;j++)
        sd[ncncv*j+i]=sd[ncncv*i+j]; /* Complete symmetry of second derivative matrix */
    dqrdc2_(sd,&ncncv,&ncncv,&ncncv,steptol,rank,aux,pivot,work);
    dqrcf_(sd,&ncncv,rank,aux,fd,&ny,diff,info);
    for(i=0,tol=0.;i<ncncv;i++) {
      if(fabs(diff[i])>tol)
        tol=fabs(diff[i]);
      beta[i]-=diff[i];
    }
    (*iter)++;
    if(*info!=0)
      break;
    if(*rank!=ncncv)
      break;
  }
  dqrcf_(sd,&ncncv,rank,aux,diag,&ncncv,hessian,info);
  for(i=0;i<*n;i++) {
    for(j=0,bot=1;j<(*nc-1);j++) {
      for(k=0,reg[j]=beta[(*ncv+1)*j];k<*ncv;k++)
        reg[j]+=beta[1+(*ncv+1)*j+k]*cov[*n*k+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        reg[j]=reg[j]*reg[j];
        break;
      case 3: /* log transform */
        reg[j]=exp(reg[j]);
        break;
      case 4: /* logit transform */
        reg[j]=exp(reg[j]);
        reg[j]=reg[j]/(1+reg[j]);
        break;
      case 5: /* cloglog transform */
        reg[j]=1.0-exp(-exp(reg[j]));
        break;
      case 6: /* square transform */
        reg[j]=sqrt(reg[j]);
        break;
      case 7: /* exp transform */
        reg[j]=log(reg[j]);
        break;
      case 8: /* loglog transform */
        reg[j]=exp(-exp(reg[j]));
        break;
      }
      bot+=exp(reg[j]);
    }
    for(j=0;j<(*nc-1);j++)
      p[j]=exp(reg[j])/bot;
    q=1/bot;
    if(y[i]!=0) {
      fit[i]=p[y[i]-1];
      *loglik+=freq[i]*log(p[y[i]-1]);
    }
    else {
      fit[i]=q;
      *loglik+=freq[i]*log(q);
    }
    pred[i]=1;
    cpred[i]=q;
    for(j=1;j<(*nc-1);j++) {
      if(p[pred[i]-1]<p[j])
        pred[i]=j+1;
      cpred[*n*j+i]=cpred[*n*(j-1)+i]+p[j-1];
    }
    if(p[pred[i]-1]<q)
      pred[i]=0;
  }
}

/* Minus log-likelihood for multinomial regresion. */
void mllmr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double *loglik,double top[]) {
  int i,j,k;
  double bot;

  for(i=0;i<*n;i++) {
    for(j=0,bot=1;j<(*nc-1);j++) {
      for(k=0,top[j]=beta[(*ncv+1)*j];k<*ncv;k++)
        top[j]+=beta[1+(*ncv+1)*j+k]*cov[*n*k+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        top[j]=top[j]*top[j];
        break;
      case 3: /* log transform */
        top[j]=exp(top[j]);
        break;
      case 4: /* logit transform */
        top[j]=exp(top[j]);
        top[j]=top[j]/(1+top[j]);
        break;
      case 5: /* cloglog transform */
        top[j]=1.0-exp(-exp(top[j]));
        break;
      case 6: /* square transform */
        top[j]=sqrt(top[j]);
        break;
      case 7: /* exp transform */
        top[j]=log(top[j]);
        break;
      case 8: /* loglog transform */
        top[j]=exp(-exp(top[j]));
        break;
      }
      bot+=exp(top[j]);
    }
    if(y[i]!=0)
      *loglik-=freq[i]*(top[y[i]-1]-log(bot));
    else
      *loglik+=freq[i]*log(bot);
  }
}

/* Fitted values for multinomial regresion. */
void fvmr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double p[],double fit[],int pred[],double cpred[]) {
  int i,j,k;
  double bot,q;

  for(i=0;i<*n;i++) {
    for(j=0,bot=1;j<(*nc-1);j++) {
      for(k=0,p[j]=beta[(*ncv+1)*j];k<*ncv;k++)
        p[j]+=beta[1+(*ncv+1)*j+k]*cov[*n*k+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        p[j]=p[j]*p[j];
        break;
      case 3: /* log transform */
        p[j]=exp(p[j]);
        break;
      case 4: /* logit transform */
        p[j]=exp(p[j]);
        p[j]=p[j]/(1+p[j]);
        break;
      case 5: /* cloglog transform */
        p[j]=1.0-exp(-exp(p[j]));
        break;
      case 6: /* square transform */
        p[j]=sqrt(p[j]);
        break;
      case 7: /* exp transform */
        p[j]=log(p[j]);
        break;
      case 8: /* loglog transform */
        p[j]=exp(-exp(p[j]));
        break;
      }
      bot+=exp(p[j]);
    }
    for(j=0;j<(*nc-1);j++)
      p[j]=exp(p[j])/bot;
    q=1/bot;
    if(y[i]!=0)
      fit[i]=p[y[i]-1];
    else
      fit[i]=q;
    pred[i]=1;
    cpred[i]=q;
    for(j=1;j<(*nc-1);j++) {
      if(p[pred[i]-1]<p[j])
        pred[i]=j+1;
      cpred[*n*j+i]=cpred[*n*(j-1)+i]+p[j-1];
    }
    if(p[pred[i]-1]<q)
      pred[i]=0;
  }
}

/* Linear simplified multinomial regression. */
void lsmr(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],double coef[],int *n,int totnc[],double totncv[],int *iter,int *iterlim,double *steptol,
          double sd[],double fd[],double reg[],int *transform,double p[],int *rank,double aux[],int pivot[],double work[],double diff[],int *info,double diag[],
          double hessian[],double *loglik,double fit[],int pred[],double cpred[]) {
  int ny,ncncv,i,j,k;
  double tol,bot,q;

  ny=1;
  tol=1.;
  ncncv=*nc-1+*ncv;
  for(i=0;i<ncncv;i++)
    beta[i]=coef[i];
  for(i=0;i<*n;i++)
    if(y[i]!=0)
      totnc[y[i]-1]+=freq[i];
  for(i=0;i<*ncv;i++)
    for(j=0;j<*n;j++)
      if(y[j]!=0)
        totncv[i]+=freq[j]*cov[*n*i+j];
  while(*iter<*iterlim&&tol>*steptol) {
    for(i=0;i<ncncv;i++) {
      for(j=0;j<ncncv;j++)
        sd[ncncv*i+j]=0.; /* Initialize: I_{a_ia_i}, I_{a_ia_l}, I_{a_ib_l}, and I_{b_lb_m} */
      fd[i]=0.; /* Initialize: S_{a_i} and S_{b_l} */
    }
    for(i=0;i<(*nc-1);i++)
      fd[i]=totnc[i]; /* First term of S_{a_i} */
    for(i=0;i<*ncv;i++)
      fd[*nc-1+i]=totncv[i]; /* First term of S_{b_l} */
    for(i=0;i<*n;i++) {
      for(j=0,bot=1;j<(*nc-1);j++) {
        for(k=0,reg[j]=beta[j];k<*ncv;k++)
          reg[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
        switch(*transform) {
        case 2: /* sqrt transform */
          reg[j]=reg[j]*reg[j];
          break;
        case 3: /* log transform */
          reg[j]=exp(reg[j]);
          break;
        case 4: /* logit transform */
          reg[j]=exp(reg[j]);
          reg[j]=reg[j]/(1+reg[j]);
          break;
        case 5: /* cloglog transform */
          reg[j]=1.0-exp(-exp(reg[j]));
          break;
        case 6: /* square transform */
          reg[j]=sqrt(reg[j]);
          break;
        case 7: /* exp transform */
          reg[j]=log(reg[j]);
          break;
        case 8: /* loglog transform */
          reg[j]=exp(-exp(reg[j]));
          break;
        }
        bot+=exp(reg[j]);
      }
      for(j=0,q=1;j<(*nc-1);j++)
        p[j]=exp(reg[j])/bot;
      q=1/bot;
      for(j=0;j<(*nc-1);j++) {
        fd[j]-=freq[i]*p[j]; /* Second term of S_{a_i} */
        sd[(ncncv+1)*j]-=freq[i]*p[j]*(1-p[j]); /* I_{a_ia_i} */
      }
      for(j=0;j<(*nc-2);j++)
        for(k=(j+1);k<(*nc-1);k++)
          sd[ncncv*j+k]+=freq[i]*p[j]*p[k]; /* I_{a_ia_l} */
      for(j=0;j<*ncv;j++) {
        fd[(*nc-1+j)]-=freq[i]*cov[*n*j+i]*(1-q); /* Second term of  S_{b_l} */
        for(k=0;k<(*nc-1);k++)
          sd[ncncv*k+(*nc-1+j)]-=freq[i]*cov[*n*j+i]*p[k]*q; /* I_{a_ib_l} */
        for(k=j;k<*ncv;k++)
          sd[ncncv*(*nc-1+j)+(*nc-1+k)]-=freq[i]*cov[*n*j+i]*cov[*n*k+i]*q*(1-q); /* I_{b_lb_m} */
      }
    }
    for(i=0;i<(ncncv-1);i++)
      for(j=(i+1);j<ncncv;j++)
        sd[ncncv*j+i]=sd[ncncv*i+j]; /* Complete symmetry of second derivative matrix */
    dqrdc2_(sd,&ncncv,&ncncv,&ncncv,steptol,rank,aux,pivot,work);
    dqrcf_(sd,&ncncv,rank,aux,fd,&ny,diff,info);
    for(i=0,tol=0.;i<ncncv;i++) {
      if(fabs(diff[i])>tol)
        tol=fabs(diff[i]);
      beta[i]-=diff[i];
    }
    (*iter)++;
    if(*info!=0)
      break;
    if(*rank!=ncncv)
      break;
  }
  dqrcf_(sd,&ncncv,rank,aux,diag,&ncncv,hessian,info);
  for(i=0;i<*n;i++) {
    for(j=0,bot=1;j<(*nc-1);j++) {
      for(k=0,reg[j]=beta[j];k<*ncv;k++)
        reg[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        reg[j]=reg[j]*reg[j];
        break;
      case 3: /* log transform */
        reg[j]=exp(reg[j]);
        break;
      case 4: /* logit transform */
        reg[j]=exp(reg[j]);
        reg[j]=reg[j]/(1+reg[j]);
        break;
      case 5: /* cloglog transform */
        reg[j]=1.0-exp(-exp(reg[j]));
        break;
      case 6: /* square transform */
        reg[j]=sqrt(reg[j]);
        break;
      case 7: /* exp transform */
        reg[j]=log(reg[j]);
        break;
      case 8: /* loglog transform */
        reg[j]=exp(-exp(reg[j]));
        break;
      }
      bot+=exp(reg[j]);
    }
    for(j=0,q=1;j<(*nc-1);j++)
      p[j]=exp(reg[j])/bot;
    q=1/bot;
    if(y[i]!=0) {
      fit[i]=p[y[i]-1];
      *loglik+=freq[i]*log(p[y[i]-1]);
    }
    else {
      fit[i]=q;
      *loglik+=freq[i]*log(q);
    }
    pred[i]=1;
    cpred[i]=q;
    for(j=1;j<(*nc-1);j++) {
      if(p[pred[i]-1]<p[j])
        pred[i]=j+1;
      cpred[*n*j+i]=cpred[*n*(j-1)+i]+p[j-1];
    }
    if(p[pred[i]-1]<q)
      pred[i]=0;
  }
}

/* Minus log-likelihood for simplified multinomial regresion. */
void mllsmr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double *loglik,double top[]) {
  int i,j,k;
  double bot;

  for(i=0;i<*n;i++) {
    for(j=0,bot=1;j<(*nc-1);j++) {
      for(k=0,top[j]=beta[j];k<*ncv;k++)
        top[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        top[j]=top[j]*top[j];
        break;
      case 3: /* log transform */
        top[j]=exp(top[j]);
        break;
      case 4: /* logit transform */
        top[j]=exp(top[j]);
        top[j]=top[j]/(1+top[j]);
        break;
      case 5: /* cloglog transform */
        top[j]=1.0-exp(-exp(top[j]));
        break;
      case 6: /* square transform */
        top[j]=sqrt(top[j]);
        break;
      case 7: /* exp transform */
        top[j]=log(top[j]);
        break;
      case 8: /* loglog transform */
        top[j]=exp(-exp(top[j]));
        break;
      }
      bot+=exp(top[j]);
    }
    bot=log(bot);
    if(y[i]!=0)
      *loglik-=freq[i]*(top[(y[i]-1)]-bot);
    else
      *loglik+=freq[i]*bot;
  }
}

/* Fitted values for simplified multinomial regresion. */
void fvsmr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double p[],double fit[],int pred[],double cpred[]) {
  int i,j,k;
  double bot,q;

  for(i=0;i<*n;i++) {
    for(j=0,bot=1;j<(*nc-1);j++) {
      for(k=0,p[j]=beta[j];k<*ncv;k++)
        p[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        p[j]=p[j]*p[j];
        break;
      case 3: /* log transform */
        p[j]=exp(p[j]);
        break;
      case 4: /* logit transform */
        p[j]=exp(p[j]);
        p[j]=p[j]/(1+p[j]);
        break;
      case 5: /* cloglog transform */
        p[j]=1.0-exp(-exp(p[j]));
        break;
      case 6: /* square transform */
        p[j]=sqrt(p[j]);
        break;
      case 7: /* exp transform */
        p[j]=log(p[j]);
        break;
      case 8: /* loglog transform */
        p[j]=exp(-exp(p[j]));
        break;
      }
      bot+=exp(p[j]);
    }
    for(j=0,q=1;j<(*nc-1);j++)
      p[j]=exp(p[j])/bot;
    q=1/bot;
    if(y[i]!=0)
      fit[i]=p[y[i]-1];
    else
      fit[i]=q;
    pred[i]=1;
    cpred[i]=q;
    for(j=1;j<(*nc-1);j++) {
      if(p[pred[i]-1]<p[j])
        pred[i]=j+1;
      cpred[*n*j+i]=cpred[*n*(j-1)+i]+p[j-1];
    }
    if(p[pred[i]-1]<q)
      pred[i]=0;
  }
}

/* Linear upwards continuation-ratio regresion. */
void lucrr(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],double coef[],int *n,int totnc[],double totncv[],int *iter,int *iterlim,double *steptol,
           double sd[],double fd[],double reg[],int *transform,double p[],double q[],int *rank,double aux[],int pivot[],double work[],double diff[],int *info,
           double diag[],double hessian[],double *loglik,double fit[],int pred[],double cpred[]) {
  int ny,ncncv,i,j,k,l;
  double tol;

  ny=1;
  tol=1.;
  ncncv=*nc-1+*ncv;
  for(i=0;i<ncncv;i++)
    beta[i]=coef[i];
  for(i=0;i<*n;i++)
    if(y[i]!=(*nc-1))
      totnc[y[i]]+=freq[i];
  for(i=0;i<*ncv;i++)
    for(j=0;j<*n;j++)
      if(y[j]!=(*nc-1))
        totncv[i]+=freq[j]*cov[*n*i+j]*(*nc-1-y[j]);
  while(*iter<*iterlim&&tol>*steptol) {
    for(i=0;i<ncncv;i++) {
      for(j=0;j<ncncv;j++)
        sd[ncncv*i+j]=0.; /* Initialize: I_{a_ia_i}, I_{a_ia_l}, I_{a_ib_l}, and I_{b_lb_m} */
      fd[i]=0.; /* Initialize: S_{a_i} and S_{b_l} */
    }
    for(i=0;i<(*nc-1);i++)
      for(j=0;j<(i+1);j++)
        fd[i]+=totnc[j]; /* First term of S_{a_i} */
    for(i=0;i<*ncv;i++)
      fd[*nc-1+i]=totncv[i]; /* First term of S_{b_l} */
    for(i=0;i<*n;i++) {
      for(j=0;j<(*nc-1);j++) {
        for(k=0,reg[j]=beta[j];k<*ncv;k++)
          reg[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
        switch(*transform) {
        case 2: /* sqrt transform */
          reg[j]=reg[j]*reg[j];
          break;
        case 3: /* log transform */
          reg[j]=exp(reg[j]);
          break;
        case 4: /* logit transform */
          reg[j]=exp(reg[j]);
          reg[j]=reg[j]/(1+reg[j]);
          break;
        case 5: /* cloglog transform */
          reg[j]=1.0-exp(-exp(reg[j]));
          break;
        case 6: /* square transform */
          reg[j]=sqrt(reg[j]);
          break;
        case 7: /* exp transform */
          reg[j]=log(reg[j]);
          break;
        case 8: /* loglog transform */
          reg[j]=exp(-exp(reg[j]));
          break;
        }
        p[j]=exp(reg[j]);
        p[j]=p[j]/(1+p[j]);
        q[j]=1-p[j];
      }
      for(j=y[i];j<(*nc-1);j++) {
        fd[j]-=freq[i]*p[j]; /* Second term of S_{a_i} */
        sd[(ncncv+1)*j]-=freq[i]*p[j]*q[j]; /* First term of I_{a_ia_i} */
      }
      if(y[i]!=0) {
        fd[(y[i]-1)]-=freq[i]*p[(y[i]-1)]; /* Third term of S_{a_i} */
        sd[(ncncv+1)*(y[i]-1)]-=freq[i]*p[(y[i]-1)]*q[(y[i]-1)]; /* Second term of I_{a_ia_i} */
      }
      for(j=0;j<*ncv;j++) {
        for(k=y[i];k<(*nc-1);k++) {
          fd[(*nc-1+j)]-=freq[i]*cov[*n*j+i]*p[k]; /* Second term of S_{b_l} */
          sd[ncncv*k+(*nc-1+j)]-=freq[i]*cov[*n*j+i]*p[k]*q[k]; /* First term of I_{a_ib_l} */
        }
        if(y[i]!=0) {
          fd[(*nc-1+j)]-=freq[i]*cov[*n*j+i]*p[(y[i]-1)]; /* Third term of S_{b_l} */
          sd[ncncv*(y[i]-1)+(*nc-1+j)]-=freq[i]*cov[*n*j+i]*p[(y[i]-1)]*q[(y[i]-1)]; /* Second term of I_{a_ib_l} */
        }
        for(k=j;k<*ncv;k++) {
          for(l=y[i];l<(*nc-1);l++)
            sd[ncncv*(*nc-1+j)+(*nc-1+k)]-=freq[i]*cov[*n*j+i]*cov[*n*k+i]*p[l]*q[l]; /* First term of I_{b_lb_m} */
          if(y[i]!=0)
            sd[ncncv*(*nc-1+j)+(*nc-1+k)]-=freq[i]*cov[*n*j+i]*cov[*n*k+i]*p[(y[i]-1)]*q[(y[i]-1)]; /* Second term of I_{b_lb_m} */
        }
      }
    }
    for(i=0;i<(ncncv-1);i++)
      for(j=(i+1);j<ncncv;j++)
        sd[ncncv*j+i]=sd[ncncv*i+j]; /* Complete symmetry of second derivative matrix */
    dqrdc2_(sd,&ncncv,&ncncv,&ncncv,steptol,rank,aux,pivot,work);
    dqrcf_(sd,&ncncv,rank,aux,fd,&ny,diff,info);
    for(i=0,tol=0.;i<ncncv;i++) {
      if(fabs(diff[i])>tol)
        tol=fabs(diff[i]);
      beta[i]-=diff[i];
    }
    (*iter)++;
    if(*info!=0)
      break;
    if(*rank!=ncncv)
      break;
  }
  dqrcf_(sd,&ncncv,rank,aux,diag,&ncncv,hessian,info);
  for(i=0;i<*n;i++) {
    for(j=0;j<(*nc-1);j++) {
      for(k=0,reg[j]=beta[j];k<*ncv;k++)
        reg[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        reg[j]=reg[j]*reg[j];
        break;
      case 3: /* log transform */
        reg[j]=exp(reg[j]);
        break;
      case 4: /* logit transform */
        reg[j]=exp(reg[j]);
        reg[j]=reg[j]/(1+reg[j]);
        break;
      case 5: /* cloglog transform */
        reg[j]=1.0-exp(-exp(reg[j]));
        break;
      case 6: /* square transform */
        reg[j]=sqrt(reg[j]);
        break;
      case 7: /* exp transform */
        reg[j]=log(reg[j]);
        break;
      case 8: /* loglog transform */
        reg[j]=exp(-exp(reg[j]));
        break;
      }
      p[j]=exp(reg[j]);
      p[j]=p[j]/(1+p[j]);
      q[j]=1-p[j];
    }
    fit[i]=1.;
    for(j=y[i];j<(*nc-1);j++) {
      fit[i]*=p[j];
      *loglik+=freq[i]*log(p[j]);
    }
    if(y[i]!=0) {
      fit[i]*=q[y[i]-1];
      *loglik+=freq[i]*log(q[y[i]-1]);
    }
    for(j=0;j<(*nc-1);j++) {
      cpred[*n*j+i]=1.;
      for(k=j;k<(*nc-1);k++)
        cpred[*n*j+i]*=p[k];
    }
    for(j=1;j<(*nc-1);j++)
      cpred[*n*j+i]*=q[j-1];
    pred[i]=0;
    for(j=1;j<(*nc-1);j++)
      if(cpred[*n*pred[i]+i]<cpred[*n*j+i])
        pred[i]=j;
    if(cpred[*n*pred[i]+i]<q[*nc-2])
      pred[i]=*nc-1;
    for(j=1;j<(*nc-1);j++)
      cpred[*n*j+i]+=cpred[*n*(j-1)+i];
  }
}

/* Minus log-likelihood for upwards continuation-ratio regresion. */
void mllucrr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double *loglik,double top[],double bot[]) {
  int i,j,k;

  for(i=0;i<*n;i++) {
    for(j=0;j<(*nc-1);j++) {
      for(k=0,top[j]=beta[j];k<*ncv;k++)
        top[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        top[j]=top[j]*top[j];
        break;
      case 3: /* log transform */
        top[j]=exp(top[j]);
        break;
      case 4: /* logit transform */
        top[j]=exp(top[j]);
        top[j]=top[j]/(1+top[j]);
        break;
      case 5: /* cloglog transform */
        top[j]=1.0-exp(-exp(top[j]));
        break;
      case 6: /* square transform */
        top[j]=sqrt(top[j]);
        break;
      case 7: /* exp transform */
        top[j]=log(top[j]);
        break;
      case 8: /* loglog transform */
        top[j]=exp(-exp(top[j]));
        break;
      }
      bot[j]=log(1+exp(top[j]));
    }
    for(j=y[i];j<(*nc-1);j++)
      *loglik+=freq[i]*(top[j]-bot[j]);
    if(y[i]!=0)
      *loglik-=freq[i]*bot[(y[i]-1)];
  }
}

/* Fitted values for upwards continuation-ratio regresion. */
void fvucrr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double p[],double q[],double fit[],int pred[],
            double cpred[]) {
  int i,j,k;

  for(i=0;i<*n;i++) {
    for(j=0;j<(*nc-1);j++) {
      for(k=0,p[j]=beta[j];k<*ncv;k++)
        p[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        p[j]=p[j]*p[j];
        break;
      case 3: /* log transform */
        p[j]=exp(p[j]);
        break;
      case 4: /* logit transform */
        p[j]=exp(p[j]);
        p[j]=p[j]/(1+p[j]);
        break;
      case 5: /* cloglog transform */
        p[j]=1.0-exp(-exp(p[j]));
        break;
      case 6: /* square transform */
        p[j]=sqrt(p[j]);
        break;
      case 7: /* exp transform */
        p[j]=log(p[j]);
        break;
      case 8: /* loglog transform */
        p[j]=exp(-exp(p[j]));
        break;
      }
      p[j]=exp(p[j]);
      p[j]=p[j]/(1+p[j]);
      q[j]=1-p[j];
    }
    fit[i]=1.;
    for(j=y[i];j<(*nc-1);j++)
      fit[i]*=p[j];
    if(y[i]!=0)
      fit[i]*=q[y[i]-1];
    for(j=0;j<(*nc-1);j++) {
      cpred[*n*j+i]=1.;
      for(k=j;k<(*nc-1);k++)
        cpred[*n*j+i]*=p[k];
    }
    for(j=1;j<(*nc-1);j++)
      cpred[*n*j+i]*=q[j-1];
    pred[i]=0;
    for(j=1;j<(*nc-1);j++)
      if(cpred[*n*pred[i]+i]<cpred[*n*j+i])
        pred[i]=j;
    if(cpred[*n*pred[i]+i]<q[*nc-2])
      pred[i]=*nc-1;
    for(j=1;j<(*nc-1);j++)
      cpred[*n*j+i]+=cpred[*n*(j-1)+i];
  }
}

/* Linear downwards continuation-ratio regresion. */
void ldcrr(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],double coef[],int *n,int totnc[],double totncv[],int *iter,int *iterlim,double *steptol,
           double sd[],double fd[],double reg[],int *transform,double p[],double q[],int *rank,double aux[],int pivot[],double work[],double diff[],int *info,
           double diag[],double hessian[],double *loglik,double fit[],int pred[],double cpred[]) {
  int ny,ncncv,i,j,k,l;
  double tol;

  ny=1;
  tol=1.;
  ncncv=*nc-1+*ncv;
  for(i=0;i<ncncv;i++)
    beta[i]=coef[i];
  for(i=0;i<*n;i++)
    if(y[i]!=0)
      totnc[y[i]-1]+=freq[i];
  for(i=0;i<*ncv;i++)
    for(j=0;j<*n;j++)
      if(y[j]!=0)
        totncv[i]+=freq[j]*cov[*n*i+j]*y[j];
  while(*iter<*iterlim&&tol>*steptol) {
    for(i=0;i<ncncv;i++) {
      for(j=0;j<ncncv;j++)
        sd[ncncv*i+j]=0.; /* Initialize: I_{a_ia_i}, I_{a_ia_l}, I_{a_ib_l}, and I_{b_lb_m} */
      fd[i]=0.; /* Initialize: S_{a_i} and S_{b_l} */
    }
    for(i=0;i<(*nc-1);i++)
      for(j=i;j<(*nc-1);j++)
        fd[i]+=totnc[j]; /* First term of S_{a_i} */
    for(i=0;i<*ncv;i++)
      fd[*nc-1+i]=totncv[i]; /* First term of S_{b_l} */
    for(i=0;i<*n;i++) {
      for(j=0;j<(*nc-1);j++) {
        for(k=0,reg[j]=beta[j];k<*ncv;k++)
          reg[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
        switch(*transform) {
        case 2: /* sqrt transform */
          reg[j]=reg[j]*reg[j];
          break;
        case 3: /* log transform */
          reg[j]=exp(reg[j]);
          break;
        case 4: /* logit transform */
          reg[j]=exp(reg[j]);
          reg[j]=reg[j]/(1+reg[j]);
          break;
        case 5: /* cloglog transform */
          reg[j]=1.0-exp(-exp(reg[j]));
          break;
        case 6: /* square transform */
          reg[j]=sqrt(reg[j]);
          break;
        case 7: /* exp transform */
          reg[j]=log(reg[j]);
          break;
        case 8: /* loglog transform */
          reg[j]=exp(-exp(reg[j]));
          break;
        }
        p[j]=exp(reg[j]);
        p[j]=p[j]/(1+p[j]);
        q[j]=1-p[j];
      }
      for(j=0;j<y[i];j++) {
        fd[j]-=freq[i]*p[j]; /* Second term of S_{a_i} */
        sd[(ncncv+1)*j]-=freq[i]*p[j]*q[j]; /* First term of I_{a_ia_i} */
      }
      if(y[i]!=(*nc-1)) {
        fd[y[i]]-=freq[i]*p[y[i]]; /* Third term of S_{a_i} */
        sd[(ncncv+1)*y[i]]-=freq[i]*p[y[i]]*q[y[i]]; /* Second term of I_{a_ia_i} */
      }
      for(j=0;j<*ncv;j++) {
        for(k=0;k<y[i];k++) {
          fd[(*nc-1+j)]-=freq[i]*cov[*n*j+i]*p[k]; /* Second term of S_{b_l} */
          sd[ncncv*k+(*nc-1+j)]-=freq[i]*cov[*n*j+i]*p[k]*q[k]; /* First term of I_{a_ib_l} */
        }
        if(y[i]!=(*nc-1)) {
          fd[(*nc-1+j)]-=freq[i]*cov[*n*j+i]*p[y[i]]; /* Third term of S_{b_l} */
          sd[ncncv*y[i]+(*nc-1+j)]-=freq[i]*cov[*n*j+i]*p[y[i]]*q[y[i]]; /* Second term of I_{a_ib_l} */
        }
        for(k=j;k<*ncv;k++) {
          for(l=0;l<y[i];l++)
            sd[ncncv*(*nc-1+j)+(*nc-1+k)]-=freq[i]*cov[*n*j+i]*cov[*n*k+i]*p[l]*q[l]; /* First term of I_{b_lb_m} */
          if(y[i]!=(*nc-1))
            sd[ncncv*(*nc-1+j)+(*nc-1+k)]-=freq[i]*cov[*n*j+i]*cov[*n*k+i]*p[y[i]]*q[y[i]]; /* Second term of I_{b_lb_m} */
        }
      }
    }
    for(i=0;i<(ncncv-1);i++)
      for(j=(i+1);j<ncncv;j++)
        sd[ncncv*j+i]=sd[ncncv*i+j]; /* Complete symmetry of second derivative matrix */
    dqrdc2_(sd,&ncncv,&ncncv,&ncncv,steptol,rank,aux,pivot,work);
    dqrcf_(sd,&ncncv,rank,aux,fd,&ny,diff,info);
    for(i=0,tol=0.;i<ncncv;i++) {
      if(fabs(diff[i])>tol)
        tol=fabs(diff[i]);
      beta[i]-=diff[i];
    }
    (*iter)++;
    if(*info!=0)
      break;
    if(*rank!=ncncv)
      break;
  }
  dqrcf_(sd,&ncncv,rank,aux,diag,&ncncv,hessian,info);
  for(i=0;i<*n;i++) {
    for(j=0;j<(*nc-1);j++) {
      for(k=0,reg[j]=beta[j];k<*ncv;k++)
        reg[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        reg[j]=reg[j]*reg[j];
        break;
      case 3: /* log transform */
        reg[j]=exp(reg[j]);
        break;
      case 4: /* logit transform */
        reg[j]=exp(reg[j]);
        reg[j]=reg[j]/(1+reg[j]);
        break;
      case 5: /* cloglog transform */
        reg[j]=1.0-exp(-exp(reg[j]));
        break;
      case 6: /* square transform */
        reg[j]=sqrt(reg[j]);
        break;
      case 7: /* exp transform */
        reg[j]=log(reg[j]);
        break;
      case 8: /* loglog transform */
        reg[j]=exp(-exp(reg[j]));
        break;
      }
      p[j]=exp(reg[j]);
      p[j]=p[j]/(1+p[j]);
      q[j]=1-p[j];
    }
    fit[i]=1.;
    for(j=0;j<y[i];j++) {
      fit[i]*=p[j];
      *loglik+=freq[i]*log(p[j]);
    }
    if(y[i]!=(*nc-1)) {
      fit[i]*=q[y[i]];
      *loglik+=freq[i]*log(q[y[i]]);
    }
    for(j=0;j<(*nc-1);j++) {
      cpred[*n*j+i]=1.;
      for(k=0;k<j;k++)
        cpred[*n*j+i]*=p[k];
    }
    for(j=0;j<(*nc-1);j++)
      cpred[*n*j+i]*=q[j];
    pred[i]=0;
    for(j=1;j<(*nc-1);j++)
      if(cpred[*n*pred[i]+i]<cpred[*n*j+i])
        pred[i]=j;
    if(cpred[*n*pred[i]+i]<((cpred[*n*(*nc-2)+i]/q[*nc-2])*p[*nc-2]))
      pred[i]=*nc-1;
    for(j=1;j<(*nc-1);j++)
      cpred[*n*j+i]+=cpred[*n*(j-1)+i];
  }
}

/* Minus log-likelihood for downwards continuation-ratio regresion. */
void mlldcrr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double *loglik,double top[],double bot[]) {
  int i,j,k;

  for(i=0;i<*n;i++) {
    for(j=0;j<(*nc-1);j++) {
      for(k=0,top[j]=beta[j];k<*ncv;k++)
        top[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        top[j]=top[j]*top[j];
        break;
      case 3: /* log transform */
        top[j]=exp(top[j]);
        break;
      case 4: /* logit transform */
        top[j]=exp(top[j]);
        top[j]=top[j]/(1+top[j]);
        break;
      case 5: /* cloglog transform */
        top[j]=1.0-exp(-exp(top[j]));
        break;
      case 6: /* square transform */
        top[j]=sqrt(top[j]);
        break;
      case 7: /* exp transform */
        top[j]=log(top[j]);
        break;
      case 8: /* loglog transform */
        top[j]=exp(-exp(top[j]));
        break;
      }
      bot[j]=log(1+exp(top[j]));
    }
    for(j=0;j<y[i];j++)
      *loglik+=freq[i]*(top[j]-bot[j]);
    if(y[i]!=(*nc-1))
      *loglik-=freq[i]*bot[y[i]];
  }
}

/* Fitted values for downwards continuation-ratio regresion. */
void fvdcrr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double p[],double q[],double fit[],int pred[],
            double cpred[]) {
  int i,j,k;

  for(i=0;i<*n;i++) {
    for(j=0;j<(*nc-1);j++) {
      for(k=0,p[j]=beta[j];k<*ncv;k++)
        p[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        p[j]=p[j]*p[j];
        break;
      case 3: /* log transform */
        p[j]=exp(p[j]);
        break;
      case 4: /* logit transform */
        p[j]=exp(p[j]);
        p[j]=p[j]/(1+p[j]);
        break;
      case 5: /* cloglog transform */
        p[j]=1.0-exp(-exp(p[j]));
        break;
      case 6: /* square transform */
        p[j]=sqrt(p[j]);
        break;
      case 7: /* exp transform */
        p[j]=log(p[j]);
        break;
      case 8: /* loglog transform */
        p[j]=exp(-exp(p[j]));
        break;
      }
      p[j]=exp(p[j]);
      p[j]=p[j]/(1+p[j]);
      q[j]=1-p[j];
    }
    fit[i]=1.;
    for(j=0;j<y[i];j++)
      fit[i]*=p[j];
    if(y[i]!=(*nc-1))
      fit[i]*=q[y[i]];
    for(j=0;j<(*nc-1);j++) {
      cpred[*n*j+i]=1.;
      for(k=0;k<j;k++)
        cpred[*n*j+i]*=p[k];
    }
    for(j=0;j<(*nc-1);j++)
      cpred[*n*j+i]*=q[j];
    pred[i]=0;
    for(j=1;j<(*nc-1);j++)
      if(cpred[*n*pred[i]+i]<cpred[*n*j+i])
        pred[i]=j;
    if(cpred[*n*pred[i]+i]<((cpred[*n*(*nc-2)+i]/q[*nc-2])*p[*nc-2]))
      pred[i]=*nc-1;
    for(j=1;j<(*nc-1);j++)
      cpred[*n*j+i]+=cpred[*n*(j-1)+i];
  }
}

/* Linear proportional odds regresion. */
void lpor(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],double coef[],int *n,int totnc[],double totncv[],int *iter,int *iterlim,double *steptol,
          double sd[],double fd[],double reg[],int *transform,int *link,double p[],double lambda[],double q[],int *rank,double aux[],int pivot[],double work[],
          double diff[],int *info,double diag[],double hessian[],double *loglik,double fit[],int pred[],double cpred[]) {
  int ny,ncncv,i,j,k;
  double tol;

  ny=1;
  tol=1.;
  ncncv=*nc-1+*ncv;
  for(i=0;i<ncncv;i++)
    beta[i]=coef[i];
  for(i=0;i<*n;i++)
    if(y[i]!=(*nc-1))
      totnc[y[i]]+=freq[i];
  for(i=0;i<*ncv;i++)
    for(j=0;j<*n;j++)
      if(y[j]!=(*nc-1))
        totncv[i]+=freq[j]*cov[*n*i+j];
  while(*iter<*iterlim&&tol>*steptol) {
    for(i=0;i<ncncv;i++) {
      for(j=0;j<ncncv;j++)
        sd[ncncv*i+j]=0.; /* Initialize: I_{a_1a_1}, I_{a_ia_i}, I_{a_{nc-1}a_{nc-1}}, I_{a_1a_2}, I_{a_ia_{i+1}}, I_{a_{nc-2}a_{nc-1}}, */
                          /*             I_{a_sa_r}, I_{a_1b_l}, I_{a_ib_l}, I_{a_{nc-1}b_l}, and I_{b_lb_m} */
      fd[i]=0.; /* Initialize: S_{a_1}, S_{a_i}, S_{a_{nc-1}} and S_{b_l} */
    }
    fd[0]=totnc[0]; /* First term of S_{a_1} */
    for(i=0;i<*ncv;i++)
      fd[(*nc-1+i)]=totncv[i]; /* First term of S_{b_l} */
    for(i=0;i<*n;i++) {
      for(j=0;j<(*nc-1);j++) {
        for(k=0,reg[j]=beta[j];k<*ncv;k++)
          reg[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
        switch(*transform) {
        case 2: /* sqrt transform */
          reg[j]=reg[j]*reg[j];
          break;
        case 3: /* log transform */
          reg[j]=exp(reg[j]);
          break;
        case 4: /* logit transform */
          reg[j]=exp(reg[j]);
          reg[j]=reg[j]/(1+reg[j]);
          break;
        case 5: /* cloglog transform */
          reg[j]=1.0-exp(-exp(reg[j]));
          break;
        case 6: /* square transform */
          reg[j]=sqrt(reg[j]);
          break;
        case 7: /* exp transform */
          reg[j]=log(reg[j]);
          break;
        case 8: /* loglog transform */
          reg[j]=exp(-exp(reg[j]));
          break;
        }
        switch(*link) {
        case 1: /* logistic distribution */
          p[j]=exp(reg[j]);
          p[j]=p[j]/(1+p[j]);
          break;
        case 2: /* normal distribution */
          p[j]=pnorm(reg[j],0.0,1.0,1,0);
          break;
        case 3: /* inverse log-log */
          p[j]=exp(-exp(reg[j]));
          break;
        case 4: /* inverse complementary log-log */
          p[j]=1.0-exp(-exp(reg[j]));
          break;
        case 5: /* uniform distribution */
          p[j]=punif(reg[j],0.0,1.0,1,0);
          break;
        case 6: /* log-normal distribution */
          p[j]=plnorm(reg[j],0.0,1.0,1,0);
          break;
        case 7: /* exponential distribution */
          p[j]=pexp(reg[j],1.0,1,0);
          break;
        case 8: /* Pareto distribution */
          p[j]=ppareto(reg[j],1.0,2.0);
          break;
        case 9: /* Cauchy distribution */
          p[j]=pcauchy(reg[j],0.0,1.0,1,0);
          break;
        case 10: /* Laplace distribution */
          p[j]=plaplace(reg[j],0.0,1.0);
          break;
        case 11: /* Levy distribution */
          p[j]=plevy(reg[j],0.0,1.0);
          break;
        case 12: /* simplex distribution */
          p[j]=psimplex(reg[j],0.5,1.0);
          break;
        case 13: /* gamma distribution */
          p[j]=pgamma(reg[j],lambda[0],1.0/lambda[0],1,0);
          break;
        case 14: /* Weibull distribution */
          p[j]=pweibull(reg[j],lambda[0],1.0,1,0);
          break;
        case 15: /* inverse Gauss distribution */
          p[j]=pinvgauss(reg[j],1.0,lambda[0]);
          break;
        case 16: /* t distribution */
          p[j]=pt(reg[j],lambda[0],1,0);
          break;
        case 17: /* chi-square distribution */
          p[j]=pchisq(reg[j],lambda[0],1,0);
          break;
        case 18: /* generalized logistic distribution */
          p[j]=pglogis(reg[j],0.0,1.0,lambda[0]);
          break;
        case 19: /* generalized extreme value distribution */
          p[j]=pgextval(reg[j],1.0,1.0,lambda[0]);
          break;
        case 20: /* Box-Cox distribution */
          p[j]=pboxcox(reg[j],1.0,1.0,lambda[0]);
          break;
        case 21: /* power exponential distribution */
          p[j]=ppowexp(reg[j],0.0,1.0,lambda[0]);
          break;
        case 22: /* Burr distribution */
          p[j]=pburr(reg[j],1.0,1.0,lambda[0]);
          break;
        case 23: /* Hjorth distribution */
          p[j]=phjorth(reg[j],1.0,1.0,lambda[0]);
          break;
        case 24: /* beta distribution */
          p[j]=pbeta(reg[j],lambda[0],lambda[1],1,0);
          break;
        case 25: /* stable distribution */
          p[j]=pstable(reg[j],0.0,1.0,lambda[0],lambda[1]);
          break;
        case 26: /* generalized gamma distribution */
          p[j]=pggamma(reg[j],lambda[0],1.0,lambda[1]);
          break;
        case 27: /* generalized Weibull distribution (Mudholkar et al, 1995) */
          p[j]=pgweibull(reg[j],lambda[0],1.0,lambda[1]);
          break;
        case 28: /* generalized inverse Gauss distribution */
          p[j]=pginvgauss(reg[j],1.0,lambda[0],lambda[1]);
          break;
        case 29: /* F distribution */
          p[j]=pf(reg[j],lambda[0],lambda[1],1,0);
          break;
        case 30: /* non-central t distribution */
          p[j]=pnt(reg[j],lambda[0],lambda[1],1,0);
          break;
        case 31: /* non-central chi-square distribution */
          p[j]=pnchisq(reg[j],lambda[0],lambda[1],1,0);
          break;
        case 32: /* studentized range distribution */
          p[j]=ptukey(reg[j],lambda[2],lambda[1],lambda[0],1,0);
          break;
        case 33: /* non-central beta distribution */
          p[j]=pnbeta(reg[j],lambda[0],lambda[1],lambda[2],1,0);
          break;
        case 34: /* non-central F distribution */
          p[j]=pnf(reg[j],lambda[0],lambda[1],lambda[2],1,0);
          break;
        }
        q[j]=1-p[j];
      }
      if(y[i]==0) {
        fd[0]-=freq[i]*p[0]; /* Second term of S_{a_1} */
        sd[0]-=freq[i]*p[0]*q[0]; /* First term of I_{a_1a_1} */
      }
      else {
        if(y[i]==(*nc-2)&&*nc>2) {
          fd[*nc-2]+=freq[i]*(p[*nc-2]*q[*nc-2])/(p[*nc-2]-p[*nc-3]); /* First term of S_{a_{nc-1}} */
          sd[(ncncv+1)*(*nc-2)]-=freq[i]*(p[*nc-2]*q[*nc-2]+((p[*nc-2]*q[*nc-2]*p[*nc-3]*q[*nc-3])/pow((p[*nc-2]-p[*nc-3]),2))); /* First term of I_{a_{nc-1}a_{nc-1}} */
        }
        else
          if(y[i]==(*nc-1)) {
            fd[*nc-2]-=freq[i]*p[*nc-2]; /* Second term of S_{a_{nc-1}} */
            sd[(ncncv+1)*(*nc-2)]-=freq[i]*p[*nc-2]*q[*nc-2]; /* Second term of I_{a_{nc-1}a_{nc-1}} */
          }
        if(*nc>3) {
          if(y[i]<(*nc-2)) {
            fd[y[i]]+=freq[i]*(p[y[i]]*q[y[i]])/(p[y[i]]-p[y[i]-1]); /* First term of S_{a_i} */
            sd[(ncncv+1)*y[i]]-=freq[i]*(p[y[i]]*q[y[i]]+((p[y[i]]*q[y[i]]*p[y[i]-1]*q[y[i]-1])/pow((p[y[i]]-p[y[i]-1]),2))); /* First term of I_{a_ia_i} */
          }
          if(y[i]!=1&&y[i]!=(*nc-1)) {
            fd[y[i]-1]-=freq[i]*(p[y[i]-1]*q[y[i]-1])/(p[y[i]]-p[y[i]-1]); /* Second term of S_{a_i} */
            sd[(ncncv+1)*(y[i]-1)]-=freq[i]*(p[y[i]-1]*q[y[i]-1]+((p[y[i]-1]*q[y[i]-1]*p[y[i]]*q[y[i]])/pow((p[y[i]]-p[y[i]-1]),2))); /* Second term of I_{a_ia_i} */
          }
          if(y[i]!=1&&y[i]<(*nc-2))
            sd[(ncncv+1)*(y[i]-1)+1]-=freq[i]*(p[y[i]-1]*q[y[i]-1]*p[y[i]]*q[y[i]])/pow((p[y[i]]-p[y[i]-1]),2); /* I_{a_ia_{i+1}} */
          if(y[i]==(*nc-2))
            sd[(ncncv+1)*(*nc-3)+1]+=freq[i]*(p[*nc-2]*q[*nc-2]*p[*nc-3]*q[*nc-3])/pow((p[*nc-2]-p[*nc-3]),2); /* I_{a_{nc-2}a_{nc-1}} */
        }
        if(y[i]==1&&*nc>2) {
          fd[0]-=freq[i]*(p[0]*q[0])/(p[1]-p[0]); /* Third term of S_{a_1} */
          sd[0]-=freq[i]*(p[0]*q[0]+(p[0]*q[0]*p[1]*q[1])/pow((p[1]-p[0]),2)); /* Second term of I_{a_1a_1} */
          sd[1]+=freq[i]*(p[0]*q[0]*p[1]*q[1])/pow((p[1]-p[0]),2); /* I_{a_1a_2} */
        }
      }
      for(j=0;j<*ncv;j++) {
        if(y[i]==0) {
          fd[(*nc-1+j)]-=freq[i]*cov[*n*j+i]*p[0]; /* Second term of S_{b_l} */
          sd[(*nc-1+j)]-=freq[i]*cov[*n*j+i]*p[0]*q[0]; /* First term of I_{a_1b_l} */
          sd[(ncncv+1)*(*nc-1+j)]-=freq[i]*cov[*n*j+i]*cov[*n*j+i]*p[0]*q[0]; /* First term of I_{b_lb_l} */
        }
        else {
          if(y[i]!=(*nc-1)) {
            fd[(*nc-1+j)]-=freq[i]*cov[*n*j+i]*(p[y[i]]+p[y[i]-1]); /* Third and fourth term of S_{b_l} */
            sd[(ncncv+1)*(*nc-1+j)]-=freq[i]*cov[*n*j+i]*cov[*n*j+i]*(p[y[i]]*q[y[i]]+p[y[i]-1]*q[y[i]-1]); /* Second term of I_{b_lb_l} */
          }
          else {
            fd[(*nc-1+j)]-=freq[i]*cov[*n*j+i]*p[*nc-2]; /* Fifth term of S_{b_l} */
            sd[ncncv*(*nc-2)+(*nc-1)+j]-=freq[i]*cov[*n*j+i]*p[*nc-2]*q[*nc-2]; /* Second term of I_{a_{nc-1}b_l} */
            sd[(ncncv+1)*(*nc-1+j)]-=freq[i]*cov[*n*j+i]*cov[*n*j+i]*p[y[i]-1]*q[y[i]-1]; /* Third term of I_{b_lb_l} */
          }
          if(y[i]==1&&*nc>2)
            sd[(*nc-1+j)]-=freq[i]*cov[*n*j+i]*p[0]*q[0]; /* Second term of I_{a_1b_l} */
          if(y[i]==(*nc-2))
            sd[ncncv*(*nc-2)+(*nc-1)+j]-=freq[i]*cov[*n*j+i]*p[*nc-2]*q[*nc-2]; /* First term of I_{a_{nc-1}b_l} */
          if(y[i]<(*nc-2))
            sd[ncncv*y[i]+(*nc-1)+j]-=freq[i]*cov[*n*j+i]*p[y[i]]*q[y[i]]; /* First term of I_{a_ib_l} */
          else
            if(y[i]!=1&&y[i]<(*nc-1))
              sd[ncncv*(y[i]-1)+(*nc-1)+j]-=freq[i]*cov[*n*j+i]*p[y[i]-1]*q[y[i]-1]; /* Second term of I_{a_ib_l} */
        }
        for(k=(j+1);k<*ncv;k++) {
          if(y[i]==0)
            sd[ncncv*(*nc-1+j)+(*nc-1)+k]-=freq[i]*cov[*n*j+i]*cov[*n*k+i]*p[0]*q[0]; /* First term of I_{b_lb_m} */
          else
            if(y[i]!=(*nc-1))
              sd[ncncv*(*nc-1+j)+(*nc-1)+k]-=freq[i]*cov[*n*j+i]*cov[*n*k+i]*(p[y[i]]*q[y[i]]+p[y[i]-1]*q[y[i]-1]); /* Second and third term of I_{b_lb_m} */
            else
              sd[ncncv*(*nc-1+j)+(*nc-1)+k]-=freq[i]*cov[*n*j+i]*cov[*n*k+i]*p[y[i]-1]*q[y[i]-1]; /* Fourth term of I_{b_lb_m} */
        }
      }
    }
    for(i=0;i<(ncncv-1);i++)
      for(j=(i+1);j<ncncv;j++)
        sd[ncncv*j+i]=sd[ncncv*i+j]; /* Complete symmetry of second derivative matrix */
    dqrdc2_(sd,&ncncv,&ncncv,&ncncv,steptol,rank,aux,pivot,work);
    dqrcf_(sd,&ncncv,rank,aux,fd,&ny,diff,info);
    for(i=0,tol=0.;i<ncncv;i++) {
      if(fabs(diff[i])>tol)
        tol=fabs(diff[i]);
      beta[i]-=diff[i];
    }
    (*iter)++;
    if(*info!=0)
      break;
    if(*rank!=ncncv)
      break;
  }
  dqrcf_(sd,&ncncv,rank,aux,diag,&ncncv,hessian,info);
  for(i=0;i<*n;i++) {
    for(j=0;j<(*nc-1);j++) {
      for(k=0,reg[j]=beta[j];k<*ncv;k++)
        reg[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        reg[j]=reg[j]*reg[j];
        break;
      case 3: /* log transform */
        reg[j]=exp(reg[j]);
        break;
      case 4: /* logit transform */
        reg[j]=exp(reg[j]);
        reg[j]=reg[j]/(1+reg[j]);
        break;
      case 5: /* cloglog transform */
        reg[j]=1.0-exp(-exp(reg[j]));
        break;
      case 6: /* square transform */
        reg[j]=sqrt(reg[j]);
        break;
      case 7: /* exp transform */
        reg[j]=log(reg[j]);
        break;
      case 8: /* loglog transform */
        reg[j]=exp(-exp(reg[j]));
        break;
      }
      switch(*link) {
      case 1: /* logistic distribution */
        p[j]=exp(reg[j]);
        p[j]=p[j]/(1+p[j]);
        break;
      case 2: /* normal distribution */
        p[j]=pnorm(reg[j],0.0,1.0,1,0);
        break;
      case 3: /* inverse log-log */
        p[j]=exp(-exp(reg[j]));
        break;
      case 4: /* inverse complementary log-log */
        p[j]=1.0-exp(-exp(reg[j]));
        break;
      case 5: /* uniform distribution */
        p[j]=punif(reg[j],0.0,1.0,1,0);
        break;
      case 6: /* log-normal distribution */
        p[j]=plnorm(reg[j],0.0,1.0,1,0);
        break;
      case 7: /* exponential distribution */
        p[j]=pexp(reg[j],1.0,1,0);
        break;
      case 8: /* Pareto distribution */
        p[j]=ppareto(reg[j],1.0,2.0);
        break;
      case 9: /* Cauchy distribution */
        p[j]=pcauchy(reg[j],0.0,1.0,1,0);
        break;
      case 10: /* Laplace distribution */
        p[j]=plaplace(reg[j],0.0,1.0);
        break;
      case 11: /* Levy distribution */
        p[j]=plevy(reg[j],0.0,1.0);
        break;
      case 12: /* simplex distribution */
        p[j]=psimplex(reg[j],0.5,1.0);
        break;
      case 13: /* gamma distribution */
        p[j]=pgamma(reg[j],lambda[0],1.0/lambda[0],1,0);
        break;
      case 14: /* Weibull distribution */
        p[j]=pweibull(reg[j],lambda[0],1.0,1,0);
        break;
      case 15: /* inverse Gauss distribution */
        p[j]=pinvgauss(reg[j],1.0,lambda[0]);
        break;
      case 16: /* t distribution */
        p[j]=pt(reg[j],lambda[0],1,0);
        break;
      case 17: /* chi-square distribution */
        p[j]=pchisq(reg[j],lambda[0],1,0);
        break;
      case 18: /* generalized logistic distribution */
        p[j]=pglogis(reg[j],0.0,1.0,lambda[0]);
        break;
      case 19: /* generalized extreme value distribution */
        p[j]=pgextval(reg[j],1.0,1.0,lambda[0]);
        break;
      case 20: /* Box-Cox distribution */
        p[j]=pboxcox(reg[j],1.0,1.0,lambda[0]);
        break;
      case 21: /* power exponential distribution */
        p[j]=ppowexp(reg[j],0.0,1.0,lambda[0]);
        break;
      case 22: /* Burr distribution */
        p[j]=pburr(reg[j],1.0,1.0,lambda[0]);
        break;
      case 23: /* Hjorth distribution */
        p[j]=phjorth(reg[j],1.0,1.0,lambda[0]);
        break;
      case 24: /* beta distribution */
        p[j]=pbeta(reg[j],lambda[0],lambda[1],1,0);
        break;
      case 25: /* stable distribution */
        p[j]=pstable(reg[j],0.0,1.0,lambda[0],lambda[1]);
        break;
      case 26: /* generalized gamma distribution */
        p[j]=pggamma(reg[j],lambda[0],1.0,lambda[1]);
        break;
      case 27: /* generalized Weibull distribution (Mudholkar et al, 1995) */
        p[j]=pgweibull(reg[j],lambda[0],1.0,lambda[1]);
        break;
      case 28: /* generalized inverse Gauss distribution */
        p[j]=pginvgauss(reg[j],1.0,lambda[0],lambda[1]);
        break;
      case 29: /* F distribution */
        p[j]=pf(reg[j],lambda[0],lambda[1],1,0);
        break;
      case 30: /* non-central t distribution */
        p[j]=pnt(reg[j],lambda[0],lambda[1],1,0);
        break;
      case 31: /* non-central chi-square distribution */
        p[j]=pnchisq(reg[j],lambda[0],lambda[1],1,0);
        break;
      case 32: /* studentized range distribution */
        p[j]=ptukey(reg[j],lambda[2],lambda[1],lambda[0],1,0);
        break;
      case 33: /* non-central beta distribution */
        p[j]=pnbeta(reg[j],lambda[0],lambda[1],lambda[2],1,0);
        break;
      case 34: /* non-central F distribution */
        p[j]=pnf(reg[j],lambda[0],lambda[1],lambda[2],1,0);
        break;
      }
      q[j]=1-p[j];
    }
    if(y[i]==0) {
      fit[i]=p[0];
      *loglik+=freq[i]*log(p[0]);
    }
    else
      if(y[i]==(*nc-1)) {
        fit[i]=q[*nc-2];
        *loglik+=freq[i]*log(q[*nc-2]);
      }
      else {
        fit[i]=p[y[i]]-p[y[i]-1];
        *loglik+=freq[i]*log(p[y[i]]-p[y[i]-1]);
      }
    pred[i]=1;
    for(j=2;j<(*nc-1);j++)
      if((p[pred[i]]-p[pred[i]-1])<(p[j]-p[j-1]))
        pred[i]=j;
    if((p[pred[i]]-p[pred[i]-1])<p[0]) {
      pred[i]=0;
      if(p[0]<q[*nc-2])
        pred[i]=*nc-1;
    }
    else
      if((p[pred[i]]-p[pred[i]-1])<q[*nc-2])
        pred[i]=*nc-1;
    for(j=0;j<(*nc-1);j++)
      cpred[*n*j+i]=p[j];
  }
}

/* Minus log-likelihood for proportional odds regresion. */
void mllpor(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,int *link,double lambda[],double *loglik,double p[],double q[]) {
  int i,j,k;

  for(i=0;i<*n;i++) {
    for(j=0;j<(*nc-1);j++) {
      for(k=0,p[j]=beta[j];k<*ncv;k++)
        p[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        p[j]=p[j]*p[j];
        break;
      case 3: /* log transform */
        p[j]=exp(p[j]);
        break;
      case 4: /* logit transform */
        p[j]=exp(p[j]);
        p[j]=p[j]/(1+p[j]);
        break;
      case 5: /* cloglog transform */
        p[j]=1.0-exp(-exp(p[j]));
        break;
      case 6: /* square transform */
        p[j]=sqrt(p[j]);
        break;
      case 7: /* exp transform */
        p[j]=log(p[j]);
        break;
      case 8: /* loglog transform */
        p[j]=exp(-exp(p[j]));
        break;
      }
      switch(*link) {
      case 1: /* logistic distribution */
        p[j]=exp(p[j]);
        p[j]=p[j]/(1+p[j]);
        break;
      case 2: /* normal distribution */
        p[j]=pnorm(p[j],0.0,1.0,1,0);
        break;
      case 3: /* inverse log-log */
        p[j]=exp(-exp(p[j]));
        break;
      case 4: /* inverse complementary log-log */
        p[j]=1.0-exp(-exp(p[j]));
        break;
      case 5: /* uniform distribution */
        p[j]=punif(p[j],0.0,1.0,1,0);
        break;
      case 6: /* log-normal distribution */
        p[j]=plnorm(p[j],0.0,1.0,1,0);
        break;
      case 7: /* exponential distribution */
        p[j]=pexp(p[j],1.0,1,0);
        break;
      case 8: /* Pareto distribution */
        p[j]=ppareto(p[j],1.0,2.0);
        break;
      case 9: /* Cauchy distribution */
        p[j]=pcauchy(p[j],0.0,1.0,1,0);
        break;
      case 10: /* Laplace distribution */
        p[j]=plaplace(p[j],0.0,1.0);
        break;
      case 11: /* Levy distribution */
        p[j]=plevy(p[j],0.0,1.0);
        break;
      case 12: /* simplex distribution */
        p[j]=psimplex(p[j],0.5,1.0);
        break;
      case 13: /* gamma distribution */
        p[j]=pgamma(p[j],lambda[0],1.0/lambda[0],1,0);
        break;
      case 14: /* Weibull distribution */
        p[j]=pweibull(p[j],lambda[0],1.0,1,0);
        break;
      case 15: /* inverse Gauss distribution */
        p[j]=pinvgauss(p[j],1.0,lambda[0]);
        break;
      case 16: /* t distribution */
        p[j]=pt(p[j],lambda[0],1,0);
        break;
      case 17: /* chi-square distribution */
        p[j]=pchisq(p[j],lambda[0],1,0);
        break;
      case 18: /* generalized logistic distribution */
        p[j]=pglogis(p[j],0.0,1.0,lambda[0]);
        break;
      case 19: /* generalized extreme value distribution */
        p[j]=pgextval(p[j],1.0,1.0,lambda[0]);
        break;
      case 20: /* Box-Cox distribution */
        p[j]=pboxcox(p[j],1.0,1.0,lambda[0]);
        break;
      case 21: /* power exponential distribution */
        p[j]=ppowexp(p[j],0.0,1.0,lambda[0]);
        break;
      case 22: /* Burr distribution */
        p[j]=pburr(p[j],1.0,1.0,lambda[0]);
        break;
      case 23: /* Hjorth distribution */
        p[j]=phjorth(p[j],1.0,1.0,lambda[0]);
        break;
      case 24: /* beta distribution */
        p[j]=pbeta(p[j],lambda[0],lambda[1],1,0);
        break;
      case 25: /* stable distribution */
        p[j]=pstable(p[j],0.0,1.0,lambda[0],lambda[1]);
        break;
      case 26: /* generalized gamma distribution */
        p[j]=pggamma(p[j],lambda[0],1.0,lambda[1]);
        break;
      case 27: /* generalized Weibull distribution (Mudholkar et al, 1995) */
        p[j]=pgweibull(p[j],lambda[0],1.0,lambda[1]);
        break;
      case 28: /* generalized inverse Gauss distribution */
        p[j]=pginvgauss(p[j],1.0,lambda[0],lambda[1]);
        break;
      case 29: /* F distribution */
        p[j]=pf(p[j],lambda[0],lambda[1],1,0);
        break;
      case 30: /* non-central t distribution */
        p[j]=pnt(p[j],lambda[0],lambda[1],1,0);
        break;
      case 31: /* non-central chi-square distribution */
        p[j]=pnchisq(p[j],lambda[0],lambda[1],1,0);
        break;
      case 32: /* studentized range distribution */
        p[j]=ptukey(p[j],lambda[2],lambda[1],lambda[0],1,0);
        break;
      case 33: /* non-central beta distribution */
        p[j]=pnbeta(p[j],lambda[0],lambda[1],lambda[2],1,0);
        break;
      case 34: /* non-central F distribution */
        p[j]=pnf(p[j],lambda[0],lambda[1],lambda[2],1,0);
        break;
      }
      q[j]=1-p[j];
    }
    if(y[i]==0)
      *loglik+=freq[i]*log(p[0]);
    else
      if(y[i]==(*nc-1))
        *loglik+=freq[i]*log(q[*nc-2]);
      else
        *loglik+=freq[i]*log(p[y[i]]-p[y[i]-1]);
  }
}

/* Fitted values for proportional odds regresion. */
void fvpor(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,int *link,double lambda[],double p[],double q[],double fit[],
           int pred[],double cpred[]) {
  int i,j,k;

  for(i=0;i<*n;i++) {
    for(j=0;j<(*nc-1);j++) {
      for(k=0,p[j]=beta[j];k<*ncv;k++)
        p[j]+=beta[(*nc-1)+k]*cov[*n*k+i];
      switch(*transform) {
      case 2: /* sqrt transform */
        p[j]=p[j]*p[j];
        break;
      case 3: /* log transform */
        p[j]=exp(p[j]);
        break;
      case 4: /* logit transform */
        p[j]=exp(p[j]);
        p[j]=p[j]/(1+p[j]);
        break;
      case 5: /* cloglog transform */
        p[j]=1.0-exp(-exp(p[j]));
        break;
      case 6: /* square transform */
        p[j]=sqrt(p[j]);
        break;
      case 7: /* exp transform */
        p[j]=log(p[j]);
        break;
      case 8: /* loglog transform */
        p[j]=exp(-exp(p[j]));
        break;
      }
      switch(*link) {
      case 1: /* logistic distribution */
        p[j]=exp(p[j]);
        p[j]=p[j]/(1+p[j]);
        break;
      case 2: /* normal distribution */
        p[j]=pnorm(p[j],0.0,1.0,1,0);
        break;
      case 3: /* inverse log-log */
        p[j]=exp(-exp(p[j]));
        break;
      case 4: /* inverse complementary log-log */
        p[j]=1.0-exp(-exp(p[j]));
        break;
      case 5: /* uniform distribution */
        p[j]=punif(p[j],0.0,1.0,1,0);
        break;
      case 6: /* log-normal distribution */
        p[j]=plnorm(p[j],0.0,1.0,1,0);
        break;
      case 7: /* exponential distribution */
        p[j]=pexp(p[j],1.0,1,0);
        break;
      case 8: /* Pareto distribution */
        p[j]=ppareto(p[j],1.0,2.0);
        break;
      case 9: /* Cauchy distribution */
        p[j]=pcauchy(p[j],0.0,1.0,1,0);
        break;
      case 10: /* Laplace distribution */
        p[j]=plaplace(p[j],0.0,1.0);
        break;
      case 11: /* Levy distribution */
        p[j]=plevy(p[j],0.0,1.0);
        break;
      case 12: /* simplex distribution */
        p[j]=psimplex(p[j],0.5,1.0);
        break;
      case 13: /* gamma distribution */
        p[j]=pgamma(p[j],lambda[0],1.0/lambda[0],1,0);
        break;
      case 14: /* Weibull distribution */
        p[j]=pweibull(p[j],lambda[0],1.0,1,0);
        break;
      case 15: /* inverse Gauss distribution */
        p[j]=pinvgauss(p[j],1.0,lambda[0]);
        break;
      case 16: /* t distribution */
        p[j]=pt(p[j],lambda[0],1,0);
        break;
      case 17: /* chi-square distribution */
        p[j]=pchisq(p[j],lambda[0],1,0);
        break;
      case 18: /* generalized logistic distribution */
        p[j]=pglogis(p[j],0.0,1.0,lambda[0]);
        break;
      case 19: /* generalized extreme value distribution */
        p[j]=pgextval(p[j],1.0,1.0,lambda[0]);
        break;
      case 20: /* Box-Cox distribution */
        p[j]=pboxcox(p[j],1.0,1.0,lambda[0]);
        break;
      case 21: /* power exponential distribution */
        p[j]=ppowexp(p[j],0.0,1.0,lambda[0]);
        break;
      case 22: /* Burr distribution */
        p[j]=pburr(p[j],1.0,1.0,lambda[0]);
        break;
      case 23: /* Hjorth distribution */
        p[j]=phjorth(p[j],1.0,1.0,lambda[0]);
        break;
      case 24: /* beta distribution */
        p[j]=pbeta(p[j],lambda[0],lambda[1],1,0);
        break;
      case 25: /* stable distribution */
        p[j]=pstable(p[j],0.0,1.0,lambda[0],lambda[1]);
        break;
      case 26: /* generalized gamma distribution */
        p[j]=pggamma(p[j],lambda[0],1.0,lambda[1]);
        break;
      case 27: /* generalized Weibull distribution (Mudholkar et al, 1995) */
        p[j]=pgweibull(p[j],lambda[0],1.0,lambda[1]);
        break;
      case 28: /* generalized inverse Gauss distribution */
        p[j]=pginvgauss(p[j],1.0,lambda[0],lambda[1]);
        break;
      case 29: /* F distribution */
        p[j]=pf(p[j],lambda[0],lambda[1],1,0);
        break;
      case 30: /* non-central t distribution */
        p[j]=pnt(p[j],lambda[0],lambda[1],1,0);
        break;
      case 31: /* non-central chi-square distribution */
        p[j]=pnchisq(p[j],lambda[0],lambda[1],1,0);
        break;
      case 32: /* studentized range distribution */
        p[j]=ptukey(p[j],lambda[2],lambda[1],lambda[0],1,0);
        break;
      case 33: /* non-central beta distribution */
        p[j]=pnbeta(p[j],lambda[0],lambda[1],lambda[2],1,0);
        break;
      case 34: /* non-central F distribution */
        p[j]=pnf(p[j],lambda[0],lambda[1],lambda[2],1,0);
        break;
      }
      q[j]=1-p[j];
    }
    if(y[i]==0)
      fit[i]=p[0];
    else
      if(y[i]==(*nc-1))
        fit[i]=q[*nc-2];
      else
        fit[i]=p[y[i]]-p[y[i]-1];
    pred[i]=1;
    for(j=2;j<(*nc-1);j++)
      if((p[pred[i]]-p[pred[i]-1])<(p[j]-p[j-1]))
        pred[i]=j;
    if((p[pred[i]]-p[pred[i]-1])<p[0]) {
      pred[i]=0;
      if(p[0]<q[*nc-2])
        pred[i]=*nc-1;
    }
    else
      if((p[pred[i]]-p[pred[i]-1])<q[*nc-2])
        pred[i]=*nc-1;
    for(j=0;j<(*nc-1);j++)
      cpred[*n*j+i]=p[j];
  }
}

/* Linear adjacent categories regresion. */
void lacr(int y[],int freq[],double cov[],int *nc,int *ncv,double beta[],double coef[],int *n,int totnc[],double totncv[],int *iter,int *iterlim,double *steptol,
          double sd[],double fd[],double reg[],int *transform,double p[],int *rank,double aux[],int pivot[],double work[],double diff[],int *info,double diag[],
          double hessian[],double *loglik,double fit[],int pred[],double cpred[]) {
  int ny,ncncv,i,j,k,l;
  double tol,bot,q,tmp,tmp2,tmp3;

  ny=1;
  tol=1.;
  ncncv=*nc-1+*ncv;
  for(i=0;i<ncncv;i++)
    beta[i]=coef[i];
  for(i=0;i<*n;i++)
    if(y[i]!=0)
      totnc[y[i]-1]+=freq[i];
  for(i=0;i<*ncv;i++)
    for(j=0;j<*n;j++)
      if(y[j]!=0)
        totncv[i]-=freq[j]*cov[*n*i+j]*y[j];
  while(*iter<*iterlim&&tol>*steptol) {
    for(i=0;i<ncncv;i++) {
      for(j=0;j<ncncv;j++)
        sd[ncncv*i+j]=0.; /* Initialize: I_{a_ia_l}, I_{a_ib_l}, and I_{b_lb_m} */
      fd[i]=0.; /* Initialize: S_{a_i} and S_{b_l} */
    }
    for(i=0;i<(*nc-1);i++)
      for(j=i;j<(*nc-1);j++)
        fd[i]-=totnc[j]; /* First term of S_{a_i} */
    for(i=0;i<*ncv;i++)
      fd[*nc-1+i]=totncv[i]; /* First term of S_{b_l} */
    for(i=0;i<*n;i++) {
      for(j=0;j<(*nc-1);j++)
        for(k=0,reg[j]=-beta[j];k<*ncv;k++)
          reg[j]-=beta[(*nc-1)+k]*cov[*n*k+i];
      for(j=1;j<(*nc-1);j++)
        reg[j]+=reg[j-1];
      for(j=0;j<(*nc-1);j++)
        switch(*transform) {
        case 2: /* sqrt transform */
          reg[j]=reg[j]*reg[j];
          break;
        case 3: /* log transform */
          reg[j]=exp(reg[j]);
          break;
        case 4: /* logit transform */
          reg[j]=exp(reg[j]);
          reg[j]=reg[j]/(1+reg[j]);
          break;
        case 5: /* cloglog transform */
          reg[j]=1.0-exp(-exp(reg[j]));
          break;
        case 6: /* square transform */
          reg[j]=sqrt(reg[j]);
          break;
        case 7: /* exp transform */
          reg[j]=log(reg[j]);
          break;
        case 8: /* loglog transform */
          reg[j]=exp(-exp(reg[j]));
          break;
        }
      for(j=0,bot=1;j<(*nc-1);j++) {
        reg[j]=exp(reg[j]);
        bot+=reg[j];
      }
      for(j=0;j<(*nc-1);j++)
        p[j]=reg[j]/bot;
      q=1/bot;
      for(j=0,tmp=0;j<(*nc-1);j++) {
        for(k=j,tmp2=0;k<(*nc-1);k++)
          tmp2+=p[k];
        fd[j]+=freq[i]*tmp2; /* Second term of S_{a_i} */
        for(k=j;k<(*nc-1);k++) {
          for(l=k,tmp3=0;l<(*nc-1);l++)
            tmp3+=p[l];
          sd[ncncv*j+k]-=freq[i]*tmp3*(1-tmp2); /* I_{a_ia_l} */
        }
        tmp+=(j+1)*p[j];
      }
      for(j=0;j<*ncv;j++) {
        fd[(*nc-1+j)]+=freq[i]*cov[*n*j+i]*tmp; /* Second term of S_{b_l} */
        for(k=0;k<(*nc-1);k++) {
          for(l=k,tmp2=0;l<(*nc-1);l++)
            tmp2+=p[l]*(l+1-tmp);
          sd[ncncv*k+(*nc-1+j)]-=freq[i]*cov[*n*j+i]*tmp2; /* I_{a_ib_l} */
        }
        for(k=j;k<*ncv;k++) {
          for(l=0,tmp2=0;l<(*nc-1);l++)
            tmp2+=p[l]*(l+1)*((l+1)-tmp);
          sd[ncncv*(*nc-1+j)+(*nc-1+k)]-=freq[i]*cov[*n*j+i]*cov[*n*k+i]*tmp2; /* I_{b_lb_m} */
        }
      }
    }
    for(i=0;i<(ncncv-1);i++)
      for(j=(i+1);j<ncncv;j++)
        sd[ncncv*j+i]=sd[ncncv*i+j]; /* Complete symmetry of second derivative matrix */
    dqrdc2_(sd,&ncncv,&ncncv,&ncncv,steptol,rank,aux,pivot,work);
    dqrcf_(sd,&ncncv,rank,aux,fd,&ny,diff,info);
    for(i=0,tol=0.;i<ncncv;i++) {
      if(fabs(diff[i])>tol)
        tol=fabs(diff[i]);
      beta[i]-=diff[i];
    }
    (*iter)++;
    if(*info!=0)
      break;
    if(*rank!=ncncv)
      break;
  }
  dqrcf_(sd,&ncncv,rank,aux,diag,&ncncv,hessian,info);
  for(i=0;i<*n;i++) {
    for(j=0;j<(*nc-1);j++)
      for(k=0,reg[j]=-beta[j];k<*ncv;k++)
        reg[j]-=beta[(*nc-1)+k]*cov[*n*k+i];
    for(j=1;j<(*nc-1);j++)
      reg[j]+=reg[j-1];
    for(j=0;j<(*nc-1);j++)
      switch(*transform) {
      case 2: /* sqrt transform */
        reg[j]=reg[j]*reg[j];
        break;
      case 3: /* log transform */
        reg[j]=exp(reg[j]);
        break;
      case 4: /* logit transform */
        reg[j]=exp(reg[j]);
        reg[j]=reg[j]/(1+reg[j]);
        break;
      case 5: /* cloglog transform */
        reg[j]=1.0-exp(-exp(reg[j]));
        break;
      case 6: /* square transform */
        reg[j]=sqrt(reg[j]);
        break;
      case 7: /* exp transform */
        reg[j]=log(reg[j]);
        break;
      case 8: /* loglog transform */
        reg[j]=exp(-exp(reg[j]));
        break;
      }
    for(j=0,bot=1;j<(*nc-1);j++) {
      reg[j]=exp(reg[j]);
      bot+=reg[j];
    }
    for(j=0;j<(*nc-1);j++)
      p[j]=reg[j]/bot;
    q=1/bot;
    if(y[i]!=0) {
      fit[i]=p[y[i]-1];
      *loglik+=freq[i]*log(p[y[i]-1]);
    }
    else {
      fit[i]=q;
      *loglik+=freq[i]*log(q);
    }
    pred[i]=1;
    cpred[i]=q;
    for(j=1;j<(*nc-1);j++) {
      if(p[pred[i]-1]<p[j])
        pred[i]=j+1;
      cpred[*n*j+i]=cpred[*n*(j-1)+i]+p[j-1];
    }
    if(p[pred[i]-1]<q)
      pred[i]=0;
  }
}

/* Minus log-likelihood for adjacent categories regresion. */
void mllacr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double *loglik,double top[]) {
  int i,j,k;
  double bot;

  for(i=0;i<*n;i++) {
    for(j=0;j<(*nc-1);j++)
      for(k=0,top[j]=-beta[j];k<*ncv;k++)
        top[j]-=beta[(*nc-1)+k]*cov[*n*k+i];
    for(j=1;j<(*nc-1);j++)
      top[j]+=top[j-1];
    for(j=0;j<(*nc-1);j++)
      switch(*transform) {
      case 2: /* sqrt transform */
        top[j]=top[j]*top[j];
        break;
      case 3: /* log transform */
        top[j]=exp(top[j]);
        break;
      case 4: /* logit transform */
        top[j]=exp(top[j]);
        top[j]=top[j]/(1+top[j]);
        break;
      case 5: /* cloglog transform */
        top[j]=1.0-exp(-exp(top[j]));
        break;
      case 6: /* square transform */
        top[j]=sqrt(top[j]);
        break;
      case 7: /* exp transform */
        top[j]=log(top[j]);
        break;
      case 8: /* loglog transform */
        top[j]=exp(-exp(top[j]));
        break;
      }
    for(j=0,bot=1;j<(*nc-1);j++)
      bot+=exp(top[j]);
    bot=log(bot);
    if(y[i]!=0)
      *loglik+=freq[i]*(top[y[i]-1]-bot);
    else
      *loglik-=freq[i]*bot;
  }
}

/* Fitted values for adjacent categories regresion. */
void fvacr(int y[],int freq[],double cov[],int *nc,int *ncv,int *n,double beta[],int *transform,double p[],double fit[],int pred[],double cpred[]) {
  int i,j,k;
  double q;

  for(i=0;i<*n;i++) {
    for(j=0;j<(*nc-1);j++)
      for(k=0,p[j]=-beta[j];k<*ncv;k++)
        p[j]-=beta[(*nc-1)+k]*cov[*n*k+i];
    for(j=1;j<(*nc-1);j++)
      p[j]+=p[j-1];
    for(j=0;j<(*nc-1);j++)
      switch(*transform) {
      case 2: /* sqrt transform */
        p[j]=p[j]*p[j];
        break;
      case 3: /* log transform */
        p[j]=exp(p[j]);
        break;
      case 4: /* logit transform */
        p[j]=exp(p[j]);
        p[j]=p[j]/(1+p[j]);
        break;
      case 5: /* cloglog transform */
        p[j]=1.0-exp(-exp(p[j]));
        break;
      case 6: /* square transform */
        p[j]=sqrt(p[j]);
        break;
      case 7: /* exp transform */
        p[j]=log(p[j]);
        break;
      case 8: /* loglog transform */
        p[j]=exp(-exp(p[j]));
        break;
      }
    for(j=0,q=1;j<(*nc-1);j++) {
      p[j]=exp(p[j]);
      q+=p[j];
    }
    for(j=0;j<(*nc-1);j++)
      p[j]=p[j]/q;
    q=1/q;
    if(y[i]!=0)
      fit[i]=p[y[i]-1];
    else
      fit[i]=q;
    pred[i]=1;
    cpred[i]=q;
    for(j=1;j<(*nc-1);j++) {
      if(p[pred[i]-1]<p[j])
        pred[i]=j+1;
      cpred[*n*j+i]=cpred[*n*(j-1)+i]+p[j-1];
    }
    if(p[pred[i]-1]<q)
      pred[i]=0;
  }
}

Generated by  Doxygen 1.6.0   Back to index