using MathNet.Numerics;
|
using MathNet.Numerics.LinearAlgebra.Double;
|
using NPOI.SS.Formula.Functions;
|
using System;
|
using System.Collections.Generic;
|
using System.Drawing;
|
using System.Linq;
|
using System.Text;
|
using System.Threading.Tasks;
|
|
namespace Hydro.CommonBase
|
{
|
public class CurveFitHelper
|
{
|
|
|
public DenseVector coefficients;
|
int degree;
|
double xMax;
|
double xMin;
|
List<PointF> data;
|
bool isFitted = false;
|
|
public CurveFitHelper(List<PointF> data, int degree = 2)
|
{
|
this.degree = degree;
|
this.data = data;
|
|
double[] xlist = data.Select(p => (double)p.X).ToArray();
|
coefficients = Fit.Polynomial(xlist, data.Select(p => (double)p.Y).ToArray(), 2);
|
isFitted = true;
|
var list = xlist.ToList();
|
xMax = list.Max();
|
xMin = list.Min();
|
}
|
|
public double Evaluate(double x)
|
{
|
while (!isFitted) ;
|
double y = 0;
|
for (int i = 0; i < coefficients.Count; i++)
|
{
|
y += (float)(coefficients[i] * Math.Pow(x, i));
|
}
|
return y;
|
}
|
|
|
public List<PointF> GetFitCurve(int count)
|
{
|
List<PointF> curve = new List<PointF>();
|
for (double x = xMin; x < xMax; x += (xMax - xMin) / count)
|
{
|
double y = Evaluate(x);
|
curve.Add(new PointF((float)x, (float)y));
|
}
|
return curve;
|
|
}
|
|
public List<List<PointF>> GetConfidenceCurve(int count)
|
{
|
List<PointF> up_list = new List<PointF>();
|
List<PointF> down_list = new List<PointF>();
|
List<PointF> avg_list = new List<PointF>();
|
|
double X =0;
|
double Y=0;
|
for (double x = xMin; x < xMax; x += (xMax - xMin) / count)
|
{
|
var list = data.FindAll(p => x - (xMax - xMin) / count / 2 <= p.X && p.X < x + (xMax - xMin) / count/2).Select(p => (double)p.Y).ToList();
|
if (list.Count >= 3)
|
{
|
var Bounds = ConfidenceRangeHelper.GetConfidenceRange(list, 95);
|
var x_avg = x;
|
if (Bounds[0] == double.NaN || Bounds[1] == double.NaN) continue;
|
down_list.Add(new PointF((float)x_avg, (float)Bounds[0]));
|
up_list.Add(new PointF((float)x_avg, (float)Bounds[1]));
|
X = x;
|
//平均值
|
Y = (Bounds[0] + Bounds[1] )/ 2;
|
avg_list.Add(new PointF((float)x_avg, (float)Y));
|
}
|
else if (list.Count == 2)
|
{
|
var x_avg = x;
|
down_list.Add(new PointF((float)x_avg, (float)list.Min()));
|
up_list.Add(new PointF((float)x_avg, (float)list.Max()));
|
X = x;
|
Y = (list.Min() + list.Max()) / 2;
|
avg_list.Add(new PointF((float)x_avg, (float)Y));
|
}
|
else if (list.Count == 1)
|
{
|
var x_avg = x;
|
down_list.Add(new PointF((float)x_avg, (float)list[0]));
|
up_list.Add(new PointF((float)x_avg, (float)list[0]));
|
X = x;
|
Y = list[0];
|
avg_list.Add(new PointF((float)x_avg, (float)Y));
|
}
|
else if(list.Count == 0)
|
{
|
//使用前一个点的值,X和Y都是前一个点的值
|
down_list.Add(new PointF((float)X, (float)Y));
|
up_list.Add(new PointF((float)X, (float)Y));
|
avg_list.Add(new PointF((float)X, (float)Y));
|
}
|
|
}
|
|
//CurveFitHelper h0 = new CurveFitHelper(down_list);
|
//var curve0 = h0.GetFitCurve(count);
|
//CurveFitHelper h1 = new CurveFitHelper(up_list);
|
//var curve1 = h1.GetFitCurve(count);
|
return new List<List<PointF>> { down_list, up_list, avg_list };
|
}
|
}
|
|
public static class ConfidenceRangeHelper
|
{
|
public static double[] GetConfidenceRange(List<double> data, double confidenceLevel = 90,bool mode=true)//true是中位数，置信区间是false
|
{
|
//如果model==true,这计算中位数,比例为confidenceLevel
|
if (mode)
|
{
|
|
data.Sort();
|
|
double lowerBound = 0;
|
double upperBound = 0;
|
//取百分位confidenceLevel/100和1-confidenceLevel/100.0的数
|
int upperIndex = (int)Math.Round( data.Count*confidenceLevel/100.0);
|
int lowerIndex = (int)Math.Round(data.Count * (1-confidenceLevel / 100.0));
|
if (upperIndex == data.Count)
|
{
|
upperIndex = data.Count - 1;
|
}
|
if (lowerIndex <0)
|
{
|
lowerIndex = 0;
|
}
|
return new double[] { data[lowerIndex], data[upperIndex] };
|
|
|
}
|
|
else
|
{
|
// 使用循环计算数据的最小值和最大值
|
double minValue = double.MaxValue;
|
double maxValue = double.MinValue;
|
foreach (double value in data)
|
{
|
if (value < minValue)
|
{
|
minValue = value;
|
}
|
if (value > maxValue)
|
{
|
maxValue = value;
|
}
|
}
|
|
// 使用循环计算样本均值
|
double sum = 0;
|
foreach (double value in data)
|
{
|
sum += value;
|
}
|
double mean = sum / data.Count;
|
|
// 使用循环计算标准差
|
double sumOfSquares = 0;
|
foreach (double value in data)
|
{
|
sumOfSquares += Math.Pow(value - mean, 2);
|
}
|
double stdDev = Math.Sqrt(sumOfSquares / (data.Count - 1));
|
|
double zValue = 0;
|
double alpha = 1 - confidenceLevel / 100.0;
|
|
double zValueOneTail = GetZValueOneTail(alpha / 2);
|
double marginOfError = zValueOneTail * stdDev;
|
double lowerBound = mean - marginOfError;
|
double upperBound = mean + marginOfError;
|
|
return new double[] { lowerBound, upperBound };
|
}
|
|
|
|
}
|
static double[,] zTable = {
|
{0.0000, 0.0039, 0.0078, 0.0117, 0.0156, 0.0195, 0.0234, 0.0274, 0.0314, 0.0353},
|
{0.0392, 0.0432, 0.0471, 0.0510, 0.0550, 0.0589, 0.0628, 0.0668, 0.0708, 0.0749},
|
{0.0789, 0.0830, 0.0871, 0.0912, 0.0953, 0.0994, 0.1035, 0.1076, 0.1117, 0.1158},
|
{0.1190, 0.1231, 0.1272, 0.1314, 0.1357, 0.1398, 0.1441, 0.1483, 0.1525, 0.1568},
|
{0.1611, 0.1654, 0.1697, 0.1740, 0.1783, 0.1826, 0.1869, 0.1912, 0.1955, 0.1999},
|
{0.2042, 0.2085, 0.2129, 0.2172, 0.2216, 0.2260, 0.2304, 0.2348, 0.2392, 0.2436},
|
{0.2480, 0.2524, 0.2568, 0.2612, 0.2658, 0.2703, 0.2747, 0.2793, 0.2838, 0.2883},
|
{0.2929, 0.2974, 0.3020, 0.3066, 0.3112, 0.3159, 0.3205, 0.3251, 0.3299, 0.3346},
|
{0.3393, 0.3441, 0.3487, 0.3535, 0.3583, 0.3632, 0.3679, 0.3727, 0.3776, 0.3823},
|
{0.3872, 0.3920, 0.3968, 0.4017, 0.4066, 0.4115, 0.4164, 0.4214, 0.4264, 0.4313}
|
};
|
static double GetZValueOneTail(double alpha)
|
{
|
if (alpha == 0.5)
|
{
|
return 0;
|
}
|
else if (alpha > 0.5)
|
{
|
// Use symmetry property of the normal distribution
|
return -GetZValueOneTail(1 - alpha);
|
}
|
double[] alphaValues = { 0.1, 0.05, 0.025, 0.01, 0.005, 0.001 };
|
double[] zValues = { 1.28, 1.645, 1.96, 2.33, 2.575, 3.09 };
|
int n = alphaValues.Length;
|
int i = 0;
|
while (i < n && alpha < alphaValues[i])
|
{
|
i++;
|
}
|
|
if (i == 0)
|
{
|
return 0;
|
}
|
else if (i == n)
|
{
|
return zValues[n - 1];
|
}
|
|
double x0 = alphaValues[i - 1];
|
double x1 = alphaValues[i];
|
double y0 = zValues[i - 1];
|
double y1 = zValues[i];
|
|
double zValue = ((y1 - y0) / (x1 - x0)) * (alpha - x0) + y0;
|
|
return zValue;
|
}
|
|
}
|
}
|
