Revision 616fe292ea4c66f8db35d04c9fa045fea56eea02 authored by Roger Koenker on 04 September 2016, 13:02 UTC, committed by cran-robot on 04 September 2016, 13:02 UTC
1 parent db6345a
mcmb.c
``````
/****************************************************************/
/*  Markov Chain Marginal Bootstrap Package for Quantile Regression   */
/*  Maria Kocherginsky <mkocherg@uchicago.edu>,Xuming He <he@uiuc.edu>*/
/*                                                                    */
/*  This code modified by R. Koenker (July 30, 2003, April 21 2010)   */
/*  to employ R random number generation conventions and to conform   */
/*  to the conventions of the quantreg package.  See the function     */
/*  boot.rq.mcmb for details on the interface to summary.rq           */
/**********************************************************************/
#include <stdio.h>
#include <math.h>
#include <sys/types.h>
#include <stdlib.h>
#include <time.h>
#include <R.h>

#define MAXN 200000
#define MAXP 100
#define XACCF .001
#define XACCD .001

#define MAXIT 100
#define BMAXIT 100

#define UNUSED (-1.11e30)
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
#define M 7
#define NSTACK 50

int allZero;

/****************************************************************/
/*  Error Handling      */
/****************************************************************/
void error(const char * format, ...);
/****************************************************************/
/*  Various Utilities off the Web (for use with sort2())        */
/****************************************************************/

unsigned int *lvector(int  nl, int  nh)
/* allocate a int vector with subscript range v[nl..nh] */
{
unsigned int *v = (unsigned int *)malloc((size_t)((nh-nl+2) * sizeof(int)));
if (v == NULL) error("allocation failure in lvector()");
return (v-nl+1);
}

void free_lvector(unsigned int *v, int nl, int nh)
/* free a int int vector allocated with lvector() */
{
free(v+nl-1);
}

/****************************************************************/
/*  Multiply two vectors                                        */
/****************************************************************/

double mprodx(double *x, double *c, int pp){
int i;
double sum;

sum=0;
for(i=0;i<pp;i++){
sum=sum+x[i]*c[i];}

return sum;
}

/****************************************************************/
/*  Sign Function                                               */
/****************************************************************/

double sign(double x){
double sign;

if(x>0.0)
sign=1;
else if (x<0.0)
sign=-1;
else
sign=0;

return sign;
}

/****************************************************************/
/*  Sort Function                                               */
/****************************************************************/

void sort2(unsigned int  n, double arr[], double brr[]){
unsigned int  i,ir=n,j,k,l=1,*istack;
int jstack=0;
double a,b,temp;
istack=lvector(1,NSTACK);

for (;;) {
if (ir-l < M) {
for (j=l+1;j<=ir;j++) {
a=arr[j];
b=brr[j];
for (i=j-1;i>=l;i--) {
if (arr[i] <= a) break;
arr[i+1]=arr[i];
brr[i+1]=brr[i];
}
arr[i+1]=a;
brr[i+1]=b;
}
if (!jstack) {
free_lvector(istack,1,NSTACK);
return;
}
ir=istack[jstack];
l=istack[jstack-1];
jstack -= 2;
} else {
k=(l+ir) >> 1;
SWAP(arr[k],arr[l+1])
SWAP(brr[k],brr[l+1])
if (arr[l] > arr[ir]) {
SWAP(arr[l],arr[ir])
SWAP(brr[l],brr[ir])
}
if (arr[l+1] > arr[ir]) {
SWAP(arr[l+1],arr[ir])
SWAP(brr[l+1],brr[ir])
}
if (arr[l] > arr[l+1]) {
SWAP(arr[l],arr[l+1])
SWAP(brr[l],brr[l+1])
}
i=l+1;
j=ir;
a=arr[l+1];
b=brr[l+1];
for (;;) {
do i++; while (arr[i] < a);
do j--; while (arr[j] > a);
if (j < i) break;
SWAP(arr[i],arr[j])
SWAP(brr[i],brr[j])
}
arr[l+1]=arr[j];
arr[j]=a;
brr[l+1]=brr[j];
brr[j]=b;
jstack += 2;
if (jstack > NSTACK) error("NSTACK too small in sort2.\n");
if (ir-i+1 >= j-l) {
istack[jstack]=ir;
istack[jstack-1]=i;
ir=j-1;
} else {
istack[jstack]=j-1;
istack[jstack-1]=l;
l=i;
}
}
}
}

