пятница, 10 октября 2025 г.

25.10.09, Stat, Test, Norm

25.10.09, Stat, Test, Norm, Box-Muller, Komlmorov-Smirnov, Chi-qudrat, Shapiro-Wilk,  DeepSeek,

1. Генератор нормального распределения (Box-Muller)

csharp
using System;
using System.Collections.Generic;
using System.Linq;

public class NormalDistributionGenerator
{, 
    private readonly Random _random;
    private readonly double _mean;
    private readonly double _stdDev;
    private bool _hasSpareValue;
    private double _spareValue;

    public NormalDistributionGenerator(double mean = 0, double stdDev = 1, int? seed = null)
    {
        _mean = mean;
        _stdDev = stdDev;
        _random = seed.HasValue ? new Random(seed.Value) : new Random();
        _hasSpareValue = false;
        _spareValue = 0;
    }

    // Метод Box-Muller для генерации нормального распределения
    public double Next()
    {
        if (_hasSpareValue)
        {
            _hasSpareValue = false;
            return _spareValue * _stdDev + _mean;
        }

        double u, v, s;
        do
        {
            u = _random.NextDouble() * 2 - 1; // [-1, 1]
            v = _random.NextDouble() * 2 - 1; // [-1, 1]
            s = u * u + v * v;
        } while (s >= 1 || s == 0);

        s = Math.Sqrt(-2.0 * Math.Log(s) / s);
        _spareValue = v * s;
        _hasSpareValue = true;

        return u * s * _stdDev + _mean;
    }

    public double[] GenerateSample(int size)
    {
        var sample = new double[size];
        for (int i = 0; i < size; i++)
        {
            sample[i] = Next();
        }
        return sample;
    }
}

2. Статистические функции для тестов

csharp
public static class StatisticalFunctions
{
    // Функция ошибок (Error Function)
    public static double Erf(double x)
    {
        // approximation using Abramowitz and Stegun formula 7.1.26
        double a1 = 0.254829592;
        double a2 = -0.284496736;
        double a3 = 1.421413741;
        double a4 = -1.453152027;
        double a5 = 1.061405429;
        double p = 0.3275911;

        int sign = x < 0 ? -1 : 1;
        x = Math.Abs(x);

        double t = 1.0 / (1.0 + p * x);
        double y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Math.Exp(-x * x);

        return sign * y;
    }

    // Нормальная функция распределения (CDF)
    public static double NormalCdf(double x, double mean = 0, double stdDev = 1)
    {
        return 0.5 * (1 + Erf((x - mean) / (stdDev * Math.Sqrt(2))));
    }

    // Гамма-функция
    public static double Gamma(double x)
    {
        // Простая реализация через ланчос
        if (x <= 0) throw new ArgumentException("x must be positive");
        
        double[] p = {
            676.5203681218851,
            -1259.1392167224028,
            771.32342877765313,
            -176.61502916214059,
            12.507343278686905,
            -0.13857109526572012,
            9.9843695780195716e-6,
            1.5056327351493116e-7
        };
        
        if (x < 0.5)
            return Math.PI / (Math.Sin(Math.PI * x) * Gamma(1 - x));
        
        x -= 1;
        double a = 0.99999999999980993;
        double t = x + 7.5;
        
        for (int i = 0; i < p.Length; i++)
        {
            a += p[i] / (x + i + 1);
        }
        
        return Math.Sqrt(2 * Math.PI) * Math.Pow(t, x + 0.5) * Math.Exp(-t) * a;
    }

    // Неполная гамма-функция
    public static double IncompleteGamma(double s, double x)
    {
        // Реализация через ряд
        if (x < 0) return double.NaN;
        if (x == 0) return 0;
        
        const int maxIterations = 100;
        const double epsilon = 1e-10;
        
        double term = 1.0 / s;
        double sum = term;
        
        for (int k = 1; k < maxIterations; k++)
        {
            term *= x / (s + k);
            sum += term;
            if (Math.Abs(term) < epsilon) break;
        }
        
        return Math.Pow(x, s) * Math.Exp(-x) * sum;
    }

    // Функция хи-квадрат распределения
    public static double ChiSquareCdf(double x, double degreesOfFreedom)
    {
        if (x <= 0) return 0;
        return IncompleteGamma(degreesOfFreedom / 2, x / 2) / Gamma(degreesOfFreedom / 2);
    }
}

3. Тест Колмогорова-Смирнова

csharp
public class KolmogorovSmirnovTest
{
    public static (double Statistic, double PValue) Test(double[] sample, double mean, double stdDev)
    {
        var sortedSample = sample.OrderBy(x => x).ToArray();
        int n = sortedSample.Length;
        
        double dPlus = 0;
        double dMinus = 0;
        
        for (int i = 0; i < n; i++)
        {
            double empirical = (i + 1.0) / n;
            double theoretical = StatisticalFunctions.NormalCdf(sortedSample[i], mean, stdDev);
            
            dPlus = Math.Max(dPlus, empirical - theoretical);
            dMinus = Math.Max(dMinus, theoretical - (double)i / n);
        }
        
        double d = Math.Max(dPlus, dMinus);
        double pValue = CalculatePValue(d, n);
        
        return (d, pValue);
    }
    
    private static double CalculatePValue(double d, int n)
    {
        // Аппроксимация p-value для теста Колмогорова-Смирнова
        if (n <= 0) return 1;
        
        double sqrtN = Math.Sqrt(n);
        double x = d * (sqrtN + 0.12 + 0.11 / sqrtN);
        
        if (x <= 0) return 1;
        if (x >= 3.0) return 0;
        
        double pValue = 0;
        for (int k = 1; k <= 100; k++)
        {
            pValue += Math.Exp(-2 * k * k * x * x);
        }
        
        return 2 * pValue;
    }
}

4. Тест хи-квадрат

csharp
public class ChiSquareTest
{
    public static (double Statistic, double PValue, int DegreesOfFreedom) Test(
        double[] sample, double mean, double stdDev, int bins = 10)
    {
        if (bins < 2) throw new ArgumentException("Number of bins must be at least 2");
        
        // Определяем границы бинов
        var min = sample.Min();
        var max = sample.Max();
        var binWidth = (max - min) / bins;
        
        // Наблюдаемые частоты
        var observed = new int[bins];
        foreach (var value in sample)
        {
            int binIndex = Math.Min((int)((value - min) / binWidth), bins - 1);
            observed[binIndex]++;
        }
        
        // Ожидаемые частоты для нормального распределения
        var expected = new double[bins];
        for (int i = 0; i < bins; i++)
        {
            double left = min + i * binWidth;
            double right = min + (i + 1) * binWidth;
            expected[i] = sample.Length * (
                StatisticalFunctions.NormalCdf(right, mean, stdDev) - 
                StatisticalFunctions.NormalCdf(left, mean, stdDev)
            );
        }
        
        // Вычисляем статистику хи-квадрат
        double chiSquare = 0;
        for (int i = 0; i < bins; i++)
        {
            if (expected[i] > 0)
            {
                chiSquare += Math.Pow(observed[i] - expected[i], 2) / expected[i];
            }
        }
        
        // Степени свободы: bins - 1 - 2 (два параметра оценены из данных)
        int degreesOfFreedom = bins - 3;
        double pValue = 1 - StatisticalFunctions.ChiSquareCdf(chiSquare, degreesOfFreedom);
        
        return (chiSquare, pValue, degreesOfFreedom);
    }
}

5. Тест Шапиро-Уилка (упрощенная версия)

csharp
public class ShapiroWilkTest
{
    // Коэффициенты для теста Шапиро-Уилка (для n <= 50)
    private static readonly Dictionary<int, double[]> Coefficients = new()
    {
        [3] = new double[] { 0.7071 },
        [4] = new double[] { 0.6872, 0.1677 },
        [5] = new double[] { 0.6646, 0.2413 },
        [6] = new double[] { 0.6431, 0.2806, 0.0875 },
        [7] = new double[] { 0.6233, 0.3031, 0.1401 },
        [8] = new double[] { 0.6052, 0.3164, 0.1743, 0.0561 },
        [9] = new double[] { 0.5888, 0.3244, 0.1976, 0.0947 },
        [10] = new double[] { 0.5739, 0.3291, 0.2141, 0.1224 },
        // Для демонстрации - в реальности нужны все коэффициенты до 50
    };

    public static (double Statistic, double PValue) Test(double[] sample)
    {
        var sorted = sample.OrderBy(x => x).ToArray();
        int n = sorted.Length;
        
        if (n < 3 || n > 50)
        {
            throw new ArgumentException("Shapiro-Wilk test requires sample size between 3 and 50");
        }
        
        if (!Coefficients.ContainsKey(n))
        {
            throw new ArgumentException($"Coefficients not available for sample size {n}");
        }
        
        // Вычисляем статистику W
        double w = CalculateWStatistic(sorted, n);
        
        // Аппроксимируем p-value (упрощенно)
        double pValue = ApproximatePValue(w, n);
        
        return (w, pValue);
    }
    
    private static double CalculateWStatistic(double[] sorted, int n)
    {
        var a = Coefficients[n];
        int k = n % 2 == 0 ? n / 2 : (n - 1) / 2;
        
        double numerator = 0;
        for (int i = 0; i < k; i++)
        {
            numerator += a[i] * (sorted[n - 1 - i] - sorted[i]);
        }
        
        double denominator = 0;
        double mean = sorted.Average();
        foreach (var x in sorted)
        {
            denominator += Math.Pow(x - mean, 2);
        }
        
        return Math.Pow(numerator, 2) / denominator;
    }
    
    private static double ApproximatePValue(double w, int n)
    {
        // Упрощенная аппроксимация p-value
        // В реальной реализации используются табличные значения
        if (w > 0.95) return 0.9;
        if (w > 0.90) return 0.7;
        if (w > 0.85) return 0.5;
        if (w > 0.80) return 0.3;
        if (w > 0.75) return 0.1;
        return 0.01;
    }
}

6. Полная система тестирования

csharp
public class NormalityTester
{
    private readonly double[] _sample;
    private readonly double _mean;
    private readonly double _stdDev;

    public NormalityTester(double[] sample)
    {
        _sample = sample ?? throw new ArgumentNullException(nameof(sample));
        _mean = sample.Average();
        _stdDev = CalculateStdDev(sample, _mean);
    }

    private double CalculateStdDev(double[] values, double mean)
    {
        double sumSquares = 0;
        foreach (var value in values)
        {
            sumSquares += Math.Pow(value - mean, 2);
        }
        return Math.Sqrt(sumSquares / (values.Length - 1));
    }

    public void RunAllTests(double significanceLevel = 0.05)
    {
        Console.WriteLine("=== ТЕСТИРОВАНИЕ НОРМАЛЬНОСТИ РАСПРЕДЕЛЕНИЯ ===");
        Console.WriteLine($"Размер выборки: {_sample.Length}");
        Console.WriteLine($"Выборочное среднее: {_mean:F4}");
        Console.WriteLine($"Выборочное СКО: {_stdDev:F4}");
        Console.WriteLine($"Уровень значимости: {significanceLevel}");
        Console.WriteLine();

        // Тест Колмогорова-Смирнова
        var (ksStat, ksPValue) = KolmogorovSmirnovTest.Test(_sample, _mean, _stdDev);
        Console.WriteLine("1. ТЕСТ КОЛМОГОРОВА-СМИРНОВА:");
        Console.WriteLine($"   Статистика: {ksStat:F4}");
        Console.WriteLine($"   P-value: {ksPValue:F4}");
        Console.WriteLine($"   Вывод: {(ksPValue > significanceLevel ? "Нормальность НЕ отвергается" : "Нормальность ОТВЕРГАЕТСЯ")}");
        Console.WriteLine();

        // Тест хи-квадрат
        var (chiStat, chiPValue, df) = ChiSquareTest.Test(_sample, _mean, _stdDev, bins: 10);
        Console.WriteLine("2. ТЕСТ ХИ-КВАДРАТ:");
        Console.WriteLine($"   Статистика: {chiStat:F4}");
        Console.WriteLine($"   P-value: {chiPValue:F4}");
        Console.WriteLine($"   Степени свободы: {df}");
        Console.WriteLine($"   Вывод: {(chiPValue > significanceLevel ? "Нормальность НЕ отвергается" : "Нормальность ОТВЕРГАЕТСЯ")}");
        Console.WriteLine();

        // Тест Шапиро-Уилка (для небольших выборок)
        if (_sample.Length >= 3 && _sample.Length <= 50)
        {
            try
            {
                var (swStat, swPValue) = ShapiroWilkTest.Test(_sample);
                Console.WriteLine("3. ТЕСТ ШАПИРО-УИЛКА:");
                Console.WriteLine($"   Статистика: {swStat:F4}");
                Console.WriteLine($"   P-value: {swPValue:F4}");
                Console.WriteLine($"   Вывод: {(swPValue > significanceLevel ? "Нормальность НЕ отвергается" : "Нормальность ОТВЕРГАЕТСЯ")}");
            }
            catch (Exception ex)
            {
                Console.WriteLine($"3. ТЕСТ ШАПИРО-УИЛКА: {ex.Message}");
            }
        }
        else
        {
            Console.WriteLine("3. ТЕСТ ШАПИРО-УИЛКА: Неприменим для данной выборки (требуется 3 ≤ n ≤ 50)");
        }

        Console.WriteLine();
        Console.WriteLine("=== СВОДКА ===");
        Console.WriteLine($"Выборка {(IsNormal(ksPValue, chiPValue, significanceLevel) ? "соответствует" : "НЕ соответствует")} нормальному распределению");
    }

