---

Июль 2015

ТОМ 30 ВЫПУСК 7

Язык программирования R - Введение в R для программистов на C#

Джеймс Маккафри

Продукты и технологии:

R, C#, Visual Studio

В статье рассматриваются:

• пример проверки по хи-квадрату (chi-square test) с применением R;

• линейный регрессионный анализ;

• проверка по хи-квадрату независимости;

• графы в R;

• пример проверки по хи-квадрату на C# для сравнения.

Исходный код можно скачать по ссылке

Язык R используется специалистами по обработке и анализу данных и программистами для статистический вычислений. R является одной из самых быстро развивающихся технологий среди моих коллег, использующих C#. Отчасти это вызвано нарастающими объемами данных, собираемых программными системами, и потребностью в анализе этих данных. Знакомство с R может оказаться ценным пополнением в вашем арсенале технических навыков.

Язык R — это проект по лицензии GNU и бесплатное программное обеспечение. R произошел от языка S (сокращение от статистики), который был создан в Bell Laboratories в 1970-х годах. По R есть много великолепных онлайновых учебных пособий, но в большинстве из них предполагается, что вы — студент университета, изучающий статистику. Эта статья нацелена на то, чтобы помочь программистам на C# быстро освоить R.

Лучший способ понять, о чем эта статья, — взглянуть на пример сеанса работы с R на рис. 1. В примере сеанса демонстрируются две несвязанные между собой тематики. Первый набор команд показывает проверку по хи-квадрату (chi-square test), также называемую хи-квадратичной проверкой (chi-squared test), для равномерного распределения. Второй набор команд иллюстрирует пример линейной регрессии, которая, по моему мнению, является чем-то вроде «Hello World» в мире статистических вычислений.

Рис. 1. Пример сеанса работы с R

Веб-сайт R находится по адресу r-project.org. На этом сайте есть ссылки на несколько зеркальных сайтов, где можно скачать и установить R. Установка заключается в простом запуске самораспаковывающегося исполняемого файла. R официально поддерживается в Windows XP и более поздних версиях, а также на большинстве платформ, отличных от Windows. Я устанавливал R на компьютеры с Windows 7 и Windows 8 безо всяких проблем. По умолчанию процесс установки предоставляет как 32-, так и 64-разрядную версию.

В этой статье предполагается, что вы умеете программировать на C# хотя бы на среднем уровне (чтобы понять пояснения сходств и различий между C# и R), но ничего не знаете о R. Демонстрационная программа на C# представлена во всей ее полноте; кроме того, вы найдете весь исходный код в сопутствующем этой статье пакете кода.

Проверка по хи-квадрату с помощью R

Взгляните на рис. 1 и первым делом обратите внимание на то, что использование R весьма заметно отличается от использования C#. Хотя R позволяет писать скрипты, он чаще всего применяется в интерактивном режиме в командной оболочке. Первый пример R — это анализ того, является обычная шестигранная игральная кость правильной или нет. При многократных бросках правильная кость должна давать примерно одинаковое количество каждого из шести возможных результатов.

Строка приглашения R в командной оболочке помечается маркером >. Первое выражение, введенное на рис. 1, начинается с символа #, который отмечает комментарий.

Первая реальная команда в примере выглядит так:

> observed <- c(20, 28, 12, 32, 22, 36)

Она создает вектор observed, использующий функцию c (сокращение от concatenate). Вектор хранит объекты одного и того же типа данных. Примерный эквивалент этого выражения на C# был бы таким:

var observed = new int[] { 20, 28, 12, 32, 22, 36 };

Первое значение, 20, — число, указывающее, сколько раз встречалось одно очко, 28 — два очка и т. д. Сумма шести значений составляет 20 + 28 + . . + 36 = 150. От правильной игральной кости следовало бы ожидать приблизительно 150/6 = 25 раз для каждого из возможных результатов. Но наблюдаемое число выпадений трех очков (12) выглядит подозрительно малым.

После создания вектора проверка по хи-квадрату выполняется с помощью функции chisq.test:

> chisq.test(observed)

В языке R при формировании имен переменных и функций вместо символа подчеркивания (_) часто используется точка (.) — такие имена легче читать. Результат вызова функции chisq.test выдает приличное количество текста:

Chi-squared test for given probabilities

