org.opensha.sha.earthquake.calc.recurInterval
Class BPT_DistCalc

java.lang.Object
  extended by org.opensha.sha.earthquake.calc.recurInterval.EqkProbDistCalc
      extended by org.opensha.sha.earthquake.calc.recurInterval.BPT_DistCalc
All Implemented Interfaces:
EventListener, ParameterChangeListener

public final class BPT_DistCalc
extends EqkProbDistCalc
implements ParameterChangeListener

Title: BPT_DistCalc.java

Description:

. This represents the Brownian Passage Time renewal model as described by Matthews et al. (2002, BSSA, vol 92, pp. 2223-2250). This adds several "Safe" methods are also supplied to avoid numerical artifacts; see description for each method for more details

Version:
1.0
Author:
Edward Field

Field Summary
 
Fields inherited from class org.opensha.sha.earthquake.calc.recurInterval.EqkProbDistCalc
adjustableParams, aperiodicity, APERIODICITY_PARAM_INFO, APERIODICITY_PARAM_NAME, aperiodicityParam, cdf, commonInfoString, DEFAULT_APERIODICITY_PARAM_VAL, DEFAULT_DELTAX_PARAM_VAL, DEFAULT_DURATION_PARAM_VAL, DEFAULT_HIST_OPEN_INTERVAL_PARAM_VAL, DEFAULT_MEAN_PARAM_VAL, DEFAULT_NUMPOINTS_PARAM_VAL, DELTA_X_DEFAULT, DELTA_X_PARAM_INFO, DELTA_X_PARAM_NAME, deltaX, deltaX_Param, duration, DURATION_PARAM_INFO, DURATION_PARAM_NAME, durationParam, HIST_OPEN_INTERVAL_PARAM_INFO, HIST_OPEN_INTERVAL_PARAM_NAME, histOpenInterval, histOpenIntParam, integratedCDF, integratedOneMinusCDF, mean, MEAN_PARAM_INFO, MEAN_PARAM_NAME, meanParam, NAME, NUM_POINTS_PARAM_INFO, NUM_POINTS_PARAM_NAME, numPoints, numPointsParam, pdf, upToDate
 
Constructor Summary
BPT_DistCalc()
           
 
Method Summary
protected  void computeDistributions()
           
 double getCondProb(double timeSinceLast, double duration)
          This is a version of the parent method getCondProb(*) that avoids numerical artifacts at high timeSinceLast (where 1-cdf gets too close to 0, and we therefore have division by zero).
static double getCondProb(double mean, double aperiodicity, double timeSinceLast, double duration)
          This computed the conditional probability using Trapezoidal integration (slightly more accurrate that the WGCEP-2002 code, which this method is modeled after).
 double getCondProbForUnknownTimeSinceLastEvent()
          This computes the probability of an event over the specified duration for the case where the date of last event is unknown (looping over all possible values), but where the historic open interval is applied (the latter defaults to zero if never set).
 EvenlyDiscretizedFunc getCondProbFunc()
          Overrides name and info assigned in parent
 double getSafeTimeSinceLastCutoff()
          This returns the maximum value of timeSinceLast (as discretized in the x-axis of the cdf) that is numerically safe (to avoid division by zero in the conditional probability calculations, where the denominator is 1-cdf).
static void main(String[] args)
          Main method for running tests.
 void testSafeCalcs()
          This tests the difference between values obtained using the local and parent getCondProbFunc(*) methods for cases within safeTimeSinceLast, where the fix for small normalized durations is applied in the local method.
 
Methods inherited from class org.opensha.sha.earthquake.calc.recurInterval.EqkProbDistCalc
computeMeanFromPDF, fitToThisFunction, getAdjParams, getAperiodicity, getCDF, getCondProbFunc, getCondProbGainFunc, getHazFunc, getMean, getMeanTimeSinceLastEventPDF, getName, getPDF, getTimeSinceLastEventPDF, initAdjParams, makeIntegratedCDFs, parameterChange, setAll, setAll, setAll, setAll, setAllParameters, setDuration, setDurationAndHistOpenInterval
 
Methods inherited from class java.lang.Object
clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
 
Methods inherited from interface org.opensha.commons.param.event.ParameterChangeListener
parameterChange
 

Constructor Detail

BPT_DistCalc

public BPT_DistCalc()
Method Detail

computeDistributions

protected void computeDistributions()

getCondProb

public static double getCondProb(double mean,
                                 double aperiodicity,
                                 double timeSinceLast,
                                 double duration)
