using MathNet.Numerics.Distributions;
|
using MathNet.Numerics.LinearAlgebra;
|
using MathNet.Numerics.LinearAlgebra.Double;
|
|
namespace IStation.Test
|
{
|
internal class Program
|
{
|
static void Main(string[] args)
|
{
|
var bll = new BLL.StationSignalRecordPacket();
|
var projectId = 661070185922629;
|
IStation.SettingsD.Project.ID = projectId;
|
var monitorDataSourcesId = 606203941007429;
|
var dateStart = new DateTime(2024, 1, 1);
|
var dateEnd = new DateTime(2025, 1, 1);
|
var stationDict = new Dictionary<int, long>();
|
stationDict.Add(1, 462958406303813);
|
stationDict.Add(2, 462958422204485);
|
|
|
|
var ptFilterList = new List<PointViewModel>();
|
foreach (var station in stationDict)
|
{
|
var stationIndex = station.Key;
|
var stationId = station.Value;
|
var packets = bll.Get(monitorDataSourcesId, stationId);
|
var ptList = new List<PointViewModel>();
|
|
var records = packets.SelectMany(x => x.StationSignalRecords).ToList();
|
records.ForEach(x =>
|
{
|
if (x.TotalPressure > 0 && x.TotalFlow > 0)
|
{
|
ptList.Add(new PointViewModel(x.Time, x.TotalFlow, x.TotalPressure, stationIndex));
|
}
|
});
|
|
var newPtList = DynamicThresholdProcessor.Filter(ptList);
|
ptFilterList.AddRange(newPtList);
|
}
|
|
|
var recordList = new List<RecordViewModel>();
|
var timeGroup = ptFilterList.GroupBy(x => x.Time);
|
foreach (var group in timeGroup)
|
{
|
if (group.Count() < 2)
|
continue;
|
var record = new RecordViewModel
|
{
|
Time = group.Key
|
};
|
foreach (var item in group)
|
{
|
if (item.Index == 1)
|
{
|
record.Flow1 = item.X;
|
record.Pressure1 = item.Y;
|
}
|
else
|
{
|
record.Flow2 = item.X;
|
record.Pressure2 = item.Y;
|
}
|
}
|
|
recordList.Add(record);
|
}
|
|
var fullPath = Path.Combine(AppDomain.CurrentDomain.BaseDirectory, "csv");
|
if (!Directory.Exists(fullPath))
|
{
|
Directory.CreateDirectory(fullPath);
|
}
|
|
|
{
|
|
string filePath = Path.Combine(fullPath, "all.csv");
|
using StreamWriter writer = new StreamWriter(filePath, false, System.Text.Encoding.UTF8);
|
|
writer.WriteLine($"Time,Flow1,Pressure1,Flow2,Pressure2");
|
foreach (var record in recordList)
|
{
|
writer.WriteLine($"{record.Time},{record.Flow1},{record.Pressure1},{record.Flow2},{record.Pressure2}");
|
}
|
|
}
|
|
|
|
var dataCellList = new List<DataCell>() {
|
new DataCell { Name = "Flow1", ValueList = recordList.Select(x => x.Flow1).ToArray() },
|
new DataCell { Name = "Pressure1", ValueList = recordList.Select(x => x.Pressure1).ToArray() },
|
new DataCell { Name = "Flow2", ValueList = recordList.Select(x => x.Flow2).ToArray() },
|
new DataCell { Name = "Pressure2", ValueList = recordList.Select(x => x.Pressure2).ToArray() } };
|
|
|
var combineList = IStation.Curve.PermutationAndCombination<DataCell>.GetCombination(dataCellList.ToArray(), 2);
|
foreach (var combine in combineList)
|
{
|
var x = combine[0];
|
var y = combine[1];
|
|
string filePath = Path.Combine(fullPath, $"{x.Name}-{y.Name}.csv");
|
using StreamWriter writer = new StreamWriter(filePath, false, System.Text.Encoding.UTF8);
|
var count = x.ValueList.Length;
|
writer.WriteLine($"{x.Name},{y.Name}");
|
for (int i = 0; i < count; i++)
|
{
|
var xv = x.ValueList[i];
|
var yv = y.ValueList[i];
|
writer.WriteLine($"{xv},{yv}");
|
}
|
}
|
|
Console.WriteLine("计算皮尔逊积差相关系数");
|
foreach (var combine in combineList)
|
{
|
var x = combine[0];
|
var y = combine[1];
|
var (r, pValue) = PearsonTest(x.ValueList, y.ValueList);
|
Console.WriteLine($"{x.Name}-{y.Name}: r = {r:F3}, p = {pValue}");
|
|
}
|
|
Console.WriteLine("计算斯皮尔曼等级相关系数");
|
foreach (var combine in combineList)
|
{
|
var x = combine[0];
|
var y = combine[1];
|
|
var (rho, pValue) = SpearmanTest(x.ValueList, y.ValueList);
|
Console.WriteLine($"{x.Name}-{y.Name}: ρ = {rho:F3}, p = {pValue}");
|
}
|
|
Console.WriteLine("ok");
|
Console.ReadKey();
|
}
|
|
|
/// <summary>
|
/// 计算皮尔逊相关系数及显著性p值
|
/// </summary>
|
public static (double r, double pValue) PearsonTest(double[] x, double[] y)
|
{
|
// 计算皮尔逊相关系数
|
double r = Correlation.Pearson(x, y);
|
|
// 计算显著性p值（双尾t检验）
|
int n = x.Length;
|
double t = r * Math.Sqrt((n - 2) / (1 - r * r));
|
var tDist = new StudentT(location: 0, scale: 1, freedom: n - 2);
|
double pValue = 2 * (1 - tDist.CumulativeDistribution(Math.Abs(t)));
|
|
return (r, pValue);
|
}
|
|
/// <summary>
|
/// 计算斯皮尔曼等级相关系数及显著性p值
|
/// </summary>
|
public static (double rho, double pValue) SpearmanTest(double[] x, double[] y)
|
{
|
// 计算斯皮尔曼相关系数
|
double rho = Correlation.Spearman(x, y);
|
|
// 大样本正态近似法计算p值
|
int n = x.Length;
|
double z = rho * Math.Sqrt(n - 1);
|
var normalDist = new Normal();
|
double pValue = 2 * (1 - normalDist.CumulativeDistribution(Math.Abs(z)));
|
|
return (rho, pValue);
|
}
|
}
|
|
|
|
public class DataCell
|
{
|
public string Name { get; set; }
|
public double[] ValueList { get; set; }
|
|
}
|
|
public class PointViewModel
|
{
|
public PointViewModel() { }
|
public PointViewModel(DateTime dt, double x, double y, int index)
|
{
|
this.Time = dt;
|
this.X = x;
|
this.Y = y;
|
this.Index = index;
|
}
|
|
public PointViewModel(PointViewModel rhs)
|
{
|
this.Time = rhs.Time;
|
this.X = rhs.X;
|
this.Y = rhs.Y;
|
this.Index = rhs.Index;
|
}
|
|
public DateTime Time { get; set; }
|
public int Index { get; set; }
|
public double X { get; set; }
|
public double Y { get; set; }
|
|
}
|
|
public class RecordViewModel
|
{
|
public RecordViewModel() { }
|
|
public DateTime Time { get; set; }
|
|
public double Flow1 { get; set; }
|
public double Pressure1 { get; set; }
|
|
public double Flow2 { get; set; }
|
|
public double Pressure2 { get; set; }
|
}
|
|
public static class Correlation
|
{
|
/// <summary>
|
/// 计算皮尔逊积差相关系数
|
/// </summary>
|
/// <param name="dataA">数据样本A</param>
|
/// <param name="dataB">数据样本B</param>
|
/// <returns>皮尔逊积差相关系数</returns>
|
public static double Pearson(IEnumerable<double> dataA, IEnumerable<double> dataB)
|
{
|
int n = 0;
|
double r = 0.0;
|
double meanA = 0;
|
double meanB = 0;
|
double varA = 0;
|
double varB = 0;
|
|
using (var ieA = dataA.GetEnumerator())
|
using (var ieB = dataB.GetEnumerator())
|
{
|
while (ieA.MoveNext())
|
{
|
if (!ieB.MoveNext())
|
{
|
throw new ArgumentOutOfRangeException("dataB", "数据长度不一致");
|
}
|
|
double currentA = ieA.Current;
|
double currentB = ieB.Current;
|
|
double deltaA = currentA - meanA;
|
double scaleDeltaA = deltaA / ++n;
|
|
double deltaB = currentB - meanB;
|
double scaleDeltaB = deltaB / n;
|
|
meanA += scaleDeltaA;
|
meanB += scaleDeltaB;
|
|
varA += scaleDeltaA * deltaA * (n - 1);
|
varB += scaleDeltaB * deltaB * (n - 1);
|
r += (deltaA * deltaB * (n - 1)) / n;
|
}
|
|
if (ieB.MoveNext())
|
{
|
throw new ArgumentOutOfRangeException("dataA", "数据长度不一致");
|
}
|
}
|
|
return r / Math.Sqrt(varA * varB);
|
}
|
|
/// <summary>
|
/// 计算斯皮尔曼等级相关系数
|
/// </summary>
|
public static double Spearman(IEnumerable<double> dataA, IEnumerable<double> dataB)
|
{
|
return Pearson(Rank(dataA), Rank(dataB));
|
}
|
|
private static IEnumerable<double> Rank(IEnumerable<double> sequence)
|
{
|
var sorted = new List<double>(sequence);
|
sorted.Sort();
|
|
return (IEnumerable<double>)sequence.Select(x => (double)sorted.IndexOf(x) + 1);
|
}
|
}
|
|
|
|
public class DynamicThresholdProcessor
|
{
|
|
public static List<PointViewModel> Filter(List<PointViewModel> ptList)
|
{
|
var pressures = ptList.Select(p => p.Y).ToList();
|
|
// 计算统计量
|
var (mean, stdDev) = CalculateStats(pressures);
|
double skewness = CalculateSkewness(pressures);
|
// 动态调整σ倍数
|
double sigmaMultiplier = CalculateSigmaMultiplier(skewness);
|
sigmaMultiplier = 3;
|
|
|
// 计算边界
|
double lower = mean - sigmaMultiplier * stdDev;
|
double upper = mean + sigmaMultiplier * stdDev;
|
|
return ptList.Where(p => p.Y >= lower && p.Y <= upper).ToList();
|
}
|
|
|
// 核心统计计算
|
private static (double mean, double stdDev) CalculateStats(List<double> values)
|
{
|
double mean = values.Average();
|
double stdDev = Math.Sqrt(values.Sum(v => Math.Pow(v - mean, 2)) / (values.Count - 1));
|
return (mean, stdDev);
|
}
|
|
// 偏度计算（Pearson's moment coefficient）
|
private static double CalculateSkewness(List<double> values)
|
{
|
double mean = values.Average();
|
double std = CalculateStats(values).stdDev;
|
double sum = values.Sum(v => Math.Pow((v - mean) / std, 3));
|
return (sum * values.Count) / ((values.Count - 1) * (values.Count - 2));
|
}
|
|
// 动态σ倍数计算规则
|
private static double CalculateSigmaMultiplier(double skewness)
|
{
|
|
return skewness switch
|
{
|
> 1.0 => 2.0, // 强正偏态
|
> 0.5 => 2.5,
|
> -0.5 => 3.0, // 近似正态
|
> -1.0 => 3.5,
|
_ => 4.0 // 强负偏态
|
};
|
}
|
|
}
|
|
|
|
public class NonlinearRegressionExample
|
{
|
//static void Main(string[] args)
|
//{
|
// // 示例数据：Flow1 和 Pressure1
|
// double[] flow1 = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 };
|
// double[] pressure1 = { 2.1, 4.5, 6.8, 9.2, 11.5, 13.8, 16.1, 18.4, 20.7, 23.0 };
|
|
// // 多项式回归的阶数（例如，二次多项式）
|
// int degree = 2;
|
|
// // 拟合多项式回归模型
|
// double[] coefficients = FitPolynomial(flow1, pressure1, degree);
|
|
// // 输出回归系数
|
// Console.WriteLine("多项式回归系数：");
|
// for (int i = 0; i <= degree; i++)
|
// {
|
// Console.WriteLine($"系数 {i}: {coefficients[i]}");
|
// }
|
|
// // 使用模型进行预测
|
// double xNew = 11.0;
|
// double yPredicted = PredictPolynomial(xNew, coefficients);
|
// Console.WriteLine($"\n当 Flow1 = {xNew} 时，预测的 Pressure1 = {yPredicted}");
|
|
// // 计算 R² 和 MSE
|
// double rSquared = CalculateRSquared(flow1, pressure1, coefficients);
|
// double mse = CalculateMSE(flow1, pressure1, coefficients);
|
// Console.WriteLine($"\nR² = {rSquared}");
|
// Console.WriteLine($"MSE = {mse}");
|
//}
|
|
public static void Test(double[] x, double[] y)
|
{
|
// 多项式回归的阶数（例如，二次多项式）
|
int degree = 2;
|
|
// 拟合多项式回归模型
|
double[] coefficients = FitPolynomial(x, y, degree);
|
|
// 输出回归系数
|
Console.WriteLine("多项式回归系数：");
|
for (int i = 0; i <= degree; i++)
|
{
|
Console.WriteLine($"系数 {i}: {coefficients[i]}");
|
}
|
|
// 使用模型进行预测
|
double xNew = x.Average();
|
double yPredicted = PredictPolynomial(xNew, coefficients);
|
Console.WriteLine($"\n当 X = {xNew} 时，预测的 Y = {yPredicted}");
|
|
// 计算 R² 和 MSE
|
double rSquared = CalculateRSquared(x, y, coefficients);
|
double mse = CalculateMSE(x, y, coefficients);
|
Console.WriteLine($"\nR² = {rSquared}");
|
Console.WriteLine($"MSE = {mse}");
|
|
}
|
|
// 多项式回归拟合
|
public static double[] FitPolynomial(double[] x, double[] y, int degree)
|
{
|
// 构建设计矩阵
|
Matrix<double> X = DenseMatrix.OfArray(new double[x.Length, degree + 1]);
|
for (int i = 0; i < x.Length; i++)
|
{
|
for (int j = 0; j <= degree; j++)
|
{
|
X[i, j] = Math.Pow(x[i], j);
|
}
|
}
|
|
// 构建目标向量
|
Vector<double> Y = DenseVector.OfArray(y);
|
|
// 使用最小二乘法求解
|
var qr = X.QR();
|
Vector<double> coefficients = qr.Solve(Y);
|
|
return coefficients.ToArray();
|
}
|
|
// 使用多项式模型进行预测
|
public static double PredictPolynomial(double x, double[] coefficients)
|
{
|
double prediction = 0.0;
|
for (int i = 0; i < coefficients.Length; i++)
|
{
|
prediction += coefficients[i] * Math.Pow(x, i);
|
}
|
return prediction;
|
}
|
|
// 计算 R²（决定系数）
|
public static double CalculateRSquared(double[] x, double[] y, double[] coefficients)
|
{
|
double yMean = 0.0;
|
double ssTotal = 0.0;
|
double ssResidual = 0.0;
|
|
for (int i = 0; i < y.Length; i++)
|
{
|
yMean += y[i];
|
}
|
yMean /= y.Length;
|
|
for (int i = 0; i < y.Length; i++)
|
{
|
double yPredicted = PredictPolynomial(x[i], coefficients);
|
ssTotal += Math.Pow(y[i] - yMean, 2);
|
ssResidual += Math.Pow(y[i] - yPredicted, 2);
|
}
|
|
return 1.0 - (ssResidual / ssTotal);
|
}
|
|
// 计算 MSE（均方误差）
|
public static double CalculateMSE(double[] x, double[] y, double[] coefficients)
|
{
|
double mse = 0.0;
|
for (int i = 0; i < y.Length; i++)
|
{
|
double yPredicted = PredictPolynomial(x[i], coefficients);
|
mse += Math.Pow(y[i] - yPredicted, 2);
|
}
|
return mse / y.Length;
|
}
|
}
|
|
|
|
}
|