data:  observed
X-squared = 15.28, df = 5, p-value = 0.009231

В терминах C# большинство функций R возвращает структуру данных, которую можно игнорировать, но также содержит много выражений, эквивалентных Console.WriteLine, которые предоставляют вывод. Заметьте, что расшифровка смысла вывода R возлагается на вас. Далее в этой статье я покажу, как создать эквивалент проверки по хи-квадрату на чистом (без библиотек) C#.

В этом примере значение «X-squared», равное 15.28, является вычисленной статистикой хи-квадрата (греческая буква «chi» напоминает прописную «X»). Значение 0.0 указывает, что наблюдавшиеся значения точно соответствуют тому, что следовало бы ожидать, если игральная кость была правильной. Более высокие значения хи-квадрата информируют о растущей вероятности того, что наблюдаемые числа не могли бы выпасть случайно, будь эта игральная кость правильной. Значение df, равное 5, — это «число степеней свободы», которое на единицу меньше, чем число наблюдавшихся значений. В этом анализе df не слишком важно.

P-значение (0.009231) — это вероятность того, что наблюдаемые счетчики могли бы быть случайными, если бы игральная кость была правильной. Поскольку p-значение столь мало, менее 1%, вы должны считать, что наблюдаемые значения вряд ли выпадали случайно, а значит, у вас есть статистическое доказательство того, что данная игральная кость скорее всего имеет смещение центра тяжести.

Линейный регрессионный анализ с помощью R

Второй набор выражений на рис. 1 показывает пример линейной регрессии. Линейная регрессия — статистический метод, применяемый для описания взаимосвязи между числовой переменной (называемой в статистике зависимой переменной) и одной или более пояснительными переменными (explanatory variables) (называемыми в статистике независимыми переменными), которые могут быть либо числовыми, либо категориальными. При наличии только одной независимой пояснительной/прогностической переменной метод называют простой линейной регрессией. Когда имеются две или более независимых переменных, как в демонстрационном примере, метод называют множественной линейной регрессией (multiple linear regression).

Прежде чем заняться анализом линейной регрессии, я создал текстовый файл DummyData.txt с восемью элементами, разделенными запятыми, и поместил его в каталог C:\IntroToR. Вот содержимое этого файла:

Color,Length,Width,Rate
blue, 5.4, 1.8, 0.9
blue, 4.8, 1.5, 0.7
blue, 4.9, 1.6, 0.8
pink, 5.0, 1.9, 0.4
pink, 5.2, 1.5, 0.3
pink, 4.7, 1.9, 0.4
teal, 3.7, 2.2, 1.4
teal, 4.2, 1.9, 1.2

Предполагается, что этот файл представляет данные о цвете, длине и ширине лепестков и показателе роста цветков. Идея в том, чтобы прогнозировать значения Rate (последний столбец) на основе значений Color, Length и Width. После выражения с комментарием первые три команды R в анализе линейной регрессии выглядят следующим образом:

> setwd("C:\\IntroToR")
> data <- read.table("DummyData.txt",
  header=TRUE, sep=",")
> print(data)

Первая команда устанавливает рабочий каталог, чтобы мне не приходилось указывать полный путь к файлу с исходными данными. Вместо привычного в C# символа (\\) я мог бы использовать символ (/), распространенный на платформах, отличных от Windows.

Вторая команда загружает данные в память в объект-таблицу с именем data. Заметьте, что R использует именованные параметры. Параметр header сообщает, является первая строка информацией о заголовке (TRUE, или T в сокращенной форме) или нет (FALSE, или F). R чувствителен к регистру букв, и для присваивания значений обычно можно использовать либо оператор <-, либо оператор =. Выбор является в основном делом личных предпочтений. Я использую <- для объектных присваиваний и = для присваивания значений параметрам.

Параметр sep (separator) указывает, как разделяются значения в каждой строке. Например, \t задавал бы разделение значений табуляторами, а (" ") — пробелами.

Функция print выводит таблицу данных в памяти. У этой функции много необязательных параметров. Заметьте, что вывод на рис. 1 отражает индексы элементов данных, начиная с 1. В случае индексов массивов, матриц и объектов R является языком, где отсчет ведется от 1, а не от 0, как в языке C#.

Анализ линейной регрессии выполняется следующими двумя командами R:

> model <- lm(data$Rate ~ (data$Color + data$Length + data$Width))
> summary(model)

Первую команду можно интерпретировать как «Сохранить в объект с именем model результат функции анализа lm (линейная модель), где прогнозируемая зависимая переменная является столбцом Rate в объекте-таблице (data$Rate), а независимые переменные-предикторы представляют собой Color, Length и Width». Вторая команда означает: «Отображать только базовые результаты анализа, сохраненные в объекте с именем model».

Функция lm генерирует большое количество информации. Допустим, вы хотите предсказать значение Rate, когда входными значениями являются Color = pink, Length = 5.0 и Width = 1.9. (Заметьте, что это соответствует элементу данных [4] с реальным значением Rate, равным 0.4.) Для прогнозирования вы должны использовать значения в столбце Estimate:

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)    -0.14758    0.48286  -0.306  0.77986
data$Colorpink -0.49083    0.04507 -10.891  0.00166 **
data$Colorteal  0.35672    0.09990   3.571  0.03754 *
data$Length     0.04159    0.07876   0.528  0.63406
data$Width      0.45200    0.11973   3.775  0.03255 *

Если X представляет значения независимой переменной, а Y представляет прогнозируемый Rate, тогда:

X = (blue = NA, pink = 1, teal = 0, Length = 5.0, Width = 1.9)
Y = -0.14758 + (-0.49083)(1) + (0.35672)(0) + (0.04159)(5.0) + (0.45200)(1.9)
  = -0.14758 + (-0.49083) + (0) + (0.20795) + (0.85880)
  = 0.42834

Заметьте, что предсказанный Rate (0.43) довольно близок к реальному (0.40).

Если на словах, то, чтобы сделать прогноз на основе модели, вы вычисляете линейную сумму произведений значений Estimate и соответствующим им значений X. Значение Intercept является константой, не связанной ни с какой переменной. Когда вы имеете дело с категориальными пояснительными переменными, одно из значений отбрасывается (в данном случае — blue).

Информация в нижней части вывода сообщает, насколько хорошо независимые переменные Color, Length и Width поясняют зависимую переменную Rate:

Residual standard error: 0.05179 on 3 degrees of freedom
Multiple R-squared:  0.9927,    Adjusted R-squared:  0.9829
F-statistic: 101.6 on 4 and 3 DF,  p-value: 0.00156

Значение Multiple R-squared (коэффициент детерминации) (0.9927) — это процентная доля вариации зависимой переменной, поясняемая линейной комбинацией независимых переменных. Выражаясь несколько иначе, R-squared — это значение между 0 и 1, где более высокие значения свидетельствуют о более качественной модели прогнозирования. Здесь значение R-squared чрезвычайно высокое, указывая на то, что Color, Length и Width могут очень точно предсказывать Rate. F-statistic (статистика дисперсии) с учетом значения R-squared и p-value — другие критерии точности модели.

Один из моментов, демонстрируемых этим примером, заключается в том, что при программировании с применением R самая крупная проблема — понимание статистики, стоящей за языковыми функциями. Большинство людей осваивают R постепенно, добавляя в копилку своих знаний по одному методу единовременно — по мере необходимости найти ответ на некий специфический вопрос. Аналогией в C# могло бы послужить изучение различных объектов-наборов в пространстве имен Collections.Generic, таких как классы Hashtable и Queue. Большинство разработчиков изучают по одной структуре данных, не пытаясь сходу запомнить информацию обо всех классах до их использования в каком-либо проекте.

Другая проверка по хи-квадрату

Тип проверки по хи-квадрату на рис. 1 часто называют тестом на равномерное распределение, поскольку он проверяет, все ли наблюдаемые данные имеют равные счетчики, т. е. равномерно ли распределяются данные. Есть несколько других видов проверок по хи-квадрату, включая проверку независимости по хи-квадрату.

Допустим, есть группа из 100 людей и вас интересует, является ли пол (male, female) независимым от принадлежности к политической партии (Democrat, Republican, Other). Представьте, что ваши данные находятся в матрице сопряженности (contingency matrix), как показано в табл. 1. Из нее следует, что мужчины (males), возможно, чаще являются республиканцами (Rep), чем женщины (females).

Табл. 1. Матрица сопряженности

	Dem	Rep	Other
Male	15	25	10	50
Female	30	15	5	50
	45	40	15	100

Чтобы использовать R для проверки того, являются ли два фактора (пол и принадлежность к партии) статистически независимыми, вы должны сначала поместить данные в численную матрицу этой командой:

> cm <- matrix( c(15,30,25,15,10,5), nrow=2, ncol=3 )

Обратите внимание, что в R матричные данные сохраняются по столбцам (сверху вниз и слева направо), а не по строкам (слева направо и сверху вниз), как в C#.

Команда проверки по хи-квадрату в R выглядит следующим образом:

> chisq.test(cm)

Результат p-value после проверки — 0.01022, и это указывает на то, что при уровне значимости (significance level) в 5% два фактора не являются независимыми. Иначе говоря, между полом и партийной принадлежностью существует статистическая взаимосвязь.

В первом примере с игральной костью входной параметр был вектором, а во втором примере со связью между полом и партийной принадлежностью — матрицей. Функция chisq.test имеет всего семь параметров, из которых один обязательный (вектор или матрица), а последующие шесть являются необязательными именованными.

Подобно многим C#-методам в пространствах имен Microsoft .NET Framework большинство функций R имеет множество перегруженных версий. В C# перегрузка обычно реализуется использованием нескольких методов с одинаковыми именами, но с разными параметрами. Применение обобщений также представляет собой одну и форм перегрузки. В языке R перегрузка реализуется использованием одной функции с несколькими необязательными именованными параметрами.

Построение графиков с помощью R

Как программист на C#, когда я хочу построить график на основе выходных данных от какой-то программы, я обычно запускаю эту программу, копирую выходные данные, вставляю их в блокнот нажатием клавиш Ctrl+V для удаления управляющих символов, копирую эти данные, вставляю в Excel, а затем создаю график с помощью Excel. Довольно замысловато, но отлично работает в большинстве случаев.

Одна из сильных сторон системы R — ее способность генерировать графики. Взгляните на пример на рис. 2. Этот тип графика невозможно создать в Excel без установки какой-нибудь надстройки.

Рис. 2. Трехмерный график, созданный с использованием R

В дополнение к программе оболочки, показанной на рис. 1, R также имеет некое подобие GUI-интерфейса, RGui.exe, для использования при построении графиков (я называю этот интерфейс «полу-GUI»).

График на рис. 2 отражает функцию z = f(x,y) = x * e^(-(x^2 + y^2)). Вот как выглядят первые три команды R, генерирующие график:

> rm(list=ls())
> x <- seq(-2, 2, length=25)
> y <- seq(-2, 2, length=25)

Функция rm удаляет объект из текущего рабочего пространства в памяти. Использованная команда является магическим заклинанием R для удаления всех объектов. Вторая и третья команды создают векторы из 25 равноотстоящих значений — от –2 до +2 включительно. Вместо этого можно было бы использовать функцию c.

Следующие две команды:

> f <- function(x,y) { x * exp(-(x^2
  + y^2)) }
> z <- outer(x,y,f)

Первая из них показывает, как определить функцию в R, используя ключевое слово function. Встроенная в R функция outer создает матрицу значений на основе векторов x и y и определения функции f. Результатом является матрица 25×25, где значение в каждой ячейке представляет собой значение функции f, которое соответствует x и y.

Еще две команды выглядят так:

> nrz <- nrow(z)
> ncz <- ncol(z)

Функции nrow и ncol возвращают количество строк или число столбцов в своем аргументе-матрице. Здесь оба значения будут по 25.

Следующая команда R использует функцию colorRampPalette для создания пользовательской палитры с цветовым градиентом, которой рисуется график:

> jet.colors <- colorRampPalette(c("midnightblue", "blue",
+ "cyan", "green", "yellow", "orange", "red", "darkred"))

При наборе в R длинной строки, если нажать клавишу Enter, система переведет курсор на следующую строку и поместит символ +, указывающий, что команда еще не закончена. Затем следуют две команды:

> nbcol <- 64
> color <- jet.colors(nbcol)

Их результат — вектор color, содержащий 64 разных цветовых значений в диапазоне от очень темно-синего до зеленого, желтого и очень темного красного. Другие две команды:

> zfacet <- z[-1,-1] + z[-1,-ncz] + z[-nrz,-1] + z[-nrz,-ncz]
> facetcol <- cut(zfacet,nbcol)

Если внимательно посмотреть на рис. 2, видно, что поверхность состоит из 25×25 (625) малых квадратов, или фасетов (facets). Две предыдущие команды создают матрицу 25×25, где значение в каждой ячейке является одним из 64 цветов, используемым для соответствующего фасета поверхности.

Наконец, трехмерный график отображается функцией persp (сокращение от perspective graph [график с перспективой]):

> persp(x,y,z, col=color[facetcol], phi=20, theta=-35,
+ ticktype="detailed", d=5, r=1, shade=0.1, expand=0.7)

У функции persp много необязательных именованных параметров. Параметр col — цвет (или цвета) для использования. Параметры phi и theta задают угол зрения (слева и справа, а также сверху и снизу) на график. Параметр ticktype управляет тем, как отображаются значения по осям x, y и z. Параметры r и d управляют кажущимся расстоянием от глаз до графика и кажущимся эффектом объемности. Параметр shade определяет имитацию затенения от виртуального источника света. Параметр expand указывает соотношение высоты и ширины графика. И это не все параметры. У функции persp гораздо больше параметров, но описанных и примененных здесь вполне достаточно в большинстве ситуаций.

Этот пример демонстрирует, что R имеет очень мощные и гибкие встроенные возможности в графике, но они сравнительно низкоуровневые и требуют немалого труда.

Проверка по хи-квадрату на C#

Чтобы получить представление о сходствах и различиях между анализом на языке R и программированием на языке C#, полезно изучить C#-реализацию проверки по хи-квадрату. Кроме того, этот код на C# может оказаться полезным пополнением вашей личной библиотеки кода. Взгляните на демонстрационную программу на C# (рис. 3).

Рис. 3. Проверка по хи-квадрату с использованием C#

Демонстрационная программа приблизительно соответствует проверке игральных костей по хи-квадрату на языке R, показанной на рис. 1. Заметьте, что выходные значения демонстрации на C# ровно те же самые, что и в сеансе работы с R.

Чтобы создать демонстрационную программу, я запустил Visual Studio и выбрал шаблон C# Console Application. Я назвал проект ChiSquare. После загрузки кода шаблона в редактор я переименовал в окне Solution Explorer файл Program.cs в более описательный ChiSquareProgram.cs, и Visual Studio автоматически переименовала класс Program за меня.

Полный исходный код демонстрационной программы на C# представлен на рис. 4. Вы заметите, что исходный код длиннее, чем можно было бы ожидать. Реализация статистического программирования на чистом C# не является сверхтрудной, но код получается длинный.

Рис. 4. Демонстрационная программа проверки по хи-квадрату на C#

