Ремарка:

В данной статье представлен авторский метод нахождения t-критерия, разработанный в процессе решения более крупной задачи. Этот метод не претендует на эффективность или универсальность в различных условиях применения, а представляет собой исключительно личный подход к решению конкретной задачи. Автор делится своим опытом и методологией работы, не утверждая их общепризнанной пригодности или оптимальности в других контекстах.

t-критерий:

t-критерий является ключевым инструментом статистического анализа, используемым для проверки гипотез о средних значениях двух выборок или для оценки значимости различий между группами данных. Этот критерий основан на распределении Стьюдента и позволяет определить, насколько значимы различия между выборочными средними, учитывая их стандартные ошибки. Значение t-критерия вычисляется как отношение разности между средними значениями выборок к их стандартной ошибке разности.

(t-распределение) используется в статистике для оценки доверительных интервалов и проверки гипотез о средних значениях, когда объем выборки мал или параметры генеральной совокупности неизвестны. Функция плотности вероятности t-распределения определяется следующей формулой:

f(t, n) = \frac{Г(\frac{n+1}{2})}{\sqrt{n\pi}Г({\frac{n}{2}})}{(1+\frac{t^2}{n})}^{-{\frac{n+1}{2}}}

где t - t-критерй, n - число степеней свободы

Основные факторы, от которых зависит t-критерий:

  1. Уровень значимости α: Уровень значимости определяет вероятность ошибки первого рода (отвержение верной нулевой гипотезы).

  2. Число степеней свободы n: Это параметр, зависящий от размера выборки и определяющий форму распределения Стьюдента.

то есть в ствою очередь t - это можно считать (не в строгом смысле) как функцию от α и n : t = f(α,n)

вырозить програмно эту зависимость я ставлю в данной статье своей задачей в виде функции ta.

//"интерфейс" этой функции 
def ta(a: Double, n: Int): Double = ??? 
// где a и n доверительный интервал и число степеней свободы соответственно

для нахождения t-критерия мне понадобилось два алгоритма ACM295 и ACM209

ACM209:

Алгоритм ACM209 представляет собой численное приближение для вычисления функции распределения Гаусса, также известной как нормальное распределение. Этот алгоритм был предложен для быстрого и точного вычисления значений интеграла от функции плотности вероятности нормального распределения. даный алгоритм всплывает в данном контексте поскольку t-распределение и распределения гауса "связаны" и в предельном случе при n -> ∞t-распределение стремиться к гаусовскому и в последсвии все сведется к вычислению разницы между t-распределением и гаусовым.

//ACM 209 реализация на scala
def ACM209(z: Double): Double = {
	var y: Double = 0.0
	var p: Double = 0.0
	var w: Double = 0.0

	if (z == 0.0)
		p = 0.0
	else {
		y = Math.abs(z) / 2
		if (y >= 3.0) {
			p = 1.0
		} else if (y < 1.0) {
			w = y * y
			p = ((((((((0.000124818987 * w
				- 0.001075204047) * w + 0.005198775019) * w
				- 0.019198292004) * w + 0.059054035642) * w
				- 0.151968751364) * w + 0.319152932694) * w
				- 0.5319230073) * w + 0.797884560593) * y * 2.0
		} else {
			y = y - 2.0
			p = (((((((((((((-0.000045255659 * y
				+ 0.00015252929) * y - 0.000019538132) * y
				- 0.000676904986) * y + 0.001390604284) * y
				- 0.00079462082) * y - 0.002034254874) * y
				+ 0.006549791214) * y - 0.010557625006) * y
				+ 0.011630447319) * y - 0.009279453341) * y
				+ 0.005353579108) * y - 0.002141268741) * y
				+ 0.000535310849) * y + 0.999936657524
		}
	}

	if (z > 0.0)
		(p + 1.0) / 2
	else
		(1.0 - p) / 2
}

Реализация алгоритма на языке Scala, представленная выше, включает несколько ключевых шагов.

Сначала вычисляется абсолютное значение входного параметра z, и в зависимости от его величины алгоритм выбирает одну из нескольких веток вычислений. Если y (половина абсолютного значения z) больше или равно 3, вероятность р устанавливается равной 1. В случае y меньше 1, применяется полиномиальная аппроксимация для вычисления p с использованием ряда коэффициентов. Для значений y между 1 и 3 используется другая полиномиальная аппроксимация, включающая другой набор коэффициентов. Наконец, в зависимости от знака z, полученное значение p корректируется для возв��ащения результата, представляющего собой значение функции распределения Гаусса для данного z.

ACM395:

Алгоритм ACM395 реализует численное приближение для вычисления значений функции распределения t-статистики. Эта функция возвращает параметр p :уровнем значимости или же вероятность ошибки первого рода. предстовляет собой функцию от числа степеней свободы n и t-критерия: p = f(t,n)

//ACM 395 моя реализация на scala
def ACM395(t: Double, n: Int): Double = {

	var a = 0.0
	var b = 0.0
	var y = 0.0

    val tSq = t * t
	y = tSq / n
	b = y + 1.0
	if (y > 1.0E-6) y = Math.log(b)
	a = n - 0.5
	b = 48.0 * a * a
	y = a * y

	y = (((((-0.4 * y - 3.3) * y - 24.0) * y - 85.5) /
		(0.8 * y * y + 100.0 + b) + y + 3.0) / b + 1.0) *
		Math.sqrt(y)
	2.0 * ACM209(-y)
}

