TABLE OF CONTENTS
CMetropolisHastings/MH_BDM [ Functions ]
NAME
MH_BDM --- birth-death-move steps for spline knots
FUNCTION
Reversible-Jump MCMC steps for adding, deleting, or moving knots as part of the adaptive knot selection procedure. Can work with either the hazard or frailty curve, see also mh.bdm.
SYNOPSIS
2371 void MH_BDM(char which, curveP hazard, curveP frailty, regressionP regression)
INPUTS
which character, can be either 'h' or 'f', for hazard or frailty respectively
hazard CCurve for the hazard
frailty CCurve for the frailty
regression CRegression structure
SOURCE
2375 { 2376 curveP theCurve; 2377 if(which == 'h') theCurve = hazard; 2378 if(which == 'f') theCurve = frailty; 2379 if(!theCurve->hasSpline) return; 2380 // Pointers to useful components 2381 int ord = theCurve->SplineOrd; 2382 int nknots = theCurve->SplineNknots; 2383 double * candknots = theCurve->SplineCandKnots; 2384 int ncandknots = theCurve->SplineNCandKnots; 2385 double * occ = theCurve->SplineCandOcc; // occupied candidate knots 2386 int * occind = malloc(nknots * sizeof(int)); 2387 // allocate temporary storage 2388 double * oldKnots = calloc(nknots+2*ord, sizeof(double)); 2389 double * oldPar = calloc(nknots+ord, sizeof(double)); 2390 double * knots2 = calloc(nknots+2*ord+1, sizeof(double)); 2391 double * params2 = calloc(nknots+ord+1, sizeof(double)); 2392 double * knots = theCurve->SplineKnots; 2393 double * params = theCurve->SplinePar; 2394 dcopyWrapper(nknots+2*ord,theCurve->SplineKnots,oldKnots); 2395 dcopyWrapper(nknots+ord,theCurve->SplinePar,oldPar); 2396 dcopyWrapper(nknots+2*ord,theCurve->SplineKnots,knots2); 2397 dcopyWrapper(nknots+ord,theCurve->SplinePar,params2); 2398 2399 int i=0; 2400 // compute probability of birth-death and move steps 2401 for(int j=ord;j<ncandknots+2*ord; j++) if(occ[j]==1) occind[i++]=j; 2402 double pk = EvalNknotsPrior( nknots, theCurve ); 2403 double pkp1 = (nknots < theCurve->SplineNknotsMax) ? EvalNknotsPrior(nknots+1, theCurve) : 0; 2404 double pkm1 = (nknots > 1) ? EvalNknotsPrior(nknots-1, theCurve) : 0; 2405 double pb = theCurve->SplineBDMConst[0] * dmin(1.0,pkp1/pk); // P(birth) 2406 double pd = theCurve->SplineBDMConst[0] * dmin(1.0,pkm1/pk); // P(death) 2407 double pm = dmax(0.0,1.0-pb-pd); // P(move) 2408 double u = runif(0,1); 2409 double baselik,candlik; 2410 // base likelihood 2411 if(theCurve->isHazard) baselik = LikelihoodSplineHazard(theCurve,frailty,regression); 2412 if(!theCurve->isHazard) baselik = LikelihoodSplineFrailty(hazard,theCurve,regression); 2413 if(u<pd){ // death step 2414 int j = (int) floor(runif(0,(double) nknots))+ord; // index of dying knot 2415 double x = knots[j]; // value of dying knot 2416 //Rprintf("Remove knot %d at %f (occ: %d)\n",j,knots[j],(int) occ[occind[j-ord]]); 2417 for(int i=0;i<nknots+ord;i++) if(i>=j-1) params2[i] = params2[i+1]; //remove params2[j]; 2418 for(int i=0;i<nknots+ord;i++) if(i>=j) knots2[i] = knots2[i+1]; //remove knots2[j2]; 2419 // update spline parameters for knot deletion 2420 if(ord>2) for(int j2=j-ord+1; j2<j-1; j2++){ 2421 double r2 = (x-knots2[j2])/(knots2[j2+ord-1]-knots2[j2]); 2422 double inner = 1/r2 * exp(params2[j2])-(1-r2)/r2*exp(params2[j2-1]); 2423 if(inner>0) params2[j2]=log(inner); else params2[j2] = theCurve->SplineMin[0]; 2424 } 2425 theCurve->SplineNknots--; theCurve->nj--; 2426 // Jacobian for the transformation 2427 double J = exp(params2[j-1]) / (exp(params[j-1]) - exp(params[j-2])); 2428 if(ord>2) for(int j2 = j-ord+2; j2<j-1; j2++){ 2429 double r2 = (x-knots2[j2])/(knots2[j2+ord-1]-knots2[j2]); 2430 J = J * exp(params2[j2])/(r2*exp(params[j2])); 2431 } 2432 // update the curve with the new knot set 2433 dcopyWrapper(nknots+2*ord, knots2,theCurve->SplineKnots); 2434 RemakeSplineBasis(theCurve,'d',j); 2435 // update curve with new parameters 2436 if(theCurve->isHazard) { 2437 UpdateSplinePar(theCurve,params2,-1); 2438 candlik = LikelihoodSplineHazard(theCurve,frailty,regression); 2439 } 2440 if(!theCurve->isHazard){ 2441 // for the frailty curve, make sure that the mean is 1 2442 double newmean = -exp(params2[j-2])*theCurve->SplineBasisExp[j-2]; 2443 for(int i=0;i<nknots+ord-1;i++) newmean += exp(params2[i])*theCurve->SplineBasisExp[i]; 2444 double newinner = -newmean/theCurve->SplineBasisExp[j-2]; 2445 if(newinner<0) candlik = -INFINITY; 2446 else{ 2447 params2[j-2] = log(newinner); 2448 UpdateSplinePar(theCurve,params2,-1); 2449 candlik = LikelihoodSplineFrailty(hazard,theCurve,regression); 2450 } 2451 } 2452 // prior ratio * Jacobian 2453 double ratio = sqrt(2*M_PI*theCurve->SplinePriorvar[0])*fabs(J); 2454 int acc = AcceptReject(baselik, candlik, ratio); 2455 if(acc){ 2456 // Update occupied index and candsd 2457 occ[occind[j-ord]] = 0; 2458 }else{ //undo damage 2459 theCurve->SplineNknots++; theCurve->nj++; 2460 dcopyWrapper(nknots+2*ord, oldKnots, theCurve->SplineKnots); 2461 RemakeSplineBasis(theCurve,'b',j); 2462 UpdateSplinePar(theCurve,oldPar,-1); 2463 } 2464 } 2465 if(u>pd & u<pd+pb){ 2466 // Birth 2467 int birthind = occind[0]; 2468 // choose a random unoccupied location for the knot to be born 2469 while(occ[birthind]) birthind = (int) floor(runif(0,(double) ncandknots))+ord; 2470 double x = candknots[birthind]; // value of the new knot 2471 int j = 0; while(knots[j]<x) j++; j--; // find the interval in which x lies 2472 //Rprintf("Birth at %f after knot %d\n",x,j); 2473 // update knot set and spline parameters 2474 for(int i=nknots+2*ord;i>j;i--) knots2[i]=knots2[i-1]; 2475 knots2[j+1]=x; 2476 for(int i=nknots+ord;i>j-ord+1;i--) params2[i]=params2[i-1]; 2477 for(int j2 = j-ord+2; j2<j+1; j2++){ 2478 double r2 = (x-knots[j2])/(knots[j2+ord-1]-knots[j2]); 2479 if(j2==j) r2 = runif(0,1); 2480 params2[j2] = log(r2*exp(params[j2])+(1-r2)*exp(params[j2-1])); 2481 } 2482 //Rprintf("Birthind,x,j: %d %f %d\n",birthind,x,j); 2483 theCurve->SplineNknots++; theCurve->nj++; 2484 // Jacobian of the transformation 2485 double J = (exp(params[j])-exp(params[j-1]))/exp(params2[j]); 2486 if(ord>2) for(int j2 = j-ord+2; j2<j; j2++){ 2487 double r2 = (x-knots[j2])/(knots[j2+ord-1]-knots[j2]); 2488 J = J * r2 * exp(params[j2])/exp(params2[j2]); 2489 } 2490 // update curve with new knots and parameters 2491 dcopyWrapper(nknots+2*ord+1, knots2,theCurve->SplineKnots); 2492 RemakeSplineBasis(theCurve,'b',j); 2493 if(theCurve->isHazard){ 2494 UpdateSplinePar(theCurve,params2,-1); 2495 candlik = LikelihoodSplineHazard(theCurve,frailty,regression); 2496 } 2497 if(!theCurve->isHazard){ 2498 // for frailty, make sure the mean is 1 2499 double newmean = -exp(params2[j])*theCurve->SplineBasisExp[j]; 2500 for(int i=0;i<nknots+ord+1;i++) newmean += exp(params2[i])*theCurve->SplineBasisExp[i]; 2501 double newinner = -newmean/theCurve->SplineBasisExp[j]; 2502 if(newinner<0) candlik = -INFINITY; 2503 else{ 2504 params2[j] = log(newinner); 2505 UpdateSplinePar(theCurve,params2,-1); 2506 candlik = LikelihoodSplineFrailty(hazard,theCurve,regression); 2507 } 2508 } 2509 // prior ratio * Jacobian 2510 double ratio = 1/sqrt(2*M_PI*theCurve->SplinePriorvar[0]) * fabs(J); 2511 int acc = AcceptReject(baselik, candlik, ratio); 2512 //Rprintf("Lik: %f %f %f %f %d\n",baselik,candlik, J, ratio,acc); 2513 if(acc){ 2514 // Update set of occupied candidate knots 2515 occ[birthind] = 1; 2516 }else{ //undo damage 2517 theCurve->SplineNknots--; theCurve->nj--; 2518 dcopyWrapper(nknots+2*ord, oldKnots, theCurve->SplineKnots); 2519 RemakeSplineBasis(theCurve,'d',j+1); 2520 UpdateSplinePar(theCurve,oldPar,-1); 2521 } 2522 } 2523 if(u>pb+pd) { // Move a knot 2524 // choose a random knot to move 2525 int moveind=0; 2526 moveind = (int) floor(runif(0,(double) nknots)); 2527 // find range of movement 2528 int leftknotind = (moveind == 0) ? ord : occind[moveind-1]+1; 2529 int rightknotind = (moveind == nknots-1) ? ncandknots + ord - 1 : occind[moveind+1]-1; 2530 // choose a candidate knot to place the new knot in 2531 int newknotind = (int) floor(runif((double) leftknotind,(double) rightknotind+1)); 2532 // if the knot doesn't stay the same, update the curve 2533 if(candknots[newknotind] != oldKnots[moveind+ord]){ 2534 theCurve->SplineKnots[moveind+ord] = candknots[newknotind]; 2535 // update spline basis with new knot position 2536 RemakeSplineBasis(theCurve,'m',moveind); 2537 if(theCurve->isHazard){ //hazard 2538 EvalSpline(theCurve,-1); 2539 candlik = LikelihoodSplineHazard(theCurve,frailty,regression); 2540 } 2541 if(!theCurve->isHazard){ //frailty 2542 // frailty needs to ensure the mean is 1 2543 double newmean = -exp(params2[moveind+ord])*theCurve->SplineBasisExp[moveind+ord]; 2544 for(int i=0;i<nknots+ord;i++) newmean += exp(params2[i])*theCurve->SplineBasisExp[i]; 2545 double newinner = -newmean/theCurve->SplineBasisExp[moveind+ord]; 2546 if(newinner<0) candlik = -INFINITY; 2547 else{ 2548 params2[moveind+ord] = log(newinner); 2549 UpdateSplinePar(theCurve,params2,-1); 2550 candlik = LikelihoodSplineFrailty(hazard,theCurve,regression); 2551 } 2552 } 2553 int acc = AcceptReject(baselik,candlik,1.0); 2554 //Rprintf("Lik: %f %f %d\n",baselik,candlik,acc); 2555 if(!acc){ 2556 // not accepted, undo the damage 2557 theCurve->SplineKnots[moveind+ord] = oldKnots[moveind+ord]; 2558 RemakeSplineBasis(theCurve,'m',moveind); 2559 if(theCurve->isHazard){ // hazard 2560 EvalSpline(theCurve,-1); 2561 }else{ // frailty 2562 UpdateSplinePar(theCurve,oldPar,-1); 2563 } 2564 }else{ // update occupied index 2565 occ[occind[moveind]]=0; 2566 occ[newknotind] = 1; 2567 } 2568 } 2569 } 2570 2571 free(oldKnots); 2572 free(oldPar); 2573 free(knots2); 2574 free(params2); 2575 free(occind); 2576 }