using System;
namespace ChiSquare
{
  class ChiSquareProgram
  {
    static void Main(string[] args)
    {
      try
      {
        Console.WriteLine(
          "\nBegin Chi-square test using C# demo\n");
        Console.WriteLine("Goal is to see if one die from a set
          of dice is biased or not\n");

        int[] observed = new int[] { 20, 28, 12, 32, 22, 36 };

        Console.WriteLine("\nStarting chi-square test");
        double p = ChiSquareTest(observed);
        Console.WriteLine("\nChi-square test complete");

        double crit = 0.05;
        if (p < crit)
        {
          Console.WriteLine(
            "\nBecause p-value is below critical value of " +
            crit.ToString("F2"));
          Console.WriteLine(
            "the null hypothsis is rejected and we conclude");
          Console.WriteLine("the data is unlikely to have
            happened by chance.");
        }
        else
        {
          Console.WriteLine(
            "\nBecause p-value is not below critical value of "
            + crit.ToString("F2"));
          Console.WriteLine("the null hypothsis is accepted
            (not rejected) and we conclude");
          Console.WriteLine("the observed data could have
            happened by chance.");
        }

        Console.WriteLine("\nEnd\n");
        Console.ReadLine();
      }
      catch (Exception ex)
      {
        Console.WriteLine(ex.Message);
        Console.ReadLine();
      }

    } // Main

    static void ShowVector(int[] v)
    {
      for (int i = 0; i < v.Length; ++i)
        Console.Write(v[i] + " ");
      Console.WriteLine("");
    }

    static double ChiSquareTest(int[] observed)
    {
      Console.WriteLine("Observed frequencies are: ");
      ShowVector(observed);

      double x = ChiSquareStatistic(observed);
      Console.WriteLine("\nX-squared = " + x.ToString("F2"));

      int df = observed.Length - 1;
      Console.WriteLine("\ndf = " + df);

      double p = ChiSquareProb(x, df);
      Console.WriteLine("\np-value = " + p.ToString("F6"));

      return p;
    }

    static double ChiSquareStatistic(int[] observed)
    {
      double sumObs = 0.0;
      for (int i = 0; i < observed.Length; ++i)
        sumObs += observed[i];
      double expected = (int)(sumObs / observed.Length);

      double result = 0.0;
      for (int i = 0; i < observed.Length; ++i)
      {
        result += ((observed[i] - expected) *
         (observed[i] - expected)) / expected;
      }
      return result;
    }

    public static double ChiSquareProb(double x, int df)
    {
      // x = a – вычисленное значение хи-квадрата, df = степени
      // свободы, вывод = prob, значение x получилось случайно.
      // Например, если результат < 0.05, существует лишь 5%-я
      // вероятность того, что значение x выпало случайно,
      // а значит, мы заключаем, что реальные данные,
      // на основе которых было получено x, не отвечают
      // ожидаемым. Эту функцию можно использовать для создания
      // процедуры ChiSquareTest.
      // Это ACM Algorithm 299 и обновление ACM TOMS June 1985.
      // Использует собственную функцию Exp(),
      // показанную ниже.

      if (x <= 0.0 || df < 1)
        throw new Exception("parameter x must be positive " +
          "and parameter df must be 1 or greater
          in ChiSquaredProb()");

      double a = 0.0; // 299 имен переменных
      double y = 0.0;
      double s = 0.0;
      double z = 0.0;
      double e = 0.0;
      double c;

      bool even; // значение df четное?

      a = 0.5 * x;
      if (df % 2 == 0) even = true; else even = false;

      if (df > 1) y = Exp(-a); // ACM update remark (4)

      if (even == true) s = y; else s = 2.0 *
        Gauss(-Math.Sqrt(x));

      if (df > 2)
      {
        x = 0.5 * (df - 1.0);
        if (even == true) z = 1.0; else z = 0.5;
        if (a > 40.0) // примечание ACM (5)
        {
          if (even == true) e = 0.0;
          else e = 0.5723649429247000870717135;
          c = Math.Log(a); // логарифм по основанию e
          while (z <= x)
          {
            e = Math.Log(z) + e;
            s = s + Exp(c * z - a - e); // примечание
                                        // в обновлении ACM (6)
            z = z + 1.0;
          }
          return s;
        } // a > 40.0
        else
        {
          if (even == true) e = 1.0;
          else e = 0.5641895835477562869480795 / Math.Sqrt(a);
          c = 0.0;
          while (z <= x)
          {
            e = e * (a / z); // примечание в обновлении ACM (7)
            c = c + e;
            z = z + 1.0;
          }
          return c * y + s;
        }
      } // df > 2
      else
      {
        return s;
      }
    } // ChiSquare()

    private static double Exp(double x) // примечание
                                        // в обновлении ACM (3)
    {
      if (x < -40.0) // примечание в обновлении ACM (8)
        return 0.0;
      else
        return Math.Exp(x);
    }

    public static double Gauss(double z)
    {
      // Ввод = z-value (-inf to +inf),
      // вывод = p при кривой нормального распределения
      // от -inf до z, например если z = 0.0, функция
      // возвращает 0.5000.
      // ACM Algorithm 209.
      double y; // scratch-переменная 209
      double p; // результат, называемый 'z' в 209
      double w; // scratch-переменная 209

      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.531923007300) * w + 0.797884560593) * y * 2.0;
        }
        else
        {
          y = y - 2.0;
          p = (((((((((((((-0.000045255659 * y
            + 0.000152529290) * y - 0.000019538132) * y
            - 0.000676904986) * y + 0.001390604284) * y
            - 0.000794620820) * 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)
        return (p + 1.0) / 2;
      else
        return (1.0 - p) / 2;
    } // Gauss()

  } // класс Program

} // ns

Метод Main состоит в основном из выражений WriteLine. Значимый код таков:

int[] observed = new int[] { 20, 28, 12, 32, 22, 36 };
double p = ChiSquareTest(observed);
double crit = 0.05;
if (p < crit) {
  // Сообщения
} else {
  // Сообщения
}

Проверка по хи-квадрату на C# принимает массив наблюдаемых значений, вычисляет статистическое значение хи-квадрата и отображает его, затем вычисляет и показывает количество степеней свободы, а также вычисляет и возвращает p-value.

Метод ChiSquareTest вызывает три вспомогательных метода:

ShowVector(observed);
double x = ChiSquareStatistic(observed);
int df = observed.Length - 1;
double p = ChiSquareProb(x, df);

Метод ShowVector отображает входной вектор по аналогии с тем, как это делается R-функцией chisq.test. Метод ChiSquareStatistic возвращает вычисленный хи-квадрат («X-squared» в выводе R), а метод ChiSquareProb использует хи-квадрат, возвращенный ChiSquareStatistic, для вычисления вероятности («p-value» в выводе R).

Метод ChiSquareStatistic — это простая проверка на равномерное распределение. Уравнение для статистики хи-квадрата:

chi-squared = Sum( (observed - expected)^2 / expected )

Поэтому, прежде чем вычислять chi-squared, вы должны рассчитать ожидаемые значения, связанные с наблюдаемыми. Как пояснялось ранее, для этого вы накапливаете значение счетчика в массиве observed (150 в данной демонстрации), затем делите сумму на количество значений в массиве (6 в данной демонстрации), что дает ожидаемое значение (25) для всех ячеек.

Самая трудная часть в написании проверки по хи-квадрату — вычисление p-value. Если вы просмотрите определение метода ChiSquareProb на рис. 4, вы быстро осознаете, что он требует глубоких специализированных знаний. Более того, этот метод вызывает вспомогательный метод Gauss, которые не менее сложен.

Кодирование сложных численных функций на самом деле не такое уж трудное, поскольку большинство этих функций решалось десятилетиями. Два из наиболее часто используемых мной ресурсов — репозитарий Association for Computing Machinery (ACM) Collected Algorithms и книга Абрамовича (Abramowitz) и Стегуна (Stegun) «Handbook of Mathematical Functions» (Dover Publications, 1965); эти авторы настолько известны, что их часто называют «A&S». Обе ссылки доступны бесплатно в нескольких местах сети.

Например, метод ChiSquareProb реализует ACM Algorithm #299, а вспомогательный метод Gauss — ACM Algorithm #209. Альтернативой ACM #209 является уравнение A&S #7.1.26 плюс простое выражение-оболочка.

Многие функции ACM Algorithms и A&S реализованы в существующих C#-библиотеках; однако большинство из этих библиотек довольно велики. По возможности я предпочитаю кодировать с нуля, используя только нужные мне методы и избегая внешних зависимостей.

Несколько замечаний

В этой статье показана лишь малая часть языка R, но ее должно быть достаточно, чтобы вы могли приступить к работе с этим языком. Программисты на C# обычно вступают в контакт с R двумя способами. Во-первых, вы можете использовать R для выполнения собственного статистического анализа. Как демонстрирует моя статья, работать с R довольно легко, если вы владеете знаниями в области статистики. Во-вторых, вам может понадобиться взаимодействие с кем-либо, кто использует R в качестве основного языка, а значит, вы должны понимать среду и терминологию R.

Вы можете применять несколько средств разработки, интегрируемых с языком R. Один из общеизвестных проектов называется RStudio. В целом, я предпочитаю использовать родную консоль R. Кроме того, несколько продуктов Microsoft работают с R. Например, сервис Microsoft Azure Machine Learning способен выполнять скрипты на языке R, и, как я подозреваю, когда-то поддержка R будет добавлена в Visual Studio.

---

Джеймс Маккафри (Dr. James McCaffrey) работает на Microsoft Research в Редмонде (штат Вашингтон). Принимал участие в создании нескольких продуктов Microsoft, в том числе Internet Explorer и Bing. С ним можно связаться по адресу jammc@microsoft.com.

Выражаю благодарность за рецензирование статьи эксперту Microsoft Research Дэну Либлингу (Dan Liebling).