On the 17th of April from 18:00 to 20:00, We will be upgrading GitLab. During this time the service will be unavailable, apologies for any inconvenience
"<div style=\"background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\"> <p><strong>Note on code implementation</strong>: you'll see that the functions in this assignment use a dictionary; this greatly reduces the number of input/output variables needed in a function. However, it can make the code inside the function more difficult to read due to the key syntax (e.g., <code>dict['variable_1']</code> versus <code>variable\n",
"_1</code>). To make this assignment easier for you to implement we have split these functions into three parts: 1) define variables from the dictionary, 2) perform analysis, 3) add results to the dictionary. Note that this is not the most efficient way to write this code; it is done here specifically for clarity and to help you focus on writing the equations properly and understanding the meaning of each term.</p></div>"
]
},
{
"cell_type": "code",
"execution_count": 18,
...
...
@@ -658,22 +663,9 @@
"def BLUE(d):\n",
" \"\"\"Calculate the Best Linear Unbiased Estimator\n",
" \n",
" Write a docstring here (an explanation of your function).\n",
" \n",
" Function to calculate the Best Linear Unbiased Estimator\n",
" \n",
" Input is a dictionary with the following keys/values:\n",
" A = A matrix (mxn)\n",
" y = vector with obervations (mx1)\n",
" Sigma_Y = Varaiance covariance matrix of the observations (mxm)\n",
" \n",
" Output is the same dictionary with the following keys/values added:\n",
" x_hat = vector with the estimates (nx1)\n",
" Sigma_X_hat = variance-covariance matrix of the unknown parameters (nxn)\n",
"\n",
" Note: it is not necessary to define the intermediate variables that are already\n",
" in the dictionary, however, it is done here for clarity so you can focus\n",
" on writing the equations properly.\n",
" Uses dict as input/output:\n",
" - inputs defined from existing values in dict\n",
"Complete the function below to help us compute the confidence intervals, then apply the function.\n",
"Complete the function below to help us compute the confidence intervals, then apply the function. Use a confidence interval of 96% in your analysis.\n",
"\n",
"<em>Hint: it can be used in exactly the same way as the <code>BLUE</code> function above, although it has one extra input.</em>\n",
"</p>\n",
...
...
@@ -853,18 +845,19 @@
"def get_CI(d, alpha):\n",
" \"\"\"Compute the confidence intervals.\n",
" \n",
" Input alpha is a float that defines the confidence level.\n",
"\n",
" As with BLUE, complete the equations for the CI.\n",
" Uses dict as input/output:\n",
" - inputs defined from existing values in dict\n",
" - outputs defined as new values in dict\n",
" \"\"\"\n",
"\n",
" std_e_hat = d['std_e_hat']\n",
" std_y = d['std_y']\n",
"\n",
" k99 = YOUR_CODE_HERE\n",
" k = YOUR_CODE_HERE\n",
" CI_y = YOUR_CODE_HERE\n",
" CI_res = YOUR_CODE_HERE\n",
"\n",
" d['alpha'] = alpha\n",
" d['CI_y'] = CI_y\n",
" d['CI_res'] = CI_res\n",
"\n",
...
...
@@ -914,7 +907,7 @@
"<p>\n",
"<b>Task 6.2:</b> \n",
" \n",
"Read the contents of file <code>figures.py</code> and identify what it is doing: you should be able to recognize that they use our model dictionary as an input and create three different figures. Note also that the function to create the figures have already been imported at the top of this notebook.\n",
"Read the contents of file <code>functions.py</code> and identify what it is doing: you should be able to recognize that they use our model dictionary as an input and create three different figures. Note also that the function to create the figures have already been imported at the top of this notebook.\n",
"\n",
"Use the functions provided to visualize the results of our two models.\n",
"</p>\n",
...
...
%% Cell type:markdown id:785b4e1d tags:
# GA 1.3: Modelling Road Deformation using Non-Linear Least-Squares
*[CEGM1000 MUDE](http://mude.citg.tudelft.nl/): Week 1.3. Due: Friday, September 20, 2024.*
%% Cell type:code id:181ccfd5 tags:
``` python
importnumpyasnp
importscipyassc
fromscipyimportinterpolate
fromscipy.statsimportnorm
importpandasaspd
importmatplotlib.pyplotasplt
fromdatetimeimportdatetime
fromdatetimeimporttimedelta
fromscipy.stats.distributionsimportchi2
fromscipy.statsimportnorm
fromfunctionsimport*
np.set_printoptions(precision=3)
```
%% Cell type:markdown id:5ca94e0e tags:
## Part 0: Dictionary Review
As described above, several functions in this assignment require the use of a Python dictionary to make it easier to keep track of important data, variables and results for the various _models_ we will be constructing and validating.
_It may be useful to revisit PA 1.1, where there was a brief infroduction to dictionaires. That PA contains all the dictionary info you need for GA 1.3. A [read-only copy is here](https://mude.citg.tudelft.nl/2024/files/Week_1_1/PA_1_1_Catch_Them_All.html) and [the source code (notebook) is here](https://gitlab.tudelft.nl/mude/2024-week-1-1)._
Test your knowledge by adding a new key <code>new_key</code> and then executing the function to print the value.
</p>
</div>
%% Cell type:code id:41c56f43 tags:
``` python
YOUR_CODE_HERE
function_that_uses_my_dictionary(my_dictionary)
```
%% Cell type:markdown id:160d6250 tags:
## Task 1: Preparing the data
Within this assignment you will work with two types of data: InSAR data and GNSS data. The cell below will load the data and visualize the observed displacements time. In this task we use the package `pandas`, which is really useful for handling time series. We will learn how to use it later in the quarter; for now, you only need to recognize that it imports the data as a `dataframe` object, which we then convert into a numpy array using the code below.
%% Cell type:markdown id:02b12781 tags:
<divstyle="background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px"><p>Tip: note that we have converted all observations to millimeters.</p></div>
Once you have used the cell above to import the data, investigate the data sets using the code cell below. Then provide some relevant summary information in the Markdown cell.
<em>Hint: at the least, you should be able to tell how many data points are in each data set and get an understanding of the mean and standard deviation of each. Make sure you compare the different datasets and use consistent units.</em>
You may have noticed that the groundwater data is available for different times than the GNSS and InSAR data. You will therefore have to *interpolate* the data to the same times for a further analysis. You can use the SciPy function ```interpolate.interp1d``` (read its [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html)).
The cells below do the following:
1. Define a function to convert the time unit
2. Convert the time stamps for all data
3. Use `interp1d` to interpolate the groundwater measurements at the time of the satellite measurements
%% Cell type:code id:f02ed4c4 tags:
``` python
defto_days_years(times):
'''Convert the observation times to days and years.'''
Describe the datasets based on the figure above and your observations from the previous tasks. What kind of deformation do you see? And what are the differences between both datasets? Be quantitative.
</p>
</div>
%% Cell type:markdown id:f9a7bdd4 tags:
Before we move on, it is time to do a little bit of housekeeping.
Have you found it confusing to keep track of two sets of variables---one for each data type? Let's use a dictionary to store relevant information about each model. We will use this in the plotting functions for this task (and again next week), so make sure you take the time to see what is happening. Review also Part 0 at the top of this notebook if you need a refresher on dictionaries.
Run the cell below to define a dictionary for storing information about the two (future) models.
</p>
</div>
%% Cell type:code id:2c27b4f3 tags:
``` python
model_insar={'data_type':'InSAR',
'y':y_insar,
'times':times_insar,
'groundwater':GW_at_InSAR_times
}
model_gnss={'data_type':'GNSS',
'y':y_gnss,
'times':times_gnss,
'groundwater':GW_at_GNSS_times
}
```
%% Cell type:markdown id:76c9115b tags:
## Task 2: Set-up linear functional model
We want to investigate how we could model the observed displacements of the road. Because the road is built in the Green Heart we expect that the observed displacements are related to the groundwater level. Furthermore, we assume that the displacements can be modeled using a constant velocity. The model is defined as
$$
d = d_0 + vt + k \ \textrm{GW},
$$
where $d$ is the displacement, $t$ is time and $\textrm{GW}$ is the groundwater level (that we assume to be deterministic).
Therefore, the model has 3 unknowns:
1. $d_0$, as the initial displacement at $t_0$;
2. $v$, as the displacement velocity;
3. $k$, as the 'groundwater factor', which can be seen as the response of the soil to changes in the groundwater level.
As a group you will construct the **functional model** that is defined as
Add the A matrix to the dictionaries for each model. This will be used to plot results later in the notebook.
</p>
</div>
%% Cell type:code id:396ac3a5 tags:
``` python
model_insar['A']=YOUR_CODE_HERE
model_gnss['A']=YOUR_CODE_HERE
```
%% Cell type:markdown id:9325d32b tags:
## 3. Set-up stochastic model
We will use the Best Linear Unbiased Estimator (BLUE) to solve for the unknown parameters. Therefore we also need a stochastic model, which is defined as
$$
\mathbb{D}(Y) = \Sigma_{Y}.
$$
where $\Sigma_{Y}$ is the covariance matrix of the observables' vector.
Write a function to apply BLUE in the cell below and use the function to estimate the unknowns for the model using the data.
Compute the modeled displacements ($\hat{\mathrm{y}}$), and corresponding residuals ($\hat{\mathrm{\epsilon}}$), as well as associated values (as requested by the blank code lines).
</p>
</div>
%% Cell type:markdown id:c5bd4bad tags:
<divstyle="background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px"><p><strong>Note on code implementation</strong>: you'll see that the functions in this assignment use a dictionary; this greatly reduces the number of input/output variables needed in a function. However, it can make the code inside the function more difficult to read due to the key syntax (e.g., <code>dict['variable_1']</code> versus <code>variable
_1</code>). To make this assignment easier for you to implement we have split these functions into three parts: 1) define variables from the dictionary, 2) perform analysis, 3) add results to the dictionary. Note that this is not the most efficient way to write this code; it is done here specifically for clarity and to help you focus on writing the equations properly and understanding the meaning of each term.</p></div>
%% Cell type:code id:d85b1826 tags:
``` python
defBLUE(d):
"""Calculate the Best Linear Unbiased Estimator
Write a docstring here (an explanation of your function).
Function to calculate the Best Linear Unbiased Estimator
Input is a dictionary with the following keys/values:
A = A matrix (mxn)
y = vector with obervations (mx1)
Sigma_Y = Varaiance covariance matrix of the observations (mxm)
Output is the same dictionary with the following keys/values added:
x_hat = vector with the estimates (nx1)
Sigma_X_hat = variance-covariance matrix of the unknown parameters (nxn)
Note: it is not necessary to define the intermediate variables that are already
in the dictionary, however, it is done here for clarity so you can focus
Do the values that you just estimated make sense? Explain, using quantitative results.
<em>Hint: all you need to do is use the figures created above to verify that the parameter values printed above are reasonable (e.g., order of magnitude, units, etc).</em>
Complete the function below to help us compute the confidence intervals, then apply the function.
Complete the function below to help us compute the confidence intervals, then apply the function. Use a confidence interval of 96% in your analysis.
<em>Hint: it can be used in exactly the same way as the <code>BLUE</code> function above, although it has one extra input.</em>
</p>
</div>
%% Cell type:code id:2711da12 tags:
``` python
defget_CI(d,alpha):
"""Compute the confidence intervals.
Input alpha is a float that defines the confidence level.
As with BLUE, complete the equations for the CI.
Uses dict as input/output:
- inputs defined from existing values in dict
- outputs defined as new values in dict
"""
std_e_hat=d['std_e_hat']
std_y=d['std_y']
k99=YOUR_CODE_HERE
k=YOUR_CODE_HERE
CI_y=YOUR_CODE_HERE
CI_res=YOUR_CODE_HERE
d['alpha']=alpha
d['CI_y']=CI_y
d['CI_res']=CI_res
returnd
```
%% Cell type:code id:d9a41ea5 tags:
``` python
model_insar=YOUR_CODE_HERE
model_gnss=YOUR_CODE_HERE
```
%% Cell type:markdown id:53cf3663 tags:
At this point we have all the important results entered in our dictionary and we will be able to use the plots that have been written for you in the next Tasks. In case you would like to easily see all of the key-value pairs that have been added to the dictionary, you can run the cell below:
Read the contents of file <code>figures.py</code> and identify what it is doing: you should be able to recognize that they use our model dictionary as an input and create three different figures. Note also that the function to create the figures have already been imported at the top of this notebook.
Read the contents of file <code>functions.py</code> and identify what it is doing: you should be able to recognize that they use our model dictionary as an input and create three different figures. Note also that the function to create the figures have already been imported at the top of this notebook.
Use the functions provided to visualize the results of our two models.
<divstyle="background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px"><p><strong>Note</strong>: remember that you will have to use the same function to look at <em>both</em> models when writing your interpretation in the Report.</p></div>
