Repository: nasa/prog_algs (language: python)

Issue:
Remove depreciated measurement_eqn
Measurement_eqn is being removed for v1.4. This PR removes it from the state estimators


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/pr-messages.yml b/.github/workflows/pr-messages.yml
new file mode 100644
index 00000000..3a839832
--- /dev/null
+++ b/.github/workflows/pr-messages.yml
@@ -0,0 +1,34 @@
+name: Print PR Message - Non Release
+
+on: 
+  pull_request:
+    branches:
+      - 'dev'
+    types: [opened]
+
+jobs:
+  benchmark_branch: 
+    runs-on: ubuntu-latest
+    steps:
+      - name: Auto Comment
+        uses: wow-actions/auto-comment@v1
+        with:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          pullRequestOpened: >
+            Thank you for opening this PR. Each PR into dev requires a code review. For the code review, look at the following:
+
+            - [ ] Reviewer should look for bugs, efficiency, readability, testing, and coverage in examples (if relevant).
+
+            - [ ] Ensure that each PR adding a new feature should include a test verifying that feature.
+
+            - [ ] All tests must be passing.
+
+            - [ ] All errors from static analysis must be resolved.
+            
+            - [ ] Review the test coverage reports (if there is a change) - will be added as comment on PR if there is a change
+
+            - [ ] Review the software benchmarking results (if there is a change) - will be added as comment on PR
+
+            - [ ] Any added dependencies are included in requirements.txt, setup.py, and dev_guide.rst (this document)
+
+            - [ ] All warnings from static analysis must be reviewed and resolved - if deemed appropriate.
diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index f26b70e8..79136947 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -32,25 +32,25 @@ jobs:
         python -m pip install testbook
     - name: Run tests
       run: python -m tests
-  test-prog_models-released:
-    runs-on: ubuntu-latest
-    strategy:
-      matrix:
-        python-version: ['3.7', '3.8', '3.9', '3.10']
-    steps:
-    - uses: actions/checkout@v2
-    - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v2
-      with:
-        python-version: ${{ matrix.python-version }}
-    - name: Install
-      run: |
-        python -m pip install --upgrade pip
-        python -m pip install -e .
-        python -m pip install notebook
-        python -m pip install testbook
-    - name: Run tests
-      run: python -m tests
+  # test-prog_models-released:
+  #   runs-on: ubuntu-latest
+  #   strategy:
+  #     matrix:
+  #       python-version: ['3.9']
+  #   steps:
+  #   - uses: actions/checkout@v2
+  #   - name: Set up Python ${{ matrix.python-version }}
+  #     uses: actions/setup-python@v2
+  #     with:
+  #       python-version: ${{ matrix.python-version }}
+  #   - name: Install
+  #     run: |
+  #       python -m pip install --upgrade pip
+  #       python -m pip install -e .
+  #       python -m pip install notebook
+  #       python -m pip install testbook
+  #   - name: Run tests
+  #     run: python -m tests
   copyright:
     runs-on: ubuntu-latest
     strategy:
diff --git a/.github/workflows/release-messages.yml b/.github/workflows/release-messages.yml
new file mode 100644
index 00000000..2842dc5a
--- /dev/null
+++ b/.github/workflows/release-messages.yml
@@ -0,0 +1,56 @@
+name: Print PR Message - Release
+
+on: 
+  pull_request:
+    branches:
+      - 'release/**'
+    types: [opened]
+
+jobs:
+  benchmark_branch: 
+    runs-on: ubuntu-latest
+    steps:
+      - name: Auto Comment
+        uses: wow-actions/auto-comment@v1
+        with:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          pullRequestOpened: >
+            Thank you for opening this PR. Since this is a release branch, the PR must complete the release checklist, below:
+
+            - [ ] Check that each new feature has corresponding tests
+
+            - [ ] Confirm all dependencies are in the following: requirements.txt, setup.py, the bottom of dev_guide.rst
+
+            - [ ] Confirm that all issues associated with the release have been closed (i.e., requirements have been met) or assigned to another release
+
+            - [ ] Run unit tests `python -m tests`
+
+            - [ ] If present, run manual tests `python -m tests.test_manual`
+
+            - [ ] Review the template(s)
+
+            - [ ] Review static-analysis/linter results
+
+            - [ ] Review the tutorial
+
+            - [ ] Run and review the examples
+
+            - [ ] Check that all examples are tested
+
+            - [ ] Check new files in PR for any accidentally added
+
+            - [ ] Check documents
+
+              - [ ] Check that all desired examples are in docs
+
+              - [ ] General review: see if any updates are required
+
+            - [ ] Rebuild sphinx documents: `sphinx-build sphinx_config/ docs/`
+
+            - [ ] Write release notes
+
+            - [ ] Update version number in src/\*/__init__.py and setup.py
+
+            - [ ] For releases adding new features- ensure that NASA release process has been followed.
+            
+            - [ ] Confirm that on GitHub Releases page, the next release has been started and that a schedule is present including at least Release Date, Release Review Date, and Release Branch Opening Date.`
diff --git a/README.md b/README.md
index c39cca93..37efef85 100644
--- a/README.md
+++ b/README.md
@@ -12,18 +12,12 @@ This is designed to be used with the [Prognostics Models Package](https://github
 `pip3 install prog_algs`
 
 ## Documentation
-See documentation [here](https://nasa.github.io/prog_algs/)
+See documentation [here](https://nasa.github.io/progpy/prog_algs_guide.html)
 
 ## Repository Directory Structure 
 
 `src/prog_algs/` - The prognostics algorithm python package<br />
-&nbsp;&nbsp; |- `metrics/` - Tools for analyzing the results of state estimation and prognostics
-&nbsp;&nbsp; |- `predictors/` - Algorithms for performing the prediction step of model-based prognostics<br />
-&nbsp;&nbsp; |- `state_estimators/` - Algorithms for performing the state estimation step of model-based prognostics<br />
-&nbsp;&nbsp; |- `uncertain_data/` - Tools for storing and manipulating data with uncertainty<br />
-&nbsp;&nbsp; |- `visualize/` - Tools for visualizing the results of state estimation and prognostics<br />
 `docs/` - Project documentation (see also [github.io](https://nasa.github.io/prog_algs/))<br />
-`sphinx_config/` - Configuration for automatic documentation generation<br />
 `examples/` - Example Python scripts using prog_algs<br />
 `tests/` - Tests for prog_models<br />
 `README.md` - The readme (this file)<br />
@@ -35,16 +29,20 @@ Use the following to cite this repository:
 @misc{2022_nasa_prog_algs,
     author    = {Christopher Teubert and Matteo Corbetta and Chetan Kulkarni},
     title     = {Prognostics Algorithm Python Package},
-    month     = May,
+    month     = Oct,
     year      = 2022,
-    version   = {1.3.0},
+    version   = {1.4.0},
     url       = {https://github.com/nasa/prog\_algs}
     }
 ```
 
 The corresponding reference should look like this:
 