    private bool IsNormal(double ksPValue, double chiPValue, double significanceLevel)
    {
        return ksPValue > significanceLevel && chiPValue > significanceLevel;
    }
}

7. Пример использования

csharp
class Program
{
    static void Main()
    {
        // Генерируем выборку нормального распределения
        var generator = new NormalDistributionGenerator(mean: 10, stdDev: 2, seed: 42);
        double[] sample = generator.GenerateSample(1000);

        // Выводим описательную статистику
        Console.WriteLine("Сгенерированная выборка:");
        Console.WriteLine($"Минимум: {sample.Min():F4}");
        Console.WriteLine($"Максимум: {sample.Max():F4}");
        Console.WriteLine($"Среднее: {sample.Average():F4}");
        Console.WriteLine($"Медиана: {CalculateMedian(sample):F4}");
        Console.WriteLine($"СКО: {CalculateStdDev(sample):F4}");
        Console.WriteLine();

        // Тестируем на нормальность
        var tester = new NormalityTester(sample);
        tester.RunAllTests(significanceLevel: 0.05);

        // Дополнительно: тест с заведомо ненормальной выборкой
        Console.WriteLine("\n" + new string('=', 50));
        Console.WriteLine("ТЕСТ С НЕНОРМАЛЬНОЙ ВЫБОРКОЙ (равномерное распределение):");
        
        var uniformSample = GenerateUniformSample(1000, 5, 15);
        var uniformTester = new NormalityTester(uniformSample);
        uniformTester.RunAllTests(significanceLevel: 0.05);
    }

    static double CalculateMedian(double[] values)
    {
        var sorted = values.OrderBy(x => x).ToArray();
        int n = sorted.Length;
        return n % 2 == 0 ? (sorted[n/2 - 1] + sorted[n/2]) / 2 : sorted[n/2];
    }

    static double CalculateStdDev(double[] values)
    {
        double mean = values.Average();
        double sumSquares = values.Sum(x => Math.Pow(x - mean, 2));
        return Math.Sqrt(sumSquares / (values.Length - 1));
    }

    static double[] GenerateUniformSample(int size, double min, double max)
    {
        var random = new Random(42);
        var sample = new double[size];
        for (int i = 0; i < size; i++)
        {
            sample[i] = min + random.NextDouble() * (max - min);
        }
        return sample;
    }
}

Особенности реализации:

  1. Box-Muller алгоритм - эффективная генерация нормальных чисел

  2. Три статистических теста - комплексная проверка нормальности

  3. Обработка граничных случаев - проверка размеров выборок

  4. Визуализация результатов - понятные выводы

  5. Сравнение с ненормальным распределением - демонстрация работы тестов

Это полная система для генерации и тестирования нормальных распределений! 🎯


на
Ярлыки: , , , , , , , ,

Комментариев нет:

Отправить комментарий

Подписаться на: Комментарии к сообщению (Atom)