После вычисления p, функция ACM395 возвращает 2 * ACM 209(-y), что представляет собой вычисление значения функции распределения t-статистики с использованием функции ACM209 для вычисления вероятности.

Вычисление доверительной вероятности

α собой предстовляет α = 1 - p = 1 - f(t,n) = f'(t,n)

// в общем тут нечего коментрировать
def alpfa(ta: Double, n: Int): Double = {
  1 - ACM395(ta, n)
}

Вычисляем t-критерий:

α принмает значения от 0 до 1 и при n = const α = f'(t,n = const) = α(t) выполняется следующие свойство: α(t1) < α(t2) при t1 < t2

так как мы ищем ищем значение t при определеном α0 = const то введем еще функцию α'(t,n) = α(t,n) - α0 и при фиксированом n = const вырождается в функцию одной переменой α'(t, n = const) = α'(t) и тогда α'(t) примает значения (-α0, 1 -α0) и сохраняется свойство: α'(t1) < α'(t2) при t1 < t2


//так я програмно отразила паралельный перенос
def alphaParallel(fun:(Double, Int) => Double, delta: Double): (Double, Int) => Double = {
  (t: Double, n: Int) => {
	fun(t, n) - delta
  }
}

//так я програмно вырозила вырождение функции двух переменых в функциию одной 
//n = const
//то что я и назвала фиксированным n
def alpfaConst(fun: (Double, Int) => Double, n: Int): Double => Double = {
  t: Double => fun(t, n)
}

val parAlpha = alphaParallel(alpfa, a0)
val conAlpha = alpfaConst(parAlpha, n) // собственно α'(t)

выше обозначеный свойства позволяют применить теорему Больцено-Коши

и найти t = α'^{-1}(α)

но для этого нужно знать интервал в котором находиться t-критерий для α0
и делаю это перебором интервалов по длиной в 1 пока не найду необходимый

//метод для нахождения интервала
def findInterval(fun: Double => Double): (Double, Double) = {

  @tailrec
  def finder(fun: Double => Double, min: Double, max: Double): (Double, Double) = {
    
    fun(min)*fun(max) match {
      case value if value < 0 => (min, max)
	  case value if value > 0 => finder(fun, min + 1, max + 1)
    }


  }

  finder(fun, 0, 1)

}

val interval = findInterval(conAlpha)

когда интервал найден можно применить метод половинго деления

//метод реализующий половинное деление
def bisection(tol: Double = 1e-10): (Double => Double, (Double, Double)) => Double = {

  (f: Double => Double, interval: (Double, Double)) => {

    @tailrec
    def loop(interval: (Double, Double)): Double = {
	
      val (a, b) = interval
      val midpoint = (a + b) / 2
      val fMid = f(midpoint)

      if (fMid == 0 || (b - a) / 2 < tol) midpoint		
      else if (f(a) * fMid < 0) loop(a, midpoint)
      else loop(midpoint, b)
				
    }
    
    loop(interval)
    
  }		
  
}

bisection()(conAlpha, interval)

и того суммарно код метода выглядит так

def ta(a0: Double, n: Int): Double = {

  def alphaParallel(fun:(Double, Int) => Double, delta: Double): (Double, Int) => Double = {
	(t: Double, n: Int) => {
      fun(t, n) - delta
	}	
  }

  def alpfaConst(fun: (Double, Int) => Double, n: Int): Double => Double = {
    t: Double => fun(t, n)
  }
  
  def findInterval(fun: Double => Double): (Double, Double) = {

    @tailrec
    def finder(fun: Double => Double, min: Double, max: Double): (Double, Double) = {
    
      fun(min)*fun(max) match {
        case value if value < 0 => (min, max)
    	case value if value > 0 => finder(fun, min + 1, max + 1)
      }
    }

    finder(fun, 0, 1)
  }
  
  def bisection(tol: Double = 1e-10): (Double => Double, (Double, Double)) => Double = {

  (f: Double => Double, interval: (Double, Double)) => {

    @tailrec
    def loop(interval: (Double, Double)): Double = {
	
      val (a, b) = interval
      val midpoint = (a + b) / 2
      val fMid = f(midpoint)

      if (fMid == 0 || (b - a) / 2 < tol) midpoint		
      else if (f(a) * fMid < 0) loop(a, midpoint)
      else loop(midpoint, b)
				
    }
    
    loop(interval)
    
  }		


  val parAlpha = alphaParallel(alpfa, a0)
  val conAlpha = alpfaConst(parAlpha, n)

  val interval = findInterval(conAlpha)

  bisection()(conAlpha, interval)

}

Результат:

сравним табличные значения и те что нагенерил мой метод

n \ α0

0.95

0.99

0.999

1

12,70

 63,65

 636,61

2

4,303

 9,925

 31,602

3

3,182

 5,841

 12,923

4

2,776

 4,604

 8,610

5

2,571

 4,032

 6,869

6

2,447

 3,707

 5,959

7

2,365

 3,499

 5,408

8

2,306

 3,355

 5,041

9

2,262

 3,250

 4,781

10

2,228

 3,169

 4,587

и вот то что мой метод выдал:

заметим что чем больше n и меньше α0 тем меньше разница с табличным значением.

P.S. код можно глянуть тут