/****************************************************************/
/* This function finds a weighted taustar quantile, instead of  */
/* using bisection to solve the MCMB equations. See new paper.  */
/****************************************************************/

double func(double *x, double *y,  double tau, double *tTilda, double *A,  double
sum_right, double sumxij, double sumabsxij, int j, int pp, int nn){
int i,m;
double *xj, *yj, *z, *wt, wtsum;
unsigned int mm;
double taustar, pwtsum, ans, large;

xj=(double *) calloc(nn+1, sizeof(double));
yj=(double *) calloc(nn+1, sizeof(double));
z=(double *) calloc(nn+2, sizeof(double));
wt=(double *) calloc(nn+2, sizeof(double));

for(i=0;i<nn;i++){
yj[i]=y[i];
xj[i]=x[i*pp+j];
}
xj[nn]=-sum_right/tau;
yj[nn]=10e16;
wtsum=sumabsxij + fabs(xj[nn]);

/*update the first n elements*/

/*sort function sorts arr[1...n], so set the first elements of z and wt to 0, so that the actual values in z and wt start from position z[1] and wt[1]*/
m=1;
z[0]=0; wt[0]=0;
for(i=0;i<nn;i++){
if(fabs(xj[i]) > 10e-16){
z[m]=(y[i]-mprodx(&x[i*pp],tTilda,pp)+tTilda[j]*xj[i])/xj[i];
wt[m]=fabs(xj[i])/wtsum;
m=m+1;
}
else{error("fabs(xj[i])<10e-16\n");}
}
z[m]=10e16*sign(xj[nn]);
large=z[m];
wt[m]=fabs(xj[nn])/wtsum;

/*calculate taustar*/
taustar=(tau-0.5)*(sumxij+xj[nn])/(wtsum)+0.5;

if(m==0){
error("Error: one design variable contains all 0s.\n");
allZero=1;}

mm=m;
sort2(m, z, wt);

pwtsum=0; i=1;
ans=z[1];
while((pwtsum<=taustar) & (i<nn+1)){
pwtsum=pwtsum+wt[i];
ans=z[i];
i++;
}

/*resample if pick the largest z*/
if(fabs(ans) > 10e15){
error("Picked infinity; need to resample\n");
return 1.0;
}

free(xj);
free(yj);
free(z);
free(wt);

return ans;

}

/****************************************************************/
/* RQMCMB                                                       */
/****************************************************************/

void bootnp(double *x, double *y,  double *tau, double *theta_tilda, double *A, double *zstar,
double *sumxij, double *sumabsxij, int *n, int *p, int *success, double *theta, int *MAXK, int *seed){

int i, j,  jj, k, nn, pp;
double sum, s[MAXP],tau2, tTilda[MAXP];
int rand_ind;
extern int allZero;

pp=(int) *p;
nn=(int) *n;
tau2=(double) *tau;
allZero=0;

for(i=0;i<pp;i++){
theta[i]=theta_tilda[i];
tTilda[i]=theta_tilda[i];
}

GetRNGstate();
success[0]=1;

for(k=0;k<*MAXK; k++){

/*bootstrap from x[i][j] individually for each j, return sum */
for(j=0;j<pp;j++){
sum=0;
for(i=0;i<nn;i++){
rand_ind=(int)(nn*j+nn*unif_rand());
sum=sum+zstar[rand_ind];
}
s[j]=sqrt(nn)/sqrt(nn-pp)* sum;
}
PutRNGstate();

for(j=0;j<pp;j++){
theta[(k+1)*pp+j]=func(x, y, tau2, tTilda, A, s[j], sumxij[j], sumabsxij[j], j, pp, nn);
if(allZero==1){
success[0]=0;
return;}
if(theta[(k+1)*pp+j]==1.0){
for(jj=0;jj<pp;jj++){
tTilda[jj]=theta[k*pp+jj];}
k=k-1;
break;
}
else {
tTilda[j]=theta[(k+1)*pp+j];
}
}
}

}

``````

Computing file changes ...