This example shows how to use sensitivity analysis to narrow down the number of parameters that you need to estimate when fitting a model. This example uses a model of the vestibulo-ocular reflex, which generates compensatory eye movements.
The vestibulo-ocular reflex (VOR) enables the eyes to move at the same speed and in the opposite direction as the head, so that vision is not blurred when the head moves during normal activity. For example, if the head turns to the right, the eyes turn to the left at the same speed. This happens even in the dark. In fact, the VOR is most easily characterized by measurements in the dark, to ensure that eye movements are predominantly driven by the VOR.
Head rotation is sensed by organs in the inner ears, known as semicircular canals. These detect head motion and transmit signals about head motion to the brain, which sends motor commands to the eye muscles, so that eye movements compensate for head motion. We would like to use eye movement data to estimate model parameters for these various stages. The model we will use is shown below. There are four parameters in the model:
sdoVOR_Data.mat contains uniformly sampled data of stimulation and eye movements. If the VOR were perfectly compensatory, then a plot of eye movement data, when flipped vertically, would overlay exactly on top of a plot of head motion data. Such a system would be described by a gain of 1 and a phase of 180 degrees. However, real eye movements are close, but not perfectly compensatory.
load sdoVOR_Data.mat; % Column vectors: Time HeadData EyeData
We will use Sensitivity Analysis UI to see how well the model output fits the data, and explore which model parameters have the most influence on the goodness of fit. To open Sensitivity Analysis UI, navigate to the Analysis menu in the Simulink model and click Sensitivity Analysis...
To associate the data with the model, click New Requirement and select a Signal Matching requirement. This specifies an objective function consisting of the sum of squared error between the data and model output. In the Signal Matching dialog, specify the output as
[Time EyeData], and specify the input as
To view the eye movement data, navigate to the data browser on the left side of the UI, right-click the
SignalMatching requirement, and select Plot & Simulate. The bottom plot shows the stimulation, consisting of a series of pulses. The top plot shows eye movement data, which resembles but does not exactly match the stimulation. It also shows that the model simulated output does not match the eye movement data, because model parameters need to be estimated.
The model attempts to capture the phenomena which cause the difference between head movements and eye movements. Here we will explore the design space formed by the model parameters. To specify the parameters to explore in Sensitivity Analysis UI, click Select Parameters and create a new parameter set. Select all model parameters:
Explore the design space by generating parameter values. Click Generate Values and select random values. For repeatability of the example, reset the random number generator.
Since there are 4 parameters, we will generate 40 samples.
Delay parameter models the fact that there is some delay in communicating the signals from the inner ear to the brain and eyes. This delay is due to the time needed for chemical neurotransmitters to traverse the synaptic clefts between nerve cells. Based on the number of synapses in the vestibulo-ocular reflex, this delay is expected to be around 5 ms. We will model it with a uniform distribution with a lower bound of 2 ms and an upper bound of 9 ms.
Gain parameter models the fact that in the dark, the eyes do not move quite as much as the head. We will model it with a uniform distribution with a lower bound of 0.6 and an upper bound of 1.
Tc parameter models the dynamics associated with the semicircular canals, as well as some additional neural processing. The canals are high-pass filters, because after a subject has been put into rotational motion, the neurally active membranes in the canals slowly relax back to resting position, so the canals stop sensing motion. Thus, after the stimulation undergoes transition edges, eye movements tend to depart from the stimulation over time. Based on mechanical characteristics of the canals, combined with additional neural processing which prolongs this time constant to improve the accuracy of the VOR, we will model
Tc with a normal (i.e., bell curve) distribution with a mean of 15 seconds and a standard deviation of 3 seconds.
Tp parameter models the dynamics of the oculomotor plant, i.e. the eye and the muscles and tissues attached to it. The plant can be modeled by two poles, however it is believed that the pole with the larger time constant is cancelled by precompensation in the brain, to enable the eye to make quick movements. Thus in the plot, when the stimulation undergoes transition edges, the eye movements follow with only a little delay. We will model
Tp with a uniform distribution with a lower bound of 0.005 seconds and an upper bound of 0.05 seconds.
When the sample values are generated, they appear in a table in the Sensitivity Analysis UI. To plot them, select
ParamSet in the data browser, click the Plots tab, and make a scatter plot. The sampling above used default options, and these are reflected in the scatter plot. For parameters modeled by a uniform distribution, the histograms appear approximately uniform. However, parameter
Tc was modeled by a normal distribution, and its histogram has a bell curve profile. If Statistics and Machine Learning Toolbox™ is available, many other distributions may be used, and sampling can be done using Sobol or Halton low-discrepancy sequences. The off-diagonal plots show scatter plots between pairs of different variables. Since we did not specify cross-correlations between parameters, the scatter plots appear uncorrelated. However, if parameters were believed to be correlated, this can be specified using Correlation Matrix tab in the dialog for generating random parameter values.
Now that we have generated values for the parameter set and specified a requirement (
SignalMatching), we can evaluate the model. In the Sensitivity Analysis tab, click Evaluate Model.
The model is run once for each set of parameter values, and the results scatter plot is updated as new computations become available. Evaluation could also be sped up using parallel computing. After evaluation is complete, all results are also displayed in a table.
From the scatter plot of evaluation results, the
SignalMatching requirement seems to vary systematically as a function of
Tc, but not
Tp. Something similar can be seen in a contour plot. Select the
EvalResults variable in the data browser, click the Plots tab, and make a contour plot. The requirement does not vary systematically from left to right as a function of
Delay, but it does vertically as a function of
We can use statistical analysis to quantify how much each parameter influences the requirement. Click the Statistics tab and select both correlation and standardized regression; and both linear and ranked analysis types. If Statistics and Machine Learning Toolbox is available, partial correlation and Kendall correlation can also be selected. Click Compute Statistics to carry out the calculations and show a tornado plot. The tornado plot displays results from top to bottom in order of which parameter most influences the requirement. The statistical values range from -1 to 1, where the magnitude indicates how much the parameter influences the requirement, and the sign indicates whether an increase in the parameter value corresponds to an increase or decrease in the requirement value. By most measures, this
SignalMatching requirement is more sensitive to
Tc, and less sensitive to
For parameter estimation, we need to specify starting values for the parameters. Click the evaluation results table and click the
SignalMatching column header to sort results. Select the row of parameter values that minimizes the
SignalMatching requirement. Right-click on the row and extract these parameter values. A new variable,
ParamValues, is shown in the data browser.
To transition from sensitivity analysis to parameter estimation, navigate to the Sensitivity Analysis tab, click Optimize, and open a parameter estimation session. In the dialog that appears, specify that you want to use the parameter values in
ParamValues, and the
Since we found above that parameters
Tc have the most influence on the value of
SignalMatching, we would like to estimate only these two parameters, since the time for estimation increases with the number of parameters being estimated. In Parameter Estimation UI, click Select Parameters and select only
Tc for estimation.
Since the experiment definition has been imported from
SignalMatching and the parameter values have been imported from
ParamValues, we have everything needed for estimation. Click Estimate to carry out parameter estimation for
Tc. Because we are only estimating the two most influential parameters, estimation converges quickly and the model output closely matches the data. As was the case with model evaluations in sensitivity analysis, parallel computing could be used to speed up estimation.
In summary, Sensitivity Analysis UI was used to explore the parameter design space and determine that two parameters,
Tc, were substantially more influential than the others. A start point for estimation was also determined. This start point and the requirement of obtaining a good fit to experimental data were imported into Parameter Estimation UI. Estimation completed quickly because only two parameters needed to be estimated, and the model output fit the data with very little residual error.
Close the model