Browse Source

Merge branch 'master' into integral

pull/655/head
Christoph Ruegg 6 years ago
committed by GitHub
parent
commit
a4c511bc9d
No known key found for this signature in database GPG Key ID: 4AEE18F83AFDEB23
  1. 1
      .gitignore
  2. 88
      .paket/Paket.Restore.targets
  3. 3
      src/Numerics.Tests/DistributionTests/CommonDistributionTests.cs
  4. 426
      src/Numerics.Tests/DistributionTests/Continuous/BurrTests.cs
  5. 419
      src/Numerics.Tests/DistributionTests/Continuous/InverseGaussianTests.cs
  6. 446
      src/Numerics.Tests/DistributionTests/Continuous/TruncatedParetoTests.cs
  7. 36
      src/Numerics.Tests/OptimizationTests/BfgsBMinimizerTests.cs
  8. 429
      src/Numerics/Distributions/Burr.cs
  9. 449
      src/Numerics/Distributions/InverseGaussian.cs
  10. 466
      src/Numerics/Distributions/TruncatedPareto.cs
  11. 9
      src/Numerics/Optimization/BfgsBMinimizer.cs
  12. 9
      src/Numerics/Precision.Comparison.cs
  13. 11
      src/Numerics/Properties/Resources.Designer.cs
  14. 3
      src/Numerics/Properties/Resources.resx

1
.gitignore

@ -12,6 +12,7 @@ out/
.vs
.vscode
.idea
.ionide
# Test & Analysis Results
[Tt]est[Rr]esult*/

88
.paket/Paket.Restore.targets

@ -27,38 +27,7 @@
<PaketBootStrapperExePath Condition=" '$(PaketBootStrapperExePath)' == '' AND Exists('$(PaketRootPath)paket.bootstrapper.exe')">$(PaketRootPath)paket.bootstrapper.exe</PaketBootStrapperExePath>
<PaketBootStrapperExePath Condition=" '$(PaketBootStrapperExePath)' == '' ">$(PaketToolsPath)paket.bootstrapper.exe</PaketBootStrapperExePath>
<PaketBootStrapperExeDir Condition=" Exists('$(PaketBootStrapperExePath)') " >$([System.IO.Path]::GetDirectoryName("$(PaketBootStrapperExePath)"))\</PaketBootStrapperExeDir>
<!-- Paket -->
<!-- windows, root => tool => proj style => bootstrapper => global -->
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' == 'Windows_NT' AND Exists('$(PaketRootPath)paket.exe') ">$(PaketRootPath)paket.exe</PaketExePath>
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' == 'Windows_NT' AND Exists('$(PaketToolsPath)paket.exe') ">$(PaketToolsPath)paket.exe</PaketExePath>
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' == 'Windows_NT' AND '$(PaketBootstrapperStyle)' == 'proj' ">$(PaketToolsPath)paket.exe</PaketExePath>
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' == 'Windows_NT' AND Exists('$(PaketBootStrapperExeDir)') ">$(_PaketBootStrapperExeDir)paket.exe</PaketExePath>
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' == 'Windows_NT' ">paket.exe</PaketExePath>
<!-- no windows, try native paket as default, root => tool => proj style => mono paket => bootstrpper => global -->
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' != 'Windows_NT' AND Exists('$(PaketRootPath)paket') ">$(PaketRootPath)paket</PaketExePath>
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' != 'Windows_NT' AND Exists('$(PaketToolsPath)paket') ">$(PaketToolsPath)paket</PaketExePath>
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' != 'Windows_NT' AND '$(PaketBootstrapperStyle)' == 'proj' ">$(PaketToolsPath)paket</PaketExePath>
<!-- no windows, try mono paket -->
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' != 'Windows_NT' AND Exists('$(PaketRootPath)paket.exe') ">$(PaketRootPath)paket.exe</PaketExePath>
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' != 'Windows_NT' AND Exists('$(PaketToolsPath)paket.exe') ">$(PaketToolsPath)paket.exe</PaketExePath>
<!-- no windows, try bootstrapper -->
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' != 'Windows_NT' AND Exists('$(PaketBootStrapperExeDir)') ">$(PaketBootStrapperExeDir)paket.exe</PaketExePath>
<!-- no windows, try global native paket -->
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' != 'Windows_NT' ">paket</PaketExePath>
<!-- Paket command -->
<_PaketExeExtension>$([System.IO.Path]::GetExtension("$(PaketExePath)"))</_PaketExeExtension>
<PaketCommand Condition=" '$(PaketCommand)' == '' AND '$(_PaketExeExtension)' == '.dll' ">dotnet "$(PaketExePath)"</PaketCommand>
<PaketCommand Condition=" '$(PaketCommand)' == '' AND '$(OS)' != 'Windows_NT' AND '$(_PaketExeExtension)' == '.exe' ">$(MonoPath) --runtime=v4.0.30319 "$(PaketExePath)"</PaketCommand>
<PaketCommand Condition=" '$(PaketCommand)' == '' ">"$(PaketExePath)"</PaketCommand>
<PaketBootStrapperCommand Condition=" '$(OS)' == 'Windows_NT'">"$(PaketBootStrapperExePath)"</PaketBootStrapperCommand>
<PaketBootStrapperCommand Condition=" '$(OS)' != 'Windows_NT' ">$(MonoPath) --runtime=v4.0.30319 "$(PaketBootStrapperExePath)"</PaketBootStrapperCommand>
@ -74,6 +43,57 @@
<PaketIntermediateOutputPath Condition=" '$(PaketIntermediateOutputPath)' == '' ">$(BaseIntermediateOutputPath.TrimEnd('\').TrimEnd('\/'))</PaketIntermediateOutputPath>
</PropertyGroup>
<!-- Check if paket is available as local dotnet cli tool -->
<Target Name="SetPaketCommand" >
<!-- Call 'dotnet paket' and see if it returns without an error. Mute all the output. -->
<Exec Command="dotnet paket --version" IgnoreExitCode="true" StandardOutputImportance="low" StandardErrorImportance="low" >
<Output TaskParameter="ExitCode" PropertyName="LocalPaketToolExitCode" />
</Exec>
<!-- If local paket tool is found, use that -->
<PropertyGroup Condition=" '$(LocalPaketToolExitCode)' == '0' ">
<InternalPaketCommand>dotnet paket</InternalPaketCommand>
</PropertyGroup>
<!-- If not, then we go through our normal steps of setting the Paket command. -->
<PropertyGroup Condition=" '$(LocalPaketToolExitCode)' != '0' ">
<!-- windows, root => tool => proj style => bootstrapper => global -->
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' == 'Windows_NT' AND Exists('$(PaketRootPath)paket.exe') ">$(PaketRootPath)paket.exe</PaketExePath>
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' == 'Windows_NT' AND Exists('$(PaketToolsPath)paket.exe') ">$(PaketToolsPath)paket.exe</PaketExePath>
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' == 'Windows_NT' AND '$(PaketBootstrapperStyle)' == 'proj' ">$(PaketToolsPath)paket.exe</PaketExePath>
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' == 'Windows_NT' AND Exists('$(PaketBootStrapperExeDir)') ">$(_PaketBootStrapperExeDir)paket.exe</PaketExePath>
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' == 'Windows_NT' ">paket.exe</PaketExePath>
<!-- no windows, try native paket as default, root => tool => proj style => mono paket => bootstrpper => global -->
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' != 'Windows_NT' AND Exists('$(PaketRootPath)paket') ">$(PaketRootPath)paket</PaketExePath>
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' != 'Windows_NT' AND Exists('$(PaketToolsPath)paket') ">$(PaketToolsPath)paket</PaketExePath>
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' != 'Windows_NT' AND '$(PaketBootstrapperStyle)' == 'proj' ">$(PaketToolsPath)paket</PaketExePath>
<!-- no windows, try mono paket -->
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' != 'Windows_NT' AND Exists('$(PaketRootPath)paket.exe') ">$(PaketRootPath)paket.exe</PaketExePath>
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' != 'Windows_NT' AND Exists('$(PaketToolsPath)paket.exe') ">$(PaketToolsPath)paket.exe</PaketExePath>
<!-- no windows, try bootstrapper -->
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' != 'Windows_NT' AND Exists('$(PaketBootStrapperExeDir)') ">$(PaketBootStrapperExeDir)paket.exe</PaketExePath>
<!-- no windows, try global native paket -->
<PaketExePath Condition=" '$(PaketExePath)' == '' AND '$(OS)' != 'Windows_NT' ">paket</PaketExePath>
<_PaketExeExtension>$([System.IO.Path]::GetExtension("$(PaketExePath)"))</_PaketExeExtension>
<InternalPaketCommand Condition=" '$(InternalPaketCommand)' == '' AND '$(_PaketExeExtension)' == '.dll' ">dotnet "$(PaketExePath)"</InternalPaketCommand>
<InternalPaketCommand Condition=" '$(InternalPaketCommand)' == '' AND '$(OS)' != 'Windows_NT' AND '$(_PaketExeExtension)' == '.exe' ">$(MonoPath) --runtime=v4.0.30319 "$(PaketExePath)"</InternalPaketCommand>
<InternalPaketCommand Condition=" '$(InternalPaketCommand)' == '' ">"$(PaketExePath)"</InternalPaketCommand>
</PropertyGroup>
<!-- The way to get a property to be available outside the target is to use this task. -->
<CreateProperty Value="$(InternalPaketCommand)">
<Output TaskParameter="Value" PropertyName="PaketCommand"/>
</CreateProperty>
</Target>
<Target Name="PaketBootstrapping" Condition="Exists('$(PaketToolsPath)paket.bootstrapper.proj')">
<MSBuild Projects="$(PaketToolsPath)paket.bootstrapper.proj" Targets="Restore" />
</Target>
@ -81,7 +101,7 @@
<!-- Official workaround for https://docs.microsoft.com/en-us/visualstudio/msbuild/getfilehash-task?view=vs-2019 -->
<UsingTask TaskName="Microsoft.Build.Tasks.GetFileHash" AssemblyName="Microsoft.Build.Tasks.Core, Version=15.1.0.0, Culture=neutral, PublicKeyToken=b03f5f7f11d50a3a" Condition=" '$(MSBuildSupportsHashing)' == 'true' And '$(DetectedMSBuildVersion)' &lt; '16.0.360' " />
<UsingTask TaskName="Microsoft.Build.Tasks.VerifyFileHash" AssemblyName="Microsoft.Build.Tasks.Core, Version=15.1.0.0, Culture=neutral, PublicKeyToken=b03f5f7f11d50a3a" Condition=" '$(MSBuildSupportsHashing)' == 'true' And '$(DetectedMSBuildVersion)' &lt; '16.0.360' " />
<Target Name="PaketRestore" Condition="'$(PaketRestoreDisabled)' != 'True'" BeforeTargets="_GenerateDotnetCliToolReferenceSpecs;_GenerateProjectRestoreGraphPerFramework;_GenerateRestoreGraphWalkPerFramework;CollectPackageReferences" DependsOnTargets="PaketBootstrapping">
<Target Name="PaketRestore" Condition="'$(PaketRestoreDisabled)' != 'True'" BeforeTargets="_GenerateDotnetCliToolReferenceSpecs;_GenerateProjectRestoreGraphPerFramework;_GenerateRestoreGraphWalkPerFramework;CollectPackageReferences" DependsOnTargets="SetPaketCommand;PaketBootstrapping">
<!-- Step 1 Check if lockfile is properly restored (if the hash of the lockfile and the cache-file match) -->
<PropertyGroup>
@ -246,7 +266,7 @@
</PropertyGroup>
</Target>
<Target Name="PaketOverrideNuspec" AfterTargets="GenerateNuspec" Condition="('$(IsPackable)' == '' Or '$(IsPackable)' == 'true') And Exists('$(PaketIntermediateOutputPath)/$(MSBuildProjectFile).references')" >
<Target Name="PaketOverrideNuspec" DependsOnTargets="SetPaketCommand" AfterTargets="GenerateNuspec" Condition="('$(IsPackable)' == '' Or '$(IsPackable)' == 'true') And Exists('$(PaketIntermediateOutputPath)/$(MSBuildProjectFile).references')" >
<ItemGroup>
<_NuspecFilesNewLocation Include="$(PaketIntermediateOutputPath)\$(Configuration)\*.nuspec"/>
<MSBuildMajorVersion Include="$(DetectedMSBuildVersion.Replace(`-`, `.`).Split(`.`)[0])" />

3
src/Numerics.Tests/DistributionTests/CommonDistributionTests.cs

@ -67,6 +67,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
{
new Beta(1.0, 1.0),
new BetaScaled(1.0, 1.5, 0.5, 2.0),
new Burr(2.0, 3.0, 1.0),
new Cauchy(1.0, 1.0),
new Chi(3.0),
new ChiSquared(3.0),
@ -76,6 +77,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
new FisherSnedecor(0.3, 0.4),
new Gamma(1.0, 1.0),
new InverseGamma(1.0, 1.0),
new InverseGaussian(1.0, 3.0),
new Laplace(1.0, 0.5),
new LogNormal(1.0, 1.0),
new Normal(0.0, 1.0),
@ -84,6 +86,7 @@ namespace MathNet.Numerics.UnitTests.DistributionTests
new Stable(0.5, 1.0, 0.5, 1.0),
new StudentT(0.0, 1.0, 5.0),
new Triangular(0, 1, 0.7),
new TruncatedPareto(1.0, 3.0, 1000),
new Weibull(1.0, 1.0),
};

426
src/Numerics.Tests/DistributionTests/Continuous/BurrTests.cs

@ -0,0 +1,426 @@
// <copyright file="BurrTests.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
//
// Copyright (c) 2009-2019 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using MathNet.Numerics.Distributions;
using NUnit.Framework;
using System;
using System.Linq;
namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
{
/// <summary>
/// <c>Burr</c> distribution tests.
/// </summary>
[TestFixture, Category("Distributions")]
public class BurrTests
{
private const int highPrecision = 12;
private const int lowPrecision = 8;
/// <summary>
/// Can create <c>Burr</c>.
/// </summary>
/// <param name="a">a parameter.</param>
/// <param name="c">c parameter.</param>
/// <param name="k">k parameter.</param>
[TestCase(1.0, 1.0, 0.5)]
[TestCase(10.0, 0.1, 0.2)]
[TestCase(5.0, 1.0, 3)]
[TestCase(2.0, 10.0, 2.0)]
[TestCase(10.0, 100.0, 1.0)]
public void CanCreateBurr(double a, double c, double k)
{
var n = new Burr(a, c, k);
Assert.AreEqual(a, n.a);
Assert.AreEqual(c, n.c);
Assert.AreEqual(k, n.k);
}
/// <summary>
/// Can create <c>Burr</c> with random source.
/// </summary>
[Test]
public void CanCreateBurrWithRandomSource()
{
var randomSource = new Numerics.Random.MersenneTwister(100);
var n = new Burr(1.0, 1.0, 0.5, randomSource);
Assert.AreEqual(randomSource, n.RandomSource);
}
/// <summary>
/// <c>Burr</c> create fails with bad parameters.
/// </summary>
/// <param name="a">a parameter.</param>
/// <param name="c">c parameter.</param>
/// <param name="k">k parameter.</param>
[TestCase(5.0, 1.0, -1.0)]
[TestCase(-1.0, 1.0, 2.05)]
[TestCase(1.0, -1.0, 2.05)]
[TestCase(Double.NaN, Double.NaN, Double.NaN)]
[TestCase(Double.NaN, 1.0, 2.4)]
[TestCase(1.0, 1.0, double.PositiveInfinity)]
[TestCase(1.0, -1.0, Double.NegativeInfinity)]
public void BurrCreateFailsWithBadParameters(double a, double c, double k)
{
Assert.That(() => new Burr(a, c, k), Throws.ArgumentException);
}
/// <summary>
/// Validate IsValidParameterSet.
/// </summary>
/// <param name="a">a parameter.</param>
/// <param name="c">c parameter.</param>
/// <param name="k">k parameter.</param>
/// <param name="validity">exoected validity of paremeter set.</param>
[TestCase(1.0, 1.0, 0.5, true)]
[TestCase(10.0, 0.1, 0.2, true)]
[TestCase(5.0, 1.0, -1.0, false)]
[TestCase(1.0, 1.0, double.PositiveInfinity, false)]
public void ValidateIsValidParameterSet(double a, double c, double k, bool validity)
{
Assert.AreEqual(Burr.IsValidParameterSet(a, c, k), validity);
}
/// <summary>
/// Validate ToString.
/// </summary>
[Test]
public void ValidateToString()
{
var n = new Burr(1d, 2d, 3d);
Assert.AreEqual("Burr(a = 1, c = 2, k = 3)", n.ToString());
}
/// <summary>
/// Validate mode.
/// </summary>
/// <param name="a">a parameter.</param>
/// <param name="c">c parameter.</param>
/// <param name="k">k parameter.</param>
/// <param name="mode">Expected value.</param>
[TestCase(2.000000, 2.00000, 1.0, 1.154700538379251)]
[TestCase(1.000000, 1.00000, 0.5, 0.0)]
[TestCase(3.000000, 4.00000, 1.0, 2.640335210380180)]
public void ValidateMode(double a, double c, double k, double mode)
{
var n = new Burr(a, c, k);
AssertHelpers.AlmostEqualRelative(mode, n.Mode, highPrecision);
}
/// <summary>
/// Validate median.
/// </summary>
/// <param name="a">a parameter.</param>
/// <param name="c">c parameter.</param>
/// <param name="k">k parameter.</param>
/// <param name="median">Expected value.</param>
[TestCase(1.000000, 1.00000, 0.5, 3.0)]
[TestCase(1.000000, 0.50000, 0.5, 9.0)]
[TestCase(1.000000, 1.00000, 5.0, 0.148698354997035)]
public void ValidateMedian(double a, double c, double k, double median)
{
var n = new Burr(a, c, k);
AssertHelpers.AlmostEqualRelative(median, n.Median, highPrecision);
}
/// <summary>
/// Validate mean.
/// </summary>
/// <param name="a">a parameter.</param>
/// <param name="c">c parameter.</param>
/// <param name="k">k parameter.</param>
/// <param name="mean">Expected value.</param>
[TestCase(1.000000, 1.00000, 1.5, 2)]
[TestCase(4.000000, 5.00000, 0.5, 6.198785110989412)]
[TestCase(4.000000, 5.00000, 5.0, 2.729694550490384)]
public void ValidateMean(double a, double c, double k, double mean)
{
var n = new Burr(a, c, k);
AssertHelpers.AlmostEqualRelative(mean, n.Mean, highPrecision);
}
/// <summary>
/// Validate variance.
/// </summary>
/// <param name="a">a parameter.</param>
/// <param name="c">c parameter.</param>
/// <param name="k">k parameter.</param>
/// <param name="variance">Expected value.</param>
[TestCase(4.0, 5.0, 2.0, 0.983559126843161)]
[TestCase(2.0, 3.5, 4.0, 0.207489170174404)]
public void ValidateVariance(double a, double c, double k, double variance)
{
var n = new Burr(a, c, k);
AssertHelpers.AlmostEqualRelative(variance, n.Variance, highPrecision);
}
/// <summary>
/// Validate standard deviation.
/// </summary>
/// <param name="a">a parameter.</param>
/// <param name="c">c parameter.</param>
/// <param name="k">k parameter.</param>
/// <param name="std">Expected value.</param>
[TestCase(4.0, 5.0, 2.0, 0.991745494995143)]
[TestCase(2.0, 3.5, 4.0, 0.455509791524182)]
public void ValidateStandardDeviation(double a, double c, double k, double std)
{
var n = new Burr(a, c, k);
AssertHelpers.AlmostEqualRelative(std, n.StdDev, highPrecision);
}
/// <summary>
/// Validate minimum.
/// </summary>
[Test]
public void ValidateMinimum()
{
var n = new Burr(1.0, 2.0, 2.0);
Assert.AreEqual(0.0, n.Minimum);
}
/// <summary>
/// Validate maximum.
/// </summary>
[Test]
public void ValidateMaximum()
{
var n = new Burr(1.0, 2.0, 1.0);
Assert.AreEqual(Double.PositiveInfinity, n.Maximum);
}
/// <summary>
/// Validate entropy throws not supported exception
/// </summary>
/// <param name="a">a parameter.</param>
/// <param name="c">c parameter.</param>
/// <param name="k">k parameter.</param>
[TestCase(1.0, 1.0, 0.5)]
public void ValidateEntropyFailsWithNotSupported(double a, double c, double k)
{
var n = new Burr(a, c, k);
Assert.Throws<NotSupportedException>(() => { var x = n.Entropy; });
}
/// <summary>
/// Validate GetMoments.
/// </summary>
/// <param name="a">a parameter.</param>
/// <param name="c">c parameter.</param>
/// <param name="k">k parameter.</param>
/// <param name="n"> order of the moment.</param>
/// <param name="value">Expected value.</param>
[TestCase(1.000000, 1.00000, 1.5, 1, 2)]
[TestCase(4.000000, 5.00000, 0.5, 1, 6.198785110989412)]
[TestCase(4.000000, 5.00000, 5.0, 1, 2.729694550490384)]
[TestCase(4.0, 5.0, 2.0, 3, 50.738165747621750)]
[TestCase(2.0, 3.5, 4.0, 3, 2.895046685294824)]
public void ValidateMoments(double a, double c, double k, int order, double value)
{
var n = new Burr(a, c, k);
AssertHelpers.AlmostEqualRelative(value, n.GetMoment(order), lowPrecision);
}
/// <summary>
/// Validate GetMoments throws when given bad parameters
/// </summary>
/// <param name="a">a parameter.</param>
/// <param name="c">c parameter.</param>
/// <param name="k">k parameter.</param>
/// <param name="order">the order of the moment.</param>
[TestCase(1.0, 1.0, 0.5, 1)]
[TestCase(10.0, 1.0, 1.5, 2)]
[TestCase(5.0, 2, 1.4, 3)]
public void ValidateGetMomentsFailsWithBadParameters(double a, double c, double k, double order)
{
var n = new Burr(a, c, k);
Assert.That(() => n.GetMoment(order), Throws.ArgumentException);
}
/// <summary>
/// Validate skewness.
/// </summary>
/// <param name="a">a parameter.</param>
/// <param name="c">c parameter.</param>
/// <param name="k">k parameter.</param>
/// <param name="skewness">Expected value.</param>
[TestCase(4.0, 5.0, 2.0, 0.635277200891842)]
[TestCase(2.0, 3.5, 4.0, 0.483073007212360)]
public void ValidateSkewness(double a, double c, double k, double skewness)
{
var n = new Burr(a, c, k);
AssertHelpers.AlmostEqualRelative(skewness, n.Skewness, lowPrecision);
}
/// <summary>
/// Validate density.
/// </summary>
/// <param name="a">a parameter.</param>
/// <param name="c">c parameter.</param>
/// <param name="k">k parameter.</param>
/// <param name="x">Input X value.</param>
/// <param name="p">Expected value.</param>
[TestCase(1.00000, 1.000000, 0.5, 0.100000, 0.433392086020724)]
[TestCase(1.000000, 1.000000, 5.0, 0.500000, 0.438957475994512)]
public void ValidateDensity(double a, double c, double k, double x, double p)
{
var n = new Burr(a, c, k);
AssertHelpers.AlmostEqualRelative(p, n.Density(x), highPrecision);
AssertHelpers.AlmostEqualRelative(p, Burr.PDF(a, c, k, x), highPrecision);
}
/// <summary>
/// Validate density log.
/// </summary>
/// <param name="a">a parameter.</param>
/// <param name="c">c parameter.</param>
/// <param name="k">k parameter.</param>
/// <param name="x">Input X value.</param>
/// <param name="p">Expected value.</param>
[TestCase(1.000000, 1.000000, 0.5, 0.500000, -1.301344842722192)]
[TestCase(1.000000, 1.000000, 5.0, 0.500000, -0.823352736214886)]
[TestCase(4.000000, 5.000000, 1.0, 5.000000, -1.682583872895781)]
public void ValidateDensityLn(double a, double c, double k, double x, double p)
{
var n = new Burr(a, c, k);
AssertHelpers.AlmostEqualRelative(p, n.DensityLn(x), highPrecision);
AssertHelpers.AlmostEqualRelative(p, Burr.PDFLn(a, c, k, x), highPrecision);
}
/// <summary>
/// Validate cumulative distribution.
/// </summary>
/// <param name="a">a parameter.</param>
/// <param name="c">c parameter.</param>
/// <param name="k">k parameter.</param>
/// <param name="x">The location at which to compute the density.</param>
/// <param name="f">Expected value.</param>
[TestCase(1.000000, 1.000000, 0.5, 0.500000, 0.183503419072274)]
[TestCase(4.000000, 5.000000, 0.5, 0.500000, 0.00001525843981653452)]
[TestCase(4.000000, 5.000000, 0.5, 5.00000, 0.503203804978536)]
[TestCase(4.000000, 5.000000, 5.0, 5.00000, 0.999084238019496)]
public void ValidateCumulativeDistribution(double a, double c, double k, double x, double f)
{
var n = new Burr(a, c, k);
AssertHelpers.AlmostEqualRelative(f, n.CumulativeDistribution(x), lowPrecision);
AssertHelpers.AlmostEqualRelative(f, Burr.CDF(a, c, k, x), lowPrecision);
}
/// <summary>
/// Can sample static.
/// </summary>
[Test]
public void CanSampleStatic()
{
Burr.Sample(new Numerics.Random.MersenneTwister(100), 1.0, 2.0, 2.0);
}
/// <summary>
/// Can fill array with samples static.
/// </summary>
[Test]
public void CanFillSampleArrayStatic()
{
double[] samples = new double[100];
Burr.Samples(new Numerics.Random.MersenneTwister(100), samples, 1.0, 2.0, 2.0);
Assert.IsTrue(!samples.Any(x => x == 0));
}
/// <summary>
/// Can sample sequence static.
/// </summary>
[Test]
public void CanSampleSequenceStatic()
{
var ied = Burr.Samples(new Numerics.Random.MersenneTwister(100), 1.0, 2.0, 2.0);
GC.KeepAlive(ied.Take(5).ToArray());
}
/// <summary>
/// Fail sample static with bad parameters.
/// </summary>
[Test]
public void FailSampleStatic()
{
Assert.That(() => { var d = Burr.Sample(new Numerics.Random.MersenneTwister(100), 1.0, -1.0, 2.0); }, Throws.ArgumentException);
}
/// <summary>
/// Fail filling array with samples with bad parameters static.
/// </summary>
[Test]
public void FailFillingSampleArrayStatic()
{
double[] samples = new double[100];
Assert.That(() => { Burr.Samples(new Numerics.Random.MersenneTwister(100), samples, -1.0, 2.0, 2.0); }, Throws.ArgumentException);
}
/// <summary>
/// Fail sample sequence static with bad parameters.
/// </summary>
[Test]
public void FailSampleSequenceStatic()
{
Assert.That(() => { var ied = Burr.Samples(new Numerics.Random.MersenneTwister(100), 1.0, -1.0, 2.0).First(); }, Throws.ArgumentException);
}
/// <summary>
/// Can sample.
/// </summary>
[Test]
public void CanSample()
{
var n = new Burr(1.0, 2.0, 2.0);
n.Sample();
}
/// <summary>
/// Can fill array with samples.
/// </summary>
[Test]
public void CanFillSampleArray()
{
double[] samples = new double[100];
var n = new Burr(1.0, 2.0, 2.0, new Numerics.Random.MersenneTwister(100));
n.Samples(samples);
Assert.IsTrue(!samples.Any(x => x == 0));
}
/// <summary>
/// Can sample sequence.
/// </summary>
[Test]
public void CanSampleSequence()
{
var n = new Burr(1.0, 2.0, 2.0);
var ied = n.Samples();
GC.KeepAlive(ied.Take(5).ToArray());
}
}
}

419
src/Numerics.Tests/DistributionTests/Continuous/InverseGaussianTests.cs

@ -0,0 +1,419 @@
// <copyright file="InverseGaussianTests.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
//
// Copyright (c) 2009-2019 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using MathNet.Numerics.Distributions;
using NUnit.Framework;
using System;
using System.Linq;
namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
{
/// <summary>
/// <c>InverseGaussian</c> distribution tests.
/// </summary>
[TestFixture, Category("Distributions")]
public class InverseGaussianTests
{
private const int precision = 12;
/// <summary>
/// Can create <c>Inverse Gaussian</c>.
/// </summary>
/// <param name="mu">Mu value.</param>
/// <param name="lambda">Lambda parameter.</param>
[TestCase(0.1, 0.1)]
[TestCase(1.0, 1.0)]
[TestCase(10.0, 10.0)]
public void CanCreateInverseGaussian(double mu, double lambda)
{
var n = new InverseGaussian(mu, lambda);
Assert.AreEqual(mu, n.Mu);
Assert.AreEqual(lambda, n.Lambda);
}
/// <summary>
/// Can create <c>Inverse Gaussian</c> with random source.
/// </summary>
[Test]
public void CanCreateInverseGaussianWithRandomSource()
{
var randomSource = new Numerics.Random.MersenneTwister(100);
var n = new InverseGaussian(1.0, 1.0, randomSource);
Assert.AreEqual(randomSource, n.RandomSource);
}
/// <summary>
/// <c>InverseGaussian</c> create fails with bad parameters.
/// </summary>
/// <param name="mu">Mu parameter.</param>
/// <param name="lambda">Lambda parameter.</param>
[TestCase(Double.NaN, 1.0)]
[TestCase(1.0, Double.NaN)]
[TestCase(Double.NaN, Double.NaN)]
[TestCase(-1.0, -1.0)]
public void InverseGaussianCreateFailsWithBadParameters(double mu, double lambda)
{
Assert.That(() => new LogNormal(mu, lambda), Throws.ArgumentException);
}
/// <summary>
/// Validate IsValidParameterSet.
/// </summary>
/// <param name="mu">Mu parameter.</param>
/// <param name="lambda">Lambda parameter.</param>
/// <param name="validity">exoected validity of paremeter set.</param>
[TestCase(1.0, 1.0, true)]
[TestCase(10.0, 10.0, true)]
[TestCase(-1.0, -1.0, false)]
public void ValidateIsValidParameterSet(double mu, double lambda, bool validity)
{
Assert.AreEqual(InverseGaussian.IsValidParameterSet(mu, lambda), validity);
}
/// <summary>
/// Validate ToString.
/// </summary>
[Test]
public void ValidateToString()
{
var n = new InverseGaussian(1d, 2d);
Assert.AreEqual("InverseGaussian(μ = 1, λ = 2)", n.ToString());
}
/// <summary>
/// Validate mode.
/// </summary>
/// <param name="mu">Mu parameter.</param>
/// <param name="lambda">Lambda parameter.</param>
/// <param name="mode">Expected value.</param>
[TestCase(0.100000, 0.100000, 0.030277563773199)]
[TestCase(1.500000, 5.500000, 1.007026953273370)]
[TestCase(2.500000, 1.500000, 0.481456008918130)]
[TestCase(5.500000, 2.500000, 0.815033614523337)]
[TestCase(5.500000, 5.500000, 1.665266007525970)]
public void ValidateMode(double mu, double lambda, double mode)
{
var n = new InverseGaussian(mu, lambda);
AssertHelpers.AlmostEqualRelative(mode, n.Mode, precision);
}
/// <summary>
/// Validate median.
/// </summary>
/// <param name="mu">Mu parameter.</param>
/// <param name="lambda">Lambda parameter.</param>
/// <param name="median">Expected value.</param>
[TestCase(0.100000, 0.100000, 0.067584130569524)]
[TestCase(1.500000, 0.100000, 0.190369896722453)]
[TestCase(2.500000, 0.100000, 0.201110669500521)]
[TestCase(5.500000, 5.500000, 3.717127181323815)]
public void ValidateMedian(double mu, double lambda, double median)
{
var n = new InverseGaussian(mu, lambda);
AssertHelpers.AlmostEqualRelative(median, n.Median, precision);
}
/// <summary>
/// Validate mean.
/// </summary>
/// <param name="mu">Mu parameter.</param>
/// <param name="lambda">Lambda parameter.</param>
/// <param name="mean">Expected value.</param>
[TestCase(0.100000, 5.500000, 0.1)]
[TestCase(1.500000, 0.100000, 1.5)]
[TestCase(2.500000, 5.500000, 2.5)]
[TestCase(5.500000, 0.100000, 5.5)]
public void ValidateMean(double mu, double lambda, double mean)
{
var n = new InverseGaussian(mu, lambda);
AssertHelpers.AlmostEqualRelative(mean, n.Mean, precision);
}
/// <summary>
/// Validate variance.
/// </summary>
/// <param name="mu">Mu parameter.</param>
/// <param name="lambda">Lambda parameter.</param>
/// <param name="variance">Expected value.</param>
[TestCase(0.100000, 5.500000, 1.818181818181819e-04)]
[TestCase(1.500000, 0.100000, 33.75)]
public void ValidateVariance(double mu, double lambda, double variance)
{
var n = new InverseGaussian(mu, lambda);
AssertHelpers.AlmostEqualRelative(variance, n.Variance, precision);
}
/// <summary>
/// Validate standard deviation.
/// </summary>
/// <param name="mu">Mu parameter.</param>
/// <param name="lambda">Lambda parameter.</param>
/// <param name="std">Expected value.</param>
[TestCase(2.500000, 5.500000, 1.685499656158105)]
[TestCase(5.500000, 0.100000, 40.789091679026143)]
public void ValidateStandardDeviation(double mu, double lambda, double std)
{
var n = new InverseGaussian(mu, lambda);
AssertHelpers.AlmostEqualRelative(std, n.StdDev, precision);
}
/// <summary>
/// Validate minimum.
/// </summary>
[Test]
public void ValidateMinimum()
{
var n = new InverseGaussian(1.0, 2.0);
Assert.AreEqual(0.0, n.Minimum);
}
/// <summary>
/// Validate maximum.
/// </summary>
[Test]
public void ValidateMaximum()
{
var n = new InverseGaussian(1.0, 2.0);
Assert.AreEqual(Double.PositiveInfinity, n.Maximum);
}
/// <summary>
/// Validate entropy throws not supported exception
/// </summary>
/// <param name="mu">Mu parameter.</param>
/// <param name="lambda">Lambda parameter.</param>
[TestCase(1.0, 1.0)]
public void ValidateEntropyFailsWithNotSupported(double mu, double lambda)
{
var n = new InverseGaussian(mu, lambda);
Assert.Throws<NotSupportedException>(() => { var x = n.Entropy; });
}
/// <summary>
/// Validate skewness.
/// </summary>
/// <param name="mu">Mu parameter.</param>
/// <param name="´lambda">Lambda parameter.</param>
/// <param name="skewness">Expected value.</param>
[TestCase(0.100000, 1.500000, 0.774596669241483)]
[TestCase(1.500000, 2.500000, 2.323790007724450)]
[TestCase(1.500000, 5.500000, 1.566698903601281)]
[TestCase(2.500000, 0.100000, 15.0)]
[TestCase(2.500000, 1.500000, 3.872983346207417)]
[TestCase(5.500000, 2.500000, 4.449719092257398)]
[TestCase(5.500000, 5.500000, 3.0)]
public void ValidateSkewness(double mu, double lambda, double skewness)
{
var n = new InverseGaussian(mu, lambda);
AssertHelpers.AlmostEqualRelative(skewness, n.Skewness, precision);
}
/// <summary>
/// Validate density.
/// </summary>
/// <param name="mu">Mu parameter.</param>
/// <param name="lambda">Lambda parameter.</param>
/// <param name="x">The location at which to compute the density.</param>
/// <param name="p">Expected value.</param>
[TestCase(1.500000, 2.500000, 0.500000, 0.587321148416470)]
[TestCase(1.500000, 2.500000, 0.800000, 0.627284170789435)]
[TestCase(2.500000, 2.500000, 0.500000, 0.360208446721537)]
[TestCase(2.500000, 2.500000, 0.800000, 0.428023216678204)]
public void ValidateDensity(double mu, double lambda, double x, double p)
{
var n = new InverseGaussian(mu, lambda);
AssertHelpers.AlmostEqualRelative(p, n.Density(x), precision);
AssertHelpers.AlmostEqualRelative(p, InverseGaussian.PDF(mu, lambda, x), precision);
}
/// <summary>
/// Validate density log.
/// </summary>
/// <param name="mu">Mu parameter.</param>
/// <param name="lambda">Sigma value.</param>
/// <param name="x">The location at which to compute the log density.</param>
/// <param name="p">Expected value.</param>
[TestCase(1.500000, 0.100000, 0.100000, 0.948091004233817)]
[TestCase(1.500000, 1.500000, 0.800000, -0.585657318845943)]
[TestCase(2.500000, 0.100000, 0.800000, -1.764415752730381)]
[TestCase(2.500000, 1.500000, 0.100000, -4.174328339659523)]
[TestCase(2.500000, 1.500000, 0.500000, -0.636485208310673)]
public void ValidateDensityLn(double mu, double lambda, double x, double p)
{
var n = new InverseGaussian(mu, lambda);
AssertHelpers.AlmostEqualRelative(p, n.DensityLn(x), precision);
AssertHelpers.AlmostEqualRelative(p, InverseGaussian.PDFLn(mu, lambda, x), precision);
}
/// <summary>
/// Validate cumulative distribution.
/// </summary>
/// <param name="mu">Mu parameter.</param>
/// <param name="lambda">Sigma value.</param>
/// <param name="x">The location at which to compute the density.</param>
/// <param name="f">Expected value.</param>
[TestCase(1.500000, 2.500000, 0.100000, 0.000002882116410867029)]
[TestCase(2.500000, 1.500000, 0.100000, 0.0001938001952257318)]
[TestCase(2.500000, 1.500000, 0.500000, 0.145457623130791)]
[TestCase(2.500000, 2.500000, 0.100000, 0.000001529605202470422)]
[TestCase(2.500000, 2.500000, 0.800000, 0.187168922781367)]
public void ValidateCumulativeDistribution(double mu, double lambda, double x, double f)
{
var n = new InverseGaussian(mu, lambda);
AssertHelpers.AlmostEqualRelative(f, n.CumulativeDistribution(x), precision);
AssertHelpers.AlmostEqualRelative(f, InverseGaussian.CDF(mu, lambda, x), precision);
}
/// <summary>
/// Validate inverse cumulative distribution.
/// </summary>
/// <param name="mu">Mu parameter.</param>
/// <param name="lambda">Sigma value.</param>
/// <param name="probability">The location at which to compute the inverse cumulative distribution function.</param>
/// <param name="f">Expected value.</param>
[TestCase(1.500000, 1.500000, 0.100000, 0.356437063090717)]
[TestCase(1.500000, 1.500000, 0.500000, 1.013761958542859)]
[TestCase(2.500000, 2.500000, 0.500000, 1.689603264238098)]
[TestCase(2.500000, 2.500000, 0.800000, 3.619719792074686)]
public void ValidateInverseCumulativeDistribution(double mu, double lambda, double probability, double f)
{
var n = new InverseGaussian(mu, lambda);
AssertHelpers.AlmostEqualRelative(f, InverseGaussian.ICDF(mu, lambda, probability), precision);
AssertHelpers.AlmostEqualRelative(f, n.InvCDF(probability), precision);
}
/// <summary>
/// Can sample static.
/// </summary>
[Test]
public void CanSampleStatic()
{
InverseGaussian.Sample(new Numerics.Random.MersenneTwister(100), 1.0, 1.0);
}
/// <summary>
/// Can fill array with samples static.
/// </summary>
[Test]
public void CanFillSampleArrayStatic()
{
double[] samples = new double[100];
InverseGaussian.Samples(new Numerics.Random.MersenneTwister(100), samples, 1.0, 1.0);
Assert.IsTrue(!samples.Any(x => x == 0));
}
/// <summary>
/// Can sample sequence static.
/// </summary>
[Test]
public void CanSampleSequenceStatic()
{
var ied = InverseGaussian.Samples(new Numerics.Random.MersenneTwister(100), 1.0, 1.0);
GC.KeepAlive(ied.Take(5).ToArray());
}
/// <summary>
/// Fail sample static with bad parameters.
/// </summary>
[Test]
public void FailSampleStatic()
{
Assert.That(() => { var d = InverseGaussian.Sample(new Numerics.Random.MersenneTwister(100), 0.0, -1.0); }, Throws.ArgumentException);
}
/// <summary>
/// Fail filling array with samples with bad parameters static.
/// </summary>
[Test]
public void FailFillingSampleArrayStatic()
{
double[] samples = new double[100];
Assert.That(() => { InverseGaussian.Samples(new Numerics.Random.MersenneTwister(100), samples, -1.0, 1.0); }, Throws.ArgumentException);
}
/// <summary>
/// Fail sample sequence static with bad parameters.
/// </summary>
[Test]
public void FailSampleSequenceStatic()
{
Assert.That(() => { var ied = InverseGaussian.Samples(new Numerics.Random.MersenneTwister(100), 0.0, -1.0).First(); }, Throws.ArgumentException);
}
/// <summary>
/// Can sample.
/// </summary>
[Test]
public void CanSample()
{
var n = new InverseGaussian(1.0, 2.0);
n.Sample();
}
/// <summary>
/// Can fill array with samples.
/// </summary>
[Test]
public void CanFillSampleArray()
{
double[] samples = new double[100];
var n = new InverseGaussian(1.0, 2.0, new Numerics.Random.MersenneTwister(100));
n.Samples(samples);
Assert.IsTrue(!samples.Any(x => x == 0));
}
/// <summary>
/// Can sample sequence.
/// </summary>
[Test]
public void CanSampleSequence()
{
var n = new InverseGaussian(1.0, 2.0);
var ied = n.Samples();
GC.KeepAlive(ied.Take(5).ToArray());
}
/// <summary>
/// Can estimate distribution parameters.
/// </summary>
[TestCase(10.0, 0.1)]
[TestCase(1.5, 1.5)]
[TestCase(2.5, 2.5)]
[TestCase(2.5, 5.0)]
[TestCase(10.0, 50.0)]
public void CanEstimateParameters(double mu, double lambda)
{
var original = new InverseGaussian(mu, lambda, new Numerics.Random.MersenneTwister(100));
var estimated = InverseGaussian.Estimate(original.Samples().Take(1000000));
AssertHelpers.AlmostEqualRelative(mu, estimated.Mu, 1);
AssertHelpers.AlmostEqualRelative(lambda, estimated.Lambda, 1);
}
}
}

446
src/Numerics.Tests/DistributionTests/Continuous/TruncatedParetoTests.cs

@ -0,0 +1,446 @@
// <copyright file="TruncatedParetoTests.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
//
// Copyright (c) 2009-2019 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using MathNet.Numerics.Distributions;
using NUnit.Framework;
using System;
using System.Linq;
namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
{
/// <summary>
/// <c>TruncatedPareto</c> distribution tests.
/// </summary>
[TestFixture, Category("Distributions")]
public class TruncatedParetoTests
{
private const int highPrecision = 12;
private const int lowPrecision = 6;
private const double highTruncation = 1E8;
/// <summary>
/// Can create <c>TruncatedPareto</c>.
/// </summary>
/// <param name="scale">Scale parameter.</param>
/// <param name="shape">Shape parameter.</param>
/// /// <param name="truncation">Truncation parameter.</param>
[TestCase(2.0, 2.0, 500.0)]
[TestCase(10.0, 2.0, 1000.0)]
[TestCase(100.0, 4.0, 10000.0)]
[TestCase(3.0, 10.0, 5.0)]
[TestCase(10.0, 100.0, 20.0)]
public void CanCreateTruncatedPareto(double scale, double shape, double truncation)
{
var n = new TruncatedPareto(scale, shape, truncation);
Assert.AreEqual(scale, n.Scale, highPrecision);
Assert.AreEqual(shape, n.Shape, highPrecision);
Assert.AreEqual(truncation, n.Truncation, highPrecision);
}
/// <summary>
/// Can create <c>TruncatedPareto</c> with random source.
/// </summary>
[Test]
public void CanCreateTruncatedParetoWithRandomSource()
{
var randomSource = new Numerics.Random.MersenneTwister(100);
var n = new TruncatedPareto(10.0, 10.0, 1000, randomSource);
Assert.AreEqual(randomSource, n.RandomSource);
}
/// <summary>
/// <c>TruncatedPareto</c> create fails with bad parameters.
/// </summary>
/// <param name="scale">Scale parameter.</param>
/// <param name="shape">Shape parameter.</param>
/// /// <param name="truncation">Truncation parameter.</param>
[TestCase(1.0, -1.0, 15.0)]
[TestCase(1.0, 3.0, 0.5)]
[TestCase(-1.0, 2.0, 15.0)]
[TestCase(double.NaN, 1.0, 1.0)]
[TestCase(1.0, double.PositiveInfinity, 0.0)]
[TestCase(5.0, 2.0, double.PositiveInfinity)]
public void TruncatedParetoCreateFailsWithBadParameters(double scale, double shape, double truncation)
{
Assert.That(() => new TruncatedPareto(scale, shape, truncation), Throws.ArgumentException);
}
/// <summary>
/// Validate IsValidParameterSet.
/// </summary>
/// <param name="scale">Scale parameter.</param>
/// <param name="shape">Shape parameter.</param>
/// /// <param name="truncation">Truncation parameter.</param>
[TestCase(2.0, 2.0, 500, true)]
[TestCase(10.0, 10.0, 1000, true)]
[TestCase(5.0, 3.0, 1, false)]
public void ValidateIsValidParameterSet(double scale, double shape, double truncation, bool validity)
{
Assert.AreEqual(TruncatedPareto.IsValidParameterSet(scale, shape, truncation), validity);
}
/// <summary>
/// Validate ToString.
/// </summary>
[Test]
public void ValidateToString()
{
var n = new TruncatedPareto(1d, 2d, 100d);
Assert.AreEqual("Truncated Pareto(Scale = 1, Shape = 2, Truncation = 100)", n.ToString());
}
/// <summary>
/// Validate mode throws not supported exception
/// </summary>
/// <param name="scale">Scale value.</param>
/// <param name="shape">Shape parameter.</param>
/// <param name="truncation">Truncation parameter.</param>
[TestCase(10.0, 10.0, highTruncation)]
public void ValidateModeFailsWithNotSupported(double scale, double shape, double truncation)
{
var n = new TruncatedPareto(scale, shape, truncation);
Assert.Throws<NotSupportedException>(() => { var x = n.Mode; });
}
/// <summary>
/// Validate median.
/// </summary>
/// <param name="scale">Scale value.</param>
/// <param name="shape">Shape parameter.</param>
/// <param name="truncation">Truncation parameter.</param>
/// <param name="median">Expected value.</param>
[TestCase(10.0, 10.0, highTruncation, 10.717734625362931)]
[TestCase(100, 3.5, 10000, 1.219013619375516e+02)]
[TestCase(1000, 5.5, 100000, 1.134312522193400e+03)]
public void ValidateMedian(double scale, double shape, double truncation, double median)
{
var n = new TruncatedPareto(scale, shape, truncation);
AssertHelpers.AlmostEqualRelative(median, n.Median, highPrecision);
}
/// <summary>
/// Validate mean.
/// </summary>
/// <param name="scale">Scale parameter.</param>
/// <param name="shape">Shape parameter.</param>
/// <param name="truncation">Truncation parameter.</param>
/// <param name="std">Expected value.</param>
[TestCase(10.0, 10.0, highTruncation, 11.111111111111111)]
[TestCase(100, 3.5, 10000, 1.399986139998614e+02)]
public void ValidateMean(double scale, double shape, double truncation, double mean)
{
var n = new TruncatedPareto(scale, shape, truncation);
AssertHelpers.AlmostEqualRelative(mean, n.Mean, highPrecision);
}
/// <summary>
/// Validate variance.
/// </summary>
/// <param name="scale">Scale parameter.</param>
/// <param name="shape">Shape parameter.</param>
/// <param name="truncation">Truncation parameter.</param>
/// <param name="variance">Expected value.</param>
[TestCase(10.0, 10.0, highTruncation, 1.543209876543210)]
[TestCase(100, 3.5, 10000, 3.710390409118045e+03)]
[TestCase(1000, 5.5, 100000, 7.760125676537910e+04)]
public void ValidateVariance(double scale, double shape, double truncation, double variance)
{
var n = new TruncatedPareto(scale, shape, truncation);
AssertHelpers.AlmostEqualRelative(variance, n.Variance, highPrecision);
}
/// <summary>
/// Validate standard deviation.
/// </summary>
/// <param name="scale">Scale parameter.</param>
/// <param name="shape">Shape parameter.</param>
/// <param name="truncation">Truncation parameter.</param>
/// <param name="std">Expected value.</param>
[TestCase(10.0, 10.0, highTruncation, 1.242259987499883)]
[TestCase(100, 3.5, 10000, 60.912974062329653)]
[TestCase(1000, 5.5, 100000, 2.785700212969427e+02)]
public void ValidateStandardDeviation(double scale, double shape, double truncation, double std)
{
var n = new TruncatedPareto(scale, shape, truncation);
AssertHelpers.AlmostEqualRelative(std, n.StdDev, highPrecision);
}
/// <summary>
/// Validate minimum.
/// </summary>
[Test]
public void ValidateMinimum()
{
var n = new TruncatedPareto(1.0, 2.0, 500.0);
Assert.AreEqual(1.0, n.Minimum);
}
/// <summary>
/// Validate maximum.
/// </summary>
[Test]
public void ValidateMaximum()
{
var n = new TruncatedPareto(1.0, 2.0, 500.0);
Assert.AreEqual(500.0, n.Maximum);
}
/// <summary>
/// Validate entropy throws not supported exception
/// </summary>
/// <param name="mu">Mu parameter.</param>
/// <param name="lambda">Lambda parameter.</param>
[TestCase(100, 3.5, 10000)]
public void ValidateEntropyFailsWithNotSupported(double scale, double shape, double truncation)
{
var n = new TruncatedPareto(scale, shape, truncation);
Assert.Throws<NotSupportedException>(() => { var x = n.Entropy; });
}
/// <summary>
/// Validate GetMoments.
/// </summary>
/// <param name="scale">Scale parameter.</param>
/// <param name="shape">Shape parameter.</param>
/// <param name="truncation">Truncation parameter.</param>
/// <param name="n"> order of the moment.</param>
/// <param name="value">Expected value.</param>
[TestCase(10.0, 10.0, highTruncation, 1, 11.111111111111111)]
[TestCase(1.0, 6.0, 10000, 1, 1.2)]
[TestCase(100, 3.5, 10000, 1, 1.399986139998614e+02)]
[TestCase(100, 3.5, 10000, 2, 2.331000233100023e+04)]
[TestCase(100, 3.5, 10000, 3, 6.300000630000063e+06)]
[TestCase(1000, 5.5, 100000, 1, 1.222222221012222e+03)]
[TestCase(1000, 5.5, 100000, 2, 1.571428414301428e+06)]
[TestCase(1000, 5.5, 100000, 4, 3.663000000036630e+12)]
[TestCase(1000, 2.0, 100000, 2, 9.211261498125996e+06)]
[TestCase(1000, 3.0, 100000, 3, 1.381552437348865e+10)]
public void ValidateMoments(double scale, double shape, double truncation, int order, double value)
{
var n = new TruncatedPareto(scale, shape, truncation);
AssertHelpers.AlmostEqualRelative(value, n.GetMoment(order), highPrecision);
}
/// <summary>
/// Validate skewness.
/// </summary>
/// <param name="scale">Scale parameter.</param>
/// <param name="shape">Shape parameter.</param>
/// <param name="truncation">Truncation parameter.</param>
/// <param name="skewness">Expected value.</param>
[TestCase(10.0, 10.0, highTruncation, 2.8110568859997356)]
public void ValidateSkewness(double scale, double shape, double truncation, double skewness)
{
var n = new TruncatedPareto(scale, shape, truncation);
AssertHelpers.AlmostEqualRelative(skewness, n.Skewness, highPrecision);
}
/// <summary>
/// Validate density.
/// </summary>
/// <param name="scale">Scale parameter.</param>
/// <param name="shape">Shape parameter.</param>
/// <param name="truncation">Truncation parameter.</param>
/// <param name="x">The location at which to compute the density.</param>
/// <param name="p">Expected value.</param>
[TestCase(1, 1, highTruncation, 1, 1)]
[TestCase(1, 1, highTruncation, 1.5, 4 / 9.0)]
[TestCase(1, 1, highTruncation, 5, 1 / 25.0)]
[TestCase(1, 1, highTruncation, 50, 1 / 2500.0)]
[TestCase(1, 4, highTruncation, 1, 4)]
[TestCase(1, 4, highTruncation, 1.5, 128 / 243.0)]
[TestCase(1, 4, highTruncation, 50, 1 / 78125000.0)]
[TestCase(3, 2, highTruncation, 3, 2 / 3.0)]
[TestCase(3, 2, highTruncation, 5, 18 / 125.0)]
[TestCase(25, 100, highTruncation, 50, 1.5777218104420236e-30)]
[TestCase(100, 25, highTruncation, 150, 6.6003546737276816e-6)]
[TestCase(100, 3.5, 10000, 1000, 1.106797291738662e-06)]
[TestCase(100, 3.5, 10000, 2000, 4.891399189920708e-08)]
[TestCase(1000, 5.5, 100000, 2000, 6.076698900882660e-05)]
public void ValidateDensity(double scale, double shape, double truncation, double x, double p)
{
var n = new TruncatedPareto(scale, shape, truncation);
AssertHelpers.AlmostEqualRelative(p, n.Density(x), lowPrecision);
AssertHelpers.AlmostEqualRelative(p, TruncatedPareto.PDF(scale, shape, truncation, x), lowPrecision);
}
/// <summary>
/// Validate density log.
/// </summary>
/// <param name="scale">Scale parameter.</param>
/// <param name="shape">Shape parameter.</param>
/// <param name="truncation">Truncation parameter.</param>
/// <param name="x">The location at which to compute the log density.</param>
/// <param name="p">Expected value.</param>
[TestCase(3, 2, highTruncation, 3, -0.405465108108164)]
[TestCase(100, 3.5, 10000, 1000, -13.714040035965924)]
[TestCase(100, 3.5, 10000, 2000, -16.833202348485678)]
public void ValidateDensityLn(double scale, double shape, double truncation, double x, double p)
{
var n = new TruncatedPareto(scale, shape, truncation);
AssertHelpers.AlmostEqualRelative(p, n.DensityLn(x), lowPrecision);
AssertHelpers.AlmostEqualRelative(p, TruncatedPareto.PDFLn(scale, shape, truncation, x), lowPrecision);
}
/// <summary>
/// Validate cumulative distribution.
/// </summary>
/// <param name="scale">Scale parameter.</param>
/// <param name="shape">Shape parameter.</param>
/// <param name="truncation">Truncation parameter.</param>
/// <param name="x">The location at which to compute the cumulative distribution function.</param>
/// <param name="expected">Expected value.</param>
[TestCase(1.0, 1.0, highTruncation, 1.0, 0)]
[TestCase(7.0, 7.0, highTruncation, 10.0, 0.9176457)]
[TestCase(10.0, 10.0, highTruncation, 12.0, 0.83849441711015427)]
[TestCase(100, 3.5, 10000, 102, 0.066961862452387)]
[TestCase(100, 3.5, 10000, 1000, 0.999683872202370)]
[TestCase(100, 3.5, 10000, 2000, 0.999972149147496)]
[TestCase(1000, 5.5, 100000, 2000, 0.977902913097699)]
public void ValidateCumulativeDistribution(double scale, double shape, double truncation, double x, double expected)
{
var n = new TruncatedPareto(scale, shape, truncation);
AssertHelpers.AlmostEqualRelative(expected, n.CumulativeDistribution(x), highPrecision);
AssertHelpers.AlmostEqualRelative(expected, TruncatedPareto.CDF(scale, shape, truncation, x), highPrecision);
}
/// <summary>
/// Validate inverse cumulative distribution.
/// </summary>
/// <param name="scale">Scale parameter.</param>
/// <param name="shape">Shape parameter.</param>
/// <param name="truncation">Truncation parameter.</param>
/// <param name="p">The location at which to compute the inverse cumulative distribution function.</param>
/// <param name="expected">Expected value.</param>
[TestCase(1.0, 1.0, highTruncation, 0, 1.0)]
[TestCase(7.0, 7.0, highTruncation, 0.9176457, 10.0)]
[TestCase(10.0, 10.0, highTruncation, 0.83849441711015427, 12.0)]
[TestCase(100, 3.5, 10000, 0.066961862452387, 102)]
[TestCase(100, 3.5, 10000, 0.999683872202370, 1000)]
[TestCase(100, 3.5, 10000, 0.999972149147496, 2000)]
[TestCase(1000, 5.5, 100000, 0.977902913097699, 2000)]
public void ValidateInverseCumulativeDistribution(double scale, double shape, double truncation, double p, double expected)
{
var n = new TruncatedPareto(scale, shape, truncation);
AssertHelpers.AlmostEqualRelative(expected, n.InvCDF(p), lowPrecision);
AssertHelpers.AlmostEqualRelative(expected, TruncatedPareto.ICDF(scale, shape, truncation, p), lowPrecision);
}
/// <summary>
/// Can sample static.
/// </summary>
[Test]
public void CanSampleStatic()
{
TruncatedPareto.Sample(new Numerics.Random.MersenneTwister(100), 10.0, 10.0, 1000.0);
}
/// <summary>
/// Can fill array with samples static.
/// </summary>
[Test]
public void CanFillSampleArrayStatic()
{
double[] samples = new double[100];
TruncatedPareto.Samples(new Numerics.Random.MersenneTwister(100), samples, 10.0, 10.0, 1000.0);
Assert.IsTrue(!samples.Any(x => x == 0));
}
/// <summary>
/// Can sample sequence static.
/// </summary>
[Test]
public void CanSampleSequenceStatic()
{
var ied = TruncatedPareto.Samples(new Numerics.Random.MersenneTwister(100), 10.0, 10.0, 1000.0);
GC.KeepAlive(ied.Take(5).ToArray());
}
/// <summary>
/// Fail sample static with bad parameters.
/// </summary>
[Test]
public void FailSampleStatic()
{
Assert.That(() => { var d = TruncatedPareto.Sample(new Numerics.Random.MersenneTwister(100), 10.0, 10.0, 5.0); }, Throws.ArgumentException);
}
/// <summary>
/// Fail filling array with samples with bad parameters static.
/// </summary>
[Test]
public void FailFillingSampleArrayStatic()
{
double[] samples = new double[100];
Assert.That(() => { TruncatedPareto.Samples(new Numerics.Random.MersenneTwister(100), samples, 10.0, 10.0, 5.0); }, Throws.ArgumentException);
}
/// <summary>
/// Fail sample sequence static with bad parameters.
/// </summary>
[Test]
public void FailSampleSequenceStatic()
{
Assert.That(() => { var ied = TruncatedPareto.Samples(new Numerics.Random.MersenneTwister(100), 10.0, 10.0, 5.0).First(); }, Throws.ArgumentException);
}
/// <summary>
/// Can sample.
/// </summary>
[Test]
public void CanSample()
{
var n = new TruncatedPareto(10.0, 10.0, 1000.0);
n.Sample();
}
/// <summary>
/// Can fill array with samples.
/// </summary>
[Test]
public void CanFillSampleArray()
{
double[] samples = new double[100];
var n = new TruncatedPareto(10.0, 10.0, 1000.0);
n.Samples(samples);
Assert.IsTrue(!samples.Any(x => x == 0));
}
/// <summary>
/// Can sample sequence.
/// </summary>
[Test]
public void CanSampleSequence()
{
var n = new TruncatedPareto(10.0, 10.0, 1000.0);
var ied = n.Samples();
GC.KeepAlive(ied.Take(5).ToArray());
}
}
}

36
src/Numerics.Tests/OptimizationTests/BfgsBMinimizerTests.cs

@ -121,6 +121,42 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3));
}
[Test]
public void FindMinimum_Quadratic()
{
var obj = ObjectiveFunction.Gradient(
x => x[0] * x[0] + x[1] * x[1],
x => new DenseVector(new[] {2 * x[0], 2 * x[1]})
);
var solver = new BfgsBMinimizer(1e-5, 1e-5, 1e-5, maximumIterations: 1000);
var lowerBound = new DenseVector(new[] {-1.0, -1.0});
var upperBound = new DenseVector(new[] {2.0, 2.0});
var initialGuess = new DenseVector(new[] {1.5, 1.5});
var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess);
Assert.That(Math.Abs(result.MinimizingPoint[0] - 0.0), Is.LessThan(1e-3));
Assert.That(Math.Abs(result.MinimizingPoint[1] - 0.0), Is.LessThan(1e-3));
}
[Test]
public void FindMinimum_Quadratic_TwoBoundaries()
{
var obj = ObjectiveFunction.Gradient(
x => x[0] * x[0] + x[1] * x[1],
x => new DenseVector(new[] {2 * x[0], 2 * x[1]})
);
var solver = new BfgsBMinimizer(1e-5, 1e-5, 1e-5, maximumIterations: 1000);
var lowerBound = new DenseVector(new[] {1.0, 1.0});
var upperBound = new DenseVector(new[] {2.0, 2.0});
var initialGuess = new DenseVector(new[] {1.5, 1.5});
var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess);
Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3));
Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3));
}
[Test]
public void FindMinimum_Rosenbrock_MinimumGreateerOrEqualToLowerBoundary()
{

429
src/Numerics/Distributions/Burr.cs

@ -0,0 +1,429 @@
// <copyright file="Burr.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
//
// Copyright (c) 2009-2019 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using MathNet.Numerics.Properties;
using MathNet.Numerics.Random;
using System;
using System.Collections.Generic;
namespace MathNet.Numerics.Distributions
{
public class Burr : IContinuousDistribution
{
private System.Random _random;
/// <summary>
/// Gets the scale (a) of the distribution. Range: a > 0.
/// </summary>
public double a { get; }
/// <summary>
/// Gets the first shape parameter (c) of the distribution. Range: c > 0.
/// </summary>
public double c { get; }
/// <summary>
/// Gets the second shape parameter (k) of the distribution. Range: k > 0.
/// </summary>
public double k { get; }
/// <summary>
/// Initializes a new instance of the Burr Type XII class.
/// </summary>
/// <param name="a">The scale parameter a of the Burr distribution. Range: a > 0.</param>
/// <param name="c">The first shape parameter c of the Burr distribution. Range: c > 0.</param>
/// <param name="k">The second shape parameter k of the Burr distribution. Range: k > 0.</param>
/// <param name="randomSource">The random number generator which is used to draw random samples. Optional, can be null.</param>
public Burr(double a, double c, double k, System.Random randomSource = null)
{
if (!IsValidParameterSet(a, c, k))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
_random = randomSource ?? SystemRandomSource.Default;
this.a = a;
this.c = c;
this.k = k;
}
/// <summary>
/// A string representation of the distribution.
/// </summary>
/// <returns>a string representation of the distribution.</returns>
public override string ToString()
{
return "Burr(a = " + a + ", c = " + c + ", k = " + k + ")";
}
/// <summary>
/// Tests whether the provided values are valid parameters for this distribution.
/// </summary>
/// <param name="a">The scale parameter a of the Burr distribution. Range: a > 0.</param>
/// <param name="c">The first shape parameter c of the Burr distribution. Range: c > 0.</param>
/// <param name="k">The second shape parameter k of the Burr distribution. Range: k > 0.</param>
public static bool IsValidParameterSet(double a, double c, double k)
{
var allFinite = a.IsFinite() && c.IsFinite() && k.IsFinite();
return allFinite && a > 0.0 && c > 0.0 && k > 0.0;
}
/// <summary>
/// Gets the random number generator which is used to draw random samples.
/// </summary>
public System.Random RandomSource
{
get { return _random; }
set { _random = value ?? SystemRandomSource.Default; }
}
/// <summary>
/// Gets the mean of the Burr distribution.
/// </summary>
public double Mean
{
get
{
return (1 / SpecialFunctions.Gamma(k)) * a * SpecialFunctions.Gamma(1 + 1 / c) * SpecialFunctions.Gamma(k - 1 / c);
}
}
/// <summary>
/// Gets the variance of the Burr distribution.
/// </summary>
public double Variance
{
get
{
return (1 / SpecialFunctions.Gamma(k)) * Math.Pow(a, 2) * SpecialFunctions.Gamma(1 + 2 / c) * SpecialFunctions.Gamma(k - 2 / c)
- Math.Pow((1 / SpecialFunctions.Gamma(k)) * a * SpecialFunctions.Gamma(1 + 1 / c) * SpecialFunctions.Gamma(k - 1 / c), 2);
}
}
/// <summary>
/// Gets the standard deviation of the Burr distribution.
/// </summary>
public double StdDev
{
get
{
return Math.Sqrt(Variance);
}
}
/// <summary>
/// Gets the mode of the Burr distribution.
/// </summary>
public double Mode
{
get
{
return a * Math.Pow((c - 1) / (c * k + 1), 1 / c);
}
}
/// <summary>
/// Gets the minimum of the Burr distribution.
/// </summary>
public double Minimum
{
get { return 0.0; }
}
/// <summary>
/// Gets the maximum of the Burr distribution.
/// </summary>
public double Maximum
{
get { return double.PositiveInfinity; }
}
/// <summary>
/// Gets the entropy of the Burr distribution (currently not supported).
/// </summary>
public double Entropy
{
get { throw new NotSupportedException(); }
}
/// <summary>
/// Gets the skewness of the Burr distribution.
/// </summary>
public double Skewness
{
get
{
var mean = Mean;
var variance = Variance;
var std = StdDev;
return (GetMoment(3) - 3 * mean * variance - mean * mean * mean) / (std * std * std);
}
}
/// <summary>
/// Gets the median of the Burr distribution.
/// </summary>
public double Median
{
get
{
return a * Math.Pow(Math.Pow(2, 1 / k) - 1, 1 / c);
}
}
/// <summary>
/// Generates a sample from the Burr distribution.
/// </summary>
/// <returns>a sample from the distribution.</returns>
public double Sample()
{
return SampleUnchecked(_random, a, c, k);
}
/// <summary>
/// Fills an array with samples generated from the distribution.
/// </summary>
/// <param name="values">The array to fill with the samples.</param>
public void Samples(double[] values)
{
SamplesUnchecked(_random, values, a, c, k);
}
/// <summary>
/// Generates a sequence of samples from the Burr distribution.
/// </summary>
/// <returns>a sequence of samples from the distribution.</returns>
public IEnumerable<double> Samples()
{
return SamplesUnchecked(_random, a, c, k);
}
/// <summary>
/// Generates a sample from the Burr distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="a">The scale parameter a of the Burr distribution. Range: a > 0.</param>
/// <param name="c">The first shape parameter c of the Burr distribution. Range: c > 0.</param>
/// <param name="k">The second shape parameter k of the Burr distribution. Range: k > 0.</param>
/// <returns>a sample from the distribution.</returns>
public static double Sample(System.Random rnd, double a, double c, double k)
{
if (!IsValidParameterSet(a, c, k))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return SampleUnchecked(rnd, a, c, k);
}
/// <summary>
/// Fills an array with samples generated from the distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="values">The array to fill with the samples.</param>
/// <param name="a">The scale parameter a of the Burr distribution. Range: a > 0.</param>
/// <param name="c">The first shape parameter c of the Burr distribution. Range: c > 0.</param>
/// <param name="k">The second shape parameter k of the Burr distribution. Range: k > 0.</param>
public static void Samples(System.Random rnd, double[] values, double a, double c, double k)
{
if (!IsValidParameterSet(a, c, k))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
SamplesUnchecked(rnd, values, a, c, k);
}
/// <summary>
/// Generates a sequence of samples from the Burr distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="a">The scale parameter a of the Burr distribution. Range: a > 0.</param>
/// <param name="c">The first shape parameter c of the Burr distribution. Range: c > 0.</param>
/// <param name="k">The second shape parameter k of the Burr distribution. Range: k > 0.</param>
/// <returns>a sequence of samples from the distribution.</returns>
public static IEnumerable<double> Samples(System.Random rnd, double a, double c, double k)
{
if (!IsValidParameterSet(a, c, k))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return SamplesUnchecked(rnd, a, c, k);
}
internal static double SampleUnchecked(System.Random rnd, double a, double c, double k)
{
var k_inv = 1 / k;
var c_inv = 1 / c;
double u = rnd.NextDouble();
return a * Math.Pow(Math.Pow(1 - u, -k_inv) - 1, c_inv);
}
internal static void SamplesUnchecked(System.Random rnd, double[] values, double a, double c, double k)
{
if (values.Length == 0)
{
return;
}
var k_inv = 1 / k;
var c_inv = 1 / c;
double[] u = rnd.NextDoubles(values.Length);
for (var j = 0; j < values.Length; ++j)
{
values[j] = a * Math.Pow(Math.Pow(1 - u[j], -k_inv) - 1, c_inv);
}
}
internal static IEnumerable<double> SamplesUnchecked(System.Random rnd, double a, double c, double k)
{
while (true)
{
yield return SampleUnchecked(rnd, a, c, k);
}
}
/// <summary>
/// Gets the n-th raw moment of the distribution.
/// </summary>
/// <param name="n">The order (n) of the moment. Range: n ≥ 1.</param>
/// <returns>the n-th moment of the distribution.</returns>
public double GetMoment(double n)
{
if (n > k * c)
{
throw new ArgumentException(Resources.ArgumentParameterSetInvalid);
}
var lambda_n = (n / c) * SpecialFunctions.Gamma(n / c) * SpecialFunctions.Gamma(k - n / c);
return Math.Pow(a, n) * lambda_n / SpecialFunctions.Gamma(k);
}
/// <summary>
/// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
/// </summary>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
/// <seealso cref="PDF"/>
public double Density(double x)
{
return DensityImpl(a, c, k, x);
}
/// <summary>
/// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x).
/// </summary>
/// <param name="x">The location at which to compute the log density.</param>
/// <returns>the log density at <paramref name="x"/>.</returns>
/// <seealso cref="PDFLn"/>
public double DensityLn(double x)
{
return DensityLnImpl(a, c, k, x);
}
/// <summary>
/// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x).
/// </summary>
/// <param name="x">The location at which to compute the cumulative distribution function.</param>
/// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
/// <seealso cref="CDF"/>
public double CumulativeDistribution(double x)
{
return CumulativeDistributionImpl(a, c, k, x);
}
/// <summary>
/// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
/// </summary>
/// <param name="a">The scale parameter a of the Burr distribution. Range: a > 0.</param>
/// <param name="c">The first shape parameter c of the Burr distribution. Range: c > 0.</param>
/// <param name="k">The second shape parameter k of the Burr distribution. Range: k > 0.</param>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
/// <seealso cref="Density"/>
public static double PDF(double a, double c, double k, double x)
{
if (!IsValidParameterSet(a, c, k))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return DensityImpl(a, c, k, x);
}
/// <summary>
/// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x).
/// </summary>
/// <param name="a">The scale parameter a of the Burr distribution. Range: a > 0.</param>
/// <param name="c">The first shape parameter c of the Burr distribution. Range: c > 0.</param>
/// <param name="k">The second shape parameter k of the Burr distribution. Range: k > 0.</param>
/// <param name="x">The location at which to compute the log density.</param>
/// <returns>the log density at <paramref name="x"/>.</returns>
/// <seealso cref="DensityLn"/>
public static double PDFLn(double a, double c, double k, double x)
{
if (!IsValidParameterSet(a, c, k))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return DensityLnImpl(a, c, k, x);
}
/// <summary>
/// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x).
/// </summary>
/// <param name="a">The scale parameter a of the Burr distribution. Range: a > 0.</param>
/// <param name="c">The first shape parameter c of the Burr distribution. Range: c > 0.</param>
/// <param name="k">The second shape parameter k of the Burr distribution. Range: k > 0.</param>
/// <param name="x">The location at which to compute the cumulative distribution function.</param>
/// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
/// <seealso cref="CumulativeDistribution"/>
public static double CDF(double a, double c, double k, double x)
{
if (!IsValidParameterSet(a, c, k))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return CumulativeDistributionImpl(a, c, k, x);
}
internal static double DensityImpl(double a, double c, double k, double x)
{
var numerator = (k * c / a) * Math.Pow(x / a, c - 1);
var denominator = Math.Pow(1 + Math.Pow(x / a, c), k + 1);
return numerator / denominator;
}
internal static double DensityLnImpl(double a, double c, double k, double x)
{
return Math.Log(DensityImpl(a, c, k, x));
}
internal static double CumulativeDistributionImpl(double a, double c, double k, double x)
{
var denominator = Math.Pow(1 + Math.Pow(x / a, c), k);
return 1 - 1 / denominator;
}
}
}

449
src/Numerics/Distributions/InverseGaussian.cs

@ -0,0 +1,449 @@
// <copyright file="InverseGaussian.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
//
// Copyright (c) 2009-2019 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using MathNet.Numerics.Properties;
using MathNet.Numerics.Random;
using MathNet.Numerics.Statistics;
using System;
using System.Collections.Generic;
using System.Linq;
namespace MathNet.Numerics.Distributions
{
public class InverseGaussian : IContinuousDistribution
{
private System.Random _random;
/// <summary>
/// Gets the mean (μ) of the distribution. Range: μ > 0.
/// </summary>
public double Mu { get; }
/// <summary>
/// Gets the shape (λ) of the distribution. Range: λ > 0.
/// </summary>
public double Lambda { get; }
/// <summary>
/// Initializes a new instance of the InverseGaussian class.
/// </summary>
/// <param name="mu">The mean (μ) of the distribution. Range: μ > 0.</param>
/// <param name="lambda">The shape (λ) of the distribution. Range: λ > 0.</param>
/// <param name="randomSource">The random number generator which is used to draw random samples.</param>
public InverseGaussian(double mu, double lambda, System.Random randomSource = null)
{
if (!IsValidParameterSet(mu, lambda))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
_random = randomSource ?? SystemRandomSource.Default;
Mu = mu;
Lambda = lambda;
}
/// <summary>
/// A string representation of the distribution.
/// </summary>
/// <returns>a string representation of the distribution.</returns>
public override string ToString()
{
return "InverseGaussian(μ = " + Mu + ", λ = " + Lambda + ")";
}
/// <summary>
/// Tests whether the provided values are valid parameters for this distribution.
/// </summary>
/// <param name="mu">The mean (μ) of the distribution. Range: μ > 0.</param>
/// <param name="lambda">The shape (λ) of the distribution. Range: λ > 0.</param>
public static bool IsValidParameterSet(double mu, double lambda)
{
var allFinite = mu.IsFinite() && lambda.IsFinite();
return allFinite && mu > 0.0 && lambda > 0.0;
}
/// <summary>
/// Gets the random number generator which is used to draw random samples.
/// </summary>
public System.Random RandomSource
{
get { return _random; }
set { _random = value ?? SystemRandomSource.Default; }
}
/// <summary>
/// Gets the mean of the Inverse Gaussian distribution.
/// </summary>
public double Mean
{
get
{
return Mu;
}
}
/// <summary>
/// Gets the variance of the Inverse Gaussian distribution.
/// </summary>
public double Variance
{
get
{
return Math.Pow(Mu, 3) / Lambda;
}
}
/// <summary>
/// Gets the standard deviation of the Inverse Gaussian distribution.
/// </summary>
public double StdDev
{
get
{
return Math.Sqrt(Variance);
}
}
/// <summary>
/// Gets the median of the Inverse Gaussian distribution.
/// No closed form analytical expression exists, so this value is approximated numerically and can throw an exception.
/// </summary>
public double Median
{
get { return InvCDF(0.5); }
}
/// <summary>
/// Gets the minimum of the Inverse Gaussian distribution.
/// </summary>
public double Minimum
{
get { return 0.0; }
}
/// <summary>
/// Gets the maximum of the Inverse Gaussian distribution.
/// </summary>
public double Maximum
{
get { return double.PositiveInfinity; }
}
/// <summary>
/// Gets the skewness of the Inverse Gaussian distribution.
/// </summary>
public double Skewness
{
get { return 3 * Math.Sqrt(Mu / Lambda); }
}
/// <summary>
/// Gets the kurtosis of the Inverse Gaussian distribution.
/// </summary>
public double Kurtosis
{
get { return 15 * Mu / Lambda; }
}
/// <summary>
/// Gets the mode of the Inverse Gaussian distribution.
/// </summary>
public double Mode
{
get { return Mu * (Math.Sqrt(1 + (9 * Mu * Mu) / (4 * Lambda * Lambda)) - (3 * Mu) / (2 * Lambda)); }
}
/// <summary>
/// Gets the entropy of the Inverse Gaussian distribution (currently not supported).
/// </summary>
public double Entropy
{
get { throw new NotSupportedException(); }
}
/// <summary>
/// Generates a sample from the inverse Gaussian distribution.
/// </summary>
/// <returns>a sample from the distribution.</returns>
public double Sample()
{
return SampleUnchecked(_random, Mu, Lambda);
}
/// <summary>
/// Fills an array with samples generated from the distribution.
/// </summary>
/// <param name="values">The array to fill with the samples.</param>
public void Samples(double[] values)
{
SamplesUnchecked(_random, values, Mu, Lambda);
}
/// <summary>
/// Generates a sequence of samples from the inverse Gaussian distribution.
/// </summary>
/// <returns>a sequence of samples from the distribution.</returns>
public IEnumerable<double> Samples()
{
return SamplesUnchecked(_random, Mu, Lambda);
}
/// <summary>
/// Generates a sample from the inverse Gaussian distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="mu">The mean (μ) of the distribution. Range: μ > 0.</param>
/// <param name="lambda">The shape (λ) of the distribution. Range: λ > 0.</param>
/// <returns>a sample from the distribution.</returns>
public static double Sample(System.Random rnd, double mu, double lambda)
{
if (!IsValidParameterSet(mu, lambda))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return SampleUnchecked(rnd, mu, lambda);
}
/// <summary>
/// Fills an array with samples generated from the distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="values">The array to fill with the samples.</param>
/// <param name="mu">The mean (μ) of the distribution. Range: μ > 0.</param>
/// <param name="lambda">The shape (λ) of the distribution. Range: λ > 0.</param>
public static void Samples(System.Random rnd, double[] values, double mu, double lambda)
{
if (!IsValidParameterSet(mu, lambda))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
SamplesUnchecked(rnd, values, mu, lambda);
}
/// <summary>
/// Generates a sequence of samples from the Burr distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="mu">The mean (μ) of the distribution. Range: μ > 0.</param>
/// <param name="lambda">The shape (λ) of the distribution. Range: λ > 0.</param>
/// <returns>a sequence of samples from the distribution.</returns>
public static IEnumerable<double> Samples(System.Random rnd, double mu, double lambda)
{
if (!IsValidParameterSet(mu, lambda))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return SamplesUnchecked(rnd, mu, lambda);
}
internal static double SampleUnchecked(System.Random rnd, double mu, double lambda)
{
double v = MathNet.Numerics.Distributions.Normal.Sample(rnd, 0, 1);
double test = rnd.NextDouble();
return InverseGaussianSampleImpl(mu, lambda, v, test);
}
internal static void SamplesUnchecked(System.Random rnd, double[] values, double mu, double lambda)
{
if (values.Length == 0)
{
return;
}
double[] v = new double[values.Length];
MathNet.Numerics.Distributions.Normal.Samples(rnd, v, 0, 1);
double[] test = rnd.NextDoubles(values.Length);
for (var j = 0; j < values.Length; ++j)
{
values[j] = InverseGaussianSampleImpl(mu, lambda, v[j], test[j]);
}
}
internal static IEnumerable<double> SamplesUnchecked(System.Random rnd, double mu, double lambda)
{
while (true)
{
yield return SampleUnchecked(rnd, mu, lambda);
}
}
internal static double InverseGaussianSampleImpl(double mu, double lambda, double normalSample, double uniformSample)
{
double y = normalSample * normalSample;
double x = mu + (mu * mu * y) / (2 * lambda) - (mu / (2 * lambda)) * Math.Sqrt(4 * mu * lambda * y + mu * mu * y * y);
if (uniformSample <= mu / (mu + x))
return x;
else
return mu * mu / x;
}
/// <summary>
/// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
/// </summary>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
/// <seealso cref="PDF"/>
public double Density(double x)
{
return DensityImpl(Mu, Lambda, x);
}
/// <summary>
/// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x).
/// </summary>
/// <param name="x">The location at which to compute the log density.</param>
/// <returns>the log density at <paramref name="x"/>.</returns>
/// <seealso cref="PDFLn"/>
public double DensityLn(double x)
{
return DensityLnImpl(Mu, Lambda, x);
}
/// <summary>
/// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x).
/// </summary>
/// <param name="x">The location at which to compute the cumulative distribution function.</param>
/// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
/// <seealso cref="CDF"/>
public double CumulativeDistribution(double x)
{
return CumulativeDistributionImpl(Mu, Lambda, x);
}
/// <summary>
/// Computes the inverse cumulative distribution (CDF) of the distribution at p, i.e. solving for P(X ≤ x) = p.
/// </summary>
/// <param name="p">The location at which to compute the inverse cumulative distribution function.</param>
/// <returns>the inverse cumulative distribution at location <paramref name="p"/>.</returns>
public double InvCDF(double p)
{
Func<double, double> equationToSolve = (x) => CumulativeDistribution(x) - p;
if (RootFinding.NewtonRaphson.TryFindRoot(equationToSolve, Density, Mode, 0, double.PositiveInfinity, 1e-8, 100, out double quantile))
return quantile;
else
throw new NonConvergenceException(Resources.NumericalEstimationFailed);
}
/// <summary>
/// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
/// </summary>
/// <param name="mu">The mean (μ) of the distribution. Range: μ > 0.</param>
/// <param name="lambda">The shape (λ) of the distribution. Range: λ > 0.</param>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
/// <seealso cref="Density"/>
public static double PDF(double mu, double lambda, double x)
{
if (!IsValidParameterSet(mu, lambda))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return DensityImpl(mu, lambda, x);
}
/// <summary>
/// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x).
/// </summary>
/// <param name="mu">The mean (μ) of the distribution. Range: μ > 0.</param>
/// <param name="lambda">The shape (λ) of the distribution. Range: λ > 0.</param>
/// <param name="x">The location at which to compute the log density.</param>
/// <returns>the log density at <paramref name="x"/>.</returns>
/// <seealso cref="DensityLn"/>
public static double PDFLn(double mu, double lambda, double x)
{
if (!IsValidParameterSet(mu, lambda))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return DensityLnImpl(mu, lambda, x);
}
/// <summary>
/// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x).
/// </summary>
/// <param name="mu">The mean (μ) of the distribution. Range: μ > 0.</param>
/// <param name="lambda">The shape (λ) of the distribution. Range: λ > 0.</param>
/// <param name="x">The location at which to compute the cumulative distribution function.</param>
/// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
/// <seealso cref="CumulativeDistribution"/>
public static double CDF(double mu, double lambda, double x)
{
if (!IsValidParameterSet(mu, lambda))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return CumulativeDistributionImpl(mu, lambda, x);
}
/// <summary>
/// Computes the inverse cumulative distribution (CDF) of the distribution at p, i.e. solving for P(X ≤ x) = p.
/// </summary>
/// <param name="mu">The mean (μ) of the distribution. Range: μ > 0.</param>
/// <param name="lambda">The shape (λ) of the distribution. Range: λ > 0.</param>
/// <param name="p">The location at which to compute the inverse cumulative distribution function.</param>
/// <returns>the inverse cumulative distribution at location <paramref name="p"/>.</returns>
/// <seealso cref="CumulativeDistribution"/>
public static double ICDF(double mu, double lambda, double p)
{
if (!IsValidParameterSet(mu, lambda))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
var igDist = new InverseGaussian(mu, lambda);
return igDist.InvCDF(p);
}
/// <summary>
/// Estimates the Inverse Gaussian parameters from sample data with maximum-likelihood.
/// </summary>
/// <param name="samples">The samples to estimate the distribution parameters from.</param>
/// <param name="randomSource">The random number generator which is used to draw random samples. Optional, can be null.</param>
/// <returns>An Inverse Gaussian distribution.</returns>
public static InverseGaussian Estimate(IEnumerable<double> samples, System.Random randomSource = null)
{
var sampleVec = samples.ToArray();
var muHat = sampleVec.Mean();
var lambdahat = 1 / (1 / samples.HarmonicMean() - 1 / muHat);
return new InverseGaussian(muHat, lambdahat, randomSource);
}
internal static double DensityImpl(double mu, double lambda, double x)
{
return Math.Sqrt(lambda / (2 * Math.PI * Math.Pow(x, 3))) * Math.Exp(-((lambda * Math.Pow(x - mu, 2)) / (2 * mu * mu * x)));
}
internal static double DensityLnImpl(double mu, double lambda, double x)
{
return Math.Log(Math.Sqrt(lambda / (2 * Math.PI * Math.Pow(x, 3)))) - ((lambda * Math.Pow(x - mu, 2)) / (2 * mu * mu * x));
}
internal static double CumulativeDistributionImpl(double mu, double lambda, double x)
{
return Normal.CDF(0, 1, Math.Sqrt(lambda / x) * (x / mu - 1)) + Math.Exp(2 * lambda / mu) * Normal.CDF(0, 1, -Math.Sqrt(lambda / x) * (x / mu + 1));
}
}
}

466
src/Numerics/Distributions/TruncatedPareto.cs

@ -0,0 +1,466 @@
// <copyright file="TruncatedPareto.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
//
// Copyright (c) 2009-2019 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using MathNet.Numerics.Properties;
using MathNet.Numerics.Random;
using System;
using System.Collections.Generic;
namespace MathNet.Numerics.Distributions
{
public class TruncatedPareto : IContinuousDistribution
{
private System.Random _random;
/// <summary>
/// Initializes a new instance of the TruncatedPareto class.
/// </summary>
/// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
/// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
/// <param name="truncation">The truncation (T) of the distribution. Range: T > xm.</param>
/// <param name="randomSource">The random number generator which is used to draw random samples.</param>
/// <exception cref="ArgumentException">If <paramref name="scale"/> or <paramref name="shape"/> are non-positive or if T ≤ xm.</exception>
public TruncatedPareto(double scale, double shape, double truncation, System.Random randomSource = null)
{
if (!IsValidParameterSet(scale, shape, truncation))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
_random = randomSource ?? SystemRandomSource.Default;
Scale = scale;
Shape = shape;
Truncation = truncation;
}
/// <summary>
/// A string representation of the distribution.
/// </summary>
/// <returns>a string representation of the distribution.</returns>
public override string ToString()
{
return "Truncated Pareto(Scale = " + Scale + ", Shape = " + Shape + ", Truncation = " + Truncation + ")";
}
/// <summary>
/// Tests whether the provided values are valid parameters for this distribution.
/// </summary>
/// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
/// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
/// <param name="truncation">The truncation (T) of the distribution. Range: T > xm.</param>
public static bool IsValidParameterSet(double scale, double shape, double truncation)
{
var allFinite = scale.IsFinite() && shape.IsFinite() && truncation.IsFinite();
return allFinite && scale > 0.0 && shape > 0.0 && truncation > scale;
}
/// <summary>
/// Gets the random number generator which is used to draw random samples.
/// </summary>
public System.Random RandomSource
{
get { return _random; }
set { _random = value ?? SystemRandomSource.Default; }
}
/// <summary>
/// Gets the scale (xm) of the distribution. Range: xm > 0.
/// </summary>
public double Scale { get; }
/// <summary>
/// Gets the shape (α) of the distribution. Range: α > 0.
/// </summary>
public double Shape { get; }
/// <summary>
/// Gets the truncation (T) of the distribution. Range: T > 0.
/// </summary>
public double Truncation { get; }
/// <summary>
/// Gets the n-th raw moment of the distribution.
/// </summary>
/// <param name="n">The order (n) of the moment. Range: n ≥ 1.</param>
/// <returns>the n-th moment of the distribution.</returns>
public double GetMoment(int n)
{
double moment;
if (Shape.AlmostEqual(n))
{
moment = ((Shape * Math.Pow(Scale, n)) / (1 - Math.Pow(Scale / Truncation, Shape))) * Math.Log(Truncation / Scale);
}
else
{
moment = ((Shape * Math.Pow(Scale, n)) / (Shape - n)) * ((1 - Math.Pow((Scale / Truncation), (Shape - n))) / (1 - Math.Pow(Scale / Truncation, Shape)));
}
return moment;
}
/// <summary>
/// Gets the mean of the truncated Pareto distribution.
/// </summary>
public double Mean
{
get
{
return GetMoment(1);
}
}
/// <summary>
/// Gets the variance of the truncated Pareto distribution.
/// </summary>
public double Variance
{
get
{
return GetMoment(2) - Math.Pow(GetMoment(1), 2);
}
}
/// <summary>
/// Gets the standard deviation of the truncated Pareto distribution.
/// </summary>
public double StdDev
{
get
{
return Math.Sqrt(Variance);
}
}
/// <summary>
/// Gets the mode of the truncated Pareto distribution (not supported).
/// </summary>
public double Mode
{
get { throw new NotSupportedException(); }
}
/// <summary>
/// Gets the minimum of the truncated Pareto distribution.
/// </summary>
public double Minimum
{
get { return Scale; }
}
/// <summary>
/// Gets the maximum of the truncated Pareto distribution.
/// </summary>
public double Maximum
{
get { return Truncation; }
}
/// <summary>
/// Gets the entropy of the truncated Pareto distribution (not supported).
/// </summary>
public double Entropy
{
get { throw new NotSupportedException(); }
}
/// <summary>
/// Gets the skewness of the truncated Pareto distribution.
/// </summary>
public double Skewness
{
get
{
var mean = Mean;
var variance = Variance;
var std = StdDev;
return (GetMoment(3) - 3.0 * mean * variance - mean * mean * mean) / (std * std * std);
}
}
/// <summary>
/// Gets the median of the truncated Pareto distribution.
/// </summary>
public double Median
{
get
{
return Scale * Math.Pow(1.0 - (1.0 / 2.0) * (1.0 - Math.Pow(Scale / Truncation, Shape)), -(1.0 / Shape));
}
}
/// <summary>
/// Generates a sample from the truncated Pareto distribution.
/// </summary>
/// <returns>a sample from the distribution.</returns>
public double Sample()
{
return SampleUnchecked(_random, Scale, Shape, Truncation);
}
/// <summary>
/// Fills an array with samples generated from the distribution.
/// </summary>
/// <param name="values">The array to fill with the samples.</param>
public void Samples(double[] values)
{
SamplesUnchecked(_random, values, Scale, Shape, Truncation);
}
/// <summary>
/// Generates a sequence of samples from the truncated Pareto distribution.
/// </summary>
/// <returns>a sequence of samples from the distribution.</returns>
public IEnumerable<double> Samples()
{
return SamplesUnchecked(_random, Scale, Shape, Truncation);
}
/// <summary>
/// Generates a sample from the truncated Pareto distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
/// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
/// <param name="truncation">The truncation (T) of the distribution. Range: T > xm.</param>
/// <returns>a sample from the distribution.</returns>
public static double Sample(System.Random rnd, double scale, double shape, double truncation)
{
if (!IsValidParameterSet(scale, shape, truncation))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return SampleUnchecked(rnd, scale, shape, truncation);
}
/// <summary>
/// Fills an array with samples generated from the distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="values">The array to fill with the samples.</param>
/// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
/// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
/// <param name="truncation">The truncation (T) of the distribution. Range: T > xm.</param>
public static void Samples(System.Random rnd, double[] values, double scale, double shape, double truncation)
{
if (!IsValidParameterSet(scale, shape, truncation))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
SamplesUnchecked(rnd, values, scale, shape, truncation);
}
/// <summary>
/// Generates a sequence of samples from the truncated Pareto distribution.
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
/// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
/// <param name="truncation">The truncation (T) of the distribution. Range: T > xm.</param>
/// <returns>a sequence of samples from the distribution.</returns>
public static IEnumerable<double> Samples(System.Random rnd, double scale, double shape, double truncation)
{
if (!IsValidParameterSet(scale, shape, truncation))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return SamplesUnchecked(rnd, scale, shape, truncation);
}
internal static double SampleUnchecked(System.Random rnd, double scale, double shape, double truncation)
{
double uniform = rnd.NextDouble();
return InvCDFUncheckedImpl(scale, shape, truncation, uniform);
}
internal static void SamplesUnchecked(System.Random rnd, double[] values, double scale, double shape, double truncation)
{
if (values.Length == 0)
{
return;
}
double[] uniforms = rnd.NextDoubles(values.Length);
for (var j = 0; j < values.Length; ++j)
{
values[j] = InvCDFUncheckedImpl(scale, shape, truncation, uniforms[j]);
}
}
internal static IEnumerable<double> SamplesUnchecked(System.Random rnd, double scale, double shape, double truncation)
{
while (true)
{
yield return SampleUnchecked(rnd, scale, shape, truncation);
}
}
/// <summary>
/// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
/// </summary>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
/// <seealso cref="PDF"/>
public double Density(double x)
{
return DensityImpl(Scale, Shape, Truncation, x);
}
/// <summary>
/// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x).
/// </summary>
/// <param name="x">The location at which to compute the log density.</param>
/// <returns>the log density at <paramref name="x"/>.</returns>
/// <seealso cref="PDFLn"/>
public double DensityLn(double x)
{
return DensityLnImpl(Scale, Shape, Truncation, x);
}
/// <summary>
/// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x).
/// </summary>
/// <param name="x">The location at which to compute the cumulative distribution function.</param>
/// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
/// <seealso cref="CDF"/>
public double CumulativeDistribution(double x)
{
return CumulativeDistributionImpl(Scale, Shape, Truncation, x);
}
/// <summary>
/// Computes the inverse cumulative distribution (CDF) of the distribution at p, i.e. solving for P(X ≤ x) = p.
/// </summary>
/// <param name="p">The location at which to compute the inverse cumulative distribution function.</param>
/// <returns>the inverse cumulative distribution at location <paramref name="p"/>.</returns>
public double InvCDF(double p)
{
return InvCDFUncheckedImpl(Scale, Shape, Truncation, p);
}
/// <summary>
/// Computes the inverse cumulative distribution (CDF) of the distribution at p, i.e. solving for P(X ≤ x) = p.
/// </summary>
/// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
/// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
/// <param name="truncation">The truncation (T) of the distribution. Range: T > xm.</param>
/// <param name="p">The location at which to compute the inverse cumulative distribution function.</param>
/// <returns>the inverse cumulative distribution at location <paramref name="p"/>.</returns>
/// <seealso cref="CumulativeDistribution"/>
public static double ICDF(double scale, double shape, double truncation, double p)
{
if (!IsValidParameterSet(scale, shape, truncation))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return InvCDFUncheckedImpl(scale, shape, truncation, p);
}
/// <summary>
/// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
/// </summary>
/// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
/// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
/// <param name="truncation">The truncation (T) of the distribution. Range: T > xm.</param>
/// <param name="x">The location at which to compute the density.</param>
/// <returns>the density at <paramref name="x"/>.</returns>
/// <seealso cref="Density"/>
public static double PDF(double scale, double shape, double truncation, double x)
{
if (!IsValidParameterSet(scale, shape, truncation))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return DensityImpl(scale, shape, truncation, x);
}
/// <summary>
/// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x).
/// </summary>
/// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
/// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
/// <param name="truncation">The truncation (T) of the distribution. Range: T > xm.</param>
/// <param name="x">The location at which to compute the log density.</param>
/// <returns>the log density at <paramref name="x"/>.</returns>
/// <seealso cref="DensityLn"/>
public static double PDFLn(double scale, double shape, double truncation, double x)
{
if (!IsValidParameterSet(scale, shape, truncation))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return DensityLnImpl(scale, shape, truncation, x);
}
/// <summary>
/// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x).
/// </summary>
/// <param name="scale">The scale (xm) of the distribution. Range: xm > 0.</param>
/// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
/// <param name="truncation">The truncation (T) of the distribution. Range: T > xm.</param>
/// <param name="x">The location at which to compute the cumulative distribution function.</param>
/// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
/// <seealso cref="CumulativeDistribution"/>
public static double CDF(double scale, double shape, double truncation, double x)
{
if (!IsValidParameterSet(scale, shape, truncation))
{
throw new ArgumentException(Resources.InvalidDistributionParameters);
}
return CumulativeDistributionImpl(scale, shape, truncation, x);
}
internal static double DensityImpl(double scale, double shape, double truncation, double x)
{
if (x < scale || x > truncation)
return 0;
else
return (shape * Math.Pow(scale, shape) * Math.Pow(x, -shape - 1)) / (1 - Math.Pow(scale / truncation, shape));
}
internal static double DensityLnImpl(double scale, double shape, double truncation, double x)
{
return Math.Log(DensityImpl(scale, shape, truncation, x));
}
internal static double CumulativeDistributionImpl(double scale, double shape, double truncation, double x)
{
if (x <= scale)
return 0;
else if (x >= truncation)
return 1;
else
return (1 - Math.Pow(scale, shape) * Math.Pow(x, -shape)) / (1 - Math.Pow(scale / truncation, shape));
}
internal static double InvCDFUncheckedImpl(double scale, double shape, double truncation, double p)
{
var numerator = p * Math.Pow(truncation, shape) - p * Math.Pow(scale, shape) - Math.Pow(truncation, shape);
var denominator = Math.Pow(truncation, shape) * Math.Pow(scale, shape);
return Math.Pow(-numerator / denominator, -(1 / shape));
}
}
}

9
src/Numerics/Optimization/BfgsBMinimizer.cs

@ -74,7 +74,7 @@ namespace MathNet.Numerics.Optimization
ValidateGradientAndObjective(objective);
// Check that we're not already done
ExitCondition currentExitCondition = ExitCriteriaSatisfied(objective, null, 0);
var currentExitCondition = ExitCriteriaSatisfied(objective, null, 0);
if (currentExitCondition != ExitCondition.None)
return new MinimizationResult(objective, 0, currentExitCondition);
@ -140,6 +140,13 @@ namespace MathNet.Numerics.Optimization
var previousPoint = objective.Fork();
var candidatePoint = lineSearchResult.FunctionInfoAtMinimum;
ValidateGradientAndObjective(candidatePoint);
// Check that we're not done
currentExitCondition = ExitCriteriaSatisfied(candidatePoint, previousPoint, 0);
if (currentExitCondition != ExitCondition.None)
return new MinimizationResult(candidatePoint, 0, currentExitCondition);
var gradient = candidatePoint.Gradient;
var step = candidatePoint.Point - initialGuess;

9
src/Numerics/Precision.Comparison.cs

@ -674,5 +674,14 @@ namespace MathNet.Numerics
return CompareToNumbersBetween(a, b, maxNumbersBetween) < 0;
}
/// <summary>
/// Checks if a given double values is finite, i.e. neither NaN nor inifnity
/// </summary>
/// <param name="value">The value to be checked fo finitenes.</param>
public static bool IsFinite(this double value)
{
return !double.IsNaN(value) && !double.IsInfinity(value);
}
}
}

11
src/Numerics/Properties/Resources.Designer.cs

@ -853,6 +853,17 @@ namespace MathNet.Numerics.Properties
}
}
/// <summary>
/// Looks up a localized string similar to: Numerical estimation of the statistic has failed.
/// </summary>
public static string NumericalEstimationFailed
{
get
{
return ResourceManager.GetString("NumericalEstimationFailed", resourceCulture);
}
}
/// <summary>
/// Looks up a localized string similar to The two arguments can&apos;t be compared (maybe they are part of a partial ordering?).
/// </summary>

3
src/Numerics/Properties/Resources.resx

@ -448,4 +448,7 @@
<data name="DegreesOfFreedomMustBeLessThanSampleSize" xml:space="preserve">
<value>The sample size must be larger than the given degrees of freedom.</value>
</data>
<data name="NumericalEstimationFailed" xml:space="preserve">
<value>Numerical estimation of the statistic has failed. The used solver did not succeed in finding a root.</value>
</data>
</root>
Loading…
Cancel
Save