Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "JohnsonSB.h"
      11             : #include "math.h"
      12             : #include "libmesh/utility.h"
      13             : 
      14             : registerMooseObject("StochasticToolsApp", JohnsonSB);
      15             : 
      16             : InputParameters
      17          32 : JohnsonSB::validParams()
      18             : {
      19          32 :   InputParameters params = Normal::validParams();
      20          32 :   params.addClassDescription("Johnson Special Bounded (SB) distribution.");
      21             : 
      22          32 :   params.set<Real>("mean") = 0.0;
      23          32 :   params.set<Real>("standard_deviation") = 1.0;
      24          32 :   params.suppressParameter<Real>("mean");
      25          32 :   params.suppressParameter<Real>("standard_deviation");
      26             : 
      27          64 :   params.addRequiredParam<Real>("a", "Lower location parameter");
      28          64 :   params.addRequiredParam<Real>("b", "Upper location parameter");
      29          64 :   params.addRequiredParam<Real>("alpha_1", "Shape parameter (sometimes called a)");
      30          64 :   params.addRequiredParam<Real>("alpha_2", "Shape parameter (sometimes called b)");
      31             : 
      32          32 :   return params;
      33           0 : }
      34             : 
      35          16 : JohnsonSB::JohnsonSB(const InputParameters & parameters)
      36             :   : Normal(parameters),
      37          16 :     _lower(getParam<Real>("a")),
      38          32 :     _upper(getParam<Real>("b")),
      39          32 :     _alpha_1(getParam<Real>("alpha_1")),
      40          48 :     _alpha_2(getParam<Real>("alpha_2"))
      41             : {
      42          16 : }
      43             : 
      44             : Real
      45          16 : JohnsonSB::pdf(
      46             :     const Real & x, const Real & a, const Real & b, const Real & alpha_1, const Real & alpha_2)
      47             : {
      48          16 :   if (x <= a)
      49             :     return 0.0;
      50          16 :   else if (x < b)
      51             :   {
      52          16 :     return (alpha_2 * (b - a)) / ((x - a) * (b - x) * std::sqrt(2.0 * M_PI)) *
      53          16 :            std::exp(-0.5 * Utility::pow<2>(alpha_1 + alpha_2 * std::log((x - a) / (b - x))));
      54             :   }
      55             :   else
      56             :     return 0.0;
      57             : }
      58             : 
      59             : Real
      60          16 : JohnsonSB::cdf(
      61             :     const Real & x, const Real & a, const Real & b, const Real & alpha_1, const Real & alpha_2)
      62             : {
      63          16 :   if (x <= a)
      64             :     return 0.0;
      65          16 :   else if (x < b)
      66             :   {
      67          16 :     return Normal::cdf(alpha_1 + alpha_2 * std::log((x - a) / (b - x)), 0.0, 1.0);
      68             :   }
      69             :   else
      70             :     return 0.0;
      71             : }
      72             : 
      73             : Real
      74          16 : JohnsonSB::quantile(
      75             :     const Real & p, const Real & a, const Real & b, const Real & alpha_1, const Real & alpha_2)
      76             : {
      77          16 :   const Real Z = Normal::quantile(p, 0.0, 1.0);
      78          16 :   return (a + b * std::exp((Z - alpha_1) / alpha_2)) / (1.0 + std::exp((Z - alpha_1) / alpha_2));
      79             : }
      80             : 
      81             : Real
      82          16 : JohnsonSB::pdf(const Real & x) const
      83             : {
      84          16 :   return pdf(x, _lower, _upper, _alpha_1, _alpha_2);
      85             : }
      86             : 
      87             : Real
      88          16 : JohnsonSB::cdf(const Real & x) const
      89             : {
      90          16 :   return cdf(x, _lower, _upper, _alpha_1, _alpha_2);
      91             : }
      92             : 
      93             : Real
      94          16 : JohnsonSB::quantile(const Real & p) const
      95             : {
      96          16 :   return quantile(p, _lower, _upper, _alpha_1, _alpha_2);
      97             : }