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 "Normal.h"
      11             : #include "math.h"
      12             : #include "libmesh/utility.h"
      13             : 
      14             : registerMooseObject("StochasticToolsApp", Normal);
      15             : 
      16             : const std::array<Real, 6> Normal::_a = {
      17             :     {-0.322232431088, -1.0, -0.342242088547, -0.0204231210245, -0.0000453642210148}};
      18             : 
      19             : const std::array<Real, 6> Normal::_b = {
      20             :     {0.099348462606, 0.588581570495, 0.531103462366, 0.10353775285, 0.0038560700634}};
      21             : 
      22             : InputParameters
      23        3780 : Normal::validParams()
      24             : {
      25        3780 :   InputParameters params = Distribution::validParams();
      26        3780 :   params.addClassDescription("Normal distribution");
      27        7560 :   params.addRequiredParam<Real>("mean", "Mean (or expectation) of the distribution.");
      28        7560 :   params.addRequiredRangeCheckedParam<Real>(
      29             :       "standard_deviation", "standard_deviation > 0", "Standard deviation of the distribution ");
      30        3780 :   return params;
      31           0 : }
      32             : 
      33        1888 : Normal::Normal(const InputParameters & parameters)
      34             :   : Distribution(parameters),
      35        1888 :     _mean(getParam<Real>("mean")),
      36        5664 :     _standard_deviation(getParam<Real>("standard_deviation"))
      37             : {
      38        1888 : }
      39             : 
      40             : Real
      41       30388 : Normal::pdf(const Real & x, const Real & mean, const Real & std_dev)
      42             : {
      43       30388 :   return 1.0 / (std_dev * std::sqrt(2.0 * M_PI)) *
      44       30388 :          std::exp(-0.5 * Utility::pow<2>((x - mean) / std_dev));
      45             : }
      46             : 
      47             : Real
      48       13554 : Normal::cdf(const Real & x, const Real & mean, const Real & std_dev)
      49             : {
      50       13554 :   return 0.5 * (1.0 + std::erf((x - mean) / (std_dev * std::sqrt(2.0))));
      51             : }
      52             : 
      53             : Real
      54       47426 : Normal::quantile(const Real & p, const Real & mean, const Real & std_dev)
      55             : {
      56       47426 :   const Real x = (p < 0.5 ? p : 1.0 - p);
      57       47426 :   const Real y = std::sqrt(-2.0 * std::log(x));
      58       47426 :   const Real y2 = y * y;
      59       47426 :   const Real y3 = y2 * y;
      60       47426 :   const Real y4 = y3 * y;
      61       47426 :   const Real sgn = (p - 0.5 < 0.0 ? -1.0 : 1.0);
      62       47426 :   const Real Zp = sgn * (y + (_a[0] + _a[1] * y + _a[2] * y2 + _a[3] * y3 + _a[4] * y4) /
      63       47426 :                                  (_b[0] + _b[1] * y + _b[2] * y2 + _b[3] * y3 + _b[4] * y4));
      64       47426 :   return Zp * std_dev + mean;
      65             : }
      66             : 
      67             : Real
      68       12176 : Normal::pdf(const Real & x) const
      69             : {
      70       12176 :   return pdf(x, _mean, _standard_deviation);
      71             : }
      72             : 
      73             : Real
      74        1984 : Normal::cdf(const Real & x) const
      75             : {
      76        1984 :   return cdf(x, _mean, _standard_deviation);
      77             : }
      78             : 
      79             : Real
      80       32816 : Normal::quantile(const Real & p) const
      81             : {
      82       32816 :   return quantile(p, _mean, _standard_deviation);
      83             : }