TABLE OF CONTENTS

CMetropolisHastings/MH_Weight [ Functions ]

NAME

    MH_Weight --- MH for weight of spline component

FUNCTION

Update the relative weight of the spline component in a curve with both spline and parametric components by Metropolis-Hastings. See also mh.weight.

SYNOPSIS

2316 void MH_Weight(curveP theCurve, curveP hazard, curveP frailty, regressionP regression)

INPUTS

    theCurve      CCurve to be updated, can be either hazard or frailty
    hazard        CCurve for the hazard
    frailty       CCurve for the frailty
    regression    CRegression structure

SOURCE

2320 {
2321     if(!theCurve->hasPar | !theCurve->hasSpline) return;
2322     // get the likelihood function, depending on whether the curve is the hazard or frailty
2323     double ( *likfun )(curveP, curveP, regressionP);
2324     if(theCurve->isHazard) likfun = &LikelihoodWeightHazard;
2325     else likfun = &LikelihoodWeightFrailty;  
2326     double oldW = theCurve->Weight[0];
2327     double w = dmin(dmax(oldW, .01), 0.99);
2328     double v = theCurve->WeightTun[0];
2329     double alpha = w*(w*(1-w)/v-1);
2330     double beta = (1-w)/w*alpha;
2331     // generate candidate as beta
2332     double cand = rbeta(alpha, beta);
2333     if(isnan(cand)){
2334         theCurve->WeightAccept[0]=0;
2335         return;
2336     }
2337     double alphac = cand*(cand*(1-cand)/v-1);
2338     double betac = (1-cand)/cand*alphac;
2339     //base likelihood
2340     double baselik = likfun(hazard, frailty, regression);
2341     theCurve->Weight[0] = cand;
2342     ReweightCurve(theCurve, -1);
2343     // candidate likelihood
2344     double candlik = likfun(hazard, frailty, regression);
2345     // transition rates
2346     double puc = dbeta(cand, alpha, beta, 0);
2347     double pcu = dbeta(w, alphac, betac, 0);
2348     int acc = AcceptReject(baselik, candlik, pcu/puc);
2349     if(!acc){
2350         theCurve->Weight[0] = oldW;
2351         ReweightCurve(theCurve, -1);
2352     }
2353     theCurve->WeightAccept[0] = (double) acc;
2354 }