Обратная функция нормального распределения

Опубликовано Aug 8, 2012 в Математика и логика, Наш опыт | Нет комментариев

, , , ,

Обратная функция нормального распределения

 Мы не скрываем, что на собеседованиях иногда задаем элементарные вопросы по математике. Нередко эти вопросы вызывают недоумение или даже обиду, как это, например, случилось совсем недавно. Разумеется,  мы не собираемся отступать от данного подхода. А чтобы обосновать его ниже приводится задача, которую пришлось решать нашим сотрудникам в рамках реального проекта по статистической обработке и визуализации результатов УЗИ исследований. Разумеется, по понятным причинам, связанным с NDA, мы не можем здесь описать полную техническую постановку задачи в контексте привязки к предметной области. Но в результате декомпозиции была выделена небольшая самостоятельная подзадача чисто математического плана. Причем ее постановка имеет настолько классическую формулировку, что ее наверняка решали все студенты технических вузов как минимум на листике, глядя в таблицы из учебников или методичек . И вот теперь пришло время программировать.

Требуется реализовать классы, которые позволяют вычислять значения обратной функции нормального распределения. Чтобы было совсем понятно – необходимо реализовать полный аналог функции NORMINV(y,0,1) из EXCEL. Ее график приведен на рисунке. Сама кривая  – это развернутый на 90 градусов и сдвинутый график функции Лапласа (см. формулу).

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

Язык решения C# – на нем  все это и работает в реальном проекте.  Интересной особенностью решения является то, как на этой с первого взгляда безнадежно старой задаче,  применяются новые веяния синтаксиса C#, в частности, LINQ.

О чем нужно подумать и что нужно сделать? Раз речь идет об интерполировании, начинать нужно с загрузки таблицы, содержащей значения узлов. Их можно вычислить прямо в  EXCELе, раз уж он упомянут в условии, и выгрузить в файл в доступном для простого чтения формате. Предварительно, нужно определиться с точностью вычислений. Этот вопрос можно озвучить напрямую, но, скорее всего разумным предложением будет использование 2-х знаков после запятой для аргумента x и 4-х знаков после запятой для функции Ф(x), как это принято во всех учебниках по теории вероятности (и до компьютеров была жизнь :) ), где обязательно приведена эта таблица.

Дальше следует принять во внимание, что не любой аргумент можно задавать (допустимый диапазон (0, 1) ) и нужно подумать о том, что возвращать, если в этот диапазон не попали, или как вызывать функцию, чтобы гарантировать попадание в нужный диапазон.

Классы, реализующие функцию нормального распределения

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
 
namespace VolusonChartXml2Dita.Xml
{
  [System.Xml.Serialization.XmlTypeAttribute(Namespace = "")]
  [System.Xml.Serialization.XmlRootAttribute
                ("root", Namespace = "", IsNullable = false)]
  public class CxStdNormalTable
  {
    [System.Xml.Serialization.XmlElement("Norm")]
    //Laplace function table values (interpolation nodes)
    public CxNorm[] Items; 
 
    public double GetZValue(double y)
    {
      // y value should be between (0.5 and 1.0)
      // and Laplace function table should be not empty
      if (y >= 0.5 && y  z.val == y).FirstOrDefault();
        if(exact != null)
        {
         //return if y is exact value from Laplace function table
            return exact.x;
        }
        //calculation of y boundaries
        var low =
          Items.Where(z => (z.val < y)).OrderBy(z =>
             (z.val*-1)).FirstOrDefault();
        var up =
          Items.Where(z => (z.val > y)).OrderBy(z => z.val).FirstOrDefault();
        if (low != null && up != null)
        {
          double delta = up.val - low.val;
          if (delta != 0.0)
          {
            //inverse interpolation
            return low.x + (y - low.val) * (up.x - low.x)/delta;
          }
          return low.x;
        }
      }
      //the value is not a number (if argument not in (0,1))
      return Double.NaN;
    }
  }
  //----------------------------------------------------------------
  [System.Xml.Serialization.XmlTypeAttribute(Namespace = "")]
  public class CxNorm
  {
    [System.Xml.Serialization.XmlAttribute("x")]
    public double x;
 
    [System.Xml.Serialization.XmlAttribute("val")]
    public double val;
  }
  //----------------------------------------------------------------
}

Пример вызова

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Xml;
using System.Xml.Serialization;
using VolusonChartXml2Dita.Xml;
using System.IO;
 
namespace LINQ_NORMAL_TABLE
{
    class Program
    {
        //-------------------------------------------------------------
        ///
<summary> /// Loads and Deserializes Xml File into Object /// </summary>
 
        static public object LoadXmlFile(string fileFullPath, Type type)
        {
            if (fileFullPath =="")
                return null;
 
            using (var fileStream =
                      new FileStream(fileFullPath, FileMode.Open))
            {
                using (var reader = new StreamReader(fileStream))
                {
                    string fileContent = reader.ReadToEnd();
 
                    var myDeserializer = new XmlSerializer(type);
                    using (var myReader = new StringReader(fileContent))
                    {
                        return myDeserializer.Deserialize(myReader);
                    }
                }
            }
        }
 
        static void Main(string[] args)
        {
 
            CxStdNormalTable NormTable;
 
            object obj = LoadXmlFile("..\\..\\Norm_table.xml",
                         typeof(CxStdNormalTable));
            if (obj != null)
            {
                NormTable = (CxStdNormalTable)obj;
 
                double val = 190;//test value 
 
                bool islower = val < 50.0;
                val = (islower) ? 100.0 - val : val;
                val /= 100.0;
                double retVal;
                retVal = NormTable.GetZValue(val) * ((islower) ? -1 : 1);
 
            }
        }
    }
}

XML с таблицей функции Лапласа


Автор публикации:

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

Ваш адрес email не будет опубликован.


*