ENH: distributionModels: replace input variable names with de facto conventions

Previous naming conventions for input variables in the
normal, multiNormal, RosinRammler, and massRosinRammler
distributions were heuristic and did not reflect
the de facto conventions being used in statistics.
This commit is contained in:
Kutalmis Bercin
2021-08-24 10:39:35 +01:00
parent 89289b0716
commit 07cd88abd0
8 changed files with 115 additions and 75 deletions

View File

@ -49,7 +49,7 @@ Foam::distributionModels::RosinRammler::RosinRammler
) )
: :
distributionModel(typeName, dict, rndGen), distributionModel(typeName, dict, rndGen),
d_(distributionModelDict_.get<scalar>("d")), lambda_(distributionModelDict_.getCompat<scalar>("lambda", {{"d", 2112}})),
n_(distributionModelDict_.get<scalar>("n")) n_(distributionModelDict_.get<scalar>("n"))
{ {
const word parcelBasisType = const word parcelBasisType =
@ -63,11 +63,11 @@ Foam::distributionModels::RosinRammler::RosinRammler
<< endl; << endl;
} }
if (d_ < VSMALL || n_ < VSMALL) if (lambda_ < VSMALL || n_ < VSMALL)
{ {
FatalErrorInFunction FatalErrorInFunction
<< "Scale/Shape parameter cannot be equal to or less than zero:" << "Scale/Shape parameter cannot be equal to or less than zero:"
<< " d = " << d_ << " lambda = " << lambda_
<< " n = " << n_ << " n = " << n_
<< exit(FatalError); << exit(FatalError);
} }
@ -79,7 +79,7 @@ Foam::distributionModels::RosinRammler::RosinRammler
Foam::distributionModels::RosinRammler::RosinRammler(const RosinRammler& p) Foam::distributionModels::RosinRammler::RosinRammler(const RosinRammler& p)
: :
distributionModel(p), distributionModel(p),
d_(p.d_), lambda_(p.lambda_),
n_(p.n_) n_(p.n_)
{} {}
@ -88,16 +88,16 @@ Foam::distributionModels::RosinRammler::RosinRammler(const RosinRammler& p)
Foam::scalar Foam::distributionModels::RosinRammler::sample() const Foam::scalar Foam::distributionModels::RosinRammler::sample() const
{ {
const scalar minValueByDPowN = pow(minValue_/d_, n_); const scalar minValueByDPowN = pow(minValue_/lambda_, n_);
const scalar K = 1 - exp(- pow(maxValue_/d_, n_) + minValueByDPowN); const scalar K = 1 - exp(- pow(maxValue_/lambda_, n_) + minValueByDPowN);
const scalar y = rndGen_.sample01<scalar>(); const scalar y = rndGen_.sample01<scalar>();
return d_*pow(minValueByDPowN - log(1 - K*y), 1/n_); return lambda_*pow(minValueByDPowN - log(1 - K*y), 1/n_);
} }
Foam::scalar Foam::distributionModels::RosinRammler::meanValue() const Foam::scalar Foam::distributionModels::RosinRammler::meanValue() const
{ {
return d_; return lambda_;
} }

View File

@ -63,13 +63,11 @@ class RosinRammler
{ {
// Private Data // Private Data
// Model coefficients //- Scale parameter
scalar lambda_;
//- Scale parameter //- Shape parameter
scalar d_; scalar n_;
//- Shape parameter
scalar n_;
public: public:

View File

@ -50,14 +50,14 @@ Foam::distributionModels::massRosinRammler::massRosinRammler
) )
: :
distributionModel(typeName, dict, rndGen), distributionModel(typeName, dict, rndGen),
d_(distributionModelDict_.get<scalar>("d")), lambda_(distributionModelDict_.getCompat<scalar>("lambda", {{"d", 2112}})),
n_(distributionModelDict_.get<scalar>("n")) n_(distributionModelDict_.get<scalar>("n"))
{ {
if (d_ < VSMALL || n_ < VSMALL) if (lambda_ < VSMALL || n_ < VSMALL)
{ {
FatalErrorInFunction FatalErrorInFunction
<< "Scale/Shape parameter cannot be equal to or less than zero:" << "Scale/Shape parameter cannot be equal to or less than zero:"
<< " d = " << d_ << " lambda = " << lambda_
<< " n = " << n_ << " n = " << n_
<< exit(FatalError); << exit(FatalError);
} }
@ -72,7 +72,7 @@ Foam::distributionModels::massRosinRammler::massRosinRammler
) )
: :
distributionModel(p), distributionModel(p),
d_(p.d_), lambda_(p.lambda_),
n_(p.n_) n_(p.n_)
{} {}
@ -89,7 +89,7 @@ Foam::scalar Foam::distributionModels::massRosinRammler::sample() const
const scalar a = 3/n_ + 1; const scalar a = 3/n_ + 1;
const scalar P = rndGen_.sample01<scalar>(); const scalar P = rndGen_.sample01<scalar>();
const scalar x = Math::invIncGamma(a, P); const scalar x = Math::invIncGamma(a, P);
d = d_*pow(x, 1/n_); d = lambda_*pow(x, 1/n_);
} while (d < minValue_ || d > maxValue_); } while (d < minValue_ || d > maxValue_);
return d; return d;
@ -98,7 +98,7 @@ Foam::scalar Foam::distributionModels::massRosinRammler::sample() const
Foam::scalar Foam::distributionModels::massRosinRammler::meanValue() const Foam::scalar Foam::distributionModels::massRosinRammler::meanValue() const
{ {
return d_; return lambda_;
} }

View File

@ -74,7 +74,7 @@ class massRosinRammler
// Private Data // Private Data
//- Scale parameter //- Scale parameter
scalar d_; scalar lambda_;
//- Shape parameter //- Shape parameter
scalar n_; scalar n_;

View File

@ -49,35 +49,67 @@ Foam::distributionModels::multiNormal::multiNormal
) )
: :
distributionModel(typeName, dict, rndGen), distributionModel(typeName, dict, rndGen),
range_(maxValue_ - minValue_), mu_
expectation_(distributionModelDict_.lookup("expectation")), (
variance_(distributionModelDict_.lookup("variance")), distributionModelDict_.lookupCompat
strength_(distributionModelDict_.lookup("strength")) (
"mu",
{{"expectation", 2112}}
)
),
sigma_
(
distributionModelDict_.lookupCompat
(
"sigma",
{{"variance", 2112}}
)
),
weight_
(
distributionModelDict_.lookupCompat
(
"weight",
{{"strength", 2112}}
)
)
{ {
check(); check();
scalar sMax = 0; scalar sum = 0;
label n = strength_.size(); for (label i = 0; i < weight_.size(); ++i)
for (label i=0; i<n; i++)
{ {
scalar x = expectation_[i]; if (i > 0 && weight_[i] < weight_[i-1])
scalar s = strength_[i];
for (label j=0; j<n; j++)
{ {
if (i!=j) FatalErrorInFunction
{ << type() << "distribution: "
scalar x2 = (x - expectation_[j])/variance_[j]; << "Weights must be specified in a monotonic order." << nl
scalar y = exp(-0.5*x2*x2); << "Please see the row i = " << i << nl
s += strength_[j]*y; << "weight[i-1] = " << weight_[i-1] << nl
} << "weight[i] = " << weight_[i]
<< exit(FatalError);
} }
sMax = max(sMax, s); sum += weight_[i];
} }
for (label i=0; i<n; i++) if (sum < VSMALL)
{ {
strength_[i] /= sMax; FatalErrorInFunction
<< type() << "distribution: "
<< "The sum of weights cannot be zero." << nl
<< "weight = " << weight_
<< exit(FatalError);
}
for (label i = 1; i < weight_.size(); ++i)
{
weight_[i] += weight_[i-1];
}
for (auto& w : weight_)
{
w /= sum;
} }
} }
@ -85,10 +117,9 @@ Foam::distributionModels::multiNormal::multiNormal
Foam::distributionModels::multiNormal::multiNormal(const multiNormal& p) Foam::distributionModels::multiNormal::multiNormal(const multiNormal& p)
: :
distributionModel(p), distributionModel(p),
range_(p.range_), mu_(p.mu_),
expectation_(p.expectation_), sigma_(p.sigma_),
variance_(p.variance_), weight_(p.weight_)
strength_(p.strength_)
{} {}
@ -98,7 +129,7 @@ Foam::scalar Foam::distributionModels::multiNormal::sample() const
{ {
scalar y = 0; scalar y = 0;
scalar x = 0; scalar x = 0;
label n = expectation_.size(); label n = mu_.size();
bool success = false; bool success = false;
while (!success) while (!success)
@ -109,9 +140,9 @@ Foam::scalar Foam::distributionModels::multiNormal::sample() const
for (label i=0; i<n; i++) for (label i=0; i<n; i++)
{ {
scalar nu = expectation_[i]; scalar nu = mu_[i];
scalar sigma = variance_[i]; scalar sigma = sigma_[i];
scalar s = strength_[i]; scalar s = weight_[i];
scalar v = (x - nu)/sigma; scalar v = (x - nu)/sigma;
p += s*exp(-0.5*v*v); p += s*exp(-0.5*v*v);
} }
@ -128,10 +159,10 @@ Foam::scalar Foam::distributionModels::multiNormal::sample() const
Foam::scalar Foam::distributionModels::multiNormal::meanValue() const Foam::scalar Foam::distributionModels::multiNormal::meanValue() const
{ {
scalar mean = 0.0; scalar mean = 0;
forAll(strength_, i) forAll(weight_, i)
{ {
mean += strength_[i]*expectation_[i]; mean += weight_[i]*mu_[i];
} }
return mean; return mean;

View File

@ -62,15 +62,15 @@ class multiNormal
{ {
// Private Data // Private Data
//- Distribution range //- List of means of the parent general normal distributions
scalar range_; List<scalar> mu_;
//- List of standard deviations of
//- the parent general normal distributions
List<scalar> sigma_;
// Model coefficients //- List of weights of a given distribution in the mixture
List<scalar> weight_;
List<scalar> expectation_;
List<scalar> variance_;
List<scalar> strength_;
public: public:

View File

@ -50,15 +50,28 @@ Foam::distributionModels::normal::normal
) )
: :
distributionModel(typeName, dict, rndGen), distributionModel(typeName, dict, rndGen),
expectation_(distributionModelDict_.get<scalar>("expectation")), mu_
variance_(distributionModelDict_.get<scalar>("variance")), (
a_(0.147) distributionModelDict_.getCompat<scalar>
(
"mu",
{{"expectation", 2112}}
)
),
sigma_
(
distributionModelDict_.getCompat<scalar>
(
"sigma",
{{"variance", 2112}}
)
)
{ {
if (mag(variance_) == 0) if (mag(sigma_) == 0)
{ {
FatalErrorInFunction FatalErrorInFunction
<< "Standard deviation cannot be zero." << nl << "Standard deviation cannot be zero." << nl
<< " variance = " << variance_ << nl << " sigma = " << sigma_ << nl
<< exit(FatalError); << exit(FatalError);
} }
@ -69,9 +82,8 @@ Foam::distributionModels::normal::normal
Foam::distributionModels::normal::normal(const normal& p) Foam::distributionModels::normal::normal(const normal& p)
: :
distributionModel(p), distributionModel(p),
expectation_(p.expectation_), mu_(p.mu_),
variance_(p.variance_), sigma_(p.sigma_)
a_(p.a_)
{} {}
@ -80,11 +92,11 @@ Foam::distributionModels::normal::normal(const normal& p)
Foam::scalar Foam::distributionModels::normal::sample() const Foam::scalar Foam::distributionModels::normal::sample() const
{ {
scalar a = erf((minValue_ - expectation_)/variance_); scalar a = erf((minValue_ - mu_)/sigma_);
scalar b = erf((maxValue_ - expectation_)/variance_); scalar b = erf((maxValue_ - mu_)/sigma_);
scalar y = rndGen_.sample01<scalar>(); scalar y = rndGen_.sample01<scalar>();
scalar x = Math::erfInv(y*(b - a) + a)*variance_ + expectation_; scalar x = Math::erfInv(y*(b - a) + a)*sigma_ + mu_;
// Note: numerical approximation of the inverse function yields slight // Note: numerical approximation of the inverse function yields slight
// inaccuracies // inaccuracies
@ -97,7 +109,7 @@ Foam::scalar Foam::distributionModels::normal::sample() const
Foam::scalar Foam::distributionModels::normal::meanValue() const Foam::scalar Foam::distributionModels::normal::meanValue() const
{ {
return expectation_; return mu_;
} }

View File

@ -63,12 +63,11 @@ class normal
{ {
// Private Data // Private Data
// Model coefficients //- Mean of the parent general normal distribution
scalar mu_;
scalar expectation_; //- Standard deviation of the parent general normal distribution
scalar variance_; scalar sigma_;
scalar a_;
public: public: