Hi, As the bugs #275892 and #302873 are preventing r-noncran-lindsey to be build, and as they are open for a very long time, I am preparing an NMU of this package.
You will find attached the patch of this NMU. I'll upload it in 3 days unless you oppose. Bye, Aurelien -- .''`. Aurelien Jarno | GPG: 1024D/F1BCDB73 : :' : Debian GNU/Linux developer | Electrical Engineer `. `' [EMAIL PROTECTED] | [EMAIL PROTECTED] `- people.debian.org/~aurel32 | www.aurel32.net
diff -u r-noncran-lindsey-1.0.20041008/debian/changelog r-noncran-lindsey-1.0.20041008/debian/changelog --- r-noncran-lindsey-1.0.20041008/debian/changelog +++ r-noncran-lindsey-1.0.20041008/debian/changelog @@ -1,3 +1,12 @@ +r-noncran-lindsey (1.0.20041008-1.1) unstable; urgency=low + + * NMU. + * Don't build-depends on libreadline4-dev (closes: bug#302873). + * Applied patch from Kåre Hviid to make build possible with gcc-3.4 and + gcc-4.0 (closes: bug#275892). + + -- Aurelien Jarno <[EMAIL PROTECTED]> Sat, 20 Aug 2005 11:34:06 +0200 + r-noncran-lindsey (1.0.20041008-1) unstable; urgency=low * New upstream releases. diff -u r-noncran-lindsey-1.0.20041008/debian/control r-noncran-lindsey-1.0.20041008/debian/control --- r-noncran-lindsey-1.0.20041008/debian/control +++ r-noncran-lindsey-1.0.20041008/debian/control @@ -3,7 +3,7 @@ Priority: optional Maintainer: Chris Lawrence <[EMAIL PROTECTED]> Standards-Version: 3.6.1 -Build-depends: r-base-dev (>> 2.0.0), debhelper (>> 4), libncurses5-dev, libreadline4-dev, g77 [!m68k], f2c [m68k] +Build-Depends: r-base-dev (>> 2.0.0), debhelper (>> 4), libncurses5-dev, g77 [!m68k], f2c [m68k] Package: r-noncran-lindsey Architecture: any diff -u r-noncran-lindsey-1.0.20041008/debian/rules r-noncran-lindsey-1.0.20041008/debian/rules --- r-noncran-lindsey-1.0.20041008/debian/rules +++ r-noncran-lindsey-1.0.20041008/debian/rules @@ -46,7 +46,9 @@ dh_installdirs (mkdir build; cd build; for lib in ../rmutil.tgz ../event.tgz ../growth.tgz ../repeated.tgz ../stable.tgz ../gnlm.tgz ../ordinal03.tgz; \ do tar -zxf $$lib; done; \ - cp ../fixed/*.c ./ordinal/src; cp ../fixed/Makefile ./ordinal/src; cp ../fixed/DESCRIPTION ./ordinal; \ + cp ../fixed/dist.c ../fixed/lcr.c ./ordinal/src; cp ../fixed/Makefile ./ordinal/src; cp ../fixed/DESCRIPTION ./ordinal; \ + cp ../fixed/ksurvb.c ../fixed/ksurvg.c ./event/src; \ + cp ../fixed/kcountb.c ./repeated/src; \ R CMD INSTALL -l $(debRlib) --clean rmutil event growth repeated stable gnlm ordinal) rm -vf $(debRlib)/R.css (find $(debRlib) -name 'COPYING' -print0 | xargs -0 rm -vf) only in patch2: unchanged: --- r-noncran-lindsey-1.0.20041008.orig/fixed/ksurvg.c +++ r-noncran-lindsey-1.0.20041008/fixed/ksurvg.c @@ -0,0 +1,223 @@ +/* + * event : A Library of Special Functions for Event Histories + * Copyright (C) 1998, 1999, 2000, 2001 J.K. 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 ksurvg(double p[],double y[],double x[],int cens[],int *nind, + * int nobs[],int *nbs,int *nccov,int *model,int *dist,int *density, + * int *dep,int *birth,int *tvc,double tvcov[],int *fit,double pred[], + * double rpred[],int *renewal,int *rf,double bb[], + * int *sf,double vv[],double *like) + * + * DESCRIPTION + * + * Function to compute the likelihood function for various distributions + * inserted in a gamma or Weibull distribution with serial dependence using + * Kalman-type update for event histories. + * + */ + +#include <math.h> +#include <stddef.h> +#include "R.h" +#include "Rmath.h" + +void ksurvg(double p[],double y[],double x[],int cens[],int *nind,int nobs[], + int *nbs,int *nccov,int *model,int *dist,int *density,int *dep, + int *birth,int *tvc,double tvcov[],int *fit,double pred[], + double rpred[],int *renewal,int *rf,double bb[], + int *sf,double vv[],double *like){ + int i,j,k,nm,c; + double a,a1,b,H,delta,lambda,omega,om,beta,bt,l1,yy,kk,ly,plap,tmp, + intercept; + + *like=0; + nm=0; + delta=exp(-p[*nccov+*birth+*tvc+1]); + if(*dep>0)omega=exp(p[*nccov+*birth+*tvc+2])/(1+exp(p[*nccov+*birth+*tvc+2])); + if(*model>1&&!*sf){ + if(*model<5)lambda=exp(p[*nccov+*birth+*tvc+2+(*dep>0)]); + else lambda=exp(p[*nccov+*birth+*tvc+2+(*dep>0)]/2);} + if(*model==4)intercept=exp(p[*nccov+*birth+*tvc+3+(*dep>0)]); + for(i=0;i<*nind;i++){ + a=b=delta; + if(!*rf){ + beta=p[0]; + for(k=0;k<*nccov;k++)beta+=p[k+1]*x[i+k**nind]; + if(*model<4){ + if(beta>40) beta=40; + if(beta<-40)beta=-40; + beta=exp(beta);}} + else if(!*tvc)bt=bb[i]; + c=0; + yy=0; + for(j=0;j<nobs[i];j++){ + if(*model>1&&*sf)lambda=vv[nm]; + /* store value and check if ties to follow */ + if(y[nm]>0){ + if(*renewal)yy=y[nm]; + else yy+=y[nm];} + c+=cens[nm]; + if(j>=nobs[i]-1||y[nm+1]!=0){ + if(*model>=5)ly=log(yy); + a1=a+c; + /* add in birth and time-varying covariates */ + if(!*rf){ + if(*tvc){ + bt=0; + for(k=0;k<*tvc;k++)bt+=p[*nccov+*birth+k+1]*tvcov[nm+*nbs*k]; + if(*model<4){ + if(bt>40) bt=40; + if(bt<-40)bt=-40; + bt=exp(bt)*beta;} + else bt+=beta;} + else bt=beta; + if(*birth){ + if(*model<4)bt*=pow(j+1.,p[*nccov+1]); + else bt+=p[*nccov+1]*log(j+1);}} + else if(*tvc)bt=bb[nm]; + if(!*density){ + /* intensity models */ + switch(*model){ + case 1: /* exponential distribution */ + H=yy/bt; + l1=1/bt; + break; + case 2: /* Weibull distribution */ + H=pow(yy/bt,lambda); + l1=lambda*pow(yy/bt,lambda-1)/bt; + break; + case 3: /* gamma distribution */ + H=-pgamma(yy,lambda,bt,0,1); + l1=dgamma(yy,lambda,bt,0)/(1-pgamma(yy,lambda,bt,1,0)); + break; + case 4: /* generalized logistic distribution */ + H=(yy+log(lambda+intercept*exp(-bt*yy))/bt)/lambda; + l1=1/(lambda+intercept*exp(-bt*yy)); + break; + case 5: /* log normal distribution */ + H=-pnorm(ly,bt,lambda,0,1); + l1=dnorm(ly,bt,lambda,0)/yy/(1-pnorm(ly,bt,lambda,1,0)); + break; + case 6: /* log logistic distribution */ + H=-plogis(ly,bt,lambda,0,1); + l1=dlogis(ly,bt,lambda,0)/yy/(1-plogis(ly,bt,lambda,1,0)); + break; + case 7: /* log Cauchy distribution */ + H=-pcauchy(ly,bt,lambda,0,1); + l1=dcauchy(ly,bt,lambda,0)/yy/(1-pcauchy(ly,bt,lambda,1,0)); + break; + case 8: /* log Laplace distribution */ + tmp=exp(-fabs(ly-bt)/lambda)/2; + plap=ly<bt?tmp:1-tmp; + H=-log(1-plap); + l1=tmp/(lambda*yy*(1-plap)); + break;}} + else{ + /* density models */ + switch(*model){ + case 1: /* exponential distribution */ + H=pexp(yy,bt,1,0); + l1=dexp(yy,bt,0); + break; + case 2: /* Weibull distribution */ + H=pweibull(yy,lambda,bt,1,0); + l1=dweibull(yy,lambda,bt,0); + break; + case 3: /* gamma distribution */ + H=pgamma(yy,lambda,bt,1,0); + l1=dgamma(yy,lambda,bt,0); + break; + case 4: /* generalized logistic distribution */ + H=exp(-yy/lambda)*pow((lambda+intercept)/(lambda+intercept*exp(-bt*yy)),1/(lambda*bt)); + l1=exp(-yy/lambda)*pow((lambda+intercept)/(lambda+intercept*exp(-bt*yy)),1/(lambda*bt)+1); + break; + case 5: /* log normal distribution */ + H=pnorm(ly,bt,lambda,1,0); + l1=dnorm(ly,bt,lambda,0)/yy; + break; + case 6: /* log logistic distribution */ + H=plogis(ly,bt,lambda,1,0); + l1=dlogis(ly,bt,lambda,0)/yy; + break; + case 7: /* log Cauchy distribution */ + H=pcauchy(ly,bt,lambda,1,0); + l1=dcauchy(ly,bt,lambda,0)/yy; + break; + case 8: /* log Laplace distribution */ + tmp=exp(-fabs(ly-bt)/lambda)/2; + H=ly<bt?tmp:1-tmp; + l1=tmp/lambda/yy; + break;}} + /* calculate likelihood */ + if(l1<1e-20)l1=1e-20; + if(H<1e-20)H=1e-20; + if(H>1e20)H=1e20; + /* is this correct when c>1 ? */ + if(*dist==2){ + if(c)*like-=c*log(l1)-b*H-lgammafn(a)+log(b)+(a-1)*log(b*H); + /* if(c)*like-=c*log(l1)+log(dgamma(H,a,1/b));*/ + else *like-=pgamma(H,a,1/b,0,1);} + else if(*dist==3){ + if(c)*like-=c*log(l1)-pow(b*H,a)+(a-1)*log(b*H)+log(b)+log(a); + else *like+=pow(b*H,a);} + if(c>1)*like+=lgammafn(c+1); + if(*fit){ + pred[nm-c+cens[nm]]=bt; + tmp=a/b; + if(!*density){ + switch(*model){ + case 1: rpred[nm-c+cens[nm]]=bt*tmp; break; + case 2: rpred[nm-c+cens[nm]]=bt*pow(tmp,1/lambda); break; + case 3: rpred[nm-c+cens[nm]]=qgamma(1-exp(-tmp),lambda,bt,1,0); break; + case 5: rpred[nm-c+cens[nm]]=exp(qnorm(1-exp(-tmp),bt,lambda,1,0)); break; + case 6: rpred[nm-c+cens[nm]]=exp(qlogis(1-exp(-tmp),bt,lambda,1,0)); break; + case 7: rpred[nm-c+cens[nm]]=exp(qcauchy(1-exp(-tmp),bt,lambda,1,0)); break; + case 8: rpred[nm-c+cens[nm]]=exp(bt+lambda*log(2*(ly<bt?exp(-tmp):1-exp(-tmp)))); break;}} + else{ + switch(*model){ + case 1: rpred[nm-c+cens[nm]]=qexp(tmp,bt,1,0); break; + case 2: rpred[nm-c+cens[nm]]=qweibull(tmp,lambda,bt,1,0); break; + case 3: rpred[nm-c+cens[nm]]=qgamma(tmp,lambda,bt,1,0); break; + case 5: rpred[nm-c+cens[nm]]=exp(qnorm(tmp,bt,lambda,1,0)); break; + case 6: rpred[nm-c+cens[nm]]=exp(qlogis(tmp,bt,lambda,1,0)); break; + case 7: rpred[nm-c+cens[nm]]=exp(qcauchy(tmp,bt,lambda,1,0)); break; + case 8: rpred[nm-c+cens[nm]]=exp(bt+lambda*log(2*(ly<bt?tmp:1-tmp))); break;}}} + /* update parameters */ + switch(*dep){ + case 1: + case 2: om=pow(omega,yy); a=om*a1+(1-om)*delta; break; + case 3: + case 4: a=omega*a1+(1-omega)*(delta-1); break; + case 5: a=a1; break; + case 6: + case 7: a=omega*a1; break; + default: ; + } + switch(*dep){ + case 1: b=pow(b/(H*delta),om)*delta; break; + case 2: b=pow(H,-om)*delta; break; + case 3: b=pow(b/(H*delta),omega)*delta; break; + case 4: b=pow(H,-omega)*delta; break; + case 5: + case 7: b=pow(b/H,omega); break; + default: ; + } + c=0;} + nm++;}} + return;} only in patch2: unchanged: --- r-noncran-lindsey-1.0.20041008.orig/fixed/ksurvb.c +++ r-noncran-lindsey-1.0.20041008/fixed/ksurvb.c @@ -0,0 +1,396 @@ +/* + * event : A Library of Special Functions for Event Histories + * Copyright (C) 1998, 1999, 2000, 2001 J.K. 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 ksurvb(double p[],double y[],double x[],int cens[],int *nind, + * int nobs[],int *nbs,int *nccov,int *model,int *density, + * int *pfamily,int *dep,int *birth,int *tvc,double tvcov[], + * int *fit,double pred[],double rpred[],int *renewal,int *rf, + * double bb[],int *sf,double vv[],double *like) + * void frailb(double p[],double y[],double x[],int cens[],int *nind, + * int nobs[],int *nbs,int *nccov,int *model,int *density,int *dep, + * int *birth,int *tvc,double tvcov[],int *fit,double pred[], + * double rpred[],int *rf,double bb[],int *sf,double vv[], + * int *frser,double *like) + * + * DESCRIPTION + * + * Function to compute the likelihood function for various distributions + * inserted in a beta distribution with serial dependence or gamma frailties + * using Kalman-type update for event histories. + * + */ + +#include <math.h> +#include <stddef.h> +#include "R.h" +#include "Rmath.h" + +void ksurvb(double p[],double y[],double x[],int cens[],int *nind, + int nobs[],int *nbs,int *nccov,int *model,int *density, + int *pfamily,int *dep,int *birth,int *tvc,double tvcov[],int *fit, + double pred[],double rpred[],int *renewal,int *rf,double bb[], + int *sf,double vv[],double *like){ + int i,j,k,nm,c; + double a,a1,b,b1,delta,lambda,omega,om,beta,bt,h,yy,kk,tmp,ly,plap,intercept, + family; + + *like=0; + nm=0; + delta=exp(-p[*nccov+*birth+*tvc+1]); + if(*dep>0) + omega=exp(p[*nccov+*birth+*tvc+2])/(1+exp(p[*nccov+*birth+*tvc+2])); + if(*pfamily)family=p[*nccov+*birth+*tvc+2+(*dep>0)]; + if(*model==4)intercept=exp(p[*nccov+*birth+*tvc+2+(*dep>0)+*pfamily]);; + if(*model>1&&!*sf){ + if(*model<5) + lambda=exp(p[*nccov+*birth+*tvc+2+(*dep>0)+*pfamily+(*model==4)]); + else lambda=exp(p[*nccov+*birth+*tvc+2+(*dep>0)+*pfamily+(*model==4)]/2);} + for(i=0;i<*nind;i++){ + a=b=delta; + if(!*rf){ + beta=p[0]; + for(k=0;k<*nccov;k++)beta+=p[k+1]*x[i+k**nind]; + if(*model<4){ + if(beta>40) beta=40; + if(beta<-40)beta=-40; + beta=exp(beta);}} + else if(!*tvc)bt=bb[i]; + c=0; + yy=0; + for(j=0;j<nobs[i];j++){ + if(*model>1&&*sf)lambda=vv[nm]; + /* store value and check if ties to follow */ + if(y[nm]>0){ + if(*renewal)yy=y[nm]; + else yy+=y[nm];} + c+=cens[nm]; + if(j>=nobs[i]-1||y[nm+1]!=0){ + /* if no ties follow, update the likelihood */ + if(*model>=5)ly=log(yy); + a1=a+c; + b1=b; + /* add in birth and time-varying covariates */ + if(!*rf){ + if(*tvc){ + bt=0; + for(k=0;k<*tvc;k++)bt+=p[*nccov+*birth+k+1]*tvcov[nm+*nbs*k]; + if(*model<4){ + if(bt>40) bt=40; + if(bt<-40)bt=-40; + bt=exp(bt)*beta;} + else bt+=beta;} + else bt=beta; + if(*birth){ + if(*model<4)bt*=pow(j+1.,p[*nccov+1]); + else bt+=p[*nccov+1]*log(j+1);}} + else if(*tvc)bt=bb[nm]; + if(!*density){ + /* intensity models */ + switch(*model){ + case 1: /* exponential distribution */ + b1+=yy/bt; + h=-log(bt); + break; + case 2: /* Weibull distribution */ + b1+=pow(yy/bt,lambda); + h=log(lambda/bt)+(lambda-1)*log(yy/bt); + break; + case 3: /* gamma distribution */ + b1-=pgamma(yy,lambda,bt,0,1); + h=dgamma(yy,lambda,bt,1)-pgamma(yy,lambda,bt,0,1); + break; + case 4: /* generalized logistic distribution */ + b1+=(yy+log(lambda+intercept*exp(-bt*yy))/bt)/lambda; + h=-log(lambda+intercept*exp(-bt*yy)); + break; + case 5: /* log normal distribution */ + b1-=pnorm(ly,bt,lambda,0,1); + h=dnorm(ly,bt,lambda,1)-ly-pnorm(ly,bt,lambda,0,1); + break; + case 6: /* log logistic distribution */ + b1-=plogis(ly,bt,lambda,0,1); + h=dlogis(ly,bt,lambda,1)-ly-plogis(ly,bt,lambda,0,1); + break; + case 7: /* log Cauchy distribution */ + b1-=pcauchy(ly,bt,lambda,0,1); + h=dcauchy(ly,bt,lambda,1)-ly-pcauchy(ly,bt,lambda,0,1); + break; + case 8: /* log Laplace distribution */ + tmp=exp(-fabs(ly-bt)/lambda)/2; + plap=ly<bt?tmp:1-tmp; + b1-=log(1-plap); + h=log(tmp)-log(lambda*yy*(1-plap)); + break;}} + else{ + /* density models */ + switch(*model){ + case 1: /* exponential distribution */ + b1+=pexp(yy,bt,1,0); + h=dexp(yy,bt,1); + break; + case 2: /* Weibull distribution */ + b1+=pweibull(yy,lambda,bt,1,0); + h=dweibull(yy,lambda,bt,1); + break; + case 3: /* gamma distribution */ + b1+=pgamma(yy,lambda,bt,1,0); + h=dgamma(yy,lambda,bt,1); + break; + case 4: /* generalized logistic distribution */ + b1+=exp(-yy/lambda)*pow((lambda+intercept)/(lambda+intercept*exp(-bt*yy)),1/(lambda*bt)); + h=-yy/lambda+(1/(lambda*bt)+1)*log((lambda+intercept)/(lambda+intercept*exp(-bt*yy))); + break; + case 5: /* log normal distribution */ + b1+=pnorm(ly,bt,lambda,1,0); + h=dnorm(ly,bt,lambda,1)-ly; + break; + case 6: /* log logistic distribution */ + b1+=plogis(ly,bt,lambda,1,0); + h=dlogis(ly,bt,lambda,1)-ly; + break; + case 7: /* log Cauchy distribution */ + b1+=pcauchy(ly,bt,lambda,1,0); + h=dcauchy(ly,bt,lambda,1)-ly; + break; + case 8: /* log Laplace distribution */ + tmp=exp(-fabs(ly-bt)/lambda)/2; + b1+=ly<bt?tmp:1-tmp; + h=log(tmp/lambda/yy); + break;}} + /* calculate likelihood */ + *like-=(c>1?lgammafn(a1)-lgammafn(a)-lgammafn(c+1):c*log(a))+c*h; + /* is this correct for c>1? - gives gamma when family -> 0 */ + /* c*(family-1)*log(b1) does not work */ + if(*pfamily)*like-=(c>0)*(family-c)*log(b1)-a*(pow(b1,family)-pow(b,family))/family; + else *like-=a*log(b)-a1*log(b1); + /* calculate fitted values */ + if(*fit){ + pred[nm-c+cens[nm]]=bt; + tmp=b/a; + if(!*density){ + switch(*model){ + case 1: rpred[nm-c+cens[nm]]=bt*tmp; break; + case 2: rpred[nm-c+cens[nm]]=bt*pow(tmp,1/lambda); break; + case 3: rpred[nm-c+cens[nm]]=qgamma(1-exp(-tmp),lambda,bt,1,0); break; + case 5: rpred[nm-c+cens[nm]]=exp(qnorm(1-exp(-tmp),bt,lambda,1,0)); break; + case 6: rpred[nm-c+cens[nm]]=exp(qlogis(1-exp(-tmp),bt,lambda,1,0)); break; + case 7: rpred[nm-c+cens[nm]]=exp(qcauchy(1-exp(-tmp),bt,lambda,1,0)); break; + case 8: rpred[nm-c+cens[nm]]=exp(bt+lambda*log(2*(ly<bt?exp(-tmp):1-exp(-tmp)))); break;}} + else{ + switch(*model){ + case 1: rpred[nm-c+cens[nm]]=qexp(tmp,bt,1,0); break; + case 2: rpred[nm-c+cens[nm]]=qweibull(tmp,lambda,bt,1,0); break; + case 3: rpred[nm-c+cens[nm]]=qgamma(tmp,lambda,bt,1,0); break; + case 5: rpred[nm-c+cens[nm]]=exp(qnorm(tmp,bt,lambda,1,0)); break; + case 6: rpred[nm-c+cens[nm]]=exp(qlogis(tmp,bt,lambda,1,0)); break; + case 7: rpred[nm-c+cens[nm]]=exp(qcauchy(tmp,bt,lambda,1,0)); break; + case 8: rpred[nm-c+cens[nm]]=exp(bt+lambda*log(2*(ly<bt?tmp:1-tmp))); break;}}} + /* update parameters */ + switch(*dep){ + case 1: + case 2: om=pow(omega,yy); a=om*a1+(1-om)*delta; break; + case 3: + case 4: a=omega*a1+(1-omega)*delta; break; + case 5: a=a1; break; + case 6: + case 7: a=omega*a1; break; + default: ; + } + switch(*dep){ + case 1: b=om*(b1-b)+delta; break; + case 2: b=om*b1+(1-om)*delta; break; + case 3: b=omega*(b1-b)+delta; break; + case 4: b=omega*b1+(1-omega)*delta; break; + case 5: + case 7: b=omega*b1; break; + default: ; + } + c=0;} + nm++;}} + return;} + +void frailb(double p[],double y[],double x[],int cens[],int *nind,int nobs[], + int *nbs,int *nccov,int *model,int *density,int *dep,int *birth, + int *tvc,double tvcov[],int *fit,double pred[],double rpred[], + int *rf,double bb[],int *sf,double vv[],int *frser, + double *like){ + int i,j,k,nm,ns,c,nn; + double b1,delta,lambda,beta,bt,l1,yy,kk,nb,ly,plap,tmp,intercept,H,btp,res; + + *like=0; + nm=0; + delta=exp(p[*nccov+*birth+*tvc+1]); + if(*model>1&&!*sf){ + if(*model<5)lambda=exp(p[*nccov+*birth+*tvc+*frser+2]); + else lambda=exp(p[*nccov+*birth+*tvc+*frser+2]/2);} + if(*model==4)intercept=exp(p[*nccov+*birth+*tvc+*frser+3]); + for(nn=i=0;i<*nind;i++)nn+=nobs[i]; + for(i=0;i<*nind;i++){ + if(!*rf){ + beta=p[0]; + for(k=0;k<*nccov;k++)beta+=p[k+1]*x[i+k**nind]; + if(*model<4){ + if(beta>40) beta=40; + if(beta<-40)beta=-40; + beta=exp(beta);}} + else if(!*tvc)bt=bb[i]; + ns=b1=0; + if(!*dep)nb=1; + else { + nb=0; + for(j=0;j<nobs[i];j++)nb+=y[nm+j]; + nb/=beta;} + res=0; + for(c=0,j=0;j<nobs[i];j++){ + if(*model>1&&*sf)lambda=vv[nm]; + ns+=cens[nm]; + /* store value and check if ties to follow */ + if(y[nm]>0)yy=y[nm]; + c+=cens[nm]; + if(j>=nobs[i]-1||y[nm+1]!=0){ + if(*model>=5)ly=log(yy); + l1=log(1+delta*j/nb); + /* add in birth and time-varying covariates */ + if(!*rf){ + if(*tvc){ + bt=0; + for(k=0;k<*tvc;k++)bt+=p[*nccov+*birth+k+1]*tvcov[nm+*nbs*k]; + if(*model<4)bt=exp(bt)*beta; + else bt+=beta;} + else bt=beta; + if(*birth){ + if(*model<4)bt*=pow(j+1.,p[*nccov+1]); + else bt+=p[*nccov+1]*log(j+1);}} + else if(*tvc)bt=bb[nm]; + /* if AR, add discounted previous residual */ + btp=bt; + if(*frser&&j>0){ + if(*model<4) + bt=exp(log(bt)+exp(p[*nccov+*birth+*tvc+2]*y[nm])*res); + else bt+=exp(p[*nccov+*birth+*tvc+2]*y[nm])*res;} + if(!*density){ + /* intensity models */ + switch(*model){ + case 1: /* exponential distribution */ + H=yy/bt; + l1+=-log(bt); + break; + case 2: /* Weibull distribution */ + H=pow(yy/bt,lambda); + l1+=log(lambda/bt)+(lambda-1)*log(yy/bt); + break; + case 3: /* gamma distribution */ + H=-pgamma(yy,lambda,bt,0,1); + l1+=dgamma(yy,lambda,bt,1)-pgamma(yy,lambda,bt,0,1); + break; + case 4: /* generalized logistic distribution */ + H=(yy+log(lambda+intercept*exp(-bt*yy))/bt)/lambda; + l1+=-log(lambda+intercept*exp(-bt*yy)); + break; + case 5: /* log normal distribution */ + H=-pnorm(ly,bt,lambda,0,1); + l1+=dnorm(ly,bt,lambda,1)-ly-pnorm(ly,bt,lambda,0,1); + break; + case 6: /* log logistic distribution */ + H=-plogis(ly,bt,lambda,0,1); + l1+=dlogis(ly,bt,lambda,1)/-ly-plogis(ly,bt,lambda,0,1); + break; + case 7: /* log Cauchy distribution */ + H=-pcauchy(ly,bt,lambda,0,1); + l1+=dcauchy(ly,bt,lambda,1)-ly-pcauchy(ly,bt,lambda,0,1); + break; + case 8: /* log Laplace distribution */ + tmp=exp(-fabs(ly-bt)/lambda)/2; + plap=ly<bt?tmp:1-tmp; + H=-log(1-plap); + l1+=tmp/(lambda*yy*(1-plap)); + break;}} + else{ + /* density models */ + switch(*model){ + case 1: /* exponential distribution */ + H=pexp(yy,bt,1,0); + l1+=dexp(yy,bt,1); + break; + case 2: /* Weibull distribution */ + H=pweibull(yy,lambda,bt,1,0); + l1+=dweibull(yy,lambda,bt,1); + break; + case 3: /* gamma distribution */ + H=pgamma(yy,lambda,bt,1,0); + l1+=dgamma(yy,lambda,bt,1); + break; + case 4: /* generalized logistic distribution */ + H=exp(-yy/lambda)*pow((lambda+intercept)/(lambda+intercept*exp(-bt*yy)),1/(lambda*bt)); + l1+=-yy/lambda+log((lambda+intercept)/(lambda+intercept*exp(-bt*yy)))/((lambda*bt)+1); + break; + case 5: /* log normal distribution */ + H=pnorm(ly,bt,lambda,1,0); + l1+=dnorm(ly,bt,lambda,1)-ly; + break; + case 6: /* log logistic distribution */ + H=plogis(ly,bt,lambda,1,0); + l1+=dlogis(ly,bt,lambda,1)-ly; + break; + case 7: /* log Cauchy distribution */ + H=pcauchy(ly,bt,lambda,1,0); + l1+=dcauchy(ly,bt,lambda,1)-ly; + break; + case 8: /* log Laplace distribution */ + tmp=exp(-fabs(ly-bt)/lambda)/2; + H=ly<bt?tmp:1-tmp; + l1+=log(tmp/lambda/yy); + break;}} + /* calculate likelihood */ + *like-=c*l1; + if(c>1)*like+=lgammafn(c+1); + if(*frser)res=(*model>6?ly:yy)-btp; + if(*fit){ + pred[nm-c+cens[nm]]=btp; + tmp=(b1+nb/(nn*delta))/(nb/(nn*delta)+ns); + if(!*density){ + switch(*model){ + case 1: rpred[nm-c+cens[nm]]=bt*tmp; break; + case 2: rpred[nm-c+cens[nm]]=bt*pow(tmp,1/lambda); break; + case 3: rpred[nm]=qgamma(1-exp(-tmp),lambda,bt,1,0); break; + case 5: rpred[nm-c+cens[nm]]=exp(qnorm(1-exp(-tmp),bt,lambda,1,0)); + break; + case 6: rpred[nm-c+cens[nm]]=exp(qlogis(1-exp(-tmp),bt,lambda,1,0)); + break; + case 7: rpred[nm-c+cens[nm]]=exp(qcauchy(1-exp(-tmp),bt,lambda,1,0)); + break; + case 8: rpred[nm-c+cens[nm]]=exp(bt+lambda*log(2*(ly<bt?exp(-tmp):1-exp(-tmp)))); + break;}} + else{ + switch(*model){ + case 1: rpred[nm-c+cens[nm]]=qexp(tmp,bt,1,0); break; + case 2: rpred[nm-c+cens[nm]]=qweibull(tmp,lambda,bt,1,0); break; + case 3: rpred[nm-c+cens[nm]]=qgamma(tmp,lambda,bt,1,0); break; + case 5: rpred[nm-c+cens[nm]]=exp(qnorm(tmp,bt,lambda,1,0)); break; + case 6: rpred[nm-c+cens[nm]]=exp(qlogis(tmp,bt,lambda,1,0)); break; + case 7: rpred[nm-c+cens[nm]]=exp(qcauchy(tmp,bt,lambda,1,0)); + break; + case 8: rpred[nm-c+cens[nm]]=exp(bt+lambda*log(2*(ly<bt?tmp:1-tmp))); + break;}}} + b1+=H; + c=0;} + nm++;} + *like+=(nb/delta+ns)*log(1+delta*b1/nb);} + return;} only in patch2: unchanged: --- r-noncran-lindsey-1.0.20041008.orig/fixed/kcountb.c +++ r-noncran-lindsey-1.0.20041008/fixed/kcountb.c @@ -0,0 +1,443 @@ +/* + * repeated : A Library of Repeated Measurements Models + * Copyright (C) 1998, 1999, 2000, 2001 J.K. 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 kcountb(double p[],double y[],double *origin,int c[],double x[], + * int *nind,int nobs[],int *nbs,int *nccov,int *model, + * int *density,int *pfamily,int *dep,int *birth,int *tvc,double tvcov[], + * int *fit,double pred[],double rpred[],int *rf, double bbb[], + * int *sf, double vv[], double *like) + * void countfb(double p[],double y[],int c[], double x[],int *nind, + * int nobs[],int *nbs,int *nccov,int *model,int *density, + * int *tvc,double tvcov[],int *fit,double pred[],double rpred[],int *rf, + * double bbb[],int *sf,double vv[],int *frser,double *like) + * + * DESCRIPTION + * + * Functions to compute the likelihood function for various distributions + * inserted in a beta distribution with serial or frailty dependence using + * Kalman-type update for longitudinal count data. + * + */ + +#include <math.h> +#include <stddef.h> +#include "R.h" +#include "Rmath.h" + +static double dpvfp2(int y, double d, double s, double s1, double f); + +void kcountb(double p[],double y[],double *origin,int c[],double x[], + int *nind,int nobs[],int *nbs,int *nccov,int *model, + int *density,int *pfamily,int *dep,int *birth,int *tvc,double tvcov[], + int *fit,double pred[],double rpred[],int *rf, double bbb[], + int *sf, double vv[], double *like){ + int i,j,j0,k,nm; + double a,a1,b,b1,bb,bb0,sc,delta,lambda,omega,om,beta,bt,H,yy,yy0, + kk,tmp,ly,ly0,plap,intercept,family; + + *like=0; + nm=0; + delta=exp(-p[*nccov+*birth+*tvc+1]); + if(*dep>0)omega=exp(p[*nccov+*birth+*tvc+2])/(1+exp(p[*nccov+*birth+*tvc+2])); + if(*pfamily)family=p[*nccov+*birth+*tvc+2+(*dep>0)]; + if(*model==4)intercept=exp(p[*nccov+*birth+*tvc+2+(*dep>0)+*pfamily]); + if(*model>1&&!*sf){ + if(*model<5) + lambda=exp(p[*nccov+*birth+*tvc+2+(*dep>0)+*pfamily+(*model==4)]); + else lambda=exp(p[*nccov+*birth+*tvc+2+(*dep>0)+*pfamily+(*model==4)]/2);} + for(i=0;i<*nind;i++){ + a=b=delta; + sc=bb0=ly0=0; + yy0=*origin; + if(!*rf){ + beta=p[0]; + for(k=0;k<*nccov;k++)beta+=p[k+1]*x[i+k**nind]; + if(*model<4){ + if(beta>40) beta=40; + if(beta<-40)beta=-40; + beta=exp(beta);}} + else if(!*tvc)bt=bbb[i]; + j0=origin==0||*tvc?0:-1; + for(j=j0;j<nobs[i];j++){ + yy=*origin+(j>-1?y[nm]:0); + if(*model>=5)ly=log(yy); + if(*model>1&&*sf)lambda=vv[nm]; + a1=a+(j>-1?c[nm]:0); + b1=b; + /* add in birth and time-varying covariates */ + if(!*rf){ + if(*tvc){ + bt=0; + for(k=0;k<*tvc;k++)bt+=p[*nccov+*birth+k+1]*tvcov[nm+*nbs*k]; + if(*model<4){ + if(bt>40) bt=40; + if(bt<-40)bt=-40; + bt=exp(bt)*beta;} + else bt+=beta;} + else bt=beta; + if(j>-1&&*birth){ + sc+=c[nm]; + if(*model<5)bt*=pow(sc,p[*nccov+1]); + else bt+=p[*nccov+1]*log(sc);}} + else if(*tvc)bt=bbb[nm]; + if(!*tvc){ + if(!*density){ + /* intensity models */ + switch(*model){ + case 1: /* exponential distribution */ + H=yy/bt; break; + case 2: /* Weibull distribution */ + H=pow(yy/bt,lambda); break; + case 3: /* gamma distribution */ + H=-pgamma(yy,lambda,bt,0,1); break; + case 4: /* generalized logistic distribution */ + H=(yy+log(lambda+intercept*exp(-bt*yy))/bt)/lambda; break; + case 5: /* log normal distribution */ + H=-pnorm(ly,bt,lambda,0,1); break; + case 6: /* log logistic distribution */ + H=-plogis(ly,bt,lambda,0,1); break; + case 7: /* log Cauchy distribution */ + H=-pcauchy(ly,bt,lambda,0,1); break; + case 8: /* log Laplace distribution */ + tmp=exp(-fabs(ly-bt)/lambda)/2; + plap=ly<bt?tmp:1-tmp; + H=-log(1-plap); + break;}} + else{ + /* density models */ + switch(*model){ + case 1: /* exponential distribution */ + H=pexp(yy,bt,1,0); break; + case 2: /* Weibull distribution */ + H=pweibull(yy,lambda,bt,1,0); break; + case 3: /* gamma distribution */ + H=pgamma(yy,lambda,bt,1,0); break; + case 4: /* generalized logistic distribution */ + H=exp(-yy/lambda)*pow((lambda+intercept)/(lambda+intercept*exp(-bt*yy)),1/(lambda*bt)); + break; + case 5: /* log normal distribution */ + H=pnorm(ly,bt,lambda,1,0); break; + case 6: /* log logistic distribution */ + H=plogis(ly,bt,lambda,1,0); break; + case 7: /* log Cauchy distribution */ + H=pcauchy(ly,bt,lambda,1,0); break; + case 8: /* log Laplace distribution */ + tmp=exp(-fabs(ly-bt)/lambda)/2; + H=ly<bt?tmp:1-tmp; + break;}} + bb=H; + H-=bb0; + bb0=bb;} + else { + /* if there are time-varying covariates, finesse the problem + by integrating over the interval, but with covariate values + from the end of the interval */ + if(!*density){ + /* intensity models */ + switch(*model){ + case 1: /* exponential distribution */ + H=(yy-yy0)/bt; break; + case 2: /* Weibull distribution */ + H=pow(yy/bt,lambda)-pow(yy0/bt,lambda); break; + case 3: /* gamma distribution */ + H=-pgamma(yy,lambda,bt,0,1)+pgamma(yy0,lambda,bt,0,1); break; + case 4: /* generalized logistic distribution */ + H=(yy-yy0+(log(lambda+intercept*exp(-bt*yy))-log(lambda+intercept*exp(-bt*yy0)))/bt)/lambda; + break; + case 5: /* log normal distribution */ + H=-pnorm(ly,bt,lambda,0,1)+(yy0>0?pnorm(ly0,bt,lambda,0,1):0); + break; + case 6: /* log logistic distribution */ + H=-plogis(ly,bt,lambda,0,1)+(yy0>0?plogis(ly0,bt,lambda,0,1):0); + break; + case 7: /* log Cauchy distribution */ + H=-pcauchy(ly,bt,lambda,0,1)+(yy0>0?pcauchy(ly0,bt,lambda,0,1):0); + break; + case 8: /* log Laplace distribution */ + tmp=exp(-fabs(ly-bt)/lambda)/2; + plap=ly<bt?tmp:1-tmp; + H=-log(1-plap); + if(yy0>0){ + tmp=exp(-fabs(ly0-bt)/lambda)/2; + plap=ly0<bt?tmp:1-tmp; + H+=log(1-plap);} + break;}} + else{ + /* density models */ + switch(*model){ + case 1: /* exponential distribution */ + H=pexp(yy-yy0,bt,1,0); break; + case 2: /* Weibull distribution */ + H=pweibull(yy,lambda,bt,1,0)-pweibull(yy0,lambda,bt,1,0); break; + case 3: /* gamma distribution */ + H=pgamma(yy,lambda,bt,1,0)-pgamma(yy0,lambda,bt,1,0); break; + case 4: /* generalized logistic distribution */ + H=exp(-yy/lambda)*pow((lambda+intercept)/(lambda+intercept*exp(-bt*yy)),1/(lambda*bt))-exp(-yy0/lambda)*pow((lambda+intercept)/(lambda+intercept*exp(-bt*yy0)),1/(lambda*bt)); + break; + case 5: /* log normal distribution */ + H=pnorm(ly,bt,lambda,1,0)-(yy0>0?pnorm(ly0,bt,lambda,1,0):0); break; + case 6: /* log logistic distribution */ + H=plogis(ly,bt,lambda,1,0)-(yy0>0?plogis(ly0,bt,lambda,1,0):0); + break; + case 7: /* log Cauchy distribution */ + H=pcauchy(ly,bt,lambda,1,0)-(yy0>0?pcauchy(ly0,bt,lambda,1,0):0); + break; + case 8: /* log Laplace distribution */ + tmp=exp(-fabs(ly-bt)/lambda)/2; + H=ly<bt?tmp:1-tmp; + if(yy0>0)H-=exp(-fabs(ly0-bt)/lambda)/2; + break;}} + ly0=ly;} + if(j>-1){ + b1+=H; + /* calculate likelihood */ + if(*pfamily)*like-=dpvfp2(c[nm],a,b,b1,family); + else *like-=lgammafn(a1)-lgammafn(a)+a*log(b)-a1*log(b1); + if(c[nm]>0)*like-=c[nm]*log(H); + if(c[nm]>1)*like+=lgammafn(c[nm]+1); + /* calculate fitted values */ + if(*fit){ + if(!*density){ + switch(*model){ + case 1: pred[nm]=1/bt; break; + case 2: pred[nm]=lambda*pow(yy/bt,lambda-1)/bt; break; + case 3: pred[nm]=dgamma(yy,lambda,bt,0)/(1-pgamma(yy,lambda,bt,1,0)); break; + case 4: pred[nm]=1/(lambda+intercept*exp(-bt*yy)); break; + case 5: pred[nm]=dnorm(ly,bt,lambda,0)/(y[nm]*(1-pnorm(ly,bt,lambda,1,0))); + break; + case 6: + pred[nm]=dlogis(ly,bt,lambda,0)/(y[nm]*(1-plogis(ly,bt,lambda,1,0))); + break; + case 7: + pred[nm]=dcauchy(ly,bt,lambda,0)/(y[nm]*(1-pcauchy(ly,bt,lambda,1,0))); + break; + case 8: + tmp=exp(-fabs(ly-bt)/lambda)/2; + plap=ly<bt?tmp:1-tmp; + pred[nm]=tmp/(lambda*y[nm]*(1-plap)); + break;}} + else{ + switch(*model){ + case 1: pred[nm]=dexp(yy,bt,0); break; + case 2: pred[nm]=dweibull(yy,lambda,bt,0); break; + case 3: pred[nm]=dgamma(yy,lambda,bt,0); break; + case 4: + pred[nm]=exp(-yy/lambda)*pow((lambda+intercept)/(lambda+intercept*exp(-bt*yy)),1/(lambda*bt)+1); + break; + case 5: pred[nm]=dnorm(ly,bt,lambda,0)/y[nm]; break; + case 6: pred[nm]=dlogis(ly,bt,lambda,0)/y[nm]; break; + case 7: pred[nm]=dcauchy(ly,bt,lambda,0)/y[nm]; break; + case 8: + tmp=exp(-fabs(ly-bt)/lambda)/2; + pred[nm]=tmp/(lambda*y[nm]); + break;}} + rpred[nm]=j==0?pred[nm]:a*H/b;} + /* update parameters */ + switch(*dep){ + case 1: a=omega*a1+(1-omega)*delta; break; + case 2: + case 6: + om=pow(omega,yy-yy0); + a=om*a1+(1-om)*delta; + break; + case 3: a=a1; break; + case 4: + case 5: a=omega*a1; break; + default: ; + } + switch(*dep){ + case 1: b=omega*b1+(1-omega)*delta; break; + case 2: b=om*b1+(1-om)*delta; break; + case 3: + case 5: b=omega*b1; break; + case 6: b=om*(b1-b)+delta; break; + default: ; + } + nm++;} + yy0=yy;}} + return;} + +void countfb(double p[],double y[],int c[], double x[],int *nind, + int nobs[],int *nbs,int *nccov,int *model,int *density, + int *tvc,double tvcov[],int *fit,double pred[],double rpred[],int *rf, + double bbb[], int *sf, double vv[], int *frser, double *like){ + int i,j,k,nm; + double a,b,bb,bb0,delta,lambda,omega,om,beta,bt,H,yy,kk,tmp,ly, + plap,intercept,res; + + *like=0; + nm=0; + delta=exp(p[*nccov+*tvc+1]); + if(*model>1&&!*sf){ + if(*model<5)lambda=exp(p[*nccov+*tvc+*frser+2]); + else lambda=exp(p[*nccov+*tvc+*frser+2]/2);} + if(*model==4)intercept=exp(p[*nccov+*tvc+*frser+3]); + for(i=0;i<*nind;i++){ + a=delta; + b=bb0=0; + if(!*rf){ + beta=p[0]; + for(k=0;k<*nccov;k++)beta+=p[k+1]*x[i+k**nind]; + if(*model<4){ + if(beta>40) beta=40; + if(beta<-40)beta=-40; + beta=exp(beta);}} + else if(!*tvc)bt=bbb[i]; + res=0; + for(j=0;j<nobs[i];j++){ + yy=y[nm]; + if(*model>=5)ly=log(yy); + if(*model>1&&*sf)lambda=vv[nm]; + a+=c[nm]; + /* add in time-varying covariates */ + if(!*rf){ + if(*tvc){ + bt=0; + for(k=0;k<*tvc;k++)bt+=p[*nccov+k+1]*tvcov[nm+*nbs*k]; + if(*model<4)bt=exp(bt)*beta; + else bt+=beta;} + else bt=beta;} + else if(*tvc)bt=bbb[nm]; + /* if AR, add discounted previous residual */ + if(*frser&&j>0){ + if(*model<4)bt=exp(log(bt)+exp(p[*nccov+*tvc+2]*(y[nm]-y[nm-1]))*res); + else bt+=exp(p[*nccov+*tvc+2]*(y[nm]-y[nm-1]))*res;} + if(!*density){ + /* intensity models */ + switch(*model){ + case 1: /* exponential distribution */ + H=yy/bt; break; + case 2: /* Weibull distribution */ + H=pow(yy/bt,lambda); break; + case 3: /* gamma distribution */ + H=-pgamma(yy,lambda,bt,0,1); break; + case 4: /* generalized logistic distribution */ + H=(yy+log(lambda+intercept*exp(-bt*yy))/bt)/lambda; break; + case 5: /* log normal distribution */ + H=-pnorm(ly,bt,lambda,0,1); break; + case 6: /* log logistic distribution */ + H=-plogis(ly,bt,lambda,0,1); break; + case 7: /* log Cauchy distribution */ + H=-pcauchy(ly,bt,lambda,0,1); break; + case 8: /* log Laplace distribution */ + tmp=exp(-fabs(ly-bt)/lambda)/2; + plap=ly<bt?tmp:1-tmp; + H=-log(1-plap); + break;}} + else{ + /* density models */ + switch(*model){ + case 1: /* exponential distribution */ + H=pexp(yy,bt,1,0); break; + case 2: /* Weibull distribution */ + H=pweibull(yy,lambda,bt,1,0); break; + case 3: /* gamma distribution */ + H=pgamma(yy,lambda,bt,1,0); break; + case 4: /* generalized logistic distribution */ + H=exp(-yy/lambda)*pow((lambda+intercept)/(lambda+intercept*exp(-bt*yy)),1/(lambda*bt)); + break; + case 5: /* log normal distribution */ + H=pnorm(ly,bt,lambda,1,0); break; + case 6: /* log logistic distribution */ + H=plogis(ly,bt,lambda,1,0); break; + case 7: /* log Cauchy distribution */ + H=pcauchy(ly,bt,lambda,1,0); break; + case 8: /* log Laplace distribution */ + tmp=exp(-fabs(ly-bt)/lambda)/2; + H=ly<bt?tmp:1-tmp; + break;}} + bb=H; + H-=bb0; + bb0=bb; + b+=H; + /* if(*tvc)b+=H; + else if(j==nobs[i]-1)b=bb; */ + /* calculate likelihood */ + if(c[nm]>0)*like-=c[nm]*log(H)-lgammafn(c[nm]+1); + /* calculate fitted values */ + if(*fit||*frser){ + if(!*density){ + switch(*model){ + case 1: pred[nm]=1/bt; break; + case 2: pred[nm]=lambda*pow(yy/bt,lambda-1)/bt; break; + case 3: pred[nm]=dgamma(yy,lambda,bt,0)/(1-pgamma(yy,lambda,bt,1,0)); break; + case 4: pred[nm]=1/(lambda+intercept*exp(-bt*yy)); break; + case 5: pred[nm]=dnorm(ly,bt,lambda,0)/(y[nm]*(1-pnorm(ly,bt,lambda,1,0))); + break; + case 6: + pred[nm]=dlogis(ly,bt,lambda,0)/(y[nm]*(1-plogis(ly,bt,lambda,1,0))); + break; + case 7: + pred[nm]=dcauchy(ly,bt,lambda,0)/(y[nm]*(1-pcauchy(ly,bt,lambda,1,0))); + break; + case 8: + tmp=exp(-fabs(ly-bt)/lambda)/2; + plap=ly<bt?tmp:1-tmp; + pred[nm]=tmp/(lambda*y[nm]*(1-plap)); + break;}} + else{ + switch(*model){ + case 1: pred[nm]=dexp(yy,bt,0); break; + case 2: pred[nm]=dweibull(yy,lambda,bt,0); break; + case 3: pred[nm]=dgamma(yy,lambda,bt,0); break; + case 4: + pred[nm]=exp(-yy/lambda)*pow((lambda+intercept)/(lambda+intercept*exp(-bt*yy)),1/(lambda*bt)+1); + break; + case 5: pred[nm]=dnorm(ly,bt,lambda,0)/y[nm]; break; + case 6: pred[nm]=dlogis(ly,bt,lambda,0)/y[nm]; break; + case 7: pred[nm]=dcauchy(ly,bt,lambda,0)/y[nm]; break; + case 8: + tmp=exp(-fabs(ly-bt)/lambda)/2; + pred[nm]=tmp/(lambda*y[nm]); + break;}} + rpred[nm]=a*H/(delta+b); + if(*frser)res=(*model>=4?ly:y[nm])-pred[nm];} + nm++;} + *like-=lgammafn(a)-a*log(delta+b);} + *like-=*nind*(-lgammafn(delta)+delta*log(delta)); + return;} + +/* power variance function Poisson */ + +static double pvfc2(int y, double d, double s1, double f){ + int i,j; + double r,*c,tmp1,tmp2,tmp3; + c=(double*)R_alloc((size_t)(y*y),sizeof(double)); + tmp1=gammafn(1.-f); + tmp2=log(d); + tmp3=log(s1); + for(i=0;i<y;i++){ + c[i*(y+1)]=1; + if(i>0){ + c[i*y]=gammafn(i+1-f)/tmp1; + if(i>1) + for(j=1;j<i;j++)c[y*i+j]=c[y*(i-1)+j-1]+c[y*(i-1)+j]*(i-(j+1)*f);}} + r=0.; + for(i=1;i<=y;i++){ + r+=c[y*(y-1)+i-1]*exp(i*tmp2+(i*f-y)*tmp3);} + return(log(r));} + +static double dpvfp2(int y, double d, double s, double s1, double f){ + double res; + if(f==0.)res=dnbinom(y,d,s/s1,1); + else { + res=-d*(pow(s1,f)-pow(s,f))/f; + if(y>0)res+=pvfc2(y,d,s1,f);} + return(res);}
signature.asc
Description: Digital signature