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

dist.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
 *
 * double plevy(double y,double m,double s)
 * double pginvgauss(double y,double m,double s,double f)
 * double ppowexp(double y,double m,double s,double f)
 * double psimplex(double y,double m,double s)
 * double pstable(double y,double loc,double disp,double skew,double tail)
 * double plaplace(double q,double m,double s)
 * double pglogis(double q,double m,double s,double f)
 * double pinvgauss(double q,double m,double s)
 * double pggamma(double q,double s,double m,double f)
 * double pgweibull(double q,double s,double m,double f)
 * double pgextval(double q,double s,double m,double f)
 * double phjorth(double q,double m,double s,double f)
 * double pburr(double q,double m,double s,double f)
 * double ppareto(double q,double m,double s)
 * double pboxcox(double q,double m,double s,double f)
 *
 *  DESCRIPTION
 *
 *    Functions to compute the cumulative probability functions for the
 * following continuous distributions:
 * Levy, generalized inverse Gaussian, power exponential, simplex, stable,
 * Laplace, generalized logistic, inverse gaussian, generalized gamma,
 * generalized Weibull, generalized extreme value, Hjorth, Burr, Pareto,
 * and Box-Cox.
 *
 */

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

extern double pcauchy(double x,double location,double scale,int lower_tail,int log_p);
extern double bessel_k(double x,double alpha,double expo);
extern double gammafn(double x);

static void interp(double x[],double fx[],int pts,double tab1[],double tab2[],double *f,double *df);
static double evalfn(int n,double a,double b,double arg1,double arg2,double arg3,double sum[],double (*fcn)(double x,double arg1,double arg2,double arg3));
static double romberg(double a,double b,double arg1,double arg2,double arg3,int pts,double tab1[],double tab2[],int iter[],int iterlim,double steptol,
                      double x[],double fx[],double (*fcn)(double x,double arg1,double arg2,double arg3));
static double dlevy(double x,double arg1,double arg2,double arg3);
double plevy(double y,double m,double s);
static double dginvgauss(double x,double arg1,double arg2,double arg3);
double pginvgauss(double y,double m,double s,double f);
static double dpowexp(double x,double arg1,double arg2,double arg3);
double ppowexp(double y,double m,double s,double f);
static double dsimplex(double x,double arg1,double arg2,double arg3);
double psimplex(double y,double m,double s);
static double fcn3(double s,double y,double alpha,double eta);
static double fcn4(double s,double y,double alpha,double eta);
double pstable(double y,double loc,double disp,double skew,double tail);
double plaplace(double q,double m,double s);
double pglogis(double q,double m,double s,double f);
double pinvgauss(double q,double m,double s);
double pggamma(double q,double s,double m,double f);
double pgweibull(double q,double s,double m,double f);
double pgextval(double q,double s,double m,double f);
double phjorth(double q,double m,double s,double f);
double pburr(double q,double m,double s,double f);
double ppareto(double q,double m,double s);
double pboxcox(double q,double m,double s,double f);

/* integration routine: stripped down version of Romberg integration from rmutil library */
static void interp(double x[],double fx[],int pts,double tab1[],double tab2[],double *f,double *df) {
  int ni,i,j;
  double tmp1,tmp2,lim1,lim2,diff1,diff2;

  ni=0;
  tmp1=fabs(x[0]);
  for(i=0;i<pts;i++) {
    tmp2=fabs(x[i]);
    if(tmp2<tmp1) {
      ni=i;
      tmp1=tmp2;
    }
    tab1[i]=tab2[i]=fx[i];
  }
  *f=fx[ni--];
  for(i=0;i<pts-1;i++) {
    for(j=0;j<pts-i-1;j++) {
      lim1=x[j];
      lim2=x[i+j+1];
      diff1=tab1[j+1]-tab2[j];
      diff2=lim1-lim2;
      if(diff2==0.0)
        return;
      diff2=diff1/diff2;
      tab2[j]=lim2*diff2;
      tab1[j]=lim1*diff2;
    }
    *df=2*ni<(pts-i-3)?tab1[ni+1]:tab2[ni--];
    *f+=*df;
  }
}