This computed the conditional probability using Trapezoidal integration (slightly more accurrate that the WGCEP-2002 code, which this method is modeled after). Although this method is static (doesn't require instantiation), it is less efficient than the non-static version here (it is also very slightly less accurate because the other interpolates the cdf). Note also that if timeSinceLast/mean > aperiodicity*10, timeSinceLast is changed to equal mean*aperiodicity*10 (to avoid numerical problems at high timeSinceLast).

Parameters:
timeSinceLast - - time since last event
rate - - average rate of events
alpha - - coefficient of variation (technically corrrect??)
duration - - forecast duration
Returns:

getCondProb

public double getCondProb(double timeSinceLast,
                          double duration)
This is a version of the parent method getCondProb(*) that avoids numerical artifacts at high timeSinceLast (where 1-cdf gets too close to 0, and we therefore have division by zero). We know that at infinite time the mean residual life (expected time to next earthquake) for BPT approaches the following asymptotically: mrl = 2*mean*aperiodicity*aperiodicity (Equation 24 of Matthews et al. (2002, BSSA, vol 92, pp. 2223-2250)), so the conditional probability becomes 1-exp(duration/mrl) at large timeSinceLast. We also define a safeTimeSinceLast as that where 1-cdf becomes close to zero (<= SAFE_ONE_MINUS_CDF). Beyond this point we assume the conditional probability varies linearly between the value at safeTimeSinceLast and 1-exp(duration/mrl), where "infinite time" is taken as timeSinceLast/mean = 10. While this choice of 10 is somewhat arbitrary, differences are only slight, and mostly only for low aperiodicity (<0.15). Also keep in mind that the probability of getting beyond safeTimeSinceLast is SAFE_ONE_MINUS_CDF (e.g., 1e-13), so we are way out on the tail, where application of the model will always be speculative anyway. The above does not solve all numerical problems, as we can get instability before safeTimeSinceLast where p2-p1 in the parent method becomes too small (e.g., 1e-14), which occurs when duration/mean is <= 1e-3 (and again, only for aperiodicities <= 0.2). To avoid this problem, the smallest duration/mean for which conditional probability is computed is MIN_NORM_DURATION, with values below this simply being the value for MIN_NORM_DURATION scaled by (duration/mean)/MIN_NORM_DURATION. The testSafeCalcs show that this produces conditional probabilities within 0.1% of true values where these numerical artifacts are lacking, and shows that discrepancies elsewhere are for very un-probable states (e.g., just below safeTimeSinceLast); nonetheless, this fixes these other numerical artifacts.

Overrides:
getCondProb in class EqkProbDistCalc
Parameters:
timeSinceLast -
duration -
Returns:

getCondProbFunc

public EvenlyDiscretizedFunc getCondProbFunc()
Overrides name and info assigned in parent

Overrides:
getCondProbFunc in class EqkProbDistCalc
Returns:

getCondProbForUnknownTimeSinceLastEvent

public double getCondProbForUnknownTimeSinceLastEvent()
This computes the probability of an event over the specified duration for the case where the date of last event is unknown (looping over all possible values), but where the historic open interval is applied (the latter defaults to zero if never set). This avoids numerical artifacts in the parent method a couple ways. First, if histOpenInterval>=safeTimeSinceLast, we simply return the conditional probability at safeTimeSinceLast (technically we should weight average the values beyond, but they are nearly constant, the weights are decaying beyond (so that at safeTimeSinceLast is highest), and the likelihood of histOpenInterval ever getting to safeTimeSinceLast is vanishingly small (as noted in the doc for the getCondProb method here)). A less efficient way of calculating this cond probability is also used if unstable results are expected (see details regarding "numer" and "denom" below).

Overrides:
getCondProbForUnknownTimeSinceLastEvent in class EqkProbDistCalc
Returns:

getSafeTimeSinceLastCutoff

public double getSafeTimeSinceLastCutoff()
This returns the maximum value of timeSinceLast (as discretized in the x-axis of the cdf) that is numerically safe (to avoid division by zero in the conditional probability calculations, where the denominator is 1-cdf). This returns Double.Nan if no x-axis values are safe (not even the first ones). The threshold for safe values was found by trial and error and checked for aperiodicity values between 0.1 and 1.0 (using the GUI).

Returns:

testSafeCalcs

public void testSafeCalcs()
This tests the difference between values obtained using the local and parent getCondProbFunc(*) methods for cases within safeTimeSinceLast, where the fix for small normalized durations is applied in the local method. This prints all cases that differ by more than 0.1%. All those printed are indeed the problems we sought to avoid, and are very rare circumstances anyway.


main

public static void main(String[] args)
Main method for running tests. Test1 compares the static getCondProb(*) method against values obtained directly from the WGCEP-2002 code; all are within 0.3%. Test2 campares the non static getCondProb(*) method against the static; all are within 0.5%. The differences is that the non-static is slightly more accurate due to interpolation of the CDF (exact same values are obtained otherwise; see commented out code). Test3 is the non-static used more efficiently (exact same values as from non-static above); this is about a factor of two faster. Test4 examines what happens if delta is changed to 0.01 in the non-static method (also about a factor of two faster).