Package com.nr.test.test_chapter15

Source Code of com.nr.test.test_chapter15.Test_mcmc

package com.nr.test.test_chapter15;

import static java.lang.Math.abs;
import static org.junit.Assert.fail;

import org.junit.After;
import org.junit.Before;
import org.junit.Test;

import com.nr.model.Plog;
import com.nr.model.Proposal;
import com.nr.model.State;
import com.nr.ran.Gammadev;

public class Test_mcmc {

  @Before
  public void setUp() throws Exception {
  }

  @After
  public void tearDown() throws Exception {
  }

  @SuppressWarnings("unused")
  @Test
  public void test() {
    int i,j,k1=1,k2=2,N=1000;
    double lambda1=3.0,lambda2=2.0,t,tc=200.0,accept;
    double llam1,llam2,kk1,kk2,ttc;
    double lam1ave,lam2ave,tcave;
    double[] times=new double[N];
    boolean localflag, globalflag=false;

   

    // Test mcmc
    System.out.println("Testing mcmc");

    Gammadev gdev1 = new Gammadev(k1,lambda1,17);
    Gammadev gdev2 = new Gammadev(k2,lambda2,17);
    t=0.0;
    for (i=0;i<N;i++) {
      if (t < tc) {
        t += gdev1.dev();
      } else {
        t += gdev2.dev();
      }
      times[i]=t;
    }

    State s=new State(1.0,3.0,100.0,1,1);
    Plog plog = new Plog(times);
    Proposal propose = new Proposal(10102,0.01);
    for (i=0;i<1000;i++) {    // Burn-in
//      if (i%10 == 0)
//        System.out.printf(s.lam1 << " %f\n", s.lam2 << " %f\n", s.k1 << " %f\n", s.k2 << " %f\n", s.tc);
      accept=Proposal.mcmcstep(1,s,plog,propose);
    }
   
    llam1=llam2=kk1=kk2=ttc=0.0;
    j=0;
    for (i=0;i<10000;i++) {   // Production
      accept = Proposal.mcmcstep(10,s,plog,propose);
      if (i%10 == 0) {
        llam1 += s.lam1;
        llam2 += s.lam2;
        kk1 = s.k1;
        kk2 = s.k2;
        ttc += s.tc;
        j++;
      }
    }
    lam1ave=llam1/j;
    lam2ave=llam2/j;
    tcave=ttc/j;

//    System.out.printf(lam1ave << " %f\n", lam2ave << " %f\n", kk1 << " %f\n", kk2 << " %f\n", tcave);

    localflag = abs(lam1ave-lambda1)/lambda1 > 0.05;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** mcmc: Model parameter lambda1 was not correctly determined");
     
    }

    localflag = abs(lam2ave-lambda2)/lambda2 > 0.05;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** mcmc: Model parameter lambda2 was not correctly determined");
     
    }

    localflag = kk1 != k1;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** mcmc: Model parameter k1 was not correctly determined");
     
    }

    localflag = kk2 != k2;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** mcmc: Model parameter k2 was not correctly determined");
     
    }
   
    localflag = abs(tcave-tc)/tc > 0.01;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** mcmc: Critical time tc was not accurately determined");
     
    }

    if (globalflag) System.out.println("Failed\n");
    else System.out.println("Passed\n");
  }

}
TOP

Related Classes of com.nr.test.test_chapter15.Test_mcmc

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.