Hi Andres,
Apologies for the late response.
I decided to make a few enhancements in Pelicun to streamline the process of adding a new distribution type. I've just released v3.1.b6 and also updated the Example notebook on DesignSafe with a lot of additional details and explanation on how Pelicun 3 works. I encourage you to take a look here:
https://www.designsafe-ci.org/data/browser/public/designsafe.storage.published/PRJ-3411v5Read on if you are still interested in adding the skew normal distribution.
First, you'll need to have the source code available locally. I also recommend linking this version of pelicun to your active Python on the local system to make testing easier. I am not sure how familiar you are with the steps to do this, so let me give you a few tips - please don't hesitate to ask me if you need more information:
- You should get pelicun from the main NHERI-SimCenter account on GitHub. The master branch always provides the latest stable release. Since you are interested in extending the code, you'll need the develop branch; you can find that here:
https://github.com/NHERI-SimCenter/pelicun/tree/develop- I suggest forking the NHERI-SimCenter/pelicun repo and then cloning the develop branch of your own fork to your local hard drive. This will give you the latest version of the code base.
- In your local system, you can store pelicun in any location and make sure Python finds it by adding it to the PYTHONPATH environment variable or by using the externals.pth file. I am happy to provide more information on either of these if you are not familiar with them.
- Once you set things up properly, you should see 3.1.b6 when you import pelicun and check pelicun.__ version__ in a python interpreter.
- Now you can make your edits, test the resulting code and, when you are happy with it, I would appreciate if you committed your contributions to the main repo so that everyone can use this new distribution. To do that, you'll need to first commit your changes to the develop branch on your fork of pelicun and then open a pull request from your develop branch to the develop branch of the main repo. When I see your pull request, I'll test your version and - if everything looks okay - I'll accept it.
That's it, that's how you can extend pelicun and share your work with the community.
Now, let's see how you would go about adding a new distribution. You'll only need to make changes in the uq module to do this. See the latest version of the script here for reference:
https://github.com/NHERI-SimCenter/pelicun/blob/develop/pelicun/uq.py . Note that if the uq module gets updated by the time you are reading this, you can always click on history in the top right and go back to today's version so that the line numbers I give below will point to the right location.
Sampling an N dimensional multivariate distribution in Pelicun follows the logic below (see generate_sample method starting at line 1465):
- Sample an N dimensional uniform distribution using Monte Carlo or Latin Hypercube Sampling
- If there are variables across the N dimensions with non-zero correlation, apply the prescribed correlations assuming a Gaussian copula function. First, we try to use a fast Cholesky transformation; if this fails because the correlation matrix is not positive semidefinite, we use another method based on Singular Value Decomposition to apply correlations that preserve as much as possible from the prescribed correlation matrix.
- Perfect correlation can be handled very efficiently with a special 'anchor' feature, but that is almost surely outside of scope for your edits; let me know if you want to know more about it.
- Finally, we take each marginal and use inverse probability integral transformation to transform the sample from uniform distribution to the desired distribution.
Adding a new distribution only affects the last step of the above (as long as you don't want to have some special, non-Gaussian copula function also included - let me know if you do). Below is a list of locations where you'd need to make edits:
- Add the name of the new distribution and a description of its parameters in the documentation at the top of the RandomVariable class - line 791-827; add more info at least under 'distribution' and 'theta' parameters, but potentially also to others, if needed. Currently, Pelicun supports up to three parameters for distributions. That should be sufficient for the skew normal, but let me know if you need more than three.
- If the new distribution requires some special checks, you can add those in the init method (line 829). See, for example, that the multinomial distribution is checked to have probabilities sum up to at most 1.0.
- The cdf() method (line 994) returns the Cumulative Distribution Function ordinates for a given set of x values. You'll need to add your distribution here in an elif clause. Take a look at how the existing distributions are handled and try to mimic the same robustness:
- You'll see that we get the parameters of the distribution in the theta vector and truncation limits in another vector. Truncation limits can be undefined which should lead to an unbounded distribution in one or both ends and not throw an error.
- As for the parameters of the distribution, they are typically mandatory, but in some cases we can have rules set up to replace missing values - e.g., see how the limits of the uniform distribution are optional and replaced with infinite values if missing.
- Make sure the input values are valid inputs to the CDF. For example, 0 and negative numbers are not part of the input domain for the lognormal CDF, so we need to make sure zeros are replaced by the smallest positive number (the computer can handle) in line 1050. This also shows that in general I try to make Pelicun work and handle issues gracefully rather than throwing error messages all the time when something out of the ordinary happens. So, someone feeding a zero to a lognormal CDF will get a probability of 0. I believe in most cases this behavior is preferred over the error message that would terminate execution.
- I suggest starting with the implementation of the non-truncated version of the distribution. After testing and making sure that it works, you shall expand it by adding the truncation option. Let me know if you need help with this.
- The next method to edit is the inverse_transform() starting at line 1070. Here, you'll need to add an elif clause; I'd add it after line 1144 because the distribution you plan to add belongs to the normal family. The task is similar to the cdf method, but you are implementing an inverse cdf function here:
- Note that the sample_size argument of the method is only used in special cases. If you implement a skew normal distribution, you should expect to have an array of values (which come from the [0,1] uniform distribution following the sampling logic I explained earlier) that your script will transform to the target distribution.
- Pay attention to handling missing inputs and truncation limits - the advice I gave for the cdf() method applies here too.
- Finally, you'll need to edit the scale_distribution() method at the top of the uq module (line 70). This method is used by various models in Pelicun to scale input parameters defined in one unit to the SI unit that is used internally. Scaling distributions so far was straightforward as you'll see in the implementation - I hope the skew normal will not be an exception. Note that I define the normal distribution with the mean and coefficient of variation so that the second parameter is unitless and does not need to be scaled. It might be worth pulling similar tricks with the skew normal.
That's it. At this point you have the skew normal distribution implemented and you can sample it by creating a RandomVariable object, adding it to a RandomVariableRegistry and then calling the generate_sample() method of that registry. This is done all over the model module:
https://github.com/NHERI-SimCenter/pelicun/blob/develop/pelicun/model.py For example, take a look at the generate_cmp_sample() method starting at line 1165. The _create_cmp_RVs() creates the RandomVariable objects and the registry (line 1139) and then at line 1180 we sample the registry.
Notice that RandomVariable objects are created by feeding the family and parameters of the distribution to the uq module automatically. Using the previous cmp example, take a look at line 1150 in the model module. As long as you feed in the correct name for the family and valid parameters, your new distribution should work immediately without making any changes to the model.py .
The only exception to this is the LossModel - currently, only normal and lognormal random variables are supported for probabilistic loss calculation. Let me know if you want to include the skewed normal there and I can help you set that up.
I hope you'll find the above helpful. Please let me know if you're working on this and don't hesitate to ask questions here if you run into any issues.