static double evalfn(int n,double a,double b,double arg1,double arg2,double arg3,double sum[],double (*fcn)(double x,double arg1,double arg2,double arg3)) {
  int i,j,nn;
  double x,tmpsum,pnt1,pnt2;

  if(n==1) {
    sum[0]=(b-a)*(*fcn)(0.5*(a+b),arg1,arg2,arg3);
    return(sum[0]);
  }
  else {
    for(i=1,j=1;i<n-1;i++)
      j*=3;
    pnt1=(b-a)/(3.0*j);
    pnt2=2.0*pnt1;
    x=a+0.5*pnt1;
    tmpsum=0.0;
    for(i=1;i<=j;i++) {
      tmpsum+=(*fcn)(x,arg1,arg2,arg3);
      x+=pnt2;
      tmpsum+=(*fcn)(x,arg1,arg2,arg3);
      x+=pnt1;
    }
    sum[0]=(sum[0]+(b-a)*tmpsum/j)/3.0;
    return(sum[0]);
  }
}

static double romberg(double a,double b,double arg1,double arg2,double arg3,int pts,double tab1[],double tab2[],int iter[],int iterlim,double steptol,
                      double x[],double fx[],double (*fcn)(double x,double arg1,double arg2,double arg3)) {
  int j,k;
  double sum[1],f,df;

  x[0]=1.0;
  for(iter[0]=0;iter[0]<iterlim;iter[0]++) {
    k=iter[0]+1;
    fx[iter[0]]=evalfn(k,a,b,arg1,arg2,arg3,sum,fcn);
    if(k>=pts) {
      interp(&x[k-pts],&fx[k-pts],pts,tab1,tab2,&f,&df);
      if(fabs(df)<steptol*fabs(f))
        return(f);
    }
    x[k]=x[iter[0]]/9.0;
    fx[k]=fx[iter[0]];
  }
  return(R_NaReal);
}

/* Levy distribution */
static double dlevy(double x,double arg1,double arg2,double arg3) {

  return(sqrt(arg2/(2.*M_PI*pow(x-arg1,3)))*exp(-arg2/(2.*(x-arg1))));
}

double plevy(double y,double m,double s) {
  int pts,iter[1],iterlim;
  double arg3,steptol,tab1[5],tab2[5],x[16],fx[16];

  if(y<=m)
    return(0.0);
  if(s<=0)
    return(R_NaReal);
  arg3=1.0;
  pts=5;
  iter[0]=0;
  iterlim=16;
  steptol=0.000001;
  return(romberg(m,y,m,s,arg3,pts,tab1,tab2,iter,iterlim,steptol,x,fx,dlevy));
}

/* generalized inverse Gaussian distribution */
static double dginvgauss(double x,double arg1,double arg2,double arg3) {

  return(pow(x,arg3-1)*exp(-(1./x+x/pow(arg1,2))/(2.*arg2))/(pow(arg1,arg3)*(2.*bessel_k(1/(arg2*arg1),fabs(arg3),1.0))));
}

double pginvgauss(double y,double m,double s,double f) {
  int pts,iter[1],iterlim;
  double t,steptol,tab1[5],tab2[5],x[16],fx[16];

  if(y<0)
    return(0.0);
  if(m<=0)
    return(R_NaReal);
  if(s<=0)
    return(R_NaReal);
  pts=5;
  iter[0]=0;
  iterlim=16;
  steptol=0.000001;
  t=0;
  return(romberg(t,y,m,s,f,pts,tab1,tab2,iter,iterlim,steptol,x,fx,dginvgauss));
}

