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 "FDistribution.h"
      11             : #include "Beta.h"
      12             : 
      13             : registerMooseObject("StochasticToolsApp", FDistribution);
      14             : 
      15             : InputParameters
      16          32 : FDistribution::validParams()
      17             : {
      18          32 :   InputParameters params = Distribution::validParams();
      19          32 :   params.addClassDescription("F-distribution or Fisher-Snedecor distribution");
      20          64 :   params.addRequiredRangeCheckedParam<unsigned int>("df1", "df1 > 0", "Degrees of freedom 1.");
      21          64 :   params.addRequiredRangeCheckedParam<unsigned int>("df2", "df2 > 0", "Degrees of freedom 2.");
      22          32 :   return params;
      23           0 : }
      24             : 
      25          16 : FDistribution::FDistribution(const InputParameters & parameters)
      26             :   : Distribution(parameters),
      27          16 :     _df1(getParam<unsigned int>("df1")),
      28          48 :     _df2(getParam<unsigned int>("df2"))
      29             : {
      30          16 : }
      31             : 
      32             : Real
      33          32 : FDistribution::pdf(const Real & x, const unsigned int & df1, const unsigned int & df2)
      34             : {
      35             :   // Handy definitions
      36          32 :   Real d1 = static_cast<Real>(df1);
      37          32 :   Real d2 = static_cast<Real>(df2);
      38          32 :   Real d1x = d1 * x;
      39          32 :   Real d2pd1x = d2 + d1x;
      40             : 
      41             :   Real a;
      42             :   Real b;
      43             :   Real y;
      44             :   Real z;
      45          32 :   if (d1 * x > d2)
      46             :   {
      47          16 :     a = d2 / 2.0;
      48          16 :     b = d1 / 2.0;
      49          16 :     y = (d2 * d1) / (d2pd1x * d2pd1x);
      50          16 :     z = d2 / d2pd1x;
      51             :   }
      52             :   else
      53             :   {
      54          16 :     a = d1 / 2.0;
      55          16 :     b = d2 / 2.0;
      56          16 :     y = (d2pd1x * d1 - d1x * d1) / (d2pd1x * d2pd1x);
      57          16 :     z = d1x / d2pd1x;
      58             :   }
      59             : 
      60          32 :   return y * Beta::pdf(z, a, b);
      61             : }
      62             : 
      63             : Real
      64          16 : FDistribution::cdf(const Real & x, const unsigned int & df1, const unsigned int & df2)
      65             : {
      66             :   // Handy definitions
      67          16 :   Real d1 = static_cast<Real>(df1);
      68          16 :   Real d2 = static_cast<Real>(df2);
      69             : 
      70          16 :   return Beta::incompleteBeta(d1 / 2.0, d2 / 2.0, d1 * x / (d1 * x + d2));
      71             : }
      72             : 
      73             : Real
      74          16 : FDistribution::quantile(const Real & p, const unsigned int & df1, const unsigned int & df2)
      75             : {
      76             :   // Handy definitions
      77          16 :   Real d1 = static_cast<Real>(df1);
      78          16 :   Real d2 = static_cast<Real>(df2);
      79             : 
      80          16 :   Real z = Beta::incompleteBetaInv(d1 / 2.0, d2 / 2.0, p);
      81          16 :   return d2 * z / d1 / (1.0 - z);
      82             : }
      83             : 
      84             : Real
      85          32 : FDistribution::pdf(const Real & x) const
      86             : {
      87          32 :   return pdf(x, _df1, _df2);
      88             : }
      89             : 
      90             : Real
      91          16 : FDistribution::cdf(const Real & x) const
      92             : {
      93          16 :   return cdf(x, _df1, _df2);
      94             : }
      95             : 
      96             : Real
      97          16 : FDistribution::quantile(const Real & p) const
      98             : {
      99          16 :   return quantile(p, _df1, _df2);
     100             : }