TABLE OF CONTENTS

CMetropolisHastings/MH_SplineHazard [ Functions ]

NAME

    MH_SplineHazard --- MH for hazard spline parameters

FUNCTION

Update the parameters of the hazard spline component by Metropolis-Hastings. See also mh.hazard.spline.

SYNOPSIS

2069 void MH_SplineHazard(curveP hazard, curveP frailty, regressionP regression)

INPUTS

    hazard        CCurve for the hazard
    frailty       CCurve for the frailty
    regression    CRegression structure

SOURCE

2073 {
2074     if(!hazard->hasSpline) return;
2075     double baselik, candlik;
2076     double sumacc=0;
2077     // base likelihood
2078     baselik = LikelihoodSplineHazard(hazard,frailty,regression); 
2079     double * cand = (double *) calloc( hazard->nj, sizeof(double));
2080     double * thiscand = (double *) calloc( hazard->nj, sizeof(double));
2081     // allocate storage for parameters of old spline
2082     // if the move is not accepted, these will be used to restore it, to make it faster
2083     double * oldPar = (double *) malloc( hazard->nj * sizeof(double));
2084     double * oldEPar = (double *) malloc( hazard->nj * sizeof(double));
2085     double * oldY = (double *) malloc( hazard->nx * sizeof(double));
2086     double * oldYcum = (double *) malloc( hazard->nx * sizeof(double));
2087     double * oldSplineY = (double *) malloc( hazard->nx * sizeof(double));
2088     double * oldSplineYcum = (double *) malloc( hazard->nx * sizeof(double));
2089     dcopyWrapper(hazard->nj, hazard->SplinePar, oldPar);
2090     dcopyWrapper(hazard->nj, hazard->SplineEPar, oldEPar);
2091     dcopyWrapper(hazard->nj, hazard->SplinePar, thiscand);
2092     dcopyWrapper(hazard->nx, hazard->Y, oldY);
2093     dcopyWrapper(hazard->nx, hazard->SplineY, oldSplineY);
2094     dcopyWrapper(hazard->nx, hazard->Ycum, oldYcum);
2095     dcopyWrapper(hazard->nx, hazard->SplineYcum, oldSplineYcum);
2096     // create candidate parameters
2097     for(int j=0;j<hazard->nj;j++)
2098         cand[j] = hazard->SplinePar[j]+hazard->SplineTun[0]*
2099             rnorm(0,hazard->SplineCandSD[j]);
2100     // update spline parameters one at a time
2101     for(int j=0; j < hazard->nj; j++){
2102         thiscand[j] = cand[j];
2103         UpdateSplinePar(hazard,thiscand,j); 
2104         // candidate likelihod
2105         candlik = LikelihoodSplineHazard(hazard,frailty,regression);
2106         int acc = AcceptReject(baselik, candlik, 1);
2107         if(acc){ // accepted, so store the new information
2108             baselik = candlik;
2109             sumacc++;
2110             dcopyWrapper(hazard->nj, hazard->SplineEPar, oldEPar);
2111             dcopyWrapper(hazard->nx, hazard->Y, oldY);
2112             dcopyWrapper(hazard->nx, hazard->SplineY, oldSplineY);
2113             dcopyWrapper(hazard->nx, hazard->Ycum, oldYcum);
2114             dcopyWrapper(hazard->nx, hazard->SplineYcum, oldSplineYcum);
2115         }else{ // rejected, so restore the old information
2116             thiscand[j]=oldPar[j];
2117             dcopyWrapper(hazard->nj, thiscand, hazard->SplinePar);
2118             dcopyWrapper(hazard->nj, oldEPar, hazard->SplineEPar);
2119             dcopyWrapper(hazard->nx, oldY, hazard->Y);
2120             dcopyWrapper(hazard->nx, oldSplineY, hazard->SplineY);
2121             dcopyWrapper(hazard->nx, oldYcum, hazard->Ycum);
2122             dcopyWrapper(hazard->nx, oldSplineYcum, hazard->SplineYcum);
2123         }
2124     }
2125     // acceptance rate
2126     hazard->SplineAccept[0] =  sumacc / ((double) hazard->nj);
2127     // free memory
2128     free(cand);
2129     free(thiscand);
2130     free(oldPar);
2131     free(oldEPar);
2132     free(oldY);
2133     free(oldYcum);
2134     free(oldSplineY);
2135     free(oldSplineYcum);
2136 }