/* power exponential distribution */
static double dpowexp(double x,double arg1,double arg2,double arg3) {
  double b,ss;

  ss=sqrt(arg2);
  b=1.+1./(2.*arg3);
  return(exp(-pow(fabs(x-arg1)/ss,2.*arg3)/2.)/(ss*gammafn(b)*pow(2.,b)));
}

double ppowexp(double y,double m,double s,double f) {
  int pts,iter[1],iterlim;
  double t,steptol,tab1[5],tab2[5],x[16],fx[16];

  if(s<=0)
    return(R_NaReal);
  if(f<=0)
    return(R_NaReal);
  pts=5;
  iter[0]=0;
  iterlim=16;
  steptol=0.000001;
  t=fabs(y-m)+m;
  if(y-m>0)
    return(0.5+romberg(m,t,m,s,f,pts,tab1,tab2,iter,iterlim,steptol,x,fx,dpowexp));
  else
    return(0.5-romberg(m,t,m,s,f,pts,tab1,tab2,iter,iterlim,steptol,x,fx,dpowexp));
}

/* simplex distribution */
static double dsimplex(double x,double arg1,double arg2,double arg3) {

  return(exp(-pow((x-arg1)/(arg1*(1-arg1)),2)/(2*x*(1-x)*arg2))/sqrt(2*M_PI*arg2*pow(x*(1-x),3)));
}

double psimplex(double y,double m,double s) {
  int pts,iter[1],iterlim;
  double arg3,t,steptol,tab1[5],tab2[5],x[16],fx[16];

  if(y<=0)
    return(0.0);
  if(y>=1)
    return(1.0);
  if(m<=0||m>=1)
    return(R_NaReal);
  if(s<=0)
    return(R_NaReal);
  arg3=1.0;
  pts=5;
  iter[0]=0;
  iterlim=16;
  steptol=0.000001;
  t=0;
  return(romberg(t,y,m,s,arg3,pts,tab1,tab2,iter,iterlim,steptol,x,fx,dsimplex));
}

/* stable distribution */
static double fcn3(double s,double y,double alpha,double eta) {
  double sa,seta,ceta;

  sa=pow(s,alpha);
  seta=sin(eta);
  ceta=cos(eta);
  return((sin(y*s-sa*seta)/s)*exp(-sa*ceta));
}

static double fcn4(double s,double y,double alpha,double eta) {
  double sa,seta,ceta;

  sa=pow(s,-alpha);
  seta=sin(eta);
  ceta=cos(eta);
  return(sin(y/s-sa*seta)*s*exp(-sa*ceta)/(s*s));
}

double pstable(double y,double loc,double disp,double skew,double tail) {
  int pts,iter[1],iterlim;
  double yy,steptol,eta,tab1[5],tab2[5],x[16],fx[16];

  if(disp<0)
    return(R_NaReal);
  if(skew<-1||skew>1)
    return(R_NaReal);
  if(tail<=0||tail>2)
    return(R_NaReal);
  if(tail==1&&skew==0)
    return(pcauchy(y,loc,disp,1,0));
  else
    if(tail==2) /* if(tail==2&&skew==0) */
      return(pnorm(y,loc,disp*sqrt(2),1,0));
    else {
      yy=(y-loc)/disp;
      pts=5;
      iter[0]=0;
      iterlim=16;
      steptol=0.000001;
      eta=skew*(1.0-fabs(1.0-tail))*M_PI/2.0;
      if(eta==0.0&&yy==0.0)
        return(0.5);
      else
        return(0.5+(romberg(0.0,1.0,yy,tail,eta,pts,tab1,tab2,iter,iterlim,steptol,x,fx,fcn3)+romberg(0.0,1.0,yy,tail,eta,pts,tab1,tab2,iter,iterlim,steptol,x,fx,fcn4))/M_PI);
    }
}

/* Laplace distribution */
double plaplace(double q,double m,double s) {
  double u,t;

  if(s<=0)
    return(R_NaReal);
  u=(q-m)/s;
  t=exp(-fabs(u))/2.0;
  if(u<0)
    return(t);
  else
    return(1.0-t);
}

/* generalized logistic distribution */
double pglogis(double q,double m,double s,double f) {

  if(s<=0)
    return(R_NaReal);
  if(f<=0)
    return(R_NaReal);
  return(pow(1.0+exp(-sqrt(3)*(q-m)/(s*M_PI)),-f));
}

/* inverse Gauss distribution */
double pinvgauss(double q,double m,double s) {
  double t,v;

  if(q<0)
    return(0.0);
  if(m<=0)
    return(R_NaReal);
  if(s<=0)
    return(R_NaReal);
  t=q/m;
  v=sqrt(q*s);
  return(pnorm((t-1.0)/v,0,1,1,0)+exp(2.0/(m*s))*pnorm(-(t+1)/v,0,1,1,0));
}

/* generalized gamma distribution */
double pggamma(double q,double s,double m,double f) {

  if(q<=0)
    return(0.0);
  if(m<=0)
    return(R_NaReal);
  if(s<=0)
    return(R_NaReal);
  if(f<=0)
    return(R_NaReal);
  return(pgamma(pow(q,f),s,pow(m/s,f),1,0));
}

/* generalized Weibull distribution (Mudholkar et al, 1995) */
double pgweibull(double q,double s,double m,double f) {

  if(q<=0)
    return(0.0);
  if(m<=0)
    return(R_NaReal);
  if(s<=0)
    return(R_NaReal);
  if(f<=0)
    return(R_NaReal);
  return(pow(1.0-exp(-pow(q/m,s)),f));
}

/* generalized extreme value distribution */
double pgextval(double q,double s,double m,double f) {
  int ind;
  double norm;

  if(q<=0)
    return(0.0);
  if(m<=0)
    return(R_NaReal);
  if(s<=0)
    return(R_NaReal);
  if(f==0) {
    return(R_NaReal); /* return(pweibull(1,s,m,1,0)); */
  }
  else {
    norm=sign(f)*exp(-pow(m,-s));
    ind=f>0;
    return((pweibull(exp(pow(q,f)/f),s,m,1,0)-ind+ind*norm)/(1-ind+norm));
  }
}

/* Hjorth distribution */
double phjorth(double q,double m,double s,double f) {

  if(q<0)
    return(0.0);
  if(m<=0)
    return(R_NaReal);
  if(s<=0)
    return(R_NaReal);
  return(1-pow(1+s*q,-f/s)*exp(-pow(q/m,2)/2));
}

/* Burr distribution */
double pburr(double q,double m,double s,double f) {

  if(q<0)
    return(0.0);
  if(m<=0)
    return(R_NaReal);
  if(s<=0)
    return(R_NaReal);
  if(f<=0)
    return(R_NaReal);
  return(1-pow(1+pow(q/m,s)/f,-f));
}

/* Pareto distribution */
double ppareto(double q,double m,double s) {

  if(q<0)
    return(0.0);
  if(m<=0)
    return(R_NaReal);
  if(s<=1)
    return(R_NaReal);
  return(1.0-pow(1.0+q/(m*(s-1.0)),-s));
}

/* Box-Cox distribution */
double pboxcox(double q,double m,double s,double f) {
  double norm;

  if(q<0)
    return(0.0);
  if(m<=0)
    return(R_NaReal);
  if(s<=0)
    return(R_NaReal);
  if(f==0) {
    return(R_NaReal); /* return(pnorm(0,m,sqrt(s),1,0)); */
  }
  else {
    norm=sign(f)*pnorm(0,m,sqrt(s),1,0);
    return((pnorm(pow(q,f)/f,m,sqrt(s),1,0)-(f>0)*norm)/(1-(f<0)-norm));
  }
}

Generated by  Doxygen 1.6.0   Back to index