#!/usr/bin/python -Wall # ================================================================ # John Kerl # kerl at math dot arizona dot edu # 2008-05-07 # See the paper for a full description. # ================================================================ from __future__ import division from math import * import random import sys # ---------------------------------------------------------------- def p(alpha, eta, t): eabt = exp(alpha - exp(eta)*t) return eabt / (1 + eabt) # ---------------------------------------------------------------- def ell(f, t, alpha, eta): n = len(f) prod = 1.0 for i in range(0, n): paet = p(alpha, eta, t[i]) if f[i] == 1: prod *= paet else: prod *= 1 - paet return prod # ---------------------------------------------------------------- def mcmc_test(f, t, alphahat, etahat, sigma_1, sigma_2, test_t, num_mcmc, mcmc_skip, mcmc_min, verbose=0): naccept_alpha = 0 naccept_eta = 0 alpha = alphahat eta = etahat paets = [] for stepi in range(0, num_mcmc): alphap = random.normalvariate(alpha, sigma_1) etap = random.normalvariate(eta, sigma_2) ell1 = ell(f, t, alphap, eta) ell2 = ell(f, t, alpha, etap) elld = ell(f, t, alpha, eta) exp1 = exp(-0.5/sigma_1**2 * \ ((alphap - alphahat)**2 - (alpha - alphahat)**2)) exp2 = exp(-0.5/sigma_2**2 * \ ((etap - etahat)**2 - (eta - etahat)**2)) r1 = ell1 / elld * exp1 r2 = ell2 / elld * exp2 if (r1 > 1): r1 = 1 if (r2 > 1): r2 = 1 paet = p(alpha, eta, test_t) # Downsample to avoid correlation if ((stepi % mcmc_skip) == 0) and (stepi >= mcmc_min): paets.append(paet) if (verbose): print "%7d %11.7f %11.7f %11.7f" % \ (stepi, alpha, eta, paet) # Metropolis update if (random.uniform(0.0, 1.0) < r1): alpha = alphap naccept_alpha += 1 if (random.uniform(0.0, 1.0) < r2): eta = etap naccept_eta += 1 if (verbose): print "#Acceptance ratios %11.7f %11.7f" % \ (naccept_alpha / num_mcmc, naccept_eta / num_mcmc) # Take the 10th percentile paets.sort() num_paets = len(paets) idx = int(num_paets * 0.1) return paets[idx] # ================================================================ # Parameters for the prior distribution alphahat = 15 etahat = -1.47 sigma_1 = 6 sigma_2 = 0.12 # Failure and temperature data. f = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, \ 1, 0, 0, 0, 0, 0] t = [53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 70, 70, 70, 70, 72, 73, 75, \ 75, 76, 76, 78, 79, 81] # MCMC parameters num_mcmc = 100000 mcmc_skip = 50 mcmc_min = 100 # ================================================================ if (len(sys.argv) == 2): test_t = int(sys.argv[1]) mcmc_test(f, t, alphahat, etahat, sigma_1, sigma_2, test_t, \ num_mcmc, mcmc_skip, mcmc_min, verbose=1) else: for test_t in [31, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100]: d = mcmc_test(f, t, alphahat, etahat, sigma_1, sigma_2, test_t, \ num_mcmc, mcmc_skip, mcmc_min, verbose=0) print "%d %11.7f" % (test_t, d)