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        1944 : Normal::validParams()
      24             : {
      25        1944 :   InputParameters params = Distribution::validParams();
      26        1944 :   params.addClassDescription("Normal distribution");
      27        3888 :   params.addRequiredParam<Real>("mean", "Mean (or expectation) of the distribution.");
      28        3888 :   params.addRequiredRangeCheckedParam<Real>(
      29             :       "standard_deviation", "standard_deviation > 0", "Standard deviation of the distribution ");
      30        1944 :   return params;
      31           0 : }
      32             : 
      33         971 : Normal::Normal(const InputParameters & parameters)
      34             :   : Distribution(parameters),
      35         971 :     _mean(getParam<Real>("mean")),
      36        2913 :     _standard_deviation(getParam<Real>("standard_deviation"))
      37             : {
      38         971 : }
      39             : 
      40             : Real
      41       40057 : Normal::pdf(const Real & x, const Real & mean, const Real & std_dev)
      42             : {
      43       40057 :   return 1.0 / (std_dev * std::sqrt(2.0 * M_PI)) *
      44       40057 :          std::exp(-0.5 * Utility::pow<2>((x - mean) / std_dev));
      45             : }
      46             : 
      47             : Real
      48       35247 : Normal::cdf(const Real & x, const Real & mean, const Real & std_dev)
      49             : {
      50       35247 :   return 0.5 * (1.0 + std::erf((x - mean) / (std_dev * std::sqrt(2.0))));
      51             : }
      52             : 
      53             : Real
      54      478087 : Normal::quantile(const Real & p, const Real & mean, const Real & std_dev)
      55             : {
      56      478087 :   const Real x = (p < 0.5 ? p : 1.0 - p);
      57      478087 :   const Real y = std::sqrt(-2.0 * std::log(x));
      58      478087 :   const Real y2 = y * y;
      59      478087 :   const Real y3 = y2 * y;
      60      478087 :   const Real y4 = y3 * y;
      61      478087 :   const Real sgn = (p - 0.5 < 0.0 ? -1.0 : 1.0);
      62      478087 :   const Real Zp = sgn * (y + (_a[0] + _a[1] * y + _a[2] * y2 + _a[3] * y3 + _a[4] * y4) /
      63      478087 :                                  (_b[0] + _b[1] * y + _b[2] * y2 + _b[3] * y3 + _b[4] * y4));
      64      478087 :   return Zp * std_dev + mean;
      65             : }
      66             : 
      67             : Real
      68        6087 : Normal::pdf(const Real & x) const
      69             : {
      70        6087 :   return pdf(x, _mean, _standard_deviation);
      71             : }
      72             : 
      73             : Real
      74        1461 : Normal::cdf(const Real & x) const
      75             : {
      76        1461 :   return cdf(x, _mean, _standard_deviation);
      77             : }
      78             : 
      79             : Real
      80      467769 : Normal::quantile(const Real & p) const
      81             : {
      82      467769 :   return quantile(p, _mean, _standard_deviation);
      83             : }