-C. Teubert, M. Corbetta, C. Kulkarni, Prognostics Algorithm Python Package, v1.3.0, May 2022. URL https://github.com/nasa/prog_algs.
+C. Teubert, M. Corbetta, C. Kulkarni, Prognostics Algorithm Python Package, v1.4.0, Oct 2022. URL https://github.com/nasa/prog_algs.
+
+Alternatively, if using both prog_models and prog_algs, you can cite the combined package as
+
+C. Teubert, C. Kulkarni, M. Corbetta, K. Jarvis, M. Daigle, ProgPy Prognostics Python Packages, v1.4, October 2022. URL https://nasa.github.io/progpy.
 
 ## Acknowledgements
 The structure and algorithms of this package are strongly inspired by the [MATLAB Prognostics Algorithm Library](https://github.com/nasa/PrognosticsAlgorithmLibrary) and the [MATLAB Prognostics Metrics Library](https://github.com/nasa/PrognosticsMetricsLibrary). We would like to recognize Matthew Daigle, Shankar Sankararaman and the rest of the team that contributed to the Prognostics Model Library for the contributions their work on the MATLAB library made to the design of prog_algs
@@ -53,7 +51,6 @@ The structure and algorithms of this package are strongly inspired by the [MATLA
 Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration.  All Rights Reserved.
 
 ## Disclaimers
-
 No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND, EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION, IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER, CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS, RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND DISTRIBUTES IT "AS IS."
 
 Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.  IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES, EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON, OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
diff --git a/examples/__init__.py b/examples/__init__.py
index 2bd4c7c0..010b98c3 100644
--- a/examples/__init__.py
+++ b/examples/__init__.py
@@ -1,3 +1,3 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration. All Rights Reserved.
 
-__all__ = ['basic_example', 'benchmarking_example', 'eol_event', 'horizon', 'kalman_filter', 'utpredictor', 'new_state_estimator_example', 'measurement_eqn_example', 'playback', 'predict_specific_event', 'thrown_object_example']
+__all__ = ['basic_example', 'basic_example_battery', 'benchmarking_example', 'eol_event', 'horizon', 'kalman_filter', 'new_state_estimator_example', 'measurement_eqn_example', 'playback', 'predict_specific_event']
diff --git a/examples/basic_example.py b/examples/basic_example.py
index 0cbe2fcd..09a19b9c 100644
--- a/examples/basic_example.py
+++ b/examples/basic_example.py
@@ -3,7 +3,7 @@
 """
 This example performs a state estimation and prediction with uncertainty given a Prognostics Model.
  
-Method: An instance of the BatteryCircuit model in prog_models is created, and the prediction process is achieved in three steps:
+Method: An instance of the ThrownObject model in prog_models is created, and the prediction process is achieved in three steps:
     1) State estimation of the current state is performed using a chosen state_estimator, and samples are drawn from this estimate
     2) Prediction of future states (with uncertainty) and the times at which the event threshold will be reached
     3) Metrics tools are used to further investigate the results of prediction
@@ -14,72 +14,60 @@
     iv) Figures illustrating results
 """
 
-from prog_models.models import BatteryCircuit as Battery
-# VVV Uncomment this to use Electro Chemistry Model VVV
-from prog_models.models import BatteryElectroChemEOD as Battery
-
+from prog_models.models import ThrownObject
 from prog_algs import *
 
 def run_example():
     # Step 1: Setup model & future loading
+    m = ThrownObject(process_noise = 1)
     def future_loading(t, x = None):
-        # Variable (piece-wise) future loading scheme 
-        if (t < 600):
-            i = 2
-        elif (t < 900):
-            i = 1
-        elif (t < 1800):
-            i = 4
-        elif (t < 3000):
-            i = 2
-        else:
-            i = 3
-        return {'i': i}
-    # Measurement noise
-    R_vars = {
-        't': 2, 
-        'v': 0.02
-    }
-    batt = Battery(measurement_noise = R_vars)
-    initial_state = batt.parameters['x0']
+        # No load for a thrown object
+        return m.InputContainer({})
+    initial_state = m.initialize()
 
     # Step 2: Demonstrating state estimator
+    # The state estimator is used to estimate the system state given sensor data. 
     print("\nPerforming State Estimation Step")
 
     # Step 2a: Setup
-    filt = state_estimators.ParticleFilter(batt, initial_state)
+    filt = state_estimators.ParticleFilter(m, initial_state)
     # VVV Uncomment this to use UKF State Estimator VVV
-    # filt = state_estimators.UnscentedKalmanFilter(batt, initial_state)
+    # filt = state_estimators.UnscentedKalmanFilter(m, initial_state)
 
     # Step 2b: Print & Plot Prior State
     print("Prior State:", filt.x.mean)
-    print('\tSOC: ', batt.event_state(filt.x.mean)['EOD'])
+    print('\nevent state: ', m.event_state(filt.x.mean))
     fig = filt.x.plot_scatter(label='prior')
 
-    # Step 2c: Perform state estimation step
-    example_measurements = {'t': 32.2, 'v': 3.915}
+    # Step 2c: Perform state estimation step, given some measurement, above what's expected
+    example_measurements = m.OutputContainer({'x': 7.5})
     t = 0.1
-    filt.estimate(t, future_loading(t), example_measurements)
+    u = future_loading(t)
+    filt.estimate(t, u, example_measurements)  # Update state, given (example) sensor data
 
     # Step 2d: Print & Plot Resulting Posterior State
+    # Note the posterior state is greater than the predicted state of 5.95
+    # This is because of the high measurement 
     print("\nPosterior State:", filt.x.mean)
-    print('\tSOC: ', batt.event_state(filt.x.mean)['EOD'])
+    # Event state for 'falling' is less, because velocity has decreased
+    print('\nEvent State: ', m.event_state(filt.x.mean))
     filt.x.plot_scatter(fig=fig, label='posterior')  # Add posterior state to figure from prior state
 
     # Note: in a prognostic application the above state estimation step would be repeated each time
     #   there is new data. Here we're doing one step to demonstrate how the state estimator is used
 
-    # Step 3: Demonstrating Predictor
+    # Step 3: Demonstrating Prediction
     print("\n\nPerforming Prediction Step")
 
     # Step 3a: Setup Predictor
-    mc = predictors.MonteCarlo(batt)
+    mc = predictors.MonteCarlo(m)
 
     # Step 3b: Perform a prediction
-    NUM_SAMPLES = 5
-    STEP_SIZE = 0.1
-    mc_results = mc.predict(filt.x, future_loading, n_samples = NUM_SAMPLES, dt=STEP_SIZE)
-    print('ToE', mc_results.time_of_event.mean)
+    NUM_SAMPLES = 50
+    STEP_SIZE = 0.01
+    mc_results = mc.predict(filt.x, future_loading, n_samples = NUM_SAMPLES, dt=STEP_SIZE, save_freq=STEP_SIZE)
+    print('Predicted time of event (ToE): ', mc_results.time_of_event.mean)
+    # Here there are 2 events predicted, when the object starts falling, and when it impacts the ground.
 
     # Step 3c: Analyze the results
 
@@ -92,19 +80,21 @@ def future_loading(t, x = None):
     # now you have all the samples corresponding to times[1]
 
     # You can also access the final state (of type UncertainData), like so:
+    # Note: to get a more accurate final state, you can decrease the step size.
     final_state = mc_results.time_of_event.final_state
-    print('Final state @EOD: ', final_state['EOD'].mean)
+    print('State when object starts falling: ', final_state['falling'].mean)
     
     # You can also use the metrics package to generate some useful metrics on the result of a prediction
     print("\nEOD Prediction Metrics")
 
     from prog_algs.metrics import prob_success
-    print('\tPortion between 3005.2 and 3005.6: ', mc_results.time_of_event.percentage_in_bounds([3005.2, 3005.6]))
-    print('\tAssuming ground truth 3002.25: ', mc_results.time_of_event.metrics(ground_truth=3005.25))
-    print('\tP(Success) if mission ends at 3002.25: ', prob_success(mc_results.time_of_event, 3005.25))
+    print('\tPortion between 3.65 and 3.8: ', mc_results.time_of_event.percentage_in_bounds([3.65, 3.8], keys='falling'))
+    print('\tAssuming ground truth 3.7: ', mc_results.time_of_event.metrics(ground_truth=3.7, keys='falling'))
+    print('\tP(Success) if mission ends at 7.6: ', prob_success(mc_results.time_of_event, 7.6, keys='impact'))
 
     # Plot state transition 
     # Here we will plot the states at t0, 25% to ToE, 50% to ToE, 75% to ToE, and ToE
+    # You should see the states move together (i.e., velocity is lowest and highest when closest to the ground (before impact, and at beginning, respectively))
     fig = mc_results.states.snapshot(0).plot_scatter(label = "t={} s".format(int(mc_results.times[0])))  # 0
     quarter_index = int(len(mc_results.times)/4)
     mc_results.states.snapshot(quarter_index).plot_scatter(fig = fig, label = "t={} s".format(int(mc_results.times[quarter_index])))  # 25%
@@ -112,7 +102,11 @@ def future_loading(t, x = None):
     mc_results.states.snapshot(quarter_index*3).plot_scatter(fig = fig, label = "t={} s".format(int(mc_results.times[quarter_index*3])))  # 75%
     mc_results.states.snapshot(-1).plot_scatter(fig = fig, label = "t={} s".format(int(mc_results.times[-1])))  # 100%
 
-    mc_results.time_of_event.plot_hist()
+    # Plot time of event for each event
+    # If you dont see many bins here, this is because there is not much variety in the estimate. 
+    # You can increase the number of bins, decrease step size, or increase the number of samples to see more of a distribution
+    mc_results.time_of_event.plot_hist(keys='impact')
+    mc_results.time_of_event.plot_hist(keys='falling')
     
     # Step 4: Show all plots
     import matplotlib.pyplot as plt  # For plotting
diff --git a/examples/basic_example_battery.py b/examples/basic_example_battery.py
new file mode 100644
index 00000000..4599d553
--- /dev/null
+++ b/examples/basic_example_battery.py
@@ -0,0 +1,148 @@
+# Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration. All Rights Reserved.
+
+"""
+This example extends the "basic example" to perform a state estimation and prediction with uncertainty given a more complicated model. Models, state estimators, and predictors can be switched out. See documentation nasa.github.io/progpy for description of options
+ 
+Method: An instance of the BatteryCircuit model in prog_models is created, and the prediction process is achieved in three steps:
+    1) State estimation of the current state is performed using a chosen state_estimator, and samples are drawn from this estimate
+    2) Prediction of future states (with uncertainty) and the times at which the event threshold will be reached
+    3) Metrics tools are used to further investigate the results of prediction
+
+Results: 
+    i) Predicted future values (inputs, states, outputs, event_states) with uncertainty from prediction
+    ii) Time event is predicted to occur (with uncertainty)
+    iii) Various prediction metrics
+    iv) Figures illustrating results
+"""
+
+from prog_models.models import BatteryCircuit as Battery
+# VVV Uncomment this to use Electro Chemistry Model VVV
+# from prog_models.models import BatteryElectroChemEOD as Battery
+
+from prog_algs.state_estimators import ParticleFilter as StateEstimator
+# VVV Uncomment this to use UKF State Estimator VVV
+# from prog_algs.state_estimators import UnscentedKalmanFilter as StateEstimator
+
+from prog_algs.predictors import MonteCarlo as Predictor
+# VVV Uncomment this to use UnscentedTransform Predictor VVV
+# from prog_algs.predictors import UnscentedTransformPredictor as Predictor
+
+def run_example():
+    # Step 1: Setup model & future loading
+    # Measurement noise
+    R_vars = {
+        't': 2, 
+        'v': 0.02
+    }
+    batt = Battery(process_noise = 0.25, measurement_noise = R_vars)
+    # Creating the input containers outside of the function accelerates prediction
+    loads = [
+        batt.InputContainer({'i': 2}),
+        batt.InputContainer({'i': 1}),
+        batt.InputContainer({'i': 4}),
+        batt.InputContainer({'i': 2}),
+        batt.InputContainer({'i': 3})
+    ]
+    def future_loading(t, x = None):
+        # Variable (piece-wise) future loading scheme 
+        if (t < 600):
+            return loads[0]
+        elif (t < 900):
+            return loads[1]
+        elif (t < 1800):
+            return loads[2]
+        elif (t < 3000):
+            return loads[3]
+        return loads[-1]
+
+    initial_state = batt.initialize()
+
+    # Step 2: Demonstrating state estimator
+    print("\nPerforming State Estimation Step")
+
+    # Step 2a: Setup
+    filt = StateEstimator(batt, initial_state)
+    
+    # Step 2b: Print & Plot Prior State
+    print("Prior State:", filt.x.mean)
+    print('\tSOC: ', batt.event_state(filt.x.mean)['EOD'])
+    fig = filt.x.plot_scatter(label='prior')
+
+    # Step 2c: Perform state estimation step
+    example_measurements = batt.OutputContainer({'t': 32.2, 'v': 3.915})
+    t = 0.1
+    u = future_loading(t)
+    filt.estimate(t, u, example_measurements)
+
+    # Step 2d: Print & Plot Resulting Posterior State
+    print("\nPosterior State:", filt.x.mean)
+    print('\tSOC: ', batt.event_state(filt.x.mean)['EOD'])
+    filt.x.plot_scatter(fig=fig, label='posterior')  # Add posterior state to figure from prior state
+
+    # Note: in a prognostic application the above state estimation step would be repeated each time
+    #   there is new data. Here we're doing one step to demonstrate how the state estimator is used
+
+    # Step 3: Demonstrating Predictor
+    print("\n\nPerforming Prediction Step")
+
+    # Step 3a: Setup Predictor
+    mc = Predictor(batt)
+
+    # Step 3b: Perform a prediction
+    NUM_SAMPLES = 25
+    STEP_SIZE = 0.1
+    SAVE_FREQ = 100  # How often to save results
+    mc_results = mc.predict(filt.x, future_loading, n_samples = NUM_SAMPLES, dt=STEP_SIZE, save_freq = SAVE_FREQ)
+    print('ToE', mc_results.time_of_event.mean)
+
+    # Step 3c: Analyze the results
+
+    # Note: The results of a sample-based prediction can be accessed by sample, e.g.,
+    from prog_algs.predictors import UnweightedSamplesPrediction
+    if isinstance(mc_results, UnweightedSamplesPrediction):
+        states_sample_1 = mc_results.states[1]
+    # now states_sample_1[n] corresponds to times[n] for the first sample
+
+    # You can also access a state distribution at a specific time using the .snapshot function
+    states_time_1 = mc_results.states.snapshot(1)
+    # now you have all the samples corresponding to times[1]
+
+    # Print Results
+    print('Results: ')
+    for i, time in enumerate(mc_results.times):
+        print('\nt = {}'.format(time))
+        print('\tu = {}'.format(mc_results.inputs.snapshot(i).mean))
+        print('\tx = {}'.format(mc_results.states.snapshot(i).mean))
+        print('\tz = {}'.format(mc_results.outputs.snapshot(i).mean))
+        print('\tevent state = {}'.format(mc_results.event_states.snapshot(i).mean))
+
+    # You can also access the final state (of type UncertainData), like so:
+    final_state = mc_results.time_of_event.final_state
+    print('Final state @EOD: ', final_state['EOD'].mean)
+    
+    # You can also use the metrics package to generate some useful metrics on the result of a prediction
+    print("\nEOD Prediction Metrics")
+
+    from prog_algs.metrics import prob_success
+    print('\tPortion between 3005.2 and 3005.6: ', mc_results.time_of_event.percentage_in_bounds([3005.2, 3005.6]))
+    print('\tAssuming ground truth 3002.25: ', mc_results.time_of_event.metrics(ground_truth=3005.25))
+    print('\tP(Success) if mission ends at 3002.25: ', prob_success(mc_results.time_of_event, 3005.25))
+
+    # Plot state transition 
+    # Here we will plot the states at t0, 25% to ToE, 50% to ToE, 75% to ToE, and ToE
+    fig = mc_results.states.snapshot(0).plot_scatter(label = "t={} s".format(int(mc_results.times[0])))  # 0
+    quarter_index = int(len(mc_results.times)/4)
+    mc_results.states.snapshot(quarter_index).plot_scatter(fig = fig, label = "t={} s".format(int(mc_results.times[quarter_index])))  # 25%
+    mc_results.states.snapshot(quarter_index*2).plot_scatter(fig = fig, label = "t={} s".format(int(mc_results.times[quarter_index*2])))  # 50%
+    mc_results.states.snapshot(quarter_index*3).plot_scatter(fig = fig, label = "t={} s".format(int(mc_results.times[quarter_index*3])))  # 75%
+    mc_results.states.snapshot(-1).plot_scatter(fig = fig, label = "t={} s".format(int(mc_results.times[-1])))  # 100%
+
+    mc_results.time_of_event.plot_hist()
+    
+    # Step 4: Show all plots
+    import matplotlib.pyplot as plt  # For plotting
+    plt.show()
+
+# This allows the module to be executed directly 
+if __name__ == '__main__':
+    run_example()
diff --git a/examples/benchmarking_example.py b/examples/benchmarking_example.py
index 1f1679fb..9831f931 100644
--- a/examples/benchmarking_example.py
+++ b/examples/benchmarking_example.py
@@ -17,25 +17,33 @@
 # VVV Uncomment this to use Electro Chemistry Model VVV
 # from prog_models.models import BatteryElectroChem as Battery
 
-from prog_algs import *
+from prog_algs import state_estimators, predictors
+from prog_algs.metrics import samples as metrics 
+import time  # For timing prediction
 
 def run_example():
     # Step 1: Setup Model and Future Loading
-    def future_loading(t, x={}):
+    batt = Battery()
+    
+    # Creating the input containers outside of the function accelerates prediction
+    loads = [
+        batt.InputContainer({'i': 2}),
+        batt.InputContainer({'i': 1}),
+        batt.InputContainer({'i': 4}),
+        batt.InputContainer({'i': 2}),
+        batt.InputContainer({'i': 3})
+    ]
+    def future_loading(t, x = None):
         # Variable (piece-wise) future loading scheme 
         if (t < 600):
-            i = 2
+            return loads[0]
         elif (t < 900):
-            i = 1
+            return loads[1]
         elif (t < 1800):
-            i = 4
+            return loads[2]
         elif (t < 3000):
-            i = 2
-        else:
-            i = 3
-        return {'i': i}
-
-    batt = Battery()   
+            return loads[3]
+        return loads[-1]
 
     # Step 2: Setup Predictor 
     pred = predictors.MonteCarlo(batt, dt= 0.05)
@@ -44,14 +52,14 @@ def future_loading(t, x={}):
     x0 = batt.initialize()
     state_estimator = state_estimators.ParticleFilter(batt, x0)
     # Send in some data to estimate state
-    state_estimator.estimate(0.1, future_loading(0.1), {'t': 32.2, 'v': 3.915})
-    state_estimator.estimate(0.2, future_loading(0.2), {'t': 32.3, 'v': 3.91})
+    z1 = batt.OutputContainer({'t': 32.2, 'v': 3.915})
+    z2 = batt.OutputContainer({'t': 32.3, 'v': 3.91})
+    state_estimator.estimate(0.1, future_loading(0.1), z1)
+    state_estimator.estimate(0.2, future_loading(0.2), z2)
 
     # Step 4: Benchmark Predictions
     # Here we're comparing the results given different numbers of samples
     print('Benchmarking...')
-    import time  # For timing prediction
-    from prog_algs.metrics import samples as metrics 
 
     # Perform benchmarking for each number of samples
     sample_counts = [1, 2, 5, 10]
diff --git a/examples/eol_event.py b/examples/eol_event.py
index 82b1eaf6..0e5c8e39 100644
--- a/examples/eol_event.py
+++ b/examples/eol_event.py
@@ -3,12 +3,15 @@
 This example demonstrates a use case where someone wants to predict the first event (i.e., End Of Life (EOL)) of a system. Many system models have multiple events that can occur. In some prognostics applications, users are not interested in predicting a specific event, and are instead interested in when the first event occurs, regardless of the event. This example demonstrates how to predict the first event of a system.
 
 Method: An instance of ThrownObject is used for this example. In this case it is trivial because the event 'falling' will always occur before 'impact', but for some other models that might not be true. The ThrownObject class is subclassed to add a new event 'EOL' which occurs if any other event occurs. The new model is then instantiated and used for prognostics like in basic_example. Prediction specifically specifies EOL as the event to be predicted.
+
 Results: 
+
     i) Predicted future values (inputs, states, outputs, event_states) with uncertainty from prediction
     ii) Time the event 'EOL' is predicted to occur (with uncertainty)
     iii) Histogram of the event 'EOL'
 """
 
+import matplotlib.pyplot as plt
 from prog_models.models import ThrownObject
 from prog_algs.predictors import MonteCarlo
 from prog_algs.uncertain_data import ScalarData
@@ -32,7 +35,7 @@ def threshold_met(self, x):
             return t_met
     
     # Step 2: Create instance of subclass
-    m = ThrownObjectWithEOL(process_noise=0.05)
+    m = ThrownObjectWithEOL(process_noise=1)
 
     # Step 3: Setup for prediction
     pred = MonteCarlo(m)
@@ -41,12 +44,12 @@ def future_loading(t=None, x=None):
     state = ScalarData(m.initialize())
 
     # Step 4: Predict to EOL event
-    pred_results = pred.predict(state, future_loading, events=['EOL'], dt=0.025)
+    pred_results = pred.predict(state, future_loading, events=['EOL'], dt=0.01, n_samples=50)
+    # In this case EOL is when the object starts falling
+    # But for some models where events aren't sequential, there might be a mixture of events in the EOL
 
     # Step 5: Plot results
     pred_results.time_of_event.plot_hist()
-
-    import matplotlib.pyplot as plt
     plt.show()
 
 # This allows the module to be executed directly 
diff --git a/examples/horizon.py b/examples/horizon.py
index 1442a22c..bd5e034b 100644
--- a/examples/horizon.py
+++ b/examples/horizon.py
@@ -6,17 +6,14 @@
 Method: An instance of the Thrown Object model in prog_models is created, and the prediction process is achieved in three steps:
     1) State estimation of the current state is performed using a chosen state_estimator, and samples are drawn from this estimate
     2) Prediction of future states (with uncertainty) and the times at which the event thresholds will be reached, within the prediction horizon. All events outside the horizon come back as None and are ignored in metrics
+
 Results: 
     i) Predicted future values (inputs, states, outputs, event_states) with uncertainty from prediction
     ii) Time event is predicted to occur (with uncertainty)
 """
 
 from prog_models.models.thrown_object import ThrownObject
-# VVV Uncomment this to use Electro Chemistry Model VVV
-# from prog_models.models import BatteryElectroChem as Battery
-
 from prog_algs import *
-
 from pprint import pprint
 
 def run_example():
@@ -24,7 +21,7 @@ def run_example():
     def future_loading(t, x = None):
         return {}
     m = ThrownObject(process_noise = 0.25, measurement_noise = 0.2)
-    initial_state = m.initialize({}, {})
+    initial_state = m.initialize()
 
     # Step 2: Demonstrating state estimator
     print("\nPerforming State Estimation Step...")
@@ -35,11 +32,13 @@ def future_loading(t, x = None):
     # VVV Uncomment this to use UKF State Estimator VVV
     # filt = state_estimators.UnscentedKalmanFilter(batt, initial_state)
 
-    # Step 2b: Print & Plot Prior State
-    filt.estimate(0.1, {}, m.output(initial_state))
+    # Step 2b: One step of state estimator
+    u = m.InputContainer({})  # No input for ThrownObject
+    filt.estimate(0.1, u, m.output(initial_state))
 
-    # Note: in a prognostic application the above state estimation step would be repeated each time
-    #   there is new data. Here we're doing one step to demonstrate how the state estimator is used
+    # Note: in a prognostic application the above state estimation 
+    # step would be repeated each time there is new data. 
+    # Here we're doing one step to demonstrate how the state estimator is used
 
     # Step 3: Demonstrating Predictor
     print("\nPerforming Prediction Step...")
diff --git a/examples/kalman_filter.py b/examples/kalman_filter.py
index eab72e2d..a914f8c2 100644
--- a/examples/kalman_filter.py
+++ b/examples/kalman_filter.py
@@ -7,10 +7,10 @@
 """
 
 import numpy as np
-
 from prog_models import LinearModel
 from prog_algs.state_estimators import KalmanFilter
 
+
 # Linear Model for an object thrown into the air
 class ThrownObject(LinearModel):
     """
@@ -68,7 +68,8 @@ class ThrownObject(LinearModel):
     C = np.array([[1, 0]])
     F = None # Will override method
 
-    # The Default parameters. Overwritten by passing parameters dictionary into constructor
+    # The Default parameters. 
+    # Overwritten by passing parameters dictionary into constructor
     default_parameters = {
         'thrower_height': 1.83,  # m
         'throwing_speed': 40,  # m/s
@@ -77,12 +78,15 @@ class ThrownObject(LinearModel):
 
     def initialize(self, u=None, z=None):
         return self.StateContainer({
-            'x': self.parameters['thrower_height'],  # Thrown, so initial altitude is height of thrower
-            'v': self.parameters['throwing_speed']  # Velocity at which the ball is thrown - this guy is a professional baseball pitcher
+            'x': self.parameters['thrower_height'],
+            # Thrown, so initial altitude is height of thrower
+            'v': self.parameters['throwing_speed']
+            # Velocity at which the ball is thrown - this guy is a professional baseball pitcher
             })
     
     # This is actually optional. Leaving thresholds_met empty will use the event state to define thresholds.
-    #  Threshold = Event State == 0. However, this implementation is more efficient, so we included it
+    #  Threshold is met when Event State == 0. 
+    # However, this implementation is more efficient, so we included it
     def threshold_met(self, x):
         return {
             'falling': x['v'] < 0,
@@ -102,11 +106,13 @@ def run_example():
 
     # Step 2: Instantiate the Kalman Filter State Estimator
     # Define the initial state to be slightly off of actual
-    x_guess = {'x': 1.75, 'v': 35} # Guess of initial state, actual is {'x': 1.83, 'v': 40}
+    x_guess = m.StateContainer({'x': 1.75, 'v': 35}) # Guess of initial state
+    # Note: actual is {'x': 1.83, 'v': 40}
     kf = KalmanFilter(m, x_guess)
 
     # Step 3: Run the Kalman Filter State Estimator
-    # Here we're using simulated data from the thrown_object. In a real application you would be using sensor data from the system
+    # Here we're using simulated data from the thrown_object. 
+    # In a real application you would be using sensor data from the system
     dt = 0.01  # Time step (s)
     print_freq = 50  # Print every print_freq'th iteration
     x = m.initialize()
diff --git a/examples/measurement_eqn_example.py b/examples/measurement_eqn_example.py
index cd2fdca4..6b9674e8 100644
--- a/examples/measurement_eqn_example.py
+++ b/examples/measurement_eqn_example.py
@@ -26,21 +26,26 @@ class MyBattery(Battery):
         outputs = ['v']
 
     # Step 2: Setup model & future loading
-    def future_loading(t, x={}):
+    batt = MyBattery()
+    loads = [  # Define loads here to accelerate prediction
+        batt.InputContainer({'i': 2}),
+        batt.InputContainer({'i': 1}),
+        batt.InputContainer({'i': 4}),
+        batt.InputContainer({'i': 2}),
+        batt.InputContainer({'i': 3})
+    ]
+    def future_loading(t, x = None):
         # Variable (piece-wise) future loading scheme 
         if (t < 600):
-            i = 2
+            return loads[0]
         elif (t < 900):
-            i = 1
+            return loads[1]
         elif (t < 1800):
-            i = 4
+            return loads[2]
         elif (t < 3000):
-            i = 2
-        else:
-            i = 3
-        return {'i': i}
+            return loads[3]
+        return loads[-1]
 
-    batt = MyBattery()
     x0 = batt.parameters['x0']
 
     # Step 3: Use the updated model
diff --git a/examples/new_state_estimator_example.py b/examples/new_state_estimator_example.py
index cc409ce0..374d2cd7 100644
--- a/examples/new_state_estimator_example.py
+++ b/examples/new_state_estimator_example.py
@@ -10,6 +10,7 @@
 from prog_algs.uncertain_data import ScalarData
 import random
 
+
 class BlindlyStumbleEstimator(StateEstimator):
     """
     A new state estimator. This is not a very effective state estimator, but one that technically works. It blindly stumbles towards the correct state by randomly generating a new state each timestep and selecting the state that's most consistant with the measurements.
diff --git a/examples/particle_filter_battery_example.py b/examples/particle_filter_battery_example.py
index 4b9a9856..68168e9a 100644
--- a/examples/particle_filter_battery_example.py
+++ b/examples/particle_filter_battery_example.py
@@ -4,16 +4,13 @@
 In this example the BatteryElectroChemEOD model is used with a particle filter to estimate the state of the battery
 """
 
-from prog_models.models import BatteryElectroChemEOD
-from prog_algs import *
-import numpy as np
 import matplotlib.pyplot as plt
+import numpy as np
+from prog_algs import *
+from prog_models.models import BatteryElectroChemEOD
 
 def run_example():
     ## Setup
-    def future_loading(t, x=None):
-        load = 1
-        return {"i": load}
     # Save battery model
     # Time increment
     dt = 1
@@ -36,6 +33,9 @@ def future_loading(t, x=None):
     battery = BatteryElectroChemEOD(process_noise= Q_vars,
                                     measurement_noise = R_vars, 
                                     dt = dt)
+    load = battery.InputContainer({"i": 1})  # Optimization
+    def future_loading(t, x=None):
+        return load
 
     # Simulate data until EOD
     start_u = future_loading(0)
diff --git a/examples/playback.py b/examples/playback.py
index bc982487..638baf78 100644
--- a/examples/playback.py
+++ b/examples/playback.py
@@ -16,6 +16,11 @@
     iv) Figures illustrating results
 """
 
+import csv
+import numpy as np
+from prog_algs.predictors import ToEPredictionProfile
+from prog_algs.uncertain_data.multivariate_normal_dist import MultivariateNormalDist
+
 from prog_models.models import BatteryCircuit as Battery
 # VVV Uncomment this to use Electro Chemistry Model VVV
 # from prog_models.models import BatteryElectroChem as Battery
@@ -24,18 +29,10 @@
 # VVV Uncomment this to use UnscentedKalmanFilter instead VVV
 # from prog_algs.state_estimators import ParticleFilter as StateEstimator
 
-from prog_algs.predictors import ToEPredictionProfile
-
 from prog_algs.predictors import UnscentedTransformPredictor as Predictor
 # VVV Uncomment this to use MonteCarloPredictor instead
 # from prog_algs.predictors import MonteCarlo as Predictor
 
-from prog_algs.uncertain_data.multivariate_normal_dist import MultivariateNormalDist
-
-import csv
-import matplotlib.pyplot as plt
-import numpy as np
-
 # Constants
 NUM_SAMPLES = 20
 NUM_PARTICLES = 1000 # For state estimator (if using ParticleFilter)
@@ -65,8 +62,9 @@ def run_example():
     filt = StateEstimator(batt, x0, num_particles = NUM_PARTICLES)
 
     # Setup Prediction
+    load = batt.InputContainer({'i': 2.35})
     def future_loading(t, x=None):
-        return {'i': 2.35}
+        return load
     Q = np.diag([batt.parameters['process_noise'][key] for key in batt.states])
     R = np.diag([batt.parameters['measurement_noise'][key] for key in batt.outputs])
     mc = Predictor(batt, Q = Q, R = R)
diff --git a/examples/predict_specific_event.py b/examples/predict_specific_event.py
index 053a179f..816a3c67 100644
--- a/examples/predict_specific_event.py
+++ b/examples/predict_specific_event.py
@@ -4,16 +4,16 @@
 In this example we are using the UTPredictor to predict a specific event, in this case impact. This will then ignore the other events which are not of interest.
 """
 
+from prog_algs import state_estimators, predictors
 from prog_models.models.thrown_object import ThrownObject
-from prog_algs import *
 
 def run_example():
     ## Setup
-    def future_loading(t, x = None):
-        return {}
-
     m = ThrownObject()
-    initial_state = m.initialize({}, {})
+    initial_state = m.initialize()
+    load = m.InputContainer({}) # Optimization - create once
+    def future_loading(t, x = None):
+        return load
 
     ## State Estimation - perform a single ukf state estimate step
     filt = state_estimators.UnscentedKalmanFilter(m, initial_state)
diff --git a/examples/thrown_object_example.py b/examples/thrown_object_example.py
deleted file mode 100644
index 54c8ac22..00000000
--- a/examples/thrown_object_example.py
+++ /dev/null
@@ -1,65 +0,0 @@
-# Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration. All Rights Reserved.
-
-"""
-This example performs a state estimation and prediction with uncertainty given a Prognostics Model. Unlike basic_example, this example uses a model with multiple events (ThrownObject). Prediction only ends when all events are met
- 
-Method: An instance of the Thrown Object model in prog_models is created, and the prediction process is achieved in three steps:
-    1) State estimation of the current state is performed using a chosen state_estimator, and samples are drawn from this estimate
-    2) Prediction of future states (with uncertainty) and the times at which the event thresholds will be reached
-Results: 
-    i) Predicted future values (inputs, states, outputs, event_states) with uncertainty from prediction
-    ii) Time event is predicted to occur (with uncertainty)
-"""
-
-from prog_models.models.thrown_object import ThrownObject
-# VVV Uncomment this to use Electro Chemistry Model VVV
-# from prog_models.models import BatteryElectroChem as Battery
-
-from prog_algs import *
-
-from pprint import pprint
-
-def run_example():
-    # Step 1: Setup model & future loading
-    def future_loading(t, x = None):
-        return {}
-    m = ThrownObject(process_noise = 0.25, measurement_noise = 0.2)
-    initial_state = m.initialize({}, {})
-
-    # Step 2: Demonstrating state estimator
-    print("\nPerforming State Estimation Step...")
-
-    # Step 2a: Setup
-    NUM_SAMPLES = 1000
-    filt = state_estimators.ParticleFilter(m, initial_state, num_particles = NUM_SAMPLES)
-    # VVV Uncomment this to use UKF State Estimator VVV
-    # filt = state_estimators.UnscentedKalmanFilter(batt, initial_state)
-
-    # Step 2b: Print & Plot Prior State
-    filt.estimate(0.1, {}, m.output(initial_state))
-
-    # Note: in a prognostic application the above state estimation step would be repeated each time
-    #   there is new data. Here we're doing one step to demonstrate how the state estimator is used
-
-    # Step 3: Demonstrating Predictor
-    print("\nPerforming Prediction Step...")
-
-    # Step 3a: Setup Predictor
-    mc = predictors.MonteCarlo(m)
-
-    # Step 3b: Perform a prediction
-    samples = filt.x  # Since we're using a particle filter, which is also sample-based, we can directly use the samples, without changes
-    STEP_SIZE = 0.01
-    mc_results = mc.predict(samples, future_loading, dt=STEP_SIZE, horizon=8)
-    print("\nPredicted Time of Event:")
-    pprint(mc_results.time_of_event.metrics())  # Note this takes some time
-    mc_results.time_of_event.plot_hist(keys = 'impact')
-    mc_results.time_of_event.plot_hist(keys = 'falling')
-    
-    # Step 4: Show all plots
-    import matplotlib.pyplot as plt  # For plotting
-    plt.show()
-
-# This allows the module to be executed directly 
-if __name__ == '__main__':
-    run_example()
diff --git a/examples/utpredictor.py b/examples/utpredictor.py
deleted file mode 100644
index 3ec1544a..00000000
--- a/examples/utpredictor.py
+++ /dev/null
@@ -1,67 +0,0 @@
-# Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration. All Rights Reserved.
-
-"""
-An example using the UnscentedTransformPredictor class to predict the degredation of a battery. 
-"""
-
-from prog_models.models import BatteryCircuit
-from prog_algs import *
-# from prog_algs.visualize import plot_hist
-# import matplotlib.pyplot as plt
-
-def run_example():
-    ## Setup
-    def future_loading(t, x = None):
-        # Variable (piece-wise) future loading scheme 
-        if (t < 600):
-            i = 2
-        elif (t < 900):
-            i = 1
-        elif (t < 1800):
-            i = 4
-        elif (t < 3000):
-            i = 2
-        else:
-            i = 3
-        return {'i': i}
-
-    batt = BatteryCircuit()
-
-    ## State Estimation - perform a single ukf state estimate step
-    filt = state_estimators.UnscentedKalmanFilter(batt, batt.parameters['x0'])
-
-    import matplotlib.pyplot as plt  # For plotting
-    print("Prior State:", filt.x.mean)
-    print('\tSOC: ', batt.event_state(filt.x.mean)['EOD'])
-    example_measurements = {'t': 32.2, 'v': 3.915}
-    t = 0.1
-    filt.estimate(t, future_loading(t), example_measurements)
-    print("Posterior State:", filt.x.mean)
-    print('\tSOC: ', batt.event_state(filt.x.mean)['EOD'])
-
-    ## Prediction - Predict EOD given current state
-    # Setup prediction
-    mc = predictors.UnscentedTransformPredictor(batt)
-
-    # Predict with a step size of 0.1
-    mc_results = mc.predict(filt.x, future_loading, dt=0.1, save_freq= 100)
-
-    # Print Results
-    for i, time in enumerate(mc_results.times):
-        print('\nt = {}'.format(time))
-        print('\tu = {}'.format(mc_results.inputs.snapshot(i).mean))
-        print('\tx = {}'.format(mc_results.states.snapshot(i).mean))
-        print('\tz = {}'.format(mc_results.outputs.snapshot(i).mean))
-        print('\tevent state = {}'.format(mc_results.event_states.snapshot(i).mean))
-
-    print('\nToE:', mc_results.time_of_event.mean)
-
-    # You can also access the final state (of type UncertainData), like so:
-    final_state = mc_results.time_of_event.final_state
-    print('Final state @EOD: ', final_state['EOD'].mean)
-    # toe.plot_hist()
-    # plt.show()
-
-# This allows the module to be executed directly 
-if __name__ == '__main__':
-    run_example()
diff --git a/predictor_template.py b/predictor_template.py
index 33cc3452..595b3b6b 100644
--- a/predictor_template.py
+++ b/predictor_template.py
@@ -1,6 +1,7 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration.  All Rights Reserved.
 
 from copy import deepcopy
+
 from prog_algs.predictors import Predictor, Prediction, PredictionResults
 # Replace the following with whatever form of UncertainData you would like to use to represent ToE 
 from prog_algs.uncertain_data import ScalarData
diff --git a/setup.py b/setup.py
index 0b9c063f..9f247a03 100644
--- a/setup.py
+++ b/setup.py
@@ -10,11 +10,11 @@
 
 setup(
     name = 'prog_algs',
-    version = '1.3.1',
+    version = '1.4.0',
     description = "The NASA Prognostics Algorithm Package is a framework for model-based prognostics (computation of remaining useful life) of engineering systems. It includes algorithms for state estimation and prediction, including uncertainty propagation. The algorithms use prognostic models (see prog_models) to perform estimation and prediction. The package enables rapid development of prognostics solutions for given models of components and systems. Algorithms can be swapped for comparative studies and evaluations",
     long_description=long_description,
     long_description_content_type='text/markdown',
-    url = 'https://github.com/nasa/prog_algs',
+    url = 'https://nasa.github.io/progpy/prog_algs_guide.html',
     author = 'Christopher Teubert',
     author_email = 'christopher.a.teubert@nasa.gov',
     classifiers = [
@@ -33,7 +33,7 @@
         'Programming Language :: Python :: 3.10',
         'Programming Language :: Python :: 3 :: Only'
     ],
-    keywords = ['prognostics', 'diagnostics', 'fault detection', 'fdir', 'prognostics and health management', 'PHM', 'health management'],
+    keywords = ['prognostics', 'diagnostics', 'fault detection', 'fdir', 'prognostics and health management', 'PHM', 'health management', 'ivhm'],
     package_dir = {"":"src"},
     packages = find_packages(where = 'src'),
     python_requires='>=3.7, <3.11',
@@ -47,7 +47,8 @@
     license = 'NOSA',
     project_urls={  # Optional
         'Bug Reports': 'https://github.com/nasa/prog_algs/issues',
-        'Organization': 'https://prognostics.nasa.gov/',
+        'Docs': 'https://nasa.github.io/progpy/prog_algs_guide.html',
+        'Organization': 'https://www.nasa.gov/content/diagnostics-prognostics',
         'Source': 'https://github.com/nasa/prog_algs',
     }
-)
\ No newline at end of file
+)
diff --git a/sphinx_config/_templates/layout.html b/sphinx_config/_templates/layout.html
deleted file mode 100644
index 8c1a51a3..00000000
--- a/sphinx_config/_templates/layout.html
+++ /dev/null
@@ -1,6 +0,0 @@
-{% extends '!layout.html' %}
-
-{% block footer %}
-
-Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration.  All Rights Reserved.
-{% endblock %}
diff --git a/sphinx_config/conf.py b/sphinx_config/conf.py
deleted file mode 100644
index e4cf31c2..00000000
--- a/sphinx_config/conf.py
+++ /dev/null
@@ -1,65 +0,0 @@
-# Copyright © 2021 United States Government as represented by the Administrator of the
-# National Aeronautics and Space Administration.  All Rights Reserved.
-
-# Configuration file for the Sphinx documentation builder.
-#
-# This file only contains a selection of the most common options. For a full
-# list see the documentation:
-# https://www.sphinx-doc.org/en/master/usage/configuration.html
-
-# -- Path setup --------------------------------------------------------------
-
-# If extensions (or modules to document with autodoc) are in another directory,
-# add these directories to sys.path here. If the directory is relative to the
-# documentation root, use os.path.abspath to make it absolute, like shown here.
-
-import os
-import sys
-sys.path.insert(0, os.path.abspath('../src'))
-sys.path.insert(0, os.path.abspath('..'))
-sys.path.insert(0, os.path.abspath('.'))
-
-
-# -- Project information -----------------------------------------------------
-
-project = 'Prognostics Algorithms Python Package'
-copyright = '2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration.  All Rights Reserved.'
-author = 'Chris Teubert, Matteo Corbetta, and Chetan Kulkarni'
-
-# The full version, including alpha/beta/rc tags
-from prog_algs import __version__ as version
-release = version
-
-# -- General configuration ---------------------------------------------------
-
-# Add any Sphinx extension module names here, as strings. They can be
-# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
-# ones.
-extensions = ['sphinx.ext.autodoc', 'sphinx.ext.coverage', "sphinx.ext.githubpages", 'sphinx.ext.napoleon', 'sphinx_reredirects']
-
-# Add any paths that contain templates here, relative to this directory.
-templates_path = ['_templates']
-
-# List of patterns, relative to source directory, that match files and
-# directories to ignore when looking for source files.
-# This pattern also affects html_static_path and html_extra_path.
-exclude_patterns = []
-
-autodoc_mock_imports = ['numpy', 'scipy', 'filterpy', 'prog_models']
-
-
-# -- Options for HTML output -------------------------------------------------
-
-# The theme to use for HTML and HTML Help pages.  See the documentation for
-# a list of builtin themes.
-#
-html_theme = 'alabaster'
-
-# Add any paths that contain custom static files (such as style sheets) here,
-# relative to this directory. They are copied after the builtin static files,
-# so a file named "default.css" will overwrite the builtin "default.css".
-html_static_path = ['_static']
-
-redirects = {
-     "*": "https://nasa.github.io/progpy"
-}
diff --git a/sphinx_config/getting_started.rst b/sphinx_config/getting_started.rst
deleted file mode 100644
index 0a95e38b..00000000
--- a/sphinx_config/getting_started.rst
+++ /dev/null
@@ -1,4 +0,0 @@
-Getting Started
-===============
-
-Documentation moved to https://nasa.github.io/progpy
diff --git a/sphinx_config/images/alpha_lambda_plt_ex.png b/sphinx_config/images/alpha_lambda_plt_ex.png
deleted file mode 100644
index c489c48f..00000000
Binary files a/sphinx_config/images/alpha_lambda_plt_ex.png and /dev/null differ
diff --git a/sphinx_config/images/package_structure.png b/sphinx_config/images/package_structure.png
deleted file mode 100644
index e5b56af1..00000000
Binary files a/sphinx_config/images/package_structure.png and /dev/null differ
diff --git a/sphinx_config/images/scatter_hist_ex.png b/sphinx_config/images/scatter_hist_ex.png
deleted file mode 100644
index b4e9907a..00000000
Binary files a/sphinx_config/images/scatter_hist_ex.png and /dev/null differ
diff --git a/sphinx_config/images/scatter_plt_ex.png b/sphinx_config/images/scatter_plt_ex.png
deleted file mode 100644
index 3352babc..00000000
Binary files a/sphinx_config/images/scatter_plt_ex.png and /dev/null differ
diff --git a/sphinx_config/index.rst b/sphinx_config/index.rst
deleted file mode 100644
index e2b86f10..00000000
--- a/sphinx_config/index.rst
+++ /dev/null
@@ -1,4 +0,0 @@
-Prognostics Algorithms Python Package
-=============================================================
-
-Documentation moved to https://nasa.github.io/progpy
diff --git a/sphinx_config/metrics.rst b/sphinx_config/metrics.rst
deleted file mode 100644
index d884a27e..00000000
--- a/sphinx_config/metrics.rst
+++ /dev/null
@@ -1,4 +0,0 @@
-Metrics and Figures
-====================
-
-Documentation moved to https://nasa.github.io/progpy
diff --git a/sphinx_config/prediction.rst b/sphinx_config/prediction.rst
deleted file mode 100644
index 22209253..00000000
--- a/sphinx_config/prediction.rst
+++ /dev/null
@@ -1,4 +0,0 @@
-Prediction
-=======================
-
-Documentation moved to https://nasa.github.io/progpy
diff --git a/sphinx_config/predictors.rst b/sphinx_config/predictors.rst
deleted file mode 100644
index c61b4be0..00000000
--- a/sphinx_config/predictors.rst
+++ /dev/null
@@ -1,4 +0,0 @@
-Predictors
-=======================
-
-Documentation moved to https://nasa.github.io/progpy
diff --git a/sphinx_config/release.rst b/sphinx_config/release.rst
deleted file mode 100644
index d0980713..00000000
--- a/sphinx_config/release.rst
+++ /dev/null
@@ -1,4 +0,0 @@
-Release Notes
-=======================
-
-Documentation moved to https://nasa.github.io/progpy
diff --git a/sphinx_config/state_estimators.rst b/sphinx_config/state_estimators.rst
deleted file mode 100644
index 2a7c2bc1..00000000
--- a/sphinx_config/state_estimators.rst
+++ /dev/null
@@ -1,4 +0,0 @@
-State Estimators
-=======================
-
-Documentation moved to https://nasa.github.io/progpy
diff --git a/sphinx_config/uncertain_data.rst b/sphinx_config/uncertain_data.rst
deleted file mode 100644
index 6f960fa0..00000000
--- a/sphinx_config/uncertain_data.rst
+++ /dev/null
@@ -1,4 +0,0 @@
-Uncertain Data
-=======================
-
-Documentation moved to https://nasa.github.io/progpy
diff --git a/src/prog_algs/__init__.py b/src/prog_algs/__init__.py
index 657faf71..0b6e1e61 100644
--- a/src/prog_algs/__init__.py
+++ b/src/prog_algs/__init__.py
@@ -5,7 +5,7 @@
 
 import warnings
 
-__version__ = '1.3.1'
+__version__ = '1.4.0'
 
 def run_prog_playback(obs, pred, future_loading, output_measurements, **kwargs):
     warnings.warn("Depreciated in 1.2.0, will be removed in a future release.", DeprecationWarning)
diff --git a/src/prog_algs/metrics/__init__.py b/src/prog_algs/metrics/__init__.py
index eb52dd7a..4c39fe28 100644
--- a/src/prog_algs/metrics/__init__.py
+++ b/src/prog_algs/metrics/__init__.py
@@ -3,6 +3,6 @@
 from . import samples
 from .toe_metrics import prob_success
 from .uncertain_data_metrics import calc_metrics
-from .toe_profile_metrics import alpha_lambda, prognostic_horizon, cumulative_relative_accuracy
+from .toe_profile_metrics import alpha_lambda, prognostic_horizon, cumulative_relative_accuracy, monotonicity
 
-__all__ = ['alpha_lambda', 'calc_metrics', 'prob_success', 'prognostic_horizon', 'cumulative_relative_accuracy']
+__all__ = ['alpha_lambda', 'calc_metrics', 'prob_success', 'prognostic_horizon', 'cumulative_relative_accuracy', 'monotonicity']
diff --git a/src/prog_algs/metrics/samples.py b/src/prog_algs/metrics/samples.py
index 6dcb5f60..d72d3e1d 100644
--- a/src/prog_algs/metrics/samples.py
+++ b/src/prog_algs/metrics/samples.py
@@ -1,13 +1,13 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration.  All Rights Reserved.
 
 # This file is kept for backwards compatability
+from numpy import mean, sqrt
+from warnings import warn
+
 from .uncertain_data_metrics import calc_metrics as eol_metrics
 from .toe_metrics import prob_success
 from .toe_profile_metrics import alpha_lambda
 
-from numpy import mean, sqrt
-from warnings import warn
-
 def mean_square_error(values : list, ground_truth : float) -> float:
     """Mean Square Error
     Args:
diff --git a/src/prog_algs/metrics/toe_metrics.py b/src/prog_algs/metrics/toe_metrics.py
index e75293c3..ec0e9603 100644
--- a/src/prog_algs/metrics/toe_metrics.py
+++ b/src/prog_algs/metrics/toe_metrics.py
@@ -3,9 +3,9 @@
 """
 This file includes functions for calculating metrics specific to Time of Event (ToE) from a single event or multiple events given the same time of prediction
 """
-
 from typing import Iterable
 from numpy import isscalar
+
 from ..uncertain_data import UncertainData, UnweightedSamples
 
 def prob_success(toe : UncertainData, time : float, **kwargs) -> float:
@@ -42,6 +42,8 @@ def prob_success(toe : UncertainData, time : float, **kwargs) -> float:
     if isinstance(toe, UncertainData):
         # Default to all keys
         keys = params.setdefault('keys', toe.keys())
+        if isinstance(keys, str):
+            keys = [keys]
 
         if isinstance(toe, UnweightedSamples):
             samples = toe
diff --git a/src/prog_algs/metrics/toe_profile_metrics.py b/src/prog_algs/metrics/toe_profile_metrics.py
index 794c9401..3f61ee28 100644
--- a/src/prog_algs/metrics/toe_profile_metrics.py
+++ b/src/prog_algs/metrics/toe_profile_metrics.py
@@ -3,10 +3,11 @@
 """
 This file includes functions for calculating metrics given a Time of Event (ToE) profile (i.e., ToE's calculated at different times of prediction resulting from running prognostics multiple times, e.g., on playback data). The metrics calculated here are specific to multiple ToE estimates (e.g. alpha-lambda metric)
 """
+from numpy import sign
 from collections import defaultdict
-from typing import Callable
+from typing import Callable, Dict
+
 from ..predictors import ToEPredictionProfile
-from typing import Dict
 
 def alpha_lambda(toe_profile : ToEPredictionProfile, ground_truth : dict, lambda_value : float, alpha : float, beta : float, **kwargs) -> dict: 
     """
@@ -107,4 +108,33 @@ def cumulative_relative_accuracy(toe_profile : ToEPredictionProfile, ground_trut
             ra_sums[event] += value
     return {event:ra_sum/len(toe_profile) for event, ra_sum in ra_sums.items()}
 
+def monotonicity(toe_profile : ToEPredictionProfile, **kwargs) -> Dict[str, float]:
+        """Calculate monotonicty for a prediction profile. 
+        Given a prediction profile, for each prediction: go through all predicted events and compare those to the next one.
+        Calculates monotonicity for each prediction key using its associated mean value in UncertainData.
+        
+        monotonoicity = |Σsign(i+1 - i) / N-1|
+        Where N is number of measurements and sign indicates sign of calculation.
+        Coble, J., et. al. (2021). Identifying Optimal Prognostic Parameters from Data: A Genetic Algorithms Approach. Annual Conference of the PHM Society.
+        http://www.papers.phmsociety.org/index.php/phmconf/article/view/1404
+        Baptistia, M., et. al. (2022). Relation between prognostics predictor evaluation metrics and local interpretability SHAP values. Aritifical Intelligence, Volume 306.
+        https://www.sciencedirect.com/science/article/pii/S0004370222000078
 
+        Args:
+            toe_profile (ToEPredictionProfile): A profile of predictions, the combination of multiple predictions
+        Returns:
+            dict (str, float): Dictionary where keys represent an event and values are float representing its respective monotonicitiy value between [0, 1].
+        """
+        result = dict()
+        by_event = defaultdict(list)
+        for time,uncertaindata in toe_profile.items():
+            # Collect and organize mean values for each event in the individual prediction v
+            for event,value in uncertaindata.mean.items():
+                by_event[event].append(value - time)
+        # For each event of this prediction v, calculate monotonicity using formula
+        for key,l in by_event.items():
+            mono_sum = []
+            for i in range(len(l)-1): 
+                mono_sum.append(sign(l[i+1] - l[i])) 
+            result[key] = abs(sum(mono_sum) / (len(l)-1))
+        return result
diff --git a/src/prog_algs/metrics/uncertain_data_metrics.py b/src/prog_algs/metrics/uncertain_data_metrics.py
index acfaeb0d..948b854c 100644
--- a/src/prog_algs/metrics/uncertain_data_metrics.py
+++ b/src/prog_algs/metrics/uncertain_data_metrics.py
@@ -7,6 +7,7 @@
 from numpy import isscalar, mean, std, array
 from scipy import stats
 from warnings import warn
+
 from ..uncertain_data import UncertainData, UnweightedSamples
 
 def calc_metrics(data : UncertainData, ground_truth : Union[float, dict] = None, **kwargs) -> dict:
@@ -30,6 +31,8 @@ def calc_metrics(data : UncertainData, ground_truth : Union[float, dict] = None,
     if isinstance(data, UncertainData):
         # Default to all keys
         keys = params.setdefault('keys', data.keys())
+        if isinstance(keys, str):
+            keys = [keys]
         
         if ground_truth and isscalar(ground_truth):
             # If ground truth is scalar, create dict (expected below)
diff --git a/src/prog_algs/predictors/monte_carlo.py b/src/prog_algs/predictors/monte_carlo.py
index a18d77c7..8111089f 100644
--- a/src/prog_algs/predictors/monte_carlo.py
+++ b/src/prog_algs/predictors/monte_carlo.py
@@ -1,12 +1,13 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration.  All Rights Reserved.
 
-from .prediction import UnweightedSamplesPrediction, PredictionResults
-from .predictor import Predictor
 from copy import deepcopy
 from typing import Callable
 from prog_models.sim_result import SimResult, LazySimResult
 from prog_algs.uncertain_data import UnweightedSamples, UncertainData
 
+from .prediction import UnweightedSamplesPrediction, PredictionResults
+from .predictor import Predictor
+
 
 class MonteCarlo(Predictor):
     """
@@ -22,7 +23,7 @@ class MonteCarlo(Predictor):
         Initial time at which prediction begins, e.g., 0
     dt : float
         Simulation step size (s), e.g., 0.1
-    events : List[string]
+    events : list[str]
         Events to predict (subset of model.events) e.g., ['event1', 'event2']
     horizon : float
         Prediction horizon (s)
@@ -30,7 +31,7 @@ class MonteCarlo(Predictor):
         Number of samples to use. If not specified, a default value is used. If state is type UnweightedSamples and n_samples is not provided, the provided unweighted samples will be used directly.
     save_freq : float
         Frequency at which results are saved (s)
-    save_pts : List[float]
+    save_pts : list[float]
         Any additional savepoints (s) e.g., [10.1, 22.5]
     """
 
@@ -49,11 +50,16 @@ def predict(self, state : UncertainData, future_loading_eqn : Callable, **kwargs
 
         params = deepcopy(self.parameters) # copy parameters
         params.update(kwargs) # update for specific run
+        params['print'] = False
+        params['progress'] = False
+
+        if len(params['events']) == 0 and 'horizon' not in params:
+            raise ValueError("If specifying no event (i.e., simulate to time), must specify horizon")
 
         # Sample from state if n_samples specified or state is not UnweightedSamples
-        state = state.sample(params['n_samples'])
+        if not isinstance(state, UnweightedSamples) or len(state) != params['n_samples']:
+            state = state.sample(params['n_samples'])
 
-        ouput_eqn = self.model.output
         es_eqn = self.model.event_state
         tm_eqn = self.model.threshold_met
         simulate_to_threshold = self.model.simulate_to_threshold
@@ -67,95 +73,97 @@ def predict(self, state : UncertainData, future_loading_eqn : Callable, **kwargs
         event_states_all = []
 
         # Perform prediction
+        t0 = params.get('t0', 0)
         for x in state:
-            events_remaining = params['events'].copy()
-            first_output = ouput_eqn(x)
+            first_output = self.model.output(x)
             
             time_of_event = {}
             last_state = {}
-            times = []
-            inputs = SimResult()
-            states = SimResult()
-            outputs = LazySimResult(fcn = ouput_eqn)
-            event_states = LazySimResult(fcn = es_eqn)
 
-            t0 = params.get('t0', 0)
+            params['t0'] = t0
+            params['x'] = x
+
             if 'save_freq' in params and not isinstance(params['save_freq'], tuple):
-                params['save_freq'] = (t0, params['save_freq'])
+                params['save_freq'] = (params['t0'], params['save_freq'])
             
-
-            # Non-vectorized prediction
-            while len(events_remaining) > 0:  # Still events to predict
-                (t, u, xi, z, es) = simulate_to_threshold(future_loading_eqn, first_output, 
-                    **{**params, 'threshold_keys': events_remaining, 'print': False, 'progress': False, 't0': t0, 'x': x}  # Merge then separate
+            if len(params['events']) == 0:  # Predict to time
+                (times, inputs, states, outputs, event_states) = simulate_to_threshold(future_loading_eqn,       
+                    first_output, 
+                    threshold_keys = [],
+                    **params
                 )
-                if len(times) != 0:
-                    # is not the first iteration, therefore the first saved state is when the last event 
-                    t.pop(0)
-                    u.pop(0)
-                    xi.pop(0)
-                    z.pop(0)
-                    es.pop(0)
-
-                # Add results
-                times.extend(t)
-                inputs.extend(u)
-                states.extend(xi)
-                outputs.extend(z)
-                event_states.extend(es)
-
-                # Get which event occurs
-                t_met = tm_eqn(states[-1])
-                t_met = {key: t_met[key] for key in events_remaining}  # Only look at remaining keys
-                try:
-                    event = list(t_met.keys())[list(t_met.values()).index(True)]
-                except ValueError:
-                    # no event has occured
-                    for event in events_remaining:
-                        time_of_event[event] = None
-                        last_state[event] = None
-                    break
-
-                # An event has occured
-                time_of_event[event] = times[-1]
-                events_remaining.remove(event)  # No longer an event to predect to
-
-                # Remove last state (event)
-                t0 = times.pop()
-                inputs.pop()
-                x = states.pop()
-                last_state[event] = x.copy()
-                outputs.pop()
-                event_states.pop()
+            else:
+                events_remaining = params['events'].copy()
+
+                times = []
+                inputs = SimResult(_copy = False)
+                states = SimResult(_copy = False)
+                outputs = LazySimResult(fcn = self.model.output, _copy = False)
+                event_states = LazySimResult(fcn = es_eqn, _copy = False)
+
+                # Non-vectorized prediction
+                while len(events_remaining) > 0:  # Still events to predict
+                    (t, u, xi, z, es) = simulate_to_threshold(future_loading_eqn,       
+                        first_output, 
+                        threshold_keys = events_remaining,
+                        **params
+                    )
+
+                    # Add results
+                    times.extend(t)
+                    inputs.extend(u)
+                    states.extend(xi)
+                    outputs.extend(z, _copy = False)
+                    event_states.extend(es, _copy = False)
+
+                    # Get which event occurs
+                    t_met = tm_eqn(states[-1])
+                    t_met = {key: t_met[key] for key in events_remaining}  # Only look at remaining keys
+
+                    try:
+                        event = list(t_met.keys())[list(t_met.values()).index(True)]
+                    except ValueError:
+                        # no event has occured - hit horizon
+                        for event in events_remaining:
+                            time_of_event[event] = None
+                            last_state[event] = None
+                        break
+
+                    # An event has occured
+                    time_of_event[event] = times[-1]
+                    events_remaining.remove(event)  # No longer an event to predect to
+
+                    # Remove last state (event)
+                    params['t0'] = times.pop()
+                    inputs.pop()
+                    params['x'] = states.pop()
+                    last_state[event] = params['x'].copy()
+                    outputs.pop()
+                    event_states.pop()
             
             # Add to "all" structures
-            times_all.append(times)
+            if len(times) > len(times_all): # Keep longest
+                times_all = times
             inputs_all.append(inputs)
             states_all.append(states)
             outputs_all.append(outputs)
             event_states_all.append(event_states)
             time_of_event_all.append(time_of_event)
             last_states.append(last_state)
-        
-        # Return longest time array
-        times_length = [len(t) for t in times_all]
-        times_max_len = max(times_length)
-        times = times_all[times_length.index(times_max_len)] 
-        
-        inputs_all = UnweightedSamplesPrediction(times, inputs_all)
-        states_all = UnweightedSamplesPrediction(times, states_all)
-        outputs_all = UnweightedSamplesPrediction(times, outputs_all)
-        event_states_all = UnweightedSamplesPrediction(times, event_states_all)
+              
+        inputs_all = UnweightedSamplesPrediction(times_all, inputs_all)
+        states_all = UnweightedSamplesPrediction(times_all, states_all)
+        outputs_all = UnweightedSamplesPrediction(times_all, outputs_all)
+        event_states_all = UnweightedSamplesPrediction(times_all, event_states_all)
         time_of_event = UnweightedSamples(time_of_event_all)
 
         # Transform final states:
-        last_states = {
+        time_of_event.final_state = {
             key: UnweightedSamples([sample[key] for sample in last_states], _type = self.model.StateContainer) for key in time_of_event.keys()
         }
-        time_of_event.final_state = last_states
 
         return PredictionResults(
-            times, 
+            times_all, 
             inputs_all, 
             states_all, 
             outputs_all, 
diff --git a/src/prog_algs/predictors/prediction.py b/src/prog_algs/predictors/prediction.py
index 7060a77a..2c124b7b 100644
--- a/src/prog_algs/predictors/prediction.py
+++ b/src/prog_algs/predictors/prediction.py
@@ -1,7 +1,7 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration. All Rights Reserved.
 
 from collections import UserList, defaultdict, namedtuple
-from typing import Dict
+from typing import Dict, List
 from numpy import sign
 from warnings import warn
 
@@ -49,7 +49,7 @@ def snapshot(self, time_index : int) -> UncertainData:
         return self.data[time_index]
 
     @property
-    def mean(self) -> list:
+    def mean(self) -> List[dict]:
         """Estimate of the mean value of the prediction at each time
 
         Returns:
@@ -72,21 +72,17 @@ def monotonicity(self) -> Dict[str, float]:
         Given a single prediction, for each event: go through all predicted states and compare those to the next one.
         Calculates monotonicity for each event key using its associated mean value in UncertainData.
         
-        monotonoicity = |Σsign(i+1 - i) / N-1|
+        :math:`monotonoicity = \| \Sigma \dfrac{sign(i+1 - i)}{N-1}\|`
 
-        Where N is number of measurements and sign indicates sign of calculation.
-        
-        Coble, J., et. al. (2021). Identifying Optimal Prognostic Parameters from D
-        ata: A Genetic Algorithms Approach. Annual Conference of the PHM Society.
-        http://www.papers.phmsociety.org/index.php/phmconf/article/view/1404
-        
-        Baptistia, M., et. al. (2022). Relation between prognostics predictor evaluation metrics and local interpretability SHAP values. Aritifical Intelligence, Volume 306.
-        https://www.sciencedirect.com/science/article/pii/S0004370222000078
+        Where N is number of measurements and sign indicates sign of calculation [0]_ [1]_.
 
-        Args:
-            None
         Returns:
-            float: Value between [0, 1] indicating monotonicity of a given event for the Prediction.
+            dict (str, float): Value between [0, 1] indicating monotonicity of a given event for the Prediction.
+
+        References:
+            .. [0] Coble, J., et. al. (2021). Identifying Optimal Prognostic Parameters from Data: A Genetic Algorithms Approach. Annual Conference of the PHM Society. http://www.papers.phmsociety.org/index.php/phmconf/article/view/1404
+            .. [1] Baptistia, M., et. al. (2022). Relation between prognostics predictor evaluation metrics and local interpretability SHAP values. Aritifical Intelligence, Volume 306. https://www.sciencedirect.com/science/article/pii/S0004370222000078
+
         """
         # Collect and organize mean values for each event
         by_event = defaultdict(list)
@@ -108,9 +104,9 @@ class UnweightedSamplesPrediction(Prediction, UserList):
     Immutable data class for the result of a prediction, where the predictions are stored as UnweightedSamples. Is returned from the predict method of a sample based prediction class (e.g., MonteCarlo). Objects of this class can be iterated and accessed like a list (e.g., prediction[0]), where prediction[n] represents a profile for sample n.
 
     Args:
-        times : array(float)
-            Times for each data point where times[n] corresponds to data[n]
-        data : array(dict)
+        times (list[float]):
+            Times for each data point where times[n] corresponds to data[:][n]
+        data (list[SimResult]):
             Data points where data[n] is a SimResult for sample n
     """
 
diff --git a/src/prog_algs/predictors/predictor.py b/src/prog_algs/predictors/predictor.py
index 53219818..4696dc81 100644
--- a/src/prog_algs/predictors/predictor.py
+++ b/src/prog_algs/predictors/predictor.py
@@ -13,12 +13,12 @@ class Predictor(ABC):
     """
     Interface class for predictors
 
-    Abstract base class for creating predictors that perform prediction. Predictor subclasses must implement this interface. Equivilant to "Observers" in NASA's Matlab Prognostics Algorithm Library
+    Abstract base class for creating predictors that perform prediction. Predictor subclasses must implement this interface. Equivilant to "Predictors" in NASA's Matlab Prognostics Algorithm Library
 
     Parameters
     ----------
     model : PrognosticsModel
-        See: Prognostics Model Package\n
+        See: :py:mod:`prog_models` package\n
         A prognostics model to be used in prediction
     kwargs : optional, keyword arguments
     """
@@ -54,26 +54,28 @@ def predict(self, state : UncertainData, future_loading_eqn : Callable, **kwargs
             Distribution representing current state of the system
         future_loading_eqn : function (t, x) -> z
             Function to generate an estimate of loading at future time t, and state x
-        options : optional, keyword arguments
-            The following configuration parameters are supported: \n
-            * dt (float): Simulation step size (s), e.g., 0.1
-            * events (list[string]): Events to predict (subset of model.events) e.g., ['event1', 'event2']
-            * horizon (float): Prediction horizon (s)
-            * save_freq (float): Frequency at which results are saved (s)
-            * save_pts (list[float]): Any additional savepoints (s) e.g., [10.1, 22.5]
 
-        Return PredictionResults namedtuple
+        Keyword Arguments
+        -------------------
+        dt : float
+            Simulation step size (s), e.g., 0.1
+        events : list[str]
+            Events to predict (subset of model.events) e.g., ['event1', 'event2']
+        horizon : float
+            Prediction horizon (s)
+        save_freq : float
+            Frequency at which results are saved (s)
+        save_pts : list[float]
+            Any additional savepoints (s) e.g., [10.1, 22.5]
+
+        Return
         ----------
-        times : List[float]
-            Times for each savepoint such that inputs.snapshot(i), states.snapshot(i), outputs.snapshot(i), and event_states.snapshot(i) are all at times[i]            
-        inputs : Prediction
-            Inputs at each savepoint such that inputs.snapshot(i) is the input distribution (type UncertainData) at times[i]
-        states : Prediction
-            States at each savepoint such that states.snapshot(i) is the state distribution (type UncertainData) at times[i]
-        outputs : Prediction
-            Outputs at each savepoint such that outputs.snapshot(i) is the output distribution (type UncertainData) at times[i]
-        event_states : Prediction
-            Event states at each savepoint such that event_states.snapshot(i) is the event state distribution (type UncertainData) at times[i]
-        time_of_event : UncertainData
-            Distribution of predicted Time of Event (ToE) for each predicted event, represented by some subclass of UncertaintData (e.g., MultivariateNormalDist)
+        NamedTuple
+            * times (list[float]): Times for each savepoint such that inputs.snapshot(i), states.snapshot(i), outputs.snapshot(i), and event_states.snapshot(i) are all at times[i]            
+            * inputs (:py:class:`Prediction`): Inputs at each savepoint such that inputs.snapshot(i) is the input distribution (type UncertainData) at times[i]
+            * states (:py:class:`Prediction`): States at each savepoint such that states.snapshot(i) is the state distribution (type UncertainData) at times[i]
+            * outputs (:py:class:`Prediction`): Outputs at each savepoint such that outputs.snapshot(i) is the output distribution (type UncertainData) at times[i]
+            * event_states (:py:class:`Prediction`): Event states at each savepoint such that event_states.snapshot(i) is the event state distribution (type UncertainData) at times[i]
+            * time_of_event (:py:class:`prog_algs.uncertain_data.UncertainData`): Distribution of predicted Time of Event (ToE) for each predicted event, represented by some subclass of UncertaintData (e.g., :py:class:`MultivariateNormalDist`)
         """
+        pass
diff --git a/src/prog_algs/predictors/toe_prediction_profile.py b/src/prog_algs/predictors/toe_prediction_profile.py
index caae3727..cfbbb2ea 100644
--- a/src/prog_algs/predictors/toe_prediction_profile.py
+++ b/src/prog_algs/predictors/toe_prediction_profile.py
@@ -45,7 +45,7 @@ def alpha_lambda(self, ground_truth : Dict[str, float], lambda_value : float, al
         """Calculate Alpha lambda metric for the prediction profile
 
         Args:
-            ground_truth (Dict[str, float]):
+            ground_truth (dict[str, float]):
                 Ground Truth time of event for each event (e.g., {'event1': 748, 'event2', 2233, ...})
             lambda_value (float):
                 Prediction time at or after which metric is evaluated. Evaluation occurs at this time (if a prediction exists) or the next prediction following.
@@ -53,13 +53,15 @@ def alpha_lambda(self, ground_truth : Dict[str, float], lambda_value : float, al
                 percentage bounds around time to event (where 0.2 allows 20% error TtE)
             beta (float):
                 portion of prediction that must be within those bounds
-            kwargs (optional, keyword arguments):
-                configuration arguments. Accepted args include: \n
-                 * keys (list[string]): list of keys to use. If not provided, all keys are used.
-                 * print (bool) : If True, print the results. Default is False.
+
+        Keyword Args:
+            keys (list[str], optional): 
+                list of keys to use. If not provided, all keys are used.
+            print (bool, optional)
+                If True, print the results to the screen. Default is False.
 
         Returns:
-            Dict[str, bool]: If alpha lambda was met for each key (e.g., {'event1': True, 'event2', False, ...})
+            dict[str, bool]: If alpha lambda was met for each key (e.g., {'event1': True, 'event2', False, ...})
         """
         from ..metrics import alpha_lambda
         return alpha_lambda(self, ground_truth, lambda_value, alpha, beta, **kwargs)
@@ -68,7 +70,7 @@ def prognostic_horizon(self, criteria_eqn, ground_truth, **kwargs) -> Dict[str,
         """
         Compute prognostic horizon metric, defined as the difference between a time ti, when the predictions meet specified performance criteria, and the time corresponding to the true Time of Event (ToE), for each event.
 
-        PH = ToE - ti
+        :math:`PH = ToE - ti`
 
         Args:
             toe_profile (ToEPredictionProfile): A profile of predictions, the combination of multiple predictions
@@ -79,8 +81,10 @@ def prognostic_horizon(self, criteria_eqn, ground_truth, **kwargs) -> Dict[str,
                 | Returns: Map of event names to boolean representing if the event has been met. 
                 |   e.g., {'event1': True, 'event2': False}
             ground_truth (dict): Dictionary containing ground truth; specified as key, value pairs for event and its value. E.g, {'event1': 47.3, 'event2': 52.1, 'event3': 46.1}
-            kwargs (optional): configuration arguments. Accepted args include:
-                * print (bool): Boolean specifying whether the prognostic horizon metric should be printed.
+        
+        Keyword Args:
+            print (bool): 
+                Boolean specifying whether the prognostic horizon metric should be printed.
 
         Returns:
             dict: Dictionary containing prognostic horizon calculations (value) for each event (key). e.g., {'event1': 12.3, 'event2': 15.1}
@@ -88,43 +92,44 @@ def prognostic_horizon(self, criteria_eqn, ground_truth, **kwargs) -> Dict[str,
         from ..metrics import prognostic_horizon
         return prognostic_horizon(self, criteria_eqn, ground_truth, **kwargs)
 
-
     def cumulative_relative_accuracy(self, ground_truth, **kwargs) -> Dict[str, float]:
-        """
+        r"""
         Compute cumulative relative accuracy for a given profile, defined as the normalized sum of relative prediction accuracies at specific time instances.
         
-        CRA = Σ(RA)/N for each event
+        :math:`CRA = \Sigma \left( \dfrac{RA}{N} \right)` for each event
 
-        Where Σ is summation of all relative accuracies for a given profile and N is the total count of profiles (Journal Prognostics Health Management, Saxena et al.)
+        Where :math:`\Sigma` is summation of all relative accuracies for a given profile and N is the total count of profiles [0]_
 
         Args:
             ground_truth (dict): Dictionary containing ground truth; specified as key, value pairs for event and its value. E.g, {'event1': 47.3, 'event2': 52.1, 'event3': 46.1}
 
         Returns:
             dict: Dictionary containing cumulative relative accuracy (value) for each event (key). e.g., {'event1': 12.3, 'event2': 15.1}
+
+        References:
+            .. [0] Journal Prognostics Health Management, Saxena et al.
         """
         from ..metrics import cumulative_relative_accuracy
         return cumulative_relative_accuracy(self, ground_truth, **kwargs)
 
     def monotonicity(self, **kwargs) -> Dict[str, float]:
-        """Calculate monotonicty for a prediction profile. 
+        r"""Calculate monotonicty for a prediction profile. 
         Given a prediction profile, for each prediction: go through all predicted states and compare those to the next one.
-        Calculates monotonicity for each prediction key using its associated mean value in UncertainData.
+        Calculates monotonicity for each prediction key using its associated mean value in :py:class:`prog_algs.uncertain_data.UncertainData`.
         
-        monotonoicity = |Σsign(i+1 - i) / N-1|
+        :math:`monotonoicity = \|\Sigma \left( \dfrac{sign(i+1 - i)}{N-1} \right) \|`
 
-        Where N is number of measurements and sign indicates sign of calculation.
-
-        Coble, J., et. al. (2021). Identifying Optimal Prognostic Parameters from Data: A Genetic Algorithms Approach. Annual Conference of the PHM Society.
-        http://www.papers.phmsociety.org/index.php/phmconf/article/view/1404
-
-        Baptistia, M., et. al. (2022). Relation between prognostics predictor evaluation metrics and local interpretability SHAP values. Aritifical Intelligence, Volume 306.
-        https://www.sciencedirect.com/science/article/pii/S0004370222000078
+        Where N is number of measurements and sign indicates sign of calculation. [0]_ [1]_
 
         Args:
             toe_profile (ToEPredictionProfile): A profile of predictions, the combination of multiple predictions
+        
         Returns:
             dict (str, dict): Dictionary where keys represent a profile and dict is a subdictionary representing an event and its respective monotonicitiy value between [0, 1].
+
+        References:
+            .. [1] Coble, J., et. al. (2021). Identifying Optimal Prognostic Parameters from Data: A Genetic Algorithms Approach. Annual Conference of the PHM Society. http://www.papers.phmsociety.org/index.php/phmconf/article/view/1404
+            .. [2] Baptistia, M., et. al. (2022). Relation between prognostics predictor evaluation metrics and local interpretability SHAP values. Aritifical Intelligence, Volume 306. https://www.sciencedirect.com/science/article/pii/S0004370222000078
         """
         from ..metrics import monotonicity
         return monotonicity(self, **kwargs)
@@ -133,15 +138,15 @@ def plot(self, ground_truth : dict = None , alpha : float = None, show : bool =
         """Produce an alpha-beta plot depicting the TtE distribution by time of prediction for each event.
 
         Args:
-            ground_truth : dict
+            ground_truth (dict):
                 Optional dictionary argument containing event and its respective ground truth value; none by default and plotted if specified
-            alpha : float
+            alpha (float):
                 Optional alpha value; none by default and plotted if specified
-            show : bool
+            show (bool):
                 Optional bool value; specify whether to display generated plots. Default is true
         
         Returns:
-            dict[str, Figure]
+            dict[str, Figure] :
                 Collection of generated matplotlib figures for each event in profile\n
                 e.g., {'event1': Fig, 'event2': Fig}
         
diff --git a/src/prog_algs/predictors/unscented_transform.py b/src/prog_algs/predictors/unscented_transform.py
index ab73fedf..91116baf 100644
--- a/src/prog_algs/predictors/unscented_transform.py
+++ b/src/prog_algs/predictors/unscented_transform.py
@@ -1,11 +1,12 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration.  All Rights Reserved.
 
+from copy import deepcopy
+from filterpy import kalman
+from numpy import diag, array, transpose, isnan
 from typing import Callable
+
 from .prediction import Prediction, UnweightedSamplesPrediction, PredictionResults
 from .predictor import Predictor
-from numpy import diag, array, transpose, isnan
-from copy import deepcopy
-from filterpy import kalman
 from prog_algs.uncertain_data import MultivariateNormalDist, UncertainData, ScalarData
 
 
@@ -69,13 +70,13 @@ class UnscentedTransformPredictor(Predictor):
         Initial time at which prediction begins, e.g., 0
     dt : float
         Simulation step size (s), e.g., 0.1
-    events : List[string]
+    events : list[str]
         Events to predict (subset of model.events) e.g., ['event1', 'event2']
     horizon : float
         Prediction horizon (s)
     save_freq : float
         Frequency at which results are saved (s)
-    save_pts : List[float]
+    save_pts : list[float]
         Any additional savepoints (s) e.g., [10.1, 22.5]
 
     Note
@@ -105,15 +106,15 @@ def __init__(self, model, **kwargs):
 
         if 'Q' not in self.parameters:
             # Default 
-            self.parameters['Q'] = diag([1.0e-1 for i in range(num_states)])
+            self.parameters['Q'] = diag([1.0e-8 for _ in range(num_states)])
         
         def measure(x):
-            x = model.StateContainer({key: value for (key, value) in zip(self.__state_keys, x)})
+            x = model.StateContainer(x)
             z = model.output(x)
-            return model.OutputContainer({array(list(z.values()))})
+            return model.OutputContainer(z)
 
         def state_transition(x, dt):
-            x = model.StateContainer({key: value for (key, value) in zip(self.__state_keys, x)})
+            x = model.StateContainer(x)
             x = model.next_state(x, self.__input, dt)
             x = model.apply_limits(x)
             return array(list(x.values()))
@@ -155,9 +156,11 @@ def predict(self, state, future_loading_eqn : Callable, **kwargs) -> PredictionR
         event_states: [[dict]]
             Estimated event state (e.g., SOH), between 1-0 where 0 is event occurance, for each sample and time in times
             where event_states[sample_id][index] corresponds to time times[sample_id][index]
-        toe: UncertainData
+        time_of_event: UncertainData
             Estimated time where a predicted event will occur for each sample. Note: Mean and Covariance Matrix will both 
             be nan if every sigma point doesnt reach threshold within horizon
+            Also, includes member final_state (time_of_event.final_state) which is the state at the last time step. 
+            time_of_event.final_state is a dict of the form {'state_name': state_value}, is equal to None if event does not occur within horizon
         """
         if isinstance(state, dict) or isinstance(state, self.model.StateContainer) or isinstance(state, ScalarData):
             raise TypeError("state must be a distribution (e.g., MultivariateNormalDist, UnweightedSamples), not scalar")
@@ -168,9 +171,12 @@ def predict(self, state, future_loading_eqn : Callable, **kwargs) -> PredictionR
 
         params = deepcopy(self.parameters) # copy parameters
         params.update(kwargs) # update for specific run
-        events_to_predict = params['events']
+
+        if len(params['events']) == 0 and 'horizon' not in params:
+            raise ValueError("If specifying no event (i.e., simulate to time), must specify horizon")
 
         # Optimizations 
+        events_to_predict = params['events']
         dt = params['dt']
         model = self.model
         filt = self.filter
@@ -188,7 +194,7 @@ def predict(self, state, future_loading_eqn : Callable, **kwargs) -> PredictionR
         t = params['t0']
         save_pt_index = 0
         ToE = {key: [float('nan') for i in range(n_points)] for key in events_to_predict}  # Keep track of final ToE values
-        last_state = {key: [float('nan') for i in range(n_points)] for key in events_to_predict}  # Keep track of final state values
+        last_state = {key: [None for i in range(n_points)] for key in events_to_predict}  # Keep track of final state values
 
         times = []
         inputs = []
@@ -225,7 +231,8 @@ def update_all():
             points = sigma_points.sigma_points(filt.x, filt.P)
             all_failed = True
             for i, point in zip(range(n_points), points):
-                x = StateContainer({key: x for (key, x) in zip(state_keys, point)})
+                # x = StateContainer({key: x for (key, x) in zip(state_keys, point)})
+                x = StateContainer(point)
                 t_met = threshold_met(x)
 
                 # Check Thresholds
@@ -249,6 +256,10 @@ def update_all():
         # Transform final state into {event_name: MultivariateNormalDist}
         final_state = {}
         for event_key in last_state.keys():
+            if any([last_state_i is None for last_state_i in last_state[event_key]]):
+                # If any sigma point has not met the event threshold
+                final_state[event_key] = None
+                continue
             last_state_pts = array([[last_state_i[state_key] for state_key in state_keys] for last_state_i in last_state[event_key]])
             # last_state_pts = transpose(last_state_pts)
             last_state_mean, last_state_cov = kalman.unscented_transform(last_state_pts, sigma_points.Wm, sigma_points.Wc)
diff --git a/src/prog_algs/state_estimators/kalman_filter.py b/src/prog_algs/state_estimators/kalman_filter.py
index 9b95661a..5a30889d 100644
--- a/src/prog_algs/state_estimators/kalman_filter.py
+++ b/src/prog_algs/state_estimators/kalman_filter.py
@@ -1,11 +1,12 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration.  All Rights Reserved.
 
 from copy import deepcopy
-from typing import Callable
+from filterpy import kalman
 import numpy as np
 from warnings import warn
-from filterpy import kalman
+
 from prog_models import LinearModel
+
 from . import state_estimator
 from ..uncertain_data import MultivariateNormalDist, UncertainData
 
@@ -18,27 +19,24 @@ class KalmanFilter(state_estimator.StateEstimator):
     The supported configuration parameters (keyword arguments) for UKF construction are described below:
 
     Args:
-        model : ProgModel
+        model (PrognosticsModel):
             A prognostics model to be used in state estimation
             See: Prognostics Model Package
-        x0 : UncertainData, model.StateContainer, or dict
+        x0 (UncertainData, model.StateContainer, or dict):
             Initial (starting) state, with keys defined by model.states \n
             e.g., x = ScalarData({'abc': 332.1, 'def': 221.003}) given states = ['abc', 'def']
 
     Keyword Args:
-        alpha: float
+        alpha (float, optional):
             KF Scaling parameter. An alpha > 1 turns this into a fading memory filter.
-        t0 : float
+        t0 (float, optional):
             Starting time (s)
-        dt : float 
+        dt (float, optional):
             time step (s)
-        Q : List[List[float]]
-            Process Noise Matrix 
-        R : List[List[float]]
-            Measurement Noise Matrix
-
-    Note:
-        The Kalman Filter does not support a custom measurement function
+        Q (list[list[float]], optional):
+            Kalman Process Noise Matrix 
+        R (list[list[float]], optional):
+            Kalman Measurement Noise Matrix
     """
     default_parameters = {
         'alpha': 1, 
@@ -46,12 +44,12 @@ class KalmanFilter(state_estimator.StateEstimator):
         'dt': 1
     } 
     
-    def __init__(self, model, x0, measurement_eqn : Callable = None, **kwargs):
+    def __init__(self, model, x0, **kwargs):
         # Note: Measurement equation kept in constructor to keep it consistent with other state estimators. This way measurement equation can be provided as an ordered argument, and will just be ignored here
         if not isinstance(model, LinearModel):
             raise Exception('Kalman Filter only supports Linear Models (i.e., models derived from prog_models.LinearModel)')
 
-        super().__init__(model, x0, None, **kwargs)
+        super().__init__(model, x0, **kwargs)
 
         self.x0 = x0
 
diff --git a/src/prog_algs/state_estimators/particle_filter.py b/src/prog_algs/state_estimators/particle_filter.py
index b218c76d..93dd31d9 100644
--- a/src/prog_algs/state_estimators/particle_filter.py
+++ b/src/prog_algs/state_estimators/particle_filter.py
@@ -1,16 +1,16 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration.  All Rights Reserved.
 
-from prog_algs.uncertain_data.uncertain_data import UncertainData
-from typing import Callable
-from . import state_estimator
-from numpy import array, empty, random, take, exp, max, take
 from filterpy.monte_carlo import residual_resample
-from numbers import Number
+from numpy import array, empty, take, exp, max, take
 from scipy.stats import norm
-from ..uncertain_data import UnweightedSamples
-from ..exceptions import ProgAlgTypeError
 from warnings import warn
 
+from prog_models.utils.containers import DictLikeMatrixWrapper
+
+from . import state_estimator
+from ..uncertain_data import UnweightedSamples, ScalarData, UncertainData
+from ..exceptions import ProgAlgTypeError
+
 
 class ParticleFilter(state_estimator.StateEstimator):
     """
@@ -21,69 +21,44 @@ class ParticleFilter(state_estimator.StateEstimator):
     The supported configuration parameters (keyword arguments) for UKF construction are described below:
 
     Args:
-        model : ProgModel
+        model (PrognosticsModel):
             A prognostics model to be used in state estimation
             See: Prognostics Model Package
-        x0 : UncertainData, model.StateContainer, or dict
+        x0 (UncertainData, model.StateContainer, or dict):
             Initial (starting) state, with keys defined by model.states \n
             e.g., x = ScalarData({'abc': 332.1, 'def': 221.003}) given states = ['abc', 'def']
 
     Keyword Args:
-        t0 : float
+        t0 (float, optional):
             Starting time (s)
-        dt : float 
+        dt (float, optional): 
             time step (s)
-        num_particles : int
+        num_particles (int, optional):
             Number of particles in particle filter
-        resample_fcn : function 
+        resample_fcn (function, optional):
             Resampling function ([weights]) -> [indexes] e.g., filterpy.monte_carlo.residual_resample
     """
     default_parameters = {
             't0': -1e-99,  # practically 0, but allowing for a 0 first estimate
             'num_particles': 20, 
             'resample_fcn': residual_resample,
-            'x0_uncertainty': 0.5
         }
 
-    def __init__(self, model, x0, measurement_eqn : Callable = None, **kwargs):
-        super().__init__(model, x0, measurement_eqn = measurement_eqn, **kwargs)
+    def __init__(self, model, x0, **kwargs):
+        super().__init__(model, x0, **kwargs)
         
-        if measurement_eqn:
-            warn("Warning: measurement_eqn depreciated as of v1.3.1, will be removed in v1.4. Use Model subclassing instead. See examples.measurement_eqn_example")
-            # update output_container
-            from prog_models.utils.containers import DictLikeMatrixWrapper
-            z0 = measurement_eqn(x0)
-            class MeasureContainer(DictLikeMatrixWrapper):
-                def __init__(self, z):
-                    super().__init__(list(z0.keys()), z)
-
-            def __measure(x):
-                return MeasureContainer(measurement_eqn(x))
-                
-            self._measure = __measure
-        else:
-            self._measure = model.output
+        self._measure = model.output
 
-        # Caching for optimization
-        paramters_x0_exist = 'x0_uncertainty' in self.parameters
-        if paramters_x0_exist: # Only create these optimizations if x0_uncertainty exists as key in self.parameters
-            parameters_x0_dict, parameters_x0_num = isinstance(self.parameters['x0_uncertainty'], dict), isinstance(self.parameters['x0_uncertainty'], Number)
         # Build array inplace
-        if isinstance(x0, UncertainData):
-            sample_gen = x0.sample(self.parameters['num_particles'])
-            samples = [array(sample_gen.key(k)) for k in x0.keys()]
-        elif paramters_x0_exist and (parameters_x0_dict or parameters_x0_num):
-            warn("Warning: x0_uncertainty depreciated as of v1.3, will be removed in v1.4. Use UncertainData type if estimating filtering with uncertain data.")
-            x = array(list(x0.values()))
-            if parameters_x0_dict:
-                sd = array([self.parameters['x0_uncertainty'][key] for key in x0.keys()])
-            elif parameters_x0_num:
-                sd = array([self.parameters['x0_uncertainty']] * len(x0))
-            samples = [random.normal(x[i], sd[i], self.parameters['num_particles']) for i in range(len(x))]
-        else:
-            raise ProgAlgTypeError("ProgAlgTypeError: x0 must be of type UncertainData or x0_uncertainty must be of type [dict, Number].")
-        self.particles = dict(zip(x0.keys(), samples))
+        if isinstance(x0, DictLikeMatrixWrapper) or isinstance(x0, dict):
+            x0 = ScalarData(x0)
+        elif not isinstance(x0, UncertainData):
+            raise ProgAlgTypeError(f"ProgAlgTypeError: x0 must be of type UncertainData or StateContainer, was {type(x0)}.")
 
+        sample_gen = x0.sample(self.parameters['num_particles'])
+        samples = [array(sample_gen.key(k)) for k in x0.keys()]
+        
+        self.particles = model.StateContainer(array(samples))
 
         if 'R' in self.parameters:
             # For backwards compatibility
@@ -100,6 +75,12 @@ def estimate(self, t : float, u, z):
         dt = t - self.t
         self.t = t
 
+        # Check Types
+        if isinstance(u, dict):
+            u = self.model.InputContainer(u)
+        if isinstance(z, dict):
+            z = self.model.OutputContainer(z)
+
         # Optimization
         particles = self.particles
         next_state = self.model.next_state
@@ -162,7 +143,7 @@ def estimate(self, t : float, u, z):
                    for state in self.particles.keys()]
 
         # Particles as a dictionary
-        self.particles = dict(zip(self.particles.keys(), samples))
+        self.particles = self.model.StateContainer(array(samples))
 
     @property
     def x(self) -> UnweightedSamples:
diff --git a/src/prog_algs/state_estimators/state_estimator.py b/src/prog_algs/state_estimators/state_estimator.py
index cb2f8d9c..fe65bc54 100644
--- a/src/prog_algs/state_estimators/state_estimator.py
+++ b/src/prog_algs/state_estimators/state_estimator.py
@@ -1,12 +1,11 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration.  All Rights Reserved.
 
 from abc import ABC, abstractmethod, abstractproperty
-from typing import Callable
+from copy import deepcopy
 
 from ..uncertain_data import UncertainData
 from ..exceptions import ProgAlgTypeError
-from copy import deepcopy
-from warnings import warn
+
 
 class StateEstimator(ABC):
     """
@@ -15,22 +14,21 @@ class StateEstimator(ABC):
     Abstract base class for creating state estimators that perform state estimation. Subclasses must implement this interface. Equivilant to "Observers" in NASA's Matlab Prognostics Algorithm Library
 
     Args:
-        model : ProgModel
+        model (PrognosticsModel):
             A prognostics model to be used in state estimation
             See: Prognostics Model Package
-        x0 : UncertainData, model.StateContainer, or dict
+        x0 (UncertainData, model.StateContainer, or dict):
             Initial (starting) state, with keys defined by model.states \n
             e.g., x = ScalarData({'abc': 332.1, 'def': 221.003}) given states = ['abc', 'def']
-        options : optional, kwargs
-            configuration options\n
-            Dictionary of any additional configuration values. See state-estimator specific documentation
+    
+    See state-estimator specific documentation for speicfic keyword arguments.
     """
 
     default_parameters = {
         't0': 0
     }
 
-    def __init__(self, model, x0, measurement_eqn : Callable = None, **kwargs):
+    def __init__(self, model, x0, **kwargs):
         # Check model
         if not hasattr(model, 'output'):
             raise ProgAlgTypeError("model must have `output` method")
@@ -42,18 +40,11 @@ def __init__(self, model, x0, measurement_eqn : Callable = None, **kwargs):
             raise ProgAlgTypeError("model must have `states` property")
         self.model = model
 
-        if measurement_eqn is not None:
-            warn('Measurement_eqn was deprecated in v1.3 in favor of model subclassing. I will remove this in v1.4. See measurement_equation example for more information', DeprecationWarning)
-
         # Check x0
         for key in model.states:
             if key not in x0:
                 raise ProgAlgTypeError("x0 missing state `{}`".format(key))
         
-        # Check measurement equation
-        if measurement_eqn and not callable(measurement_eqn):
-            raise ProgAlgTypeError("measurement_eqn must be callable")
-        
         # Process kwargs (configuration)
         self.parameters = deepcopy(self.default_parameters)
         self.parameters.update(kwargs)
@@ -70,12 +61,12 @@ def estimate(self, t : float, u, z, **kwargs) -> None:
         t : float
             Current timestamp in seconds (≥ 0.0)
             e.g., t = 3.4
-        u : dict
+        u : InputContainer
             Measured inputs, with keys defined by model.inputs.
-            e.g., u = {'i':3.2} given inputs = ['i']
-        z : dict
+            e.g., u = m.InputContainer({'i':3.2}) given inputs = ['i']
+        z : OutputContainer
             Measured outputs, with keys defined by model.outputs.
-            e.g., z = {'t':12.4, 'v':3.3} given inputs = ['t', 'v']
+            e.g., z = m.OutputContainer({'t':12.4, 'v':3.3}) given outputs = ['t', 'v']
 
         Note
         ----
diff --git a/src/prog_algs/state_estimators/unscented_kalman_filter.py b/src/prog_algs/state_estimators/unscented_kalman_filter.py
index 442e9740..6118782b 100644
--- a/src/prog_algs/state_estimators/unscented_kalman_filter.py
+++ b/src/prog_algs/state_estimators/unscented_kalman_filter.py
@@ -1,10 +1,10 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration.  All Rights Reserved.
 
-from warnings import warn
-from typing import Callable
-from . import state_estimator
 from filterpy import kalman
 from numpy import diag, array
+from warnings import warn
+
+from . import state_estimator
 from ..uncertain_data import MultivariateNormalDist, UncertainData
 
 class UnscentedKalmanFilter(state_estimator.StateEstimator):
@@ -16,23 +16,27 @@ class UnscentedKalmanFilter(state_estimator.StateEstimator):
     The supported configuration parameters (keyword arguments) for UKF construction are described below:
 
     Args:
-        model : ProgModel
+        model (PrognosticsModel):
             A prognostics model to be used in state estimation
             See: Prognostics Model Package
-        x0 : UncertainData, model.StateContainer, or dict
+        x0 (UncertainData, model.StateContainer, or dict):
             Initial (starting) state, with keys defined by model.states \n
             e.g., x = ScalarData({'abc': 332.1, 'def': 221.003}) given states = ['abc', 'def']
 
-    Keyword Args
-        alpha, beta, kappa: float
-            UKF Scaling parameters
-        t0 : float
+    Keyword Args:
+        alpha (float, optional):
+            UKF Scaling parameter
+        beta (float, optional):
+            UKF Scaling parameter
+        kappa (float, optional):
+            UKF Scaling parameter
+        t0 (float, optional):
             Starting time (s)
-        dt : float 
+        dt (float, optional):
             time step (s)
-        Q : List[List[float]]
+        Q (list[list[float]], optional):
             Process Noise Matrix 
-        R : List[List[float]]
+        R (list[list[float]], optional):
             Measurement Noise Matrix 
     """
     default_parameters = {
@@ -43,27 +47,20 @@ class UnscentedKalmanFilter(state_estimator.StateEstimator):
         'dt': 1
     } 
 
-    def __init__(self, model, x0, measurement_eqn : Callable = None, **kwargs):
-        super().__init__(model, x0, measurement_eqn, **kwargs)
+    def __init__(self, model, x0, **kwargs):
+        super().__init__(model, x0, **kwargs)
 
         self.__input = None
         self.x0 = x0
         # Saving for reduce pickling
 
-        if measurement_eqn is None: 
-            def measure(x):
-                x = model.StateContainer({key: value for (key, value) in zip(x0.keys(), x)})
-                R_err = model.parameters['measurement_noise'].copy()
-                model.parameters['measurement_noise'] = dict.fromkeys(R_err, 0)
-                z = model.output(x)
-                model.parameters['measurement_noise'] = R_err
-                return array(list(z.values())).ravel()
-        else:
-            warn("Warning: measurement_eqn depreciated as of v1.3.1, will be removed in v1.4. Use Model subclassing instead. See examples.measurement_eqn_example")
-            def measure(x):
-                x = model.StateContainer({key: value for (key, value) in zip(x0.keys(), x)})
-                z = measurement_eqn(x)
-                return array(list(z.values())).ravel()
+        def measure(x):
+            x = model.StateContainer({key: value for (key, value) in zip(x0.keys(), x)})
+            R_err = model.parameters['measurement_noise'].copy()
+            model.parameters['measurement_noise'] = dict.fromkeys(R_err, 0)
+            z = model.output(x)
+            model.parameters['measurement_noise'] = R_err
+            return array(list(z.values())).ravel()
 
         if 'Q' not in self.parameters:
             self.parameters['Q'] = diag([1.0e-3 for i in x0.keys()])
@@ -81,7 +78,7 @@ def state_transition(x, dt):
         self.filter = kalman.UnscentedKalmanFilter(num_states, num_measurements, self.parameters['dt'], measure, state_transition, points)
         
         if isinstance(x0, dict) or isinstance(x0, model.StateContainer):
-            warn("Warning: x0_uncertainty depreciated as of v1.3, will be removed in v1.4. Use UncertainData type if estimating filtering with uncertain data.")
+            warn("Use UncertainData type if estimating filtering with uncertain data.")
             self.filter.x = array(list(x0.values()))
             self.filter.P = self.parameters['Q'] / 10
         elif isinstance(x0, UncertainData):
diff --git a/src/prog_algs/uncertain_data/multivariate_normal_dist.py b/src/prog_algs/uncertain_data/multivariate_normal_dist.py
index 1010ceb4..2f8b40c7 100644
--- a/src/prog_algs/uncertain_data/multivariate_normal_dist.py
+++ b/src/prog_algs/uncertain_data/multivariate_normal_dist.py
@@ -1,9 +1,10 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration. All Rights Reserved.
 
-from . import UncertainData, UnweightedSamples
 from numpy import array
 from numpy.random import multivariate_normal
 
+from . import UncertainData, UnweightedSamples
+
 
 class MultivariateNormalDist(UncertainData):
     """
diff --git a/src/prog_algs/uncertain_data/scalar_data.py b/src/prog_algs/uncertain_data/scalar_data.py
index c680fee3..28fe8926 100644
--- a/src/prog_algs/uncertain_data/scalar_data.py
+++ b/src/prog_algs/uncertain_data/scalar_data.py
@@ -2,6 +2,7 @@
 
 from typing import Union
 from numpy import array
+
 from . import UncertainData, UnweightedSamples
 
 
@@ -10,7 +11,7 @@ class ScalarData(UncertainData):
     Data without uncertainty- single value
 
     Args:
-            state (dict or Container): Single state in the form of dict or model.*Container (InputContainer, OutputContainer, Statecontainer) representing states and respective values.
+        state (dict or Container): Single state in the form of dict or model.*Container (InputContainer, OutputContainer, Statecontainer) representing states and respective values.
     """
     def __init__(self, state, _type = dict): 
         self.__state = state
@@ -55,7 +56,7 @@ def __isub__(self, other : int) -> "UncertainData":
         return self
 
     @property
-    def median(self) -> array:
+    def median(self) -> dict:
         return self.mean
         
     @property
diff --git a/src/prog_algs/uncertain_data/uncertain_data.py b/src/prog_algs/uncertain_data/uncertain_data.py
index 18332e08..7641511a 100644
--- a/src/prog_algs/uncertain_data/uncertain_data.py
+++ b/src/prog_algs/uncertain_data/uncertain_data.py
@@ -11,12 +11,13 @@
 
 
 class UncertainData(ABC):
-    def __init__(self, _type = dict):
-        self._type = _type
-
     """
     Abstract base class for data with uncertainty. Any new uncertainty type must implement this class
     """
+
+    def __init__(self, _type = dict):
+        self._type = _type
+    
     @abstractmethod
     def sample(self, nSamples : int = 1):
         """Generate samples from data
@@ -28,31 +29,37 @@ def sample(self, nSamples : int = 1):
             samples (UnweightedSamples): Array of nSamples samples
 
         Example:
-            samples = data.samples(100)
+            ::
+            
+                samples = data.samples(100)
         """
 
     @property
     @abstractproperty
-    def median(self) -> array:
+    def median(self) -> dict:
         """The median of the UncertainData distribution or samples 
 
         Returns:
-            Dict[str, float]: Median value. e.g., {'key1': 23.2, ...}
+            dict[str, float]: Median value. e.g., {'key1': 23.2, ...}
 
         Example:
-            median_value = data.median
+            ::
+            
+                median_value = data.median
         """
 
     @property
     @abstractproperty
-    def mean(self) -> array:
+    def mean(self) -> dict:
         """The mean of the UncertainData distribution or samples 
 
         Returns:
-            Dict[str, float]: Mean value. e.g., {'key1': 23.2, ...}
+            dict[str, float]: Mean value. e.g., {'key1': 23.2, ...}
 
         Example:
-            mean_value = data.mean
+            ::
+                
+                mean_value = data.mean
         """
 
     @property
@@ -61,20 +68,31 @@ def cov(self) -> array:
         """The covariance matrix of the UncertiantyData distribution or samples in order of keys (i.e., cov[1][1] is the standard deviation for key keys()[1])
 
         Returns:
-            array[array[float]]: Covariance matrix
+            np.array[np.array[float]]: Covariance matrix
+
+        Example:
+            ::
+
+                covariance_matrix = data.cov
         """
 
     def relative_accuracy(self, ground_truth : dict) -> dict:
-        """The relative accuracy of the mean of the distribution. 
+        """The relative accuracy is how close the mean of the distribution is to the ground truth, on relative terms
         
-        RA = 1 - |r-p|/r 
+        :math:`RA = 1 - \dfrac{\| r-p \|}{r}`
 
-        Where r is ground truth and p is mean of predicted distribution
-        
-        Prognostics: The Science of Making Predictions (Goebel et al, 239)
+        Where r is ground truth and p is mean of predicted distribution [0]_
 
         Returns:
             dict[str, float]: Relative accuracy for each event where value is relative accuracy between [0,1]
+
+        Example:
+            ::
+
+                ra = data.relative_accuracy({'key1': 22, 'key2': 57})
+
+        References:
+            .. [0] Prognostics: The Science of Making Predictions (Goebel et al, 239)
         """
         # if this check isn't here, goes to divide by zero check and raises AttributeError instead of TypeError. Keep? There are unittests checking for type
         if not (isinstance(ground_truth, dict) or isinstance(ground_truth, DictLikeMatrixWrapper)):
@@ -88,7 +106,12 @@ def keys(self):
         """Get the keys for the property represented
 
         Returns:
-            [string]: keys
+            list[str]: keys
+
+        Example:
+            ::
+                
+                keys = data.keys()
         """
 
     def __contains__(self, key):
@@ -101,12 +124,18 @@ def percentage_in_bounds(self, bounds : tuple, keys : list = None, n_samples = 1
             bounds (tuple[float, float] or dict): Lower and upper bounds. \n
                 if tuple: (lower, upper)\n
                 if dict: {key: (lower, upper), ...}
-            keys (list of strings, optional): UncertainData keys to consider when calculating. Defaults to all keys.
+            keys (list[str], optional): UncertainData keys to consider when calculating. Defaults to all keys.
             n_samples (int, optional): Number of samples to use when calculating
 
         Returns:
-            dict: Percentage within bounds for each key in keys (where 0.5 = 50%).\n
-                e.g., {'key1': 1, 'key2': 0.75}
+            dict: Percentage within bounds for each key in keys (where 0.5 = 50%). e.g., {'key1': 1, 'key2': 0.75}
+
+        Example:
+            ::
+            
+                data.percentage_in_bounds((1025, 1075))
+                data.percentage_in_bounds({'key1': (1025, 1075), 'key2': (2520, 2675)})
+                data.percentage_in_bounds((1025, 1075), keys=['key1', 'key3'])
         """
         return self.sample(n_samples).percentage_in_bounds(bounds, keys=keys)
 
@@ -116,10 +145,17 @@ def metrics(self, **kwargs) -> dict:
         Keyword Args:
             ground_truth (int or dict, optional): Ground truth value. Defaults to None.
             n_samples (int, optional): Number of samples to use for calculating metrics (if not UnweightedSamples)
-            keys (List[str], optional): Keys to calculate metrics for. Defaults to all keys.
+            keys (list[str], optional): Keys to calculate metrics for. Defaults to all keys.
 
         Returns:
             dict: Dictionary of metrics
+
+        Example:
+            ::
+
+                print(data.metrics())
+                m = data.metrics(ground_truth={'key1': 200, 'key2': 350})
+                m = data.metrics(keys=['key1', 'key3'])
         """
         from ..metrics import calc_metrics
         return calc_metrics(self, **kwargs)
@@ -130,7 +166,7 @@ def plot_scatter(self, fig : Figure = None, keys : list = None, num_samples : in
 
         Args:
             fig (Figure, optional): Existing figure previously used to plot states. If passed a figure argument additional data will be added to the plot. Defaults to creating new figure
-            keys (list of strings, optional): Keys to plot. Defaults to all keys.
+            keys (list[str], optional): Keys to plot. Defaults to all keys.
             num_samples (int, optional): Number of samples to plot. Defaults to 100
             **kwargs (optional): Additional keyword arguments passed to scatter function.
 
@@ -159,7 +195,7 @@ def plot_hist(self, fig = None, keys = None, num_samples = 100, **kwargs):
         Args:
             fig (MatPlotLib Figure, optional): Existing histogram figure to be overritten. Defaults to create new figure.
             num_samples (int, optional): Number of samples to plot. Defaults to 100
-            keys (List(String), optional): Keys to be plotted. Defaults to None.
+            keys (list(String), optional): Keys to be plotted. Defaults to None.
 
         Example:
             ::
@@ -189,6 +225,11 @@ def describe(self, title : str = "UncertainData Metrics", print : bool = True) -
         Returns:
             defaultdict
                 Dictionary of lists used to print metrics.
+        
+        Example:
+            ::
+
+                data.describe()
         """
         recursive_metrics_table = print_table_recursive(self.metrics(), title, print)
         return recursive_metrics_table
diff --git a/src/prog_algs/uncertain_data/unweighted_samples.py b/src/prog_algs/uncertain_data/unweighted_samples.py
index 032deb4b..da9130f8 100644
--- a/src/prog_algs/uncertain_data/unweighted_samples.py
+++ b/src/prog_algs/uncertain_data/unweighted_samples.py
@@ -1,11 +1,14 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration. All Rights Reserved.
 
-from . import UncertainData
 from collections import UserList
 from collections.abc import Iterable
 from numpy import array, cov, random
 from warnings import warn
 
+from prog_models.utils.containers import DictLikeMatrixWrapper
+
+from . import UncertainData
+
 
 class UnweightedSamples(UncertainData, UserList):
     """
@@ -19,7 +22,7 @@ class UnweightedSamples(UncertainData, UserList):
     """
     def __init__(self, samples : list = [], _type = dict):
         super().__init__(_type)
-        if isinstance(samples, dict):
+        if isinstance(samples, dict) or isinstance(samples, DictLikeMatrixWrapper):
             # Is in form of {key: [value, ...], ...}
             # Convert to array of samples
             if len(samples.keys()) == 0:
@@ -163,6 +166,8 @@ def size(self) -> int:
     def percentage_in_bounds(self, bounds, keys : list = None) -> dict:
         if not keys:
             keys = self.keys()
+        if isinstance(keys, str):
+            keys = [keys]
         if isinstance(bounds, list):
             bounds = {key: bounds for key in self.keys()}
         if not isinstance(bounds, dict) or all([isinstance(b, list) and len(b) == 2 for b in bounds]):
diff --git a/src/prog_algs/visualize/plot_hist.py b/src/prog_algs/visualize/plot_hist.py
index 428e3fb2..267a5511 100644
--- a/src/prog_algs/visualize/plot_hist.py
+++ b/src/prog_algs/visualize/plot_hist.py
@@ -3,7 +3,6 @@
 from matplotlib.figure import Figure
 import matplotlib.pyplot as plt
 from numpy import array
-
 from typing import List
 
 def plot_hist(samples : array, fig : Figure = None, keys : List[str] = None, **kwargs) -> Figure:
diff --git a/src/prog_algs/visualize/plot_scatter.py b/src/prog_algs/visualize/plot_scatter.py
index 38c6a2b6..cab22872 100644
--- a/src/prog_algs/visualize/plot_scatter.py
+++ b/src/prog_algs/visualize/plot_scatter.py
@@ -3,8 +3,6 @@
 from matplotlib.figure import Figure
 import matplotlib.pyplot as plt
 from math import sqrt
-
-from numpy import array
 from typing import List
 
 def plot_scatter(samples : List[dict], fig : Figure = None, keys : List[str] = None, legend : str = 'auto', **kwargs) -> Figure:
diff --git a/state_estimator_template.py b/state_estimator_template.py
index 9a6ef7d1..9aec63c6 100644
--- a/state_estimator_template.py
+++ b/state_estimator_template.py
@@ -15,11 +15,11 @@ class TemplateStateEstimator(state_estimators.StateEstimator):
         't0' : 0.0  # Initial timestamp
     } 
 
-    def __init__(self, model, x0, measurement_eqn=None, **kwargs):
+    def __init__(self, model, x0, **kwargs):
         """
         Constructor (optional)
         """
-        super().__init__(model, x0, measurement_eqn, **kwargs)
+        super().__init__(model, x0, **kwargs)
         # ADD PARAMETER CHECKS HERE
         # e.g., self.parameters['some_value'] < 0
 
diff --git a/tutorial.ipynb b/tutorial.ipynb
index 1eab11d1..cb91125b 100644
--- a/tutorial.ipynb
+++ b/tutorial.ipynb
@@ -24,7 +24,7 @@
     "* `prog_algs.predictors` - Tools for performing prediction\n",
     "* `prog_algs.uncertain_data` - Tools for representing data with uncertainty\n",
     "\n",
-    "In addition to the `prog_algs` package, this repo includes examples showing how to use the package (see `examples/`), a template for implementing a new state estimator (`state_estimator_template`), a template for implementing a new predictor (`predictor_template`), documentation (<https://nasa.github.io/prog_algs>), and this tutorial (`tutorial.ipynb`).\n",
+    "In addition to the `prog_algs` package, this repo includes examples showing how to use the package (see `examples/`), a template for implementing a new state estimator (`state_estimator_template`), a template for implementing a new predictor (`predictor_template`), documentation (<https://nasa.github.io/progpy/guide>), and this tutorial (`tutorial.ipynb`).\n",
     "\n",
     "Before you start, install `prog_algs` using pip:\n",
     "\n",
@@ -192,7 +192,7 @@
     "\n",
     "In `prog_algs` a State Estimator is used to estimate the system state. \n",
     "\n",
-    "To start, import the needed packages. Here we will import the `BatteryCircuit` model and the `UnscentedKalmanFilter` state estimator. See <https://nasa.github.io/prog_algs> for more details on the available state estimators.\n"
+    "To start, import the needed packages. Here we will import the `BatteryCircuit` model and the `UnscentedKalmanFilter` state estimator. See <https://nasa.github.io/progpy/guide> for more details on the available state estimators.\n"
    ]
   },
   {
@@ -253,8 +253,8 @@
     "fig = est.x.plot_scatter(label='prior')\n",
     "\n",
     "t = 0.1\n",
-    "u = {'i': 2}\n",
-    "example_measurements = {'t': 32.2, 'v': 3.915}\n",
+    "u = m.InputContainer({'i': 2})\n",
+    "example_measurements = m.OutputContainer({'t': 32.2, 'v': 3.915})\n",
     "est.estimate(t, u, example_measurements)\n",
     "\n",
     "print(\"Posterior State:\", est.x.mean)\n",
@@ -539,8 +539,8 @@
     "se = BlindlyStumbleEstimator(m, x0)\n",
     "\n",
     "for i in range(25):\n",
-    "    u = {'i': 0}\n",
-    "    z = {'t': 18.95, 'v': 4.183}\n",
+    "    u = m.InputContainer({'i': 0})\n",
+    "    z = m.OutputContainer({'t': 18.95, 'v': 4.183})\n",
     "    se.estimate(i, u, z)\n",
     "    print(se.x.mean)\n",
     "    print(\"\\tcorrect: {'tb': 18.95, 'qb': 7856.3254, 'qcp': 0, 'qcs': 0}\")"
@@ -565,16 +565,14 @@
    "metadata": {},
    "source": [
     "## Conclusion\n",
-    "This is just the basics, there's much more to learn. Please see the documentation at <https://nasa.github.io/prog_algs> and the examples in the `examples/` folder for more details on how to use the package, including:\n",
+    "This is just the basics, there's much more to learn. Please see the documentation at <https://nasa.github.io/progpy/guide> and the examples in the `examples/` folder for more details on how to use the package, including:\n",
     "* `examples.basic_example` : A basic Example using prog_algs for Prognostics \n",
-    "* `examples.benchmarking_example` : An example benchmarking the performance of prognostics algorithms\n",
+    "* `examples.basic_example_battery` : A basic Example using prog_algs for Prognostics, using the more complex battery model\n",
     "* `examples.eol_event` : An example where a model has multiple events, but the user is only interested in predicting the time when the first event occurs (whatever it is).\n",
     "* `examples.measurement_eqn_example` : An example where not every output is measured or measurements are not in the same format as outputs, so a measurement equation is defined to translate between outputs and what's measured. \n",
     "* `examples.new_state_estimator_example` : An example of extending StateEstimator to create a new state estimator class\n",
     "* `examples.playback` : A full example performing prognostics using playback data.\n",
     "* `examples.predict_specific_event` : An example where the model has multiple events, but the user is only interested in predicting a specific event (or events).\n",
-    "* `examples.thrown_object_example` : An example performing prognostics with the simplified ThrownObject model\n",
-    "* `examples.utpredictor` : An example using the Unscented Transform Predictor for prediction.\n",
     "\n",
     "Thank you for trying out this tutorial. Open an issue on github (<https://github.com/nasa/prog_algs/issues>) or email Chris Teubert (christopher.a.teubert@nasa.gov) with any questions or issues."
    ]
@@ -588,11 +586,8 @@
   }
  ],
  "metadata": {
-  "interpreter": {
-   "hash": "ff94885aa2d97705a9dae03869c2058fa855d1acd9df351499300343e2e591a2"
-  },
   "kernelspec": {
-   "display_name": "Python 3.8.13 64-bit",
+   "display_name": "Python 3.8.9 64-bit",
    "language": "python",
    "name": "python3"
   },
@@ -606,14 +601,19 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.13"
+   "version": "3.8.9"
   },
   "metadata": {
    "interpreter": {
     "hash": "ff94885aa2d97705a9dae03869c2058fa855d1acd9df351499300343e2e591a2"
    }
   },
-  "orig_nbformat": 2
+  "orig_nbformat": 2,
+  "vscode": {
+   "interpreter": {
+    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
+   }
+  }
  },
  "nbformat": 4,
  "nbformat_minor": 2


Test patch:
diff --git a/tests/benchmark.py b/tests/benchmark.py
new file mode 100644
index 00000000..a9778489
--- /dev/null
+++ b/tests/benchmark.py
@@ -0,0 +1,126 @@
+# Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration.  All Rights Reserved.
+
+import time
+from io import StringIO
+import sys
+
+FORMAT_STR = '{:40s}'
+CLOCK = time.process_time
+
+if __name__ == "__main__":
+    print(FORMAT_STR.format('import main'), end='')
+    t = CLOCK()
+    from prog_algs import *
+    t2 = CLOCK()
+    print(t2-t)
+
+    from prog_models.models import BatteryElectroChemEOD as Battery
+    def future_loading(t, x = None):
+        # Variable (piece-wise) future loading scheme 
+        if (t < 600):
+            i = 2
+        elif (t < 900):
+            i = 1
+        elif (t < 1800):
+            i = 4
+        elif (t < 3000):
+            i = 2
+        else:
+            i = 3
+        return {'i': i}
+    R_vars = {
+        't': 2, 
+        'v': 0.02
+    }
+    batt = Battery(measurement_noise = R_vars)
+    initial_state = batt.parameters['x0']
+
+    print('PF')
+    print(FORMAT_STR.format('   Initialize'), end='')
+    temp_out = StringIO()
+    sys.stdout = temp_out
+    sys.stderr = temp_out
+    t = CLOCK()
+    filt = state_estimators.ParticleFilter(batt, initial_state)
+    t2 = CLOCK()
+    sys.stdout = sys.__stdout__
+    sys.stderr = sys.__stderr__
+    print(t2-t)
+
+    print(FORMAT_STR.format('   Step'), end='')
+    example_measurements = {'t': 32.2, 'v': 3.915}
+    t = 0.1
+    t = CLOCK()
+    filt.estimate(t, future_loading(t), example_measurements)
+    t2 = CLOCK()
+    print(t2-t)
+
+    print('UKF')
+    print(FORMAT_STR.format('   Initialize '), end='')
+    temp_out = StringIO()
+    sys.stdout = temp_out
+    sys.stderr = temp_out
+    t = CLOCK()
+    filt = state_estimators.UnscentedKalmanFilter(batt, initial_state)
+    t2 = CLOCK()
+    sys.stdout = sys.__stdout__
+    sys.stderr = sys.__stderr__
+    print(t2-t)
+
+    # print(FORMAT_STR.format('   Step'), end='')
+    # example_measurements = {'t': 32.2, 'v': 3.915}
+    # t = 0.1
+    # temp_out = StringIO()
+    # sys.stdout = temp_out
+    # sys.stderr = temp_out
+    # t = CLOCK()
+    # filt.estimate(t, future_loading(t), example_measurements)
+    # t2 = CLOCK()
+    # sys.stdout = sys.__stdout__
+    # sys.stderr = sys.__stderr__
+    # print(t2-t)
+
+    print(FORMAT_STR.format('Plot Results'), end='')
+    t = CLOCK()
+    filt.x.plot_scatter(label='prior')
+    t2 = CLOCK()
+    print(t2-t)
+
+    print('MC')
+    print(FORMAT_STR.format('   Initialize '), end='')
+    t = CLOCK()
+    mc = predictors.MonteCarlo(batt)
+    t2 = CLOCK()
+    print(t2-t)
+
+    print(FORMAT_STR.format('   Prediction'), end='')
+    NUM_SAMPLES = 5
+    STEP_SIZE = 0.1
+    t = CLOCK()
+    mc_results = mc.predict(batt.initialize(), future_loading, n_samples = NUM_SAMPLES, dt=STEP_SIZE)
+    t2 = CLOCK()
+    print(t2-t)
+
+    print(FORMAT_STR.format('Metrics'), end='')
+    t = CLOCK()
+    mc_results.time_of_event.percentage_in_bounds([3005.2, 3005.6])
+    mc_results.time_of_event.metrics(ground_truth=3005.25)
+    t2 = CLOCK()
+    print(t2-t)
+
+    print(FORMAT_STR.format('Plot Scatter'), end='')
+    t = CLOCK()
+    fig = mc_results.states.snapshot(0).plot_scatter(label = "t={} s".format(int(mc_results.times[0])))  # 0
+    quarter_index = int(len(mc_results.times)/4)
+    mc_results.states.snapshot(quarter_index).plot_scatter(fig = fig, label = "t={} s".format(int(mc_results.times[quarter_index])))  # 25%
+    mc_results.states.snapshot(quarter_index*2).plot_scatter(fig = fig, label = "t={} s".format(int(mc_results.times[quarter_index*2])))  # 50%
+    mc_results.states.snapshot(quarter_index*3).plot_scatter(fig = fig, label = "t={} s".format(int(mc_results.times[quarter_index*3])))  # 75%
+    mc_results.states.snapshot(-1).plot_scatter(fig = fig, label = "t={} s".format(int(mc_results.times[-1])))  # 100%
+    t2 = CLOCK()
+    print(t2-t)
+
+    print(FORMAT_STR.format('Plot Hist'), end='')
+    t = CLOCK()
+    mc_results.time_of_event.plot_hist()
+    t2 = CLOCK()
+    print(t2-t)
diff --git a/tests/test_metrics.py b/tests/test_metrics.py
index 2b87cc93..560576a6 100644
--- a/tests/test_metrics.py
+++ b/tests/test_metrics.py
@@ -454,10 +454,10 @@ def test_toe_profile_cumulative_relative_accuracy(self):
             )
         # Test positive floats ground truth
         GROUND_TRUTH = {'a': 9.0, 'b': 8.0, 'c': 18.0}
-        self.assertEquals(profile.cumulative_relative_accuracy(GROUND_TRUTH), {'a': 0.4444444444444445, 'b': 0.375, 'c': 0.6388888888888888})
+        self.assertEqual(profile.cumulative_relative_accuracy(GROUND_TRUTH), {'a': 0.4444444444444445, 'b': 0.375, 'c': 0.6388888888888888})
         # Test negative floats ground truth
         GROUND_TRUTH = {'a': -9.0, 'b': -8.0, 'c': -18.0}
-        self.assertEquals(profile.cumulative_relative_accuracy(GROUND_TRUTH), {'a': 3.555555555555556, 'b': 3.625, 'c': 3.305555555555556})
+        self.assertEqual(profile.cumulative_relative_accuracy(GROUND_TRUTH), {'a': 3.555555555555556, 'b': 3.625, 'c': 3.305555555555556})
         # Test ground truth values of 0; already caught by relative_accuracy
         with self.assertRaises(ZeroDivisionError):
             GROUND_TRUTH = {'a': 0, 'b': 0, 'c': 0}
@@ -466,7 +466,70 @@ def test_toe_profile_cumulative_relative_accuracy(self):
         with self.assertRaises(TypeError):
             GROUND_TRUTH = []
             raise_error = profile.cumulative_relative_accuracy(GROUND_TRUTH)
-    
+
+    def test_toe_profile_monotonicity(self):
+        from prog_algs.predictors import ToEPredictionProfile
+
+        # Test monotonically increasing and decreasing
+        profile = ToEPredictionProfile()  # Empty profile
+        for i in range(10):
+            data = [{'a': 1+i/10, 'b': 2-i/5} for i in range(10)]
+            profile.add_prediction(
+                i,  # Time (reverse so data is decreasing)
+                UnweightedSamples(data)  # ToE Prediction
+            )
+        self.assertDictEqual(profile.monotonicity(), {'a': 1.0, 'b': 1.0})
+
+        # Test no monotonicity
+        profile = ToEPredictionProfile()  # Empty profile
+        for i in range(0,11):
+            data = [{'a': i + i*(i%2-0.5), 'b': i + i*(i%2-0.5)}]*10
+            profile.add_prediction(
+                i,  # Time (reverse so data is decreasing)
+                UnweightedSamples(data)  # ToE Prediction
+            )
+        self.assertDictEqual(profile.monotonicity(), {'a': 0.0, 'b': 0.0})
+
+        # Test monotonicity between range [0,1]
+        profile = ToEPredictionProfile()  # Empty profile
+        for i in range(11):
+            data = [{'a': i+i*(i%3-0.5), 'b': i+i*(i%3-0.5)}] * 10
+            profile.add_prediction(
+                i,  # Time (reverse so data is decreasing)
+                UnweightedSamples(data)  # ToE Prediction
+            )
+        self.assertDictEqual(profile.monotonicity(), {'a': 0.4, 'b': 0.4})
+
+        # Test mixed
+        profile = ToEPredictionProfile()  # Empty profile
+        for i in range(11):
+            data = [{'a': i+i*(i%3-0.5), 'b': i + i*(i%2-0.5), 'c': 1+i/10}] * 10
+            profile.add_prediction(
+                i,  # Time (reverse so data is decreasing)
+                UnweightedSamples(data)  # ToE Prediction
+            )
+        self.assertDictEqual(profile.monotonicity(), {'a': 0.4, 'b': 0.0, 'c': 1.0})
+
+        # Test MultivariateNormalDist
+        profile = ToEPredictionProfile()  # Empty profile
+        covar = [[0.1, 0.01], [0.01, 0.1]]
+        for i in range(11):
+            data = [{'a': i, 'b': i*(i%3+5)}] * 11
+            profile.add_prediction(
+                i,  # Time (reverse so data is decreasing)
+                MultivariateNormalDist(data[i].keys(), data[i].values(), covar)  # ToE Prediction
+            )
+        self.assertDictEqual(profile.monotonicity(), {'a': 0.0, 'b': 0.5})
+
+        # Test ScalarData
+        profile = ToEPredictionProfile()  # Empty profile
+        for i in range(11):
+            data = {'a': i+i*(i%3-0.5), 'b': i+i*(i%3-0.5)}
+            profile.add_prediction(
+                i,  # Time (reverse so data is decreasing)
+                ScalarData(data)  # ToE Prediction
+            )
+        self.assertDictEqual(profile.monotonicity(), {'a': 0.4, 'b': 0.4})
 
 
 # This allows the module to be executed directly    
diff --git a/tests/test_predictors.py b/tests/test_predictors.py
index 998ff8d3..34a1d1eb 100644
--- a/tests/test_predictors.py
+++ b/tests/test_predictors.py
@@ -83,9 +83,9 @@ def future_loading(t, x={}):
             return {}
 
         mc_results = pred.predict(samples, future_loading, dt=0.01, events=['falling'], save_freq=1)
-        self.assertAlmostEqual(mc_results.time_of_event.mean['falling'], 4.15, 0)
+        self.assertAlmostEqual(mc_results.time_of_event.mean['falling'], 3.8, 0)
         self.assertTrue('impact' not in mc_results.time_of_event.mean)
-        self.assertAlmostEqual(mc_results.times[-1], 4, 1)  # Saving every second, last time should be around the nearest 1s before falling event
+        self.assertAlmostEqual(mc_results.times[-1], 3, 1)  # Saving every second, last time should be around the nearest 1s before falling event
 
     def test_UKP_Battery(self):
         def future_loading(t, x = None):
diff --git a/tests/test_state_estimators.py b/tests/test_state_estimators.py
index 37d8aad9..49c37524 100644
--- a/tests/test_state_estimators.py
+++ b/tests/test_state_estimators.py
@@ -234,10 +234,10 @@ def test_UKF_incorrect_input(self):
         self.__incorrect_input_tests(UnscentedKalmanFilter)
 
     def test_PF(self):
-        m = ThrownObject(process_noise={'x': 1, 'v': 3}, measurement_noise=1, num_particles = 1000)
+        m = ThrownObject(process_noise={'x': 0.75, 'v': 0.75}, measurement_noise=1)
         x_guess = {'x': 1.75, 'v': 38.5} # Guess of initial state, actual is {'x': 1.83, 'v': 40}
 
-        filt = ParticleFilter(m, x_guess)
+        filt = ParticleFilter(m, x_guess, num_particles = 1000)
         self.__test_state_est(filt, m)
 
         # Test ParticleFilter ScalarData
@@ -273,69 +273,15 @@ def test_PF(self):
                 self.assertAlmostEqual(filt_us.x.cov[i][j], x_us.cov[i][j], delta=0.1)
 
         # Test x0 if-else Control
-        # Case 0: Both isinstance(x0, UncertainData) and x0_uncertainty parameter provided; expect x0_uncertainty to be skipped
+        # Case 0: isinstance(x0, UncertainData) 
         x_scalar = ScalarData({'x': 1.75, 'v': 38.5}) # Testing with ScalarData
-        filt_scalar = ParticleFilter(m, x_scalar, num_particles = 20, x0_uncertainty = 0.5) # Sample count does not affect ScalarData testing
+        filt_scalar = ParticleFilter(m, x_scalar, num_particles = 20) # Sample count does not affect ScalarData testing
         mean1 = filt_scalar.x.mean
         mean2 = x_scalar.mean
         self.assertSetEqual(set(mean1.keys()), set(mean2.keys()))
         for k in mean1.keys():
             self.assertEqual(mean1[k], mean2[k])
         self.assertTrue((filt_scalar.x.cov == x_scalar.cov).all())
-        # Case 1: Only x0_uncertainty provided; expect a warning issued
-        with self.assertWarns(Warning):
-            filt_scalar = ParticleFilter(m, {'x': 1.75, 'v': 38.5}, x0_uncertainty = 0.5)
-        # Case 2: Raise ProgAlgTypeError if x0 not UncertainData or x0_uncertainty not of type {{dict, Number}}.
-        with self.assertRaises(ProgAlgTypeError):
-            filt_scalar = ParticleFilter(m, {'x': 1.75, 'v': 38.5}, num_particles = 20, x0_uncertainty = [])
-
-    def test_measurement_eq_UKF(self):
-        m = MockProgModel2()
-        x0 = m.initialize()
-
-        # Setup
-        filt = UnscentedKalmanFilter(m, x0)
-        
-        # Try using
-        filt.estimate(0.2, {'i1': 1, 'i2': 2}, {'o1': -2.0, 'o2': 7})
-
-        # Add Measurement eqn
-        def measurement_eqn(x):
-            z = m.output(x)
-            del z['o2']
-            return z
-        filt = UnscentedKalmanFilter(m, x0, measurement_eqn=measurement_eqn)
-        filt.estimate(0.1, {'i1': 1, 'i2': 2}, {'o1': -2.0})
-
-        # New Measurement eqn method
-        class MyModel(MockProgModel2):
-            outputs = ['o1']
-        filt = UnscentedKalmanFilter(MyModel(), x0)
-        filt.estimate(0.1, {'i1': 1, 'i2': 2}, {'o1': -2.0})
-
-    def test_measurement_eq_PF(self):
-        m = MockProgModel2()
-        x0 = m.initialize()
-
-        # Setup
-        filt = ParticleFilter(m, x0)
-        
-        # This one should work
-        filt.estimate(0.2, {'i1': 1, 'i2': 2}, {'o1': -2.0, 'o2': 7})
-
-        # Add Measurement eqn
-        def measurement_eqn(x):
-            z = m.output(x)
-            del z['o2']
-            return z
-        filt = ParticleFilter(m, x0, measurement_eqn=measurement_eqn)
-        filt.estimate(0.1, {'i1': 1, 'i2': 2}, {'o1': -2.0}) 
-
-        # New Measurement eqn method
-        class MyModel(MockProgModel2):
-            outputs = ['o1']
-        filt = ParticleFilter(MyModel(), x0)
-        filt.estimate(0.1, {'i1': 1, 'i2': 2}, {'o1': -2.0})
         
     def test_PF_incorrect_input(self):
         self.__incorrect_input_tests(ParticleFilter)