Commit 7b584fca authored by Andreas Tille's avatar Andreas Tille

New upstream version 2.1.10

parent fe917fc7
Package: logspline
Version: 2.1.9
Date: 2016-02-01
Title: Logspline Density Estimation Routines
Version: 2.1.10
Date: 2018-06-01
Title: Routines for Logspline Density Estimation
Author: Charles Kooperberg <clk@fredhutch.org>
Maintainer: Charles Kooperberg <clk@fredhutch.org>
Description: Routines for the logspline density estimation. oldlogspline()
uses the same algorithm as the logspline 1.0.x package - the Kooperberg
and Stone (1992) <DOI: 10.2307/1390786> algorithm (with an improved interface).
The recommended routine logspline() uses an algorithm from Stone et al (1997).
<DOI: 10.1214/aos/1031594728>
Description: Contains routines for logspline density estimation.
The function oldlogspline() uses the same algorithm as the logspline package
version 1.0.x; i.e. the Kooperberg and Stone (1992) <DOI:10.2307/1390786>
algorithm (with an improved interface). The recommended routine logspline()
uses an algorithm from Stone et al (1997) <DOI:10.1214/aos/1031594728>.
Imports: stats, graphics
License: Apache License 2.0
Packaged: 2016-02-02 17:55:13 UTC; clk
NeedsCompilation: yes
Packaged: 2018-06-01 23:09:26 UTC; clk
Repository: CRAN
Date/Publication: 2016-02-03 01:15:45
Date/Publication: 2018-06-02 04:58:19 UTC
dlogspline Logspline Density Estimation
doldlogspline Logspline Density Estimation - 1992 version
logspline Logspline Density Estimation
oldlogspline Logspline Density Estimation - 1992 version
plot.logspline Logspline Density Estimation
plot.oldlogspline Logspline Density Estimation - 1992 version
summary.logspline Logspline Density Estimation
summary.oldlogspline Logspline Density Estimation - 1992 version
unstrip Reformat data as vector or matrix
4dd8e8f0d49a03888c4135fa3883adf3 *DESCRIPTION
4a653206f0b058e399f40562e9449249 *INDEX
9ed89ae264cb6adf6563a680edbb7431 *NAMESPACE
058393f025c217ad9deb7a05afe38747 *R/logspline.R
e005c87432b1a0afdd14162adcb4e956 *DESCRIPTION
73d98ba76a12b8b7aed3bb0ac4d6183c *NAMESPACE
09b85350d31e7ea491ce12e0af07b944 *R/logspline.R
e5937390332be13e210110b2f5030d8d *man/dlogspline.Rd
1905f0e1bc2347092ee3a4d6fc806b08 *man/doldlogspline.Rd
e8f70418b42dc00732c0bba662e65537 *man/logspline-internal.Rd
0e63737114e136fef7ed9e2b893d8372 *man/logspline.Rd
85aef4792b0f000ec0cf45e3e64e700d *man/oldlogspline.Rd
08c579c346e1f7b79c687bf6dc674269 *man/oldlogspline.to.logspline.Rd
......@@ -14,6 +14,7 @@ eabaae5f996c1402466f5a54717a3121 *man/summary.logspline.Rd
d9d2d511c4cec309ecc1c2f05d87b121 *man/unstrip.Rd
8290d2e9740414e315237f0d5d4024bb *src/Makevars
157084291a6fa50c11e5d7ae2325507f *src/allpack.f
08ec5922327358575ea46e2fe28e4d3f *src/lsdall.c
baecd678233af49d9e3174f53caab65b *src/nlsd.c
51b9bf35a21a9698214eb75ee5033258 *src/lsdall.c
6df20d87970ec1b9a5eef7f70933c11e *src/nlsd.c
c625d1f0667c03581e95c09fdf044aae *src/registerDynamicSymbol.c
782c6ba6b56e9842d5854775ce3653e3 *src/x2c.h
# Remove the previous line if you edit this file
useDynLib(logspline)
useDynLib(logspline, .registration = TRUE)
# Export all names
exportPattern(".")
......
#
# Copyright [1993-2016] [Charles Kooperberg]
# Copyright [1993-2018] [Charles Kooperberg]
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
......@@ -353,6 +353,7 @@ oldlogspline <- function(uncensored, right, left, interval, lbound, ubound,
SorC <- vector(mode = "integer", length = 35)
SorC[1] <- 1 # the actual function call
nsample[6] <- nsample[6]-1
if(length(table(sample))<3)stop("Not enough unique values")
z <- .C("logcensor",
as.integer(delete),
as.integer(iautoknot),
......@@ -438,6 +439,7 @@ logspline <- function(x, lbound, ubound, maxknots=0, knots, nknots=0,
call <- match.call()
if(!missing(x))x <- unstrip(x)
data <- x
if(length(table(data))<3)stop("Not enough unique values")
ilx <- 0; iux <- 0
if(!missing(lbound)){ilx <- 1;jlx <- lbound}
if(!missing(ubound)){iux <- 1;jux <- ubound}
......
\name{logspline-internal}
\title{Internal glmnet functions}
\alias{logcensor}
\alias{nlogcensor}
\alias{nlogcensorx}
\alias{pqlsd}
\alias{rpqlsd}
\description{Internal logspline functions}
\author{Charles Koopeeberg}
\details{These are not intended for use by users.}
\keyword{internal}
/*
*
* Copyright [1993-2016] [Charles Kooperberg]
* Copyright [1993-2018] [Charles Kooperberg]
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
......@@ -41,7 +41,7 @@ static double pqexp(),pqnum(),lpqexpi(),pqdens();
/* this is the main program */
/* remove follows at the end */
int logcensor(idelete,iknotauto,sample,nsample,bound,
void logcensor(idelete,iknotauto,sample,nsample,bound,
SorC,ynknots,yknots,ycoef,alpha,wk,wk2,logl)
int *idelete,*iknotauto,nsample[],SorC[],*ynknots;
......@@ -73,11 +73,10 @@ double aic,aicmin,r1,rknots[NC],xcoef2[NC][NC];
alpha - alpha value in aic
rknots - copy of knots */
int i,j,nkstart,iremove=0,iknots2[NC],iknots[NC],xiknots[NC];
int i,j,nkstart,iremove=0,iknots[NC],xiknots[NC];
/* local integers
i,j,k - counter, utility
nkstart - number of knots at the beginning of the algorithm
iknots2 - copy of iknots
iremove - number of the knot that is removed */
......@@ -96,7 +95,7 @@ int i,j,nkstart,iremove=0,iknots2[NC],iknots[NC],xiknots[NC];
(void)Rprintf("sample is too small\n");
else
SorC[0] = 2;
return 0;
return;
}
/* determine the number of starting knots */
......@@ -129,9 +128,9 @@ int i,j,nkstart,iremove=0,iknots2[NC],iknots[NC],xiknots[NC];
for(i=0;i<nknots;i++)yknots[i]=knots[i];
/* some possible errors */
if(SorC[0] == -2)return 0;
if(SorC[0] == 647)return 0;
if(SorC[0] == 23)return 0;
if(SorC[0] == -2)return;
if(SorC[0] == 647)return;
if(SorC[0] == 23)return;
/* Compute stuff later used for computing sufficient statistics. */
suffstat1(suffcombine,sample,nsample);
......@@ -172,7 +171,7 @@ int i,j,nkstart,iremove=0,iknots2[NC],iknots[NC],xiknots[NC];
if(nknots==nkstart){
SorC[0] = -1;
SorC[27]=0;
return 0;
return;
}
if(SorC[0]==-647)
(void)Rprintf("Smallest number of knots tried is %d\n",nknots+1);
......@@ -201,22 +200,19 @@ int i,j,nkstart,iremove=0,iknots2[NC],iknots[NC],xiknots[NC];
logl[nknots-1]=loglikelihood*nsample[0]+log(qt[1])*nsample[1]; /* &&&&& */
if(aic <= aicmin){
/* Then we mark the solution iknots2 is put to 0, so we can fill it later */
/* Then we mark the solution */
aicmin = aic;
xczheta = czheta;
xnknots = nknots;
for(i=0; i<NC; i++){
xzheta[i] = zheta[i];
xiknots[i] = iknots[i];
iknots2[i] = 0;
for(j=0;j<NC;j++)
xcoef2[i][j]=coef2[i][j];
}
/* that is we fill it here, iknots contains the rank numbers of the remaining
knots */
for(i=0;i<nknots;i++)
iknots2[iknots[i]] = 1;
}
else{
r1 = -2. * loglikelihood * nsample[0] + *alpha * 2;
......@@ -270,7 +266,7 @@ int i,j,nkstart,iremove=0,iknots2[NC],iknots[NC],xiknots[NC];
ycoef[0]=ycoef[0]+ycoef[1]*qt[0]+log(qt[1]);
ycoef[1]=ycoef[1]*qt[1];
return nkstart;
return;
}
/******************************************************************************/
......@@ -891,7 +887,6 @@ double candidate[],sample[],bound[];
int nsample[],accuracy;
{
double r0,r1,likl,r3[NC+1],aa[6],bb[6];
double rtt;
int i1,i2,i3,i4,iv,iw;
r0=exp((double)(-740));
......@@ -955,7 +950,6 @@ int nsample[],accuracy;
}
/* the right censored data */
for(i1=0;i1<nsample[3];i1++){
rtt=likl;
i2=i1+nsample[1]+2*nsample[2];
/* in which interval is the datapoint */
for(i3=0;knots[i3]<sample[i2] && i3 < nknots;i3++);
......@@ -1888,7 +1882,7 @@ double rknots[],sample[],bound[],smp2[],smp3[],qt[];
/* these quantities are defined in lhead.h and the files where they originate */
{
int i,j=0,j2,k,kk,ll,ia=0,il;
int i,j=0,j2,k,kk,ll,ia=0,il,kx;
/* local integers
i k - counters
j j2 - is there an odd or an even number of knots? */
......@@ -2075,7 +2069,8 @@ double rknots[],sample[],bound[],smp2[],smp3[],qt[];
k = 0;
for(i=1;i<nknots-1;i++){
rknots[i]=(rknots[i]-0.5);
for(k=k;k<=il;k++){
kx=k;
for(k=kx;k<=il;k++){
if(smp3[k]>=rknots[i]){
knots[i]=((rknots[i]-smp3[k-1])*smp2[k]+
(smp3[k]-rknots[i])*smp2[k-1])/(smp3[k]-smp3[k-1]);
......@@ -2889,6 +2884,8 @@ int *ipq,*lk,*lp;
{
double v1[2],v2[2];
int ij;
v2[0]=0;
v2[1]=0;
if((*ipq)==1)
qtop(coef,knots,bound,pp,qq,*lp,*lk);
else{
......@@ -2907,6 +2904,7 @@ int lp,lk;
{
double l0,l1,r0,r1, s1,s2,s3,s4,x1,x2,x3,sj,xj;
int i,j,vr,vl,k;
j=0;
l0 = coef[0];
l1 = coef[1];
r0 = l0;
......
/*
*
* Copyright [1993-2016] [Charles Kooperberg]
* Copyright [1993-2018] [Charles Kooperberg]
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
......@@ -1720,6 +1720,7 @@ struct datastruct *data;
double a,b;
struct basisfunct *bn;
/* no deep meanings */
m=0;
/* get (*spc).basis.c3 */
for(i=0;i<(*spc).ndim;i++){
bn=&((*spc).basis[i]);
......@@ -2055,7 +2056,7 @@ int *ipq,*lk,*lp;
{
double *kpl,**cpl,*ppl,*pqx,r;
double *zz,cor;
int i,j,nk,fst,lst;
int i,j,nk,fst,lst,jx;
/* Gaussian quadrature coefficients */
ww6[1 ]= 0.467913934572691; yy6[1 ]= 0.238619186083197;
ww6[2 ]= 0.360761573048139; yy6[2 ]= 0.661209386466265;
......@@ -2197,7 +2198,8 @@ int *ipq,*lk,*lp;
/* before the first knot */
fst=j;
lst=j-1;
for(j=j;j<(*lp) && pq[j]<=zz[1];j++) lst=j;
jx=j;
for(j=jx;j<(*lp) && pq[j]<=zz[1];j++) lst=j;
if(lst>=fst){
if((*ipq)==0) getq0(pq,pqx,fst,lst,cpl[0],kpl[0],bnd[0],cor);
else getp0(pq,pqx,fst,lst,cpl[0],kpl[0],bnd[0],cor);
......@@ -2206,7 +2208,8 @@ int *ipq,*lk,*lp;
for(i=1;i<nk-1;i++){
fst=j;
lst=j-1;
for(j=j;j<(*lp) && pq[j]<=zz[i+1];j++) lst=j;
jx=j;
for(j=jx;j<(*lp) && pq[j]<=zz[i+1];j++) lst=j;
if(lst>=fst){
if((*ipq)==0)
getq1(pq,pqx,fst,lst,cpl[i],kpl[i],kpl[i+1],ppl[i],ppl[i+1]);
......@@ -2216,13 +2219,15 @@ int *ipq,*lk,*lp;
/* beyond the larst knot */
fst=j;
lst=j-1;
for(j=j;j<(*lp) && pq[j]<zz[nk];j++) lst=j;
jx=j;
for(j=jx;j<(*lp) && pq[j]<zz[nk];j++) lst=j;
if(lst>=fst){
if((*ipq)==0) getq2(pq,pqx,fst,lst,cpl[nk-1],kpl[nk],bnd[2],cor);
else getp2(pq,pqx,fst,lst,cpl[nk-1],kpl[nk],bnd[2],cor);
}
/* outside the range */
for(j=j;j<(*lp);j++){
jx=j;
for(j=jx;j<(*lp);j++){
if((*ipq)==0){
if(bnd[2]>0.5)pqx[j]=kpl[nk];
else pqx[j]= 1.0e100;
......@@ -2295,6 +2300,8 @@ int f,l;
for(i=1;i<=50;i++)y[i]=y[i-1]+r*(f1[2*(i-1)]+4*f1[2*i-1]+f1[2*i])/3.;
s=(p1-p0)/y[50];
for(i=0;i<=50;i++)y[i]=p0+y[i]*s;
y[0]=p0;
y[50]=p1;
i=0;
s=2.*r;
for(j=f;j<=l;j++){
......@@ -2303,7 +2310,7 @@ int f,l;
if(p[j]>=y[i] && p[j]<=y[i+1]) q[j]=k0+s*i+s*(p[j]-y[i])/(y[i+1]-y[i]);
else i++;
}
while (q[j]<k0);
while ((q[j]<k0)&(i<50));
}
}
/******************************************************************************/
......
#include <stdlib.h> // for NULL
#include <R_ext/Rdynload.h>
/* FIXME:
Check these declarations against the C/Fortran source code.
*/
/* .C calls */
extern void logcensor(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *);
extern void nlogcensor(void *, void *, void *, void *, void *, void *);
extern void nlogcensorx(void *);
extern void pqlsd(void *, void *, void *, void *, void *, void *, void *, void *);
extern void rpqlsd(void *, void *, void *, void *, void *, void *, void *);
static const R_CMethodDef CEntries[] = {
{"logcensor", (DL_FUNC) &logcensor, 13},
{"nlogcensor", (DL_FUNC) &nlogcensor, 6},
{"nlogcensorx", (DL_FUNC) &nlogcensorx, 1},
{"pqlsd", (DL_FUNC) &pqlsd, 8},
{"rpqlsd", (DL_FUNC) &rpqlsd, 7},
{NULL, NULL, 0}
};
void R_init_logspline(DllInfo *dll)
{
R_registerRoutines(dll, CEntries, NULL, NULL, NULL);
R_useDynamicSymbols(dll, FALSE);
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment