Here are the examples of the csharp api System.Math.Log(double) taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.
1089 Examples
19
Source : Mathd.cs
with MIT License
from 734843327
with MIT License
from 734843327
public static double Log(double d) {
return Math.Log(d);
}
19
Source : DataManager.cs
with MIT License
from ABTSoftware
with MIT License
from ABTSoftware
public double GetGaussianRandomNumber(double mean, double stdDev)
{
double u1 = _random.NextDouble(); //these are uniform(0,1) random doubles
double u2 = _random.NextDouble();
double randStdNormal = Math.Sqrt(-2.0 * Math.Log(u1)) *
Math.Sin(2.0 * Math.PI * u2); //random normal(0,1)
double randNormal =
mean + stdDev * randStdNormal; //random normal(mean,stdDev^2)
return randNormal;
}
19
Source : MainControl.xaml.cs
with MIT License
from Actipro
with MIT License
from Actipro
private void InitializeSampleDataContext() {
for (int i = 0; i < 1000; i++) {
int modulus = i % 2;
double xm = i / (20.0d + 3 * modulus);
double ym = 10.0d + 2 * modulus;
double x = random.NextDouble() * xm + 1;
double y = Math.Log(ym * (x - 1.0) + 1.0) * (random.NextDouble() + 0.9);
if (modulus == 0)
primaryChartPoints1.Add(new Point(x, y));
else
primaryChartPoints2.Add(new Point(x, y));
}
}
19
Source : Hazard.cs
with BSD 3-Clause "New" or "Revised" License
from ActuarialIntelligence
with BSD 3-Clause "New" or "Revised" License
from ActuarialIntelligence
public IList<Point<decimal, decimal>> GetHazardFunctionOverEachPeriod()
{
nelsonAalen = new KaplanMeier(observationsInternal);
var result = new List<Point<decimal, decimal>>();
var cnt = 0;
foreach (var set in observationsInternal)
{
var survival = nelsonAalen.GetSurvivalOverPeriod(cnt);
var v1 = survival == 0 ? 0 : (-1) * Math.Log((double)survival);
var v2 = (double)(set.unitTime);
var div = (decimal)((-1) * v1 / v2);
var point = new Point<decimal, decimal>(set.unitTime * cnt, (div));
result.Add(point);
cnt++;
}
return result;
}
19
Source : UnivariateRegressionFitting.cs
with BSD 3-Clause "New" or "Revised" License
from ActuarialIntelligence
with BSD 3-Clause "New" or "Revised" License
from ActuarialIntelligence
private static void GammaDistributionFit(IList<decimal> points)
{
var sum = 0m;
var sumxLNx = 0m;
var sumLNx = 0m;
int n = points.Count();
foreach (var d in points)
{
sum += d;
sumxLNx += d * ((decimal)Math.Log((double)d));
sumLNx += (decimal)Math.Log((double)d);
}
k = ((n * sum) / ((n * sumxLNx) - (sumLNx * sum)));
α = (decimal)((1 / (Math.Pow((double)n, 2))) * (double)((n * sumxLNx) - (sumLNx * sum)));
}
19
Source : TestSoftPlus.cs
with MIT License
from adamtiger
with MIT License
from adamtiger
private double SoftPlusFunc(double x)
{
return Math.Log(1 + Math.Exp(x));
}
19
Source : PrimeFactory.cs
with GNU General Public License v3.0
from AdamWhiteHat
with GNU General Public License v3.0
from AdamWhiteHat
public static BigInteger GetApproximateValueFromIndex(UInt64 n)
{
if (n < 6)
{
return primes[(int)n];
}
double fn = (double)n;
double flogn = Math.Log(n);
double flog2n = Math.Log(flogn);
double upper;
if (n >= 688383) /* Dusart 2010 page 2 */
{
upper = fn * (flogn + flog2n - 1.0 + ((flog2n - 2.00) / flogn));
}
else if (n >= 178974) /* Dusart 2010 page 7 */
{
upper = fn * (flogn + flog2n - 1.0 + ((flog2n - 1.95) / flogn));
}
else if (n >= 39017) /* Dusart 1999 page 14 */
{
upper = fn * (flogn + flog2n - 0.9484);
}
else /* Modified from Robin 1983 for 6-39016 _only_ */
{
upper = fn * (flogn + 0.6000 * flog2n);
}
if (upper >= (double)UInt64.MaxValue)
{
throw new OverflowException($"{upper} > {UInt64.MaxValue}");
}
return new BigInteger((UInt64)Math.Ceiling(upper));
}
19
Source : FileSize.cs
with MIT License
from Adoxio
with MIT License
from Adoxio
public string ToString(int precision, IFormatProvider formatProvider = null)
{
var pow = Math.Floor((_value > 0 ? Math.Log(_value) : 0) / Math.Log(1024));
pow = Math.Min(pow, _units.Length - 1);
var value = _value / Math.Pow(1024, pow);
var precisionString = formatProvider == null
? precision.ToString(CultureInfo.CurrentCulture)
: precision.ToString(formatProvider);
return value.ToString(Math.Abs(pow - 0) < double.Epsilon ? "F0" : "F" + precisionString) + " " + _units[(int)pow];
}
19
Source : Math.cs
with Mozilla Public License 2.0
from ahyahy
with Mozilla Public License 2.0
from ahyahy
[ContextMethod("Логарифм", "Log")]
public double Log(double p1)
{
return System.Math.Log(p1);
}
19
Source : Conversions.cs
with MIT License
from alen-smajic
with MIT License
from alen-smajic
public static Vector2d GeoToWorldPosition(double lat, double lon, Vector2d refPoint, float scale = 1)
{
var posx = lon * OriginShift / 180;
var posy = Math.Log(Math.Tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180);
posy = posy * OriginShift / 180;
return new Vector2d((posx - refPoint.x) * scale, (posy - refPoint.y) * scale);
}
19
Source : Conversions.cs
with MIT License
from alen-smajic
with MIT License
from alen-smajic
public static Vector2d LatLonToMeters(double lat, double lon)
{
var posx = lon * OriginShift / 180;
var posy = Math.Log(Math.Tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180);
posy = posy * OriginShift / 180;
return new Vector2d(posx, posy);
}
19
Source : Conversions.cs
with MIT License
from alen-smajic
with MIT License
from alen-smajic
public static UnwrappedTileId LareplacedudeLongitudeToTileId(double lareplacedude, double longitude, int zoom)
{
var x = (int)Math.Floor((longitude + 180.0) / 360.0 * Math.Pow(2.0, zoom));
var y = (int)Math.Floor((1.0 - Math.Log(Math.Tan(lareplacedude * Math.PI / 180.0)
+ 1.0 / Math.Cos(lareplacedude * Math.PI / 180.0)) / Math.PI) / 2.0 * Math.Pow(2.0, zoom));
return new UnwrappedTileId(zoom, x, y);
}
19
Source : MercatorProjection.cs
with MIT License
from alen-smajic
with MIT License
from alen-smajic
public static double latToY(double lat)
{
lat = Math.Min(89.5, Math.Max(lat, -89.5));
double phi = DegToRad(lat);
double sinphi = Math.Sin(phi);
double con = ECCENT * sinphi;
con = Math.Pow(((1.0 - con) / (1.0 + con)), COM);
double ts = Math.Tan(0.5 * ((Math.PI * 0.5) - phi)) / con;
return 0 - R_MAJOR * Math.Log(ts);
}
19
Source : TileMapAnnotation.cs
with MIT License
from AlexGyver
with MIT License
from AlexGyver
public override void Render(IRenderContext rc, PlotModel model)
{
base.Render(rc, model);
var clippingRect = this.GetClippingRect();
var lon0 = this.XAxis.ActualMinimum;
var lon1 = this.XAxis.ActualMaximum;
var lat0 = this.YAxis.ActualMinimum;
var lat1 = this.YAxis.ActualMaximum;
// the desired number of tiles horizontally
double tilesx = model.Width / this.TileSize;
// calculate the desired zoom level
var n = tilesx / (((lon1 + 180) / 360) - ((lon0 + 180) / 360));
var zoom = (int)Math.Round(Math.Log(n) / Math.Log(2));
if (zoom < this.MinZoomLevel)
{
zoom = this.MinZoomLevel;
}
if (zoom > this.MaxZoomLevel)
{
zoom = this.MaxZoomLevel;
}
// find tile coordinates for the corners
double x0, y0;
LatLonToTile(lat0, lon0, zoom, out x0, out y0);
double x1, y1;
LatLonToTile(lat1, lon1, zoom, out x1, out y1);
double xmax = Math.Max(x0, x1);
double xmin = Math.Min(x0, x1);
double ymax = Math.Max(y0, y1);
double ymin = Math.Min(y0, y1);
// Add the tiles
for (var x = (int)xmin; x < xmax; x++)
{
for (var y = (int)ymin; y < ymax; y++)
{
string uri = this.GetTileUri(x, y, zoom);
var img = this.GetImage(uri, rc.RendersToScreen);
if (img == null)
{
continue;
}
// transform from tile coordinates to lat/lon
double lareplacedude0, lareplacedude1, longitude0, longitude1;
TileToLatLon(x, y, zoom, out lareplacedude0, out longitude0);
TileToLatLon(x + 1, y + 1, zoom, out lareplacedude1, out longitude1);
// transform from lat/lon to screen coordinates
var s00 = this.Transform(longitude0, lareplacedude0);
var s11 = this.Transform(longitude1, lareplacedude1);
var r = OxyRect.Create(s00.X, s00.Y, s11.X, s11.Y);
// draw the image
rc.DrawClippedImage(clippingRect, img, r.Left, r.Top, r.Width, r.Height, this.Opacity, true);
}
}
// draw the copyright notice
var p = new ScreenPoint(clippingRect.Right - 5, clippingRect.Bottom - 5);
var textSize = rc.MeasureText(this.CopyrightNotice, null, 12);
rc.DrawRectangle(new OxyRect(p.X - textSize.Width - 2, p.Y - textSize.Height - 2, textSize.Width + 4, textSize.Height + 4), OxyColors.White.ChangeAlpha(200), null);
rc.DrawText(
p,
this.CopyrightNotice,
OxyColors.Black,
null,
12,
500,
0,
HorizontalAlignment.Right,
VerticalAlignment.Bottom);
}
19
Source : TileMapAnnotation.cs
with MIT License
from AlexGyver
with MIT License
from AlexGyver
private static void LatLonToTile(double lareplacedude, double longitude, int zoom, out double x, out double y)
{
// http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
int n = 1 << zoom;
double lat = lareplacedude / 180 * Math.PI;
x = (longitude + 180.0) / 360.0 * n;
y = (1.0 - Math.Log(Math.Tan(lat) + 1.0 / Math.Cos(lat)) / Math.PI) / 2.0 * n;
}
19
Source : LogarithmicAxis.cs
with MIT License
from AlexGyver
with MIT License
from AlexGyver
public override void GetTickValues(
out IList<double> majorLabelValues, out IList<double> majorTickValues, out IList<double> minorTickValues)
{
if (this.ActualMinimum <= 0)
{
this.ActualMinimum = 0.1;
}
double logBase = Math.Log(this.Base);
var e0 = (int)Math.Floor(Math.Log(this.ActualMinimum) / logBase);
var e1 = (int)Math.Ceiling(Math.Log(this.ActualMaximum) / logBase);
// find the min & max values for the specified base
// round to max 10 digits
double p0 = Math.Pow(this.Base, e0);
double p1 = Math.Pow(this.Base, e1);
double d0 = Math.Round(p0, 10);
double d1 = Math.Round(p1, 10);
if (d0 <= 0)
{
d0 = p0;
}
double d = d0;
majorTickValues = new List<double>();
minorTickValues = new List<double>();
double epsMin = this.ActualMinimum * 1e-6;
double epsMax = this.ActualMaximum * 1e-6;
while (d <= d1 + epsMax)
{
// d = RemoveNoiseFromDoubleMath(d);
if (d >= this.ActualMinimum - epsMin && d <= this.ActualMaximum + epsMax)
{
majorTickValues.Add(d);
}
for (int i = 1; i < this.Base; i++)
{
double d2 = d * (i + 1);
if (d2 > d1 + double.Epsilon)
{
break;
}
if (d2 > this.ActualMaximum)
{
break;
}
if (d2 >= this.ActualMinimum && d2 <= this.ActualMaximum)
{
minorTickValues.Add(d2);
}
}
d *= this.Base;
if (double.IsInfinity(d))
{
break;
}
if (d < double.Epsilon)
{
break;
}
if (double.IsNaN(d))
{
break;
}
}
if (majorTickValues.Count < 2)
{
base.GetTickValues(out majorLabelValues, out majorTickValues, out minorTickValues);
}
else
{
majorLabelValues = majorTickValues;
}
}
19
Source : LogarithmicAxis.cs
with MIT License
from AlexGyver
with MIT License
from AlexGyver
public override double Transform(double x)
{
Debug.replacedert(x > 0, "Value should be positive.");
if (x <= 0)
{
return -1;
}
return (Math.Log(x) - this.offset) * this.scale;
}
19
Source : LogarithmicAxis.cs
with MIT License
from AlexGyver
with MIT License
from AlexGyver
internal override double PreTransform(double x)
{
Debug.replacedert(x > 0, "Value should be positive.");
if (x <= 0)
{
return 0;
}
return Math.Log(x);
}
19
Source : LogarithmicAxis.cs
with MIT License
from AlexGyver
with MIT License
from AlexGyver
internal override void UpdateActualMaxMin()
{
if (this.PowerPadding)
{
double logBase = Math.Log(this.Base);
var e0 = (int)Math.Floor(Math.Log(this.ActualMinimum) / logBase);
var e1 = (int)Math.Ceiling(Math.Log(this.ActualMaximum) / logBase);
if (!double.IsNaN(this.ActualMinimum))
{
this.ActualMinimum = Math.Exp(e0 * logBase).RemoveNoiseFromDoubleMath();
}
if (!double.IsNaN(this.ActualMaximum))
{
this.ActualMaximum = Math.Exp(e1 * logBase).RemoveNoiseFromDoubleMath();
}
}
base.UpdateActualMaxMin();
}
19
Source : Log.cs
with MIT License
from alexshtf
with MIT License
from alexshtf
public override void Eval()
{
Value = Math.Log(Arg.Value);
}
19
Source : Log.cs
with MIT License
from alexshtf
with MIT License
from alexshtf
public override void Diff()
{
var arg = Arg.Value;
Value = Math.Log(arg);
Inputs.SetWeight(ArgIdx, 1 / arg);
}
19
Source : BasicDifferentiationTests.cs
with MIT License
from alexshtf
with MIT License
from alexshtf
[Fact]
public void DiffTermPower()
{
var x = new Variable();
var y = new Variable();
var func = Power(x, y);
var grad = func.Differentiate(Vec(x, y), NumVec(2, 3));
replacedert.Equal(NumVec(12, 8 * Math.Log(2)), grad);
}
19
Source : BasicDifferentiationTests.cs
with MIT License
from alexshtf
with MIT License
from alexshtf
[Fact]
public void DiffTermPowerSingleVariable()
{
var x = new Variable();
var func = Power(x, x);
var grad = func.Differentiate(Vec(x), NumVec(2.5));
var expectedGrad = NumVec(Math.Pow(2.5, 2.5) * (Math.Log(2.5) + 1));
replacedert.Equal(expectedGrad[0], grad[0], 10);
}
19
Source : TermPower.cs
with MIT License
from alexshtf
with MIT License
from alexshtf
public override void Diff()
{
var baseVal = Base.Value;
var expVal = Exponent.Value;
Value = Math.Pow(baseVal, expVal);
Inputs.SetWeight(BaseIdx, expVal * Math.Pow(baseVal, expVal - 1));
Inputs.SetWeight(ExpIdx, Value * Math.Log(baseVal));
}
19
Source : Distributions.cs
with Apache License 2.0
from alexyakunin
with Apache License 2.0
from alexyakunin
private static bool PolarTransform(double x01, double y01, out double x, out double y)
{
var xn11 = 2.0 * x01 - 1.0;
var yn11 = 2.0 * y01 - 1.0;
var r2 = xn11 * xn11 + yn11 * yn11;
if (r2 >= 1.0 || r2 == 0.0) {
x = 0.0;
y = 0.0;
return false;
}
var factor = Math.Sqrt(-2.0 * Math.Log(r2) / r2);
x = xn11 * factor;
y = yn11 * factor;
return true;
}
19
Source : GlobalMercator.cs
with MIT License
from AliFlux
with MIT License
from AliFlux
public CoordinatePair LatLonToMeters(double lat, double lon)
{
CoordinatePair retval = new CoordinatePair();
try
{
retval.X = lon * this.originShift / 180.0;
retval.Y = Math.Log(Math.Tan((90 + lat) * Math.PI / 360.0)) / (Math.PI / 180.0);
retval.Y *= this.originShift / 180.0;
return retval;
}
catch (Exception ex)
{
throw ex;
}
}
19
Source : RandomExtensions.cs
with Apache License 2.0
from allenai
with Apache License 2.0
from allenai
public static double NextGaussian(this Random r, double mu = 0, double sigma = 1) {
var u1 = r.NextDouble();
var u2 = r.NextDouble();
var rand_std_normal = Math.Sqrt(-2.0 * Math.Log(u1)) * Math.Sin(2.0 * Math.PI * u2);
var rand_normal = mu + sigma * rand_std_normal;
// Very rarely, it is possible to have underflow or an infinity if u1/u1 = 0.
// In such a case, just return the mean.
if (
Double.IsNaN(rand_normal)
|| Double.IsInfinity(rand_normal)
|| Double.IsNegativeInfinity(rand_normal)
) {
return mu;
}
return rand_normal;
}
19
Source : DoubleExtension.cs
with MIT License
from AlphaYu
with MIT License
from AlphaYu
public static double Log(this double d)
{
return Math.Log(d);
}
19
Source : RandomExtension.cs
with MIT License
from AlphaYu
with MIT License
from AlphaYu
public static double NextGauss(this Random rand, double mean, double stdDev)
{
double u1 = 1.0 - rand.NextDouble();
double u2 = 1.0 - rand.NextDouble();
double randStdNormal = Math.Sqrt(-2.0 * Math.Log(u1)) * Math.Sin(2.0 * Math.PI * u2);
return mean + stdDev * randStdNormal;
}
19
Source : Program.cs
with MIT License
from altimesh
with MIT License
from altimesh
[MethodImpl(MethodImplOptions.AggressiveInlining), IntrinsicFunction("logf")]
public static float Logf(float f)
{
return (float)Math.Log((double)f);
}
19
Source : Program.cs
with MIT License
from altimesh
with MIT License
from altimesh
[MethodImpl(MethodImplOptions.AggressiveInlining), IntrinsicFunction("__logf")]
public static float Logf(float f)
{
return (float)Math.Log((double)f);
}
19
Source : MathFunctions.cs
with MIT License
from altimesh
with MIT License
from altimesh
[IntrinsicFunction("logf")]
public static float logf(float x)
{
return (float)Math.Log(x);
}
19
Source : NumberToHumanReadableConverter.cs
with GNU General Public License v3.0
from Amebis
with GNU General Public License v3.0
from Amebis
public object Convert(object value, Type targetType, object parameter, CultureInfo culture)
{
if (value == null)
return null;
double number = System.Convert.ToDouble(value);
int b = parameter != null ? System.Convert.ToInt32(parameter) : Base;
if (number <= 0.5 && EmptyIfZero)
return "";
int n = number > 0.5 ? Math.Min((int)Math.Truncate(Math.Log(Math.Abs(number)) / Math.Log(b)), Prefixes.Length) : 0;
double x = number / Math.Pow(b, n);
return string.Format(
Views.Resources.Strings.NumberToHumanReadable,
n > 0 && Math.Abs(x) < 10 ?
(Math.Truncate(x * 10) / 10).ToString("N1") :
Math.Truncate(x).ToString(),
Prefixes[n],
Unit);
}
19
Source : Program.cs
with MIT License
from anastasios-stamoulis
with MIT License
from anastasios-stamoulis
int sample(Random random, float[] preds, double temperature=1.0) {
// step 1: apply temperature to predictions, and normalize them to create a probability distribution
float sum = 0;
for (int i=0; i<preds.Length; i++) {
var p = (float)Math.Exp((Math.Log(Math.Max(preds[i], 1e-10)) / temperature));
sum += p;
preds[i] = p;
}
for (int i = 0; i < preds.Length; i++) { preds[i] /= sum; }
// step 2: draw a random sample from this distribution
var d = random.NextDouble();
sum = 0;
for (int i=0; i<preds.Length; i++) {
sum += preds[i];
if ( d<sum ) { return i; }
}
return preds.Length - 1;
}
19
Source : EdPeltChangePointDetector.cs
with MIT License
from AndreyAkinshin
with MIT License
from AndreyAkinshin
public double GetCost(int tau0, int tau1, int tau2)
{
double sum = 0;
int offset = tau1; // offset of partialSums'[i, tau1] in the single-dimenstional `partialSums` array
int tauDiff = tau2 - tau1;
for (int i = 0; i < k; i++)
{
// actualSum is (count(data[j] < t) * 2 + count(data[j] == t) * 1) for j=tau1..tau2-1
int actualSum =
partialSums[offset + tauDiff] -
partialSums[offset]; // partialSums'[i, tau2] - partialSums'[i, tau1]
// We skip these two cases (correspond to fit = 0 or fit = 1) because of invalid Math.Log values
if (actualSum != 0 && actualSum != tauDiff * 2)
{
// Empirical CDF $\hat{F}_i(t)$ (Section 2.1 "Model" in [Haynes2017])
double fit = actualSum * 0.5 / tauDiff;
// Segment cost $\mathcal{L}_{np}$ (Section 2.2 "Nonparametric maximum likelihood" in [Haynes2017])
double lnp = tauDiff * (fit * Math.Log(fit) + (1 - fit) * Math.Log(1 - fit));
sum += lnp;
}
offset += n + 1;
}
double c = -Math.Log(2 * n - 1); // Constant from Lemma 3.1 in [Haynes2017]
return 2.0 * c / k * sum; // See Section 3.1 "Discrete approximation" in [Haynes2017]
}
19
Source : BetaDistribution.cs
with MIT License
from AndreyAkinshin
with MIT License
from AndreyAkinshin
public double Pdf(double x)
{
if (x < 0 || x > 1)
return 0;
if (x < 1e-9)
{
if (Alpha > 1)
return 0;
if (Abs(Alpha - 1) < 1e-9)
return Beta;
return double.PositiveInfinity;
}
if (x > 1 - 1e-9)
{
if (Beta > 1)
return 0;
if (Abs(Beta - 1) < 1e-9)
return Alpha;
return double.PositiveInfinity;
}
if (Alpha < 1e-9 || Beta < 1e-9)
return 0;
return Exp((Alpha - 1) * Log(x) + (Beta - 1) * Log(1 - x) - BetaFunction.CompleteLogValue(Alpha, Beta));
}
19
Source : LogNormalDistribution.cs
with MIT License
from AndreyAkinshin
with MIT License
from AndreyAkinshin
public double Pdf(double x)
{
if (x < 1e-9)
return 0;
return Exp(-(Log(x) - Mean).Sqr() / (2 * StandardDeviation.Sqr())) / (x * StandardDeviation * Constants.Sqrt2Pi);
}
19
Source : LogNormalDistribution.cs
with MIT License
from AndreyAkinshin
with MIT License
from AndreyAkinshin
public double Cdf(double x)
{
if (x < 1e-9)
return 0;
return 0.5 * (1 + ErrorFunction.Value((Log(x) - Mean) / (Constants.Sqrt2 * StandardDeviation)));
}
19
Source : NormalDistribution.cs
with MIT License
from AndreyAkinshin
with MIT License
from AndreyAkinshin
public override double Next()
{
double u = 0, v = 0;
while (u < 1e-100)
{
u = Random.NextDouble();
v = Random.NextDouble();
}
double stdDevFactor = Sqrt(-2.0 * Log(u)) * Sin(2.0 * PI * v);
return distribution.Mean + distribution.StandardDeviation * stdDevFactor;
}
19
Source : BetaFunction.cs
with MIT License
from AndreyAkinshin
with MIT License
from AndreyAkinshin
public static double IncompleteLogValue(double a, double b, double x)
{
return a * Log(x) - Log(a) + Log(HypergeometricFunction.Value(a, 1 - b, a + 1, x, (int) Round(b)));
}
19
Source : BetaFunction.cs
with MIT License
from AndreyAkinshin
with MIT License
from AndreyAkinshin
public static double RegularizedIncompleteValue(double a, double b, double x)
{
// The implementation is inspired by "Incomplete Beta Function in C" (Lewis Van Winkle, 2017)
// See https://codeplea.com/incomplete-beta-function-c for details
//
// We calculate the regularized incomplete beta function using a continued fraction (https://dlmf.nist.gov/8.17#v):
// Ix(a, b) = x^a * (1-x)^b / (a*B(a, b)) * 1 / (1 + d[1] / (1 + d[2] / (1 + d[3] / (...))))
// where
// d[2m] = m(b-m)x / (a+2m-1)(a+2m)
// d[2m+1] = -(a+m)(a+b+m)x / (a+2m)(a+2m+1)
//
// The approximated value of the continued fraction is calculated using the Lentz's algorithm
replacedertion.NonNegative(nameof(a), a);
replacedertion.NonNegative(nameof(b), b);
const double eps = 1e-8;
if (x < eps)
return 0;
if (x > 1 - eps)
return 1;
if (a < eps && b < eps)
return 0.5;
if (a < eps)
return 1;
if (b < eps)
return 0;
// According to https://dlmf.nist.gov/8.17#v, the continued fraction converges rapidly for x<(a+1)/(a+b+2)
// If x>=(a+1)/(a+b+2), we use Ix(a, b) = I{1-x}(b, a)
if (x > (a + 1) / (a + b + 2))
return 1.0 - RegularizedIncompleteValue(b, a, 1 - x);
// We use the Lentz's algorithm to calculate the continued fraction
// f = 1 + d[1] / (1 + d[2] / (1 + d[3] / (...)))
// The implementation is based on the following formulas:
// u[0] = 1, v[0] = 0, f[0] = 1
// u[i] = 1 + d[i] / u[i - 1]
// v[i] = 1 / (1 + d[i] * v[i - 1])
// f[i] = f[i - 1] * u[i] * v[i]
const int maxIterationCount = 300;
static double Normalize(double z) => Abs(z) < 1e-30 ? 1e-30 : z; // Normalization prevents getting zero values
double u = 1, v = 0, f = 1;
for (int i = 0; i <= maxIterationCount; i++)
{
double d; // d[i]
int m = i / 2;
if (i == 0)
d = 1.0; // d[0]
else if (i % 2 == 0)
d = m * (b - m) * x / ((a + 2 * m - 1) * (a + 2 * m)); // d[2m]
else
d = -((a + m) * (a + b + m) * x) / ((a + 2.0 * m) * (a + 2.0 * m + 1)); // d[2m+1]
u = Normalize(1 + d / u);
v = 1 / Normalize(1 + d * v);
double uv = u * v;
f *= uv;
if (Abs(uv - 1) < eps)
break;
}
// Ix(a, b) = x^a * (1-x)^b / (a*B(a, b)) * 1 / (1 + d[1] / (1 + d[2] / (1 + d[3] / (...))))
return Exp(Log(x) * a + Log(1.0 - x) * b - CompleteLogValue(a, b)) / a * (f - 1);
}
19
Source : BetaFunction.cs
with MIT License
from AndreyAkinshin
with MIT License
from AndreyAkinshin
public static double RegularizedIncompleteInverseValue(double a, double b, double p)
{
// The implementation is based on "Incomplete Beta Function" from "Numerical Recipes", 3rd edition, page 273
replacedertion.NonNegative(nameof(a), a);
replacedertion.NonNegative(nameof(b), b);
if (p <= 0)
return 0;
if (p >= 1)
return 1;
const double eps = 1e-8;
double t, u, x, w;
if (a >= 1 && b >= 1)
{
double pp = p < 0.5 ? p : 1.0 - p;
t = Sqrt(-2.0 * Log(pp));
x = (2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t;
if (p < 0.5)
x = -x;
double al = (x.Sqr() - 3.0) / 6.0;
double h = 2.0 / (1.0 / (2.0 * a - 1.0) + 1.0 / (2.0 * b - 1.0));
w = x * Sqrt(al + h) / h - (1.0 / (2.0 * b - 1) - 1.0 / (2.0 * a - 1.0)) * (al + 5.0 / 6.0 - 2.0 / (3.0 * h));
x = a / (a + b * Exp(2.0 * w));
}
else
{
double lna = Log(a / (a + b));
double lnb = Log(b / (a + b));
t = Exp(a * lna) / a;
u = Exp(b * lnb) / b;
w = t + u;
x = p < t / w
? Pow(a * w * p, 1.0 / a)
: 1.0 - Pow(b * w * (1.0 - p), 1.0 / b);
}
double afac = -GammaFunction.LogValue(a) - GammaFunction.LogValue(b) + GammaFunction.LogValue(a + b);
for (int iteration = 0; iteration < 10; iteration++)
{
if (x < eps || x > 1.0 - eps)
return x; // a or b are too small for accurate calculations
double error = RegularizedIncompleteValue(a, b, x) - p;
t = Exp((a - 1) * Log(x) + (b - 1) * Log(1.0 - x) + afac);
u = error / t;
t = u / (1.0 - 0.5 * Min(1.0, u * ((a - 1) / x - (b - 1) / (1.0 - x)))); // Halley's method
x -= t;
if (x <= 0.0)
x = 0.5 * (x + t);
if (x >= 1.0)
x = 0.5 * (x + t + 1.0); // Bisect if x tries to go negative or > 1
if (Abs(t) < eps * x && iteration > 0)
break;
}
return x;
}
19
Source : ErrorFunction.cs
with MIT License
from AndreyAkinshin
with MIT License
from AndreyAkinshin
public static double InverseValue(double p)
{
replacedertion.InRangeExclusive(nameof(p), p, -1, 1);
p = 1 - p;
double pp = p < 1.0 ? p : 2 - p;
double t = Sqrt(-2 * Log(pp / 2));
double x = -0.70711 * ((2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t);
for (int i = 0; i < 2; i++)
{
double err = (1 - Value(x)) - pp;
x += err / (1.12837916709551257 * Exp(-x.Sqr()) - x * err);
}
return p < 1.0 ? x : -x;
}
19
Source : FactorialFunction.cs
with MIT License
from AndreyAkinshin
with MIT License
from AndreyAkinshin
public static double LogValue(int n)
{
replacedertion.NonNegative(nameof(n), n);
if (n <= 20)
return Math.Log(Value(n));
return GammaFunction.LogValue(n + 1);
}
19
Source : GammaFunction.cs
with MIT License
from AndreyAkinshin
with MIT License
from AndreyAkinshin
public static double LogValue(double x)
{
if (x < 1e-5)
throw new ArgumentOutOfRangeException(nameof(x), "x should be positive");
if (x < 3)
return Log(Value(x));
return StirlingApproximationLog(x);
}
19
Source : GammaFunction.cs
with MIT License
from AndreyAkinshin
with MIT License
from AndreyAkinshin
private static double StirlingApproximationLog(double x)
{
return x * Log(x) - x + Log(2 * PI / x) / 2 + GetSeriesValue(x);
}
19
Source : ExponentialDecaySequence.cs
with MIT License
from AndreyAkinshin
with MIT License
from AndreyAkinshin
public static ISequence CreateFromHalfLife(int halfLife) => new ExponentialDecaySequence(1, Math.Log(2) / halfLife);
19
Source : BetaFunctionTests.cs
with MIT License
from AndreyAkinshin
with MIT License
from AndreyAkinshin
[Fact]
public void BetaCompleteLogValue()
{
var comparer = new AbsoluteEqualityComparer(0.0000001);
for (int a = 1; a <= 20; a++)
for (int b = 1; b <= 20; b++)
{
double actual = BetaFunction.CompleteLogValue(a, b);
double expected = Math.Log(Factorial(a - 1) * Factorial(b - 1) / Factorial(a + b - 1));
replacedert.Equal(expected, actual, comparer);
}
}
19
Source : GammaFunctionTests.cs
with MIT License
from AndreyAkinshin
with MIT License
from AndreyAkinshin
[Fact]
public void GammaFunctionLogValue()
{
var comparer = new AbsoluteEqualityComparer(0.0001);
for (int i = 0; i < knownX.Length; i ++)
{
double expected = Math.Log(knownY[i]);
double actual = GammaFunction.LogValue(knownX[i]);
output.WriteLine($"X = {knownX[i]}");
output.WriteLine($"Expected = {expected}");
output.WriteLine($"Actual = {actual}");
output.WriteLine("");
replacedert.Equal(expected, actual, comparer);
}
}
19
Source : TimeSeriesAndForecasting.cs
with MIT License
from AngeloCresta
with MIT License
from AngeloCresta
private void Regression( RegressionType regressionType, double [][] inputValues, out double [][] outputValues, int polynomialDegree, int forecastingPeriod )
{
if( regressionType == RegressionType.Exponential )
{
double [] oldYValues = new double[ inputValues[1].Length ];
for( int index = 0; index < inputValues[1].Length; index++ )
{
oldYValues[ index ] = inputValues[1][index];
if( inputValues[1][index] <= 0 )
{
throw new InvalidOperationException(SR.ExceptionForecastingExponentialRegressionHasZeroYValues);
}
inputValues[1][index] = Math.Log( inputValues[1][index] );
}
PolynomialRegression( regressionType, inputValues, out outputValues, 2, forecastingPeriod, 0 );
inputValues[1] = oldYValues;
}
else if( regressionType == RegressionType.Logarithmic )
{
double interval;
double first = inputValues[0][0];
// Find Interval for X values
interval = Math.Abs( inputValues[0][0] - inputValues[0][inputValues[0].Length - 1] ) / ( inputValues[0].Length - 1 );
if( interval <= 0 )
interval = 1;
for( int index = 0; index < inputValues[0].Length; index++ )
{
inputValues[0][index] = Math.Log( inputValues[0][index] );
}
PolynomialRegression( regressionType, inputValues, out outputValues, 2, forecastingPeriod, interval );
// Create Y values based on approximation.
for( int i = 0; i < outputValues[0].Length; i++ )
{
// Set X value
outputValues[0][i] = first + i * interval;
}
}
else if( regressionType == RegressionType.Power )
{
double [] oldYValues = new double[ inputValues[1].Length ];
double interval;
double first = inputValues[0][0];
for( int index = 0; index < inputValues[1].Length; index++ )
{
oldYValues[ index ] = inputValues[1][index];
if( inputValues[1][index] <= 0 )
{
throw new InvalidOperationException(SR.ExceptionForecastingPowerRegressionHasZeroYValues);
}
}
// Find Interval for X values
interval = Math.Abs( inputValues[0][0] - inputValues[0][inputValues[0].Length - 1] ) / ( inputValues[0].Length - 1 );
if( interval <= 0 )
interval = 1;
PolynomialRegression( regressionType, inputValues, out outputValues, 2, forecastingPeriod, interval );
inputValues[1] = oldYValues;
// Create Y values based on approximation.
for( int i = 0; i < outputValues[0].Length; i++ )
{
// Set X value
outputValues[0][i] = first + i * interval;
}
}
else
{
PolynomialRegression( regressionType, inputValues, out outputValues, polynomialDegree, forecastingPeriod, 0 );
}
}
See More Examples