Hi Phil,
We have A SciPyDistribution class for that purpose in OpenTURNS but unfortunately it is currently broken and will be fixed for the 1.10 version to be released soon.
Luckily it's pure Python, so here is a script that redefines SciPyDistribution with the fixed range computation used with the johnsonSU distribution.
j
De : [email protected] <[email protected]> de la part de Phil Fernandes <[email protected]>
Envoyé : samedi 11 novembre 2017 00:30:54
À : regis lebrun; [email protected]
Objet : Re: [ot-users] [External] Re: Custom distribution in FORM
Envoyé : samedi 11 novembre 2017 00:30:54
À : regis lebrun; [email protected]
Objet : Re: [ot-users] [External] Re: Custom distribution in FORM
I tried implementing via distribution algebra, but for some reason my program just hangs, so I decided to try implementing the distribution as a subclass of PythonDistribution as per the example
https://github.com/openturns/openturns/blob/master/python/test/t_Distribution_python.py. Unfortunately there seem to be inconsistencies in the argument types that are accepted by the methods of Distribution objects and the example PythonDistribution. What
are the required object types for the outputs of computeXXX(), e.g., computeQuantile() in order for the distribution to work in FORM? Would it suffice to output a list?
For example
a = ot.Normal()
x = np.linspace(0,1,5)[:,None]
a.computePDF(x)
returns
class=Sample name=Unnamed implementation=class=SampleImplementation name=Unnamed size=5 dimension=1 data="">
However
b=UniformNdPy()
b.computePDF(x)
returns 1.0.
Thank you.
-----Original Message-----
From: regis lebrun [mailto:[email protected]]
Sent: Friday, November 10, 2017 12:12 PM
To: [email protected]; Phil Fernandes
Subject: [External] Re: [ot-users] Custom distribution in FORM
Hi,
You can easily implement this distribution using OpenTURNS unique feature regarding distribution algebra (see https://en.wikipedia.org/wiki/Johnson%27s_SU-distribution):
import openturns as ot
lamb = 1.5
xi = 1.1
delta = 2.0
gamma = 1.0
distJU = ((ot.Normal() - gamma) / delta).sinh() * lamb + xi print("distJU=", distJU)
ot.Show(distJU.drawPDF())
You will get:
distJU= RandomMixture(1.1 + 1.5 * CompositeDistribution=f(RandomMixture(Normal(mu = -0.5, sigma = 0.5))) with f=[x]->[sinh(x)])
and a graph similar to the one given on the wikipedia page.
Your script contains some bugs (computeQuantile, getMean, getStandardDeviation should return the result as a float sequence of size 1) and a missing method, namely getRange(). To get a running distribution you must implement getRange() and computeCDF(), all the other methods have a generic implementation, but these generic algorithms may be slow or inaccurate in difficult situations, so the more methods you provide the most efficient your distribution is.
The online documentation (http://secure-web.cisco.com/1_QlR51NTnci244KawP-NQBpuYV0mfhSTp4JXBwTpGJAqWfkDBxUY1JGCpr_XGFfQIZiEjVXqGphin3yoL6fV0Ro5q1JtmYip6xNs_iSbp14Iqwp3nL4GrrNfOhdGLYYpkwqSHmNCB0HmQ9TV4e4AAXGci3esSntXt8UjGeLuSqsWVIAgSfMQHC-yrAm6_JuU7HLDfmDuKjh0tAEtKH4Exm1PXrpvRwVwODuNBTOUej_7Q49C1pP-1sswRmGgOGm3NLSy3q3ZTJfNSogVMSuyZ4wcTzp0YH2CW-VW6edm4x21iG0omiNXB3pYDqQmrnNqz_uUj9AvsTr72Dh6iW8jGw/http%3A%2F%2Fopenturns.github.io%2Fopenturns%2Fmaster%2Fuser_manual%2F_generated%2Fopenturns.PythonDistribution.html%3Fhighlight%3Dpythondistribution) is very poor and will be updated. You can have a look at https://github.com/openturns/openturns/blob/master/python/test/t_Distribution_python.py for an example of a custom Python distribution.
You could also have used the SciPyDistribution wrapper of scipy distributions (see http://secure-web.cisco.com/1B_pKGvPkYtC86WGBDSnrF--4exrvoUqVoGHUYGIgi2rsg8OcFilo_WAXcF13h2kI6KqhmosBqsMMIjA_MZ9-V-t5GJ0wnk2hxZCqcogsNY4Q-1P7-gw4Jaj5Q1Y4l_n0WzDuI9V0YpSEvzmwpq65VeVL9VfSWO3Ec7QcmzW3i2NaK5oYsp4X5rrLp03MjlZiJscbJaKKstp9LrMGohAul3vBWC30HVxArmix5guZPggdaAavQcloR4ZVn0HF0ZGPR1LI6wJzM4LA4ZLfRh_EnO3yzEybS_kpyPYtyPUHVZY4xpUX8ojw-cDz54hEAWAx3a4pqsF5QN57pa9gdsppCw/http%3A%2F%2Fopenturns.github.io%2Fopenturns%2Fmaster%2Fuser_manual%2F_generated%2Fopenturns.SciPyDistribution.html and https://github.com/openturns/openturns/blob/master/python/test/t_Distribution_scipy.py):
import openturns as ot
import scipy.stats as st
lamb = 1.5
xi = 1.1
delta = 2.0
gamma = 1.0
distJU = ot.Distribution(ot.ScyPyDistribution(st.johnsonsu(gamma, delta, loc=xi, scale=lamb)))
but unfortunately this wrapper has a bug for unbounded distributions, resulting in a wrong range and a boggus computeQuantile() method.
Thanks for the question, it raised a lot of problems in OT!
Best regards,
Régis
Le vendredi 10 novembre 2017 à 18:42:19 UTC+1, Phil Fernandes <[email protected]> a écrit :
Hello,
I am attempting to use a custom continuous probability distribution for a probability of failure calculation via FORM, but when I try to create a ComposedDistribution with my custom dist, the program fails with
NotImplementedError: Wrong number or type of arguments for overloaded function 'new_ComposedDistribution'.
The custom dist is defined as
class JohnsonSU(ot.PythonDistribution):
def __init__(self, gamma=1, xi=0, delta=0.5, lam=1):
super(JohnsonSU, self).__init__(1)
if np.any(delta <= 0):
raise ValueError('Delta must be >0.')
if np.any(lam <= 0):
raise ValueError('Lambda must be >0.')
self.gamma = gamma # shape 1
self.xi = xi # location
self.delta = delta # shape 2, >0
self.lam = lam # scale, >0
self.scipy_dist = st.johnsonsu(self.gamma, self.delta, loc=self.xi, scale=self.lam)
def computeCDF(self, x):
return self.scipy_dist.cdf(x)
def computePDF(self, x):
return self.scipy_dist.pdf(x)
def computeQuantile(self, p):
return self.scipy_dist.ppf(p)
def getMean(self):
return self.scipy_dist.mean()
def getStandardDeviation(self):
return self.scipy_dist.std()
Are there additional functions that must be defined in order for the PythonDistribution to be compatible with existing OpenTurns Distributions? Or, is there a straightforward way that I could add arbitrary distributions to the OpenTurns source code?
Many thanks!
Phil Fernandes P.Eng, MASc
Engineer, Reliability Assessment
-
ENBRIDGE PIPELINES INC.
TEL: 780-420-8210 | FAX: 780-420-5234
7045 Enbridge Centre, 10175 101 Street NW, Edmonton, AB, T5J 0H3
www.enbridge.com
Integrity. Safety. Respect.
_______________________________________________
OpenTURNS users mailing list
[email protected]
http://openturns.org/mailman/listinfo/users
_______________________________________________
OpenTURNS users mailing list
[email protected]
http://openturns.org/mailman/listinfo/users
For example
a = ot.Normal()
x = np.linspace(0,1,5)[:,None]
a.computePDF(x)
returns
class=Sample name=Unnamed implementation=class=SampleImplementation name=Unnamed size=5 dimension=1 data="">
However
b=UniformNdPy()
b.computePDF(x)
returns 1.0.
Thank you.
-----Original Message-----
From: regis lebrun [mailto:[email protected]]
Sent: Friday, November 10, 2017 12:12 PM
To: [email protected]; Phil Fernandes
Subject: [External] Re: [ot-users] Custom distribution in FORM
Hi,
You can easily implement this distribution using OpenTURNS unique feature regarding distribution algebra (see https://en.wikipedia.org/wiki/Johnson%27s_SU-distribution):
import openturns as ot
lamb = 1.5
xi = 1.1
delta = 2.0
gamma = 1.0
distJU = ((ot.Normal() - gamma) / delta).sinh() * lamb + xi print("distJU=", distJU)
ot.Show(distJU.drawPDF())
You will get:
distJU= RandomMixture(1.1 + 1.5 * CompositeDistribution=f(RandomMixture(Normal(mu = -0.5, sigma = 0.5))) with f=[x]->[sinh(x)])
and a graph similar to the one given on the wikipedia page.
Your script contains some bugs (computeQuantile, getMean, getStandardDeviation should return the result as a float sequence of size 1) and a missing method, namely getRange(). To get a running distribution you must implement getRange() and computeCDF(), all the other methods have a generic implementation, but these generic algorithms may be slow or inaccurate in difficult situations, so the more methods you provide the most efficient your distribution is.
The online documentation (http://secure-web.cisco.com/1_QlR51NTnci244KawP-NQBpuYV0mfhSTp4JXBwTpGJAqWfkDBxUY1JGCpr_XGFfQIZiEjVXqGphin3yoL6fV0Ro5q1JtmYip6xNs_iSbp14Iqwp3nL4GrrNfOhdGLYYpkwqSHmNCB0HmQ9TV4e4AAXGci3esSntXt8UjGeLuSqsWVIAgSfMQHC-yrAm6_JuU7HLDfmDuKjh0tAEtKH4Exm1PXrpvRwVwODuNBTOUej_7Q49C1pP-1sswRmGgOGm3NLSy3q3ZTJfNSogVMSuyZ4wcTzp0YH2CW-VW6edm4x21iG0omiNXB3pYDqQmrnNqz_uUj9AvsTr72Dh6iW8jGw/http%3A%2F%2Fopenturns.github.io%2Fopenturns%2Fmaster%2Fuser_manual%2F_generated%2Fopenturns.PythonDistribution.html%3Fhighlight%3Dpythondistribution) is very poor and will be updated. You can have a look at https://github.com/openturns/openturns/blob/master/python/test/t_Distribution_python.py for an example of a custom Python distribution.
You could also have used the SciPyDistribution wrapper of scipy distributions (see http://secure-web.cisco.com/1B_pKGvPkYtC86WGBDSnrF--4exrvoUqVoGHUYGIgi2rsg8OcFilo_WAXcF13h2kI6KqhmosBqsMMIjA_MZ9-V-t5GJ0wnk2hxZCqcogsNY4Q-1P7-gw4Jaj5Q1Y4l_n0WzDuI9V0YpSEvzmwpq65VeVL9VfSWO3Ec7QcmzW3i2NaK5oYsp4X5rrLp03MjlZiJscbJaKKstp9LrMGohAul3vBWC30HVxArmix5guZPggdaAavQcloR4ZVn0HF0ZGPR1LI6wJzM4LA4ZLfRh_EnO3yzEybS_kpyPYtyPUHVZY4xpUX8ojw-cDz54hEAWAx3a4pqsF5QN57pa9gdsppCw/http%3A%2F%2Fopenturns.github.io%2Fopenturns%2Fmaster%2Fuser_manual%2F_generated%2Fopenturns.SciPyDistribution.html and https://github.com/openturns/openturns/blob/master/python/test/t_Distribution_scipy.py):
import openturns as ot
import scipy.stats as st
lamb = 1.5
xi = 1.1
delta = 2.0
gamma = 1.0
distJU = ot.Distribution(ot.ScyPyDistribution(st.johnsonsu(gamma, delta, loc=xi, scale=lamb)))
but unfortunately this wrapper has a bug for unbounded distributions, resulting in a wrong range and a boggus computeQuantile() method.
Thanks for the question, it raised a lot of problems in OT!
Best regards,
Régis
Le vendredi 10 novembre 2017 à 18:42:19 UTC+1, Phil Fernandes <[email protected]> a écrit :
Hello,
I am attempting to use a custom continuous probability distribution for a probability of failure calculation via FORM, but when I try to create a ComposedDistribution with my custom dist, the program fails with
NotImplementedError: Wrong number or type of arguments for overloaded function 'new_ComposedDistribution'.
The custom dist is defined as
class JohnsonSU(ot.PythonDistribution):
def __init__(self, gamma=1, xi=0, delta=0.5, lam=1):
super(JohnsonSU, self).__init__(1)
if np.any(delta <= 0):
raise ValueError('Delta must be >0.')
if np.any(lam <= 0):
raise ValueError('Lambda must be >0.')
self.gamma = gamma # shape 1
self.xi = xi # location
self.delta = delta # shape 2, >0
self.lam = lam # scale, >0
self.scipy_dist = st.johnsonsu(self.gamma, self.delta, loc=self.xi, scale=self.lam)
def computeCDF(self, x):
return self.scipy_dist.cdf(x)
def computePDF(self, x):
return self.scipy_dist.pdf(x)
def computeQuantile(self, p):
return self.scipy_dist.ppf(p)
def getMean(self):
return self.scipy_dist.mean()
def getStandardDeviation(self):
return self.scipy_dist.std()
Are there additional functions that must be defined in order for the PythonDistribution to be compatible with existing OpenTurns Distributions? Or, is there a straightforward way that I could add arbitrary distributions to the OpenTurns source code?
Many thanks!
Phil Fernandes P.Eng, MASc
Engineer, Reliability Assessment
-
ENBRIDGE PIPELINES INC.
TEL: 780-420-8210 | FAX: 780-420-5234
7045 Enbridge Centre, 10175 101 Street NW, Edmonton, AB, T5J 0H3
www.enbridge.com
Integrity. Safety. Respect.
_______________________________________________
OpenTURNS users mailing list
[email protected]
http://openturns.org/mailman/listinfo/users
_______________________________________________
OpenTURNS users mailing list
[email protected]
http://openturns.org/mailman/listinfo/users
#!/usr/bin/env python from __future__ import print_function import openturns as ot import scipy as sp import scipy.stats as st
class SciPyDistribution(ot.PythonDistribution):
"""
Allow to override Distribution from a scipy distribution.
Parameters
----------
dist : a scipy.stats distribution
The distribution to wrap
Examples
--------
>>> import openturns as ot
>>> # import scipy.stats as st
>>> # scipy_dist = st.johnsonsu(2.55, 2.25)
>>> # distribution = ot.Distribution(ot.SciPyDistribution(scipy_dist))
>>> # distribution.getRealization()
"""
def __init__(self, dist):
super(SciPyDistribution, self).__init__(1)
if dist.__class__.__name__ != 'rv_frozen':
raise TypeError('Argument is not a scipy distribution')
self._dist = dist
# compute range
cdf_epsilon = ot.ResourceMap.GetAsScalar('Distribution-DefaultCDFEpsilon')
lb = dist.ppf(0.0)
ub = dist.ppf(1.0)
flb = lb != float('-inf')
fub = ub != float('+inf')
if not flb:
lb = dist.ppf(cdf_epsilon)
if not fub:
ub = dist.ppf(1.0 - cdf_epsilon)
self.__range = Interval([lb], [ub])
self.__range.setFiniteLowerBound([int(flb)])
self.__range.setFiniteUpperBound([int(fub)])
def getRange(self):
return self.__range
def getRealization(self):
rvs = self._dist.rvs()
return [rvs]
def getSample(self, size):
rvs = self._dist.rvs(size)
return rvs.reshape(size, 1)
def computePDF(self, X):
pdf = self._dist.pdf(X[0])
return pdf
def computeCDF(self, X):
cdf = self._dist.cdf(X[0])
return cdf
def getMean(self):
mean = float(self._dist.stats('m'))
return [mean]
def getStandardDeviation(self):
var = float(self._dist.stats('v'))
std = var ** 0.5
return [std]
def getSkewness(self):
skewness = float(self._dist.stats('s'))
return [skewness]
def getKurtosis(self):
kurtosis = float(self._dist.stats('k')) + 3.0
return [kurtosis]
def getMoment(self, n):
moment = self._dist.moment(n)
return [moment]
def computeScalarQuantile(self, p):
q = self._dist.ppf(p)
return q
def computeQuantile(self, p):
q = self._dist.ppf(p)
return [q]
scipy_dist = st.johnsonsu(2.55, 2.25)
sp.random.seed(42)
# create an openturns distribution
py_dist = ot.SciPyDistribution(scipy_dist)
distribution = ot.Distribution(py_dist)
print('distribution=', distribution)
print('realization=', distribution.getRealization())
sample = distribution.getSample(10000)
print('sample=', sample[0:5])
point = [0.6]
print('pdf= %.6g' % distribution.computePDF(point))
cdf = distribution.computeCDF(point)
print('cdf= %.6g' % cdf)
print('quantile=', distribution.computeQuantile(cdf))
print('mean=', distribution.getMean())
print('mean(sampling)=', sample.computeMean())
print('std=', distribution.getStandardDeviation())
print('std(sampling)=', sample.computeStandardDeviationPerComponent())
print('skewness=', distribution.getSkewness())
print('skewness(sampling)=', sample.computeSkewness())
print('kurtosis=', distribution.getKurtosis())
print('kurtosis(sampling)=', sample.computeKurtosis())
print('range=', distribution.getRange())
_______________________________________________ OpenTURNS users mailing list [email protected] http://openturns.org/mailman/listinfo/users
