Подгонка данных к распределению? - PullRequest
27 голосов
/ 27 ноября 2010

Я не статистика (скорее исследовательский веб-разработчик), но я много слышал о scipy и R в эти дни. Поэтому из любопытства я хотел задать этот вопрос (хотя это может показаться глупым для экспертов здесь), потому что я не уверен в достижениях в этой области и хочу знать, как люди без достаточного фона статистики подходят к этим проблемам.

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

Реальный сценарий: Например, допустим, я инициировал небольшой опрос и записал информацию о том, сколько человек разговаривает с каждым человеком, скажем, по 300 человек, и у меня есть следующая информация:

1 10
2 5
3 20
...
...

где X Y говорит мне, что человек X разговаривал с Y людьми во время опроса. Теперь, используя информацию от 300 человек, я хочу вписать это в модель. Вопрос сводится к тому, существуют ли какие-либо автоматизированные способы определения правильных параметров распределения и распределения для этих данных или, если нет, есть ли хорошая пошаговая процедура для достижения того же самого?

Ответы [ 6 ]

38 голосов
/ 27 ноября 2010

Это сложный вопрос, и нет идеальных ответов. Я постараюсь дать вам обзор основных понятий и укажу вам направление на полезное чтение по этой теме.

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

  1. Наименьшие квадраты
  2. Максимальное правдоподобие

По моему опыту, максимальное правдоподобие было предпочтительным в последние годы, хотя это может быть не во всех областях.

Вот конкретный пример того, как оценить параметры в R. Рассмотрим набор случайных точек, сгенерированных из распределения Гаусса со средним значением 0 и стандартным отклонением 1:

x = rnorm( n = 100, mean = 0, sd = 1 )

Предположим, вы знаете, что данные были сгенерированы с использованием гауссовского процесса, но вы забыли (или никогда не знали!) Параметры для гауссовского. Вы хотели бы использовать данные, чтобы дать вам разумные оценки среднего и стандартного отклонения. В R есть стандартная библиотека, которая делает это очень просто:

library(MASS)
params = fitdistr( x, "normal" )
print( params )

Это дало мне следующий вывод:

      mean           sd     
  -0.17922360    1.01636446 
 ( 0.10163645) ( 0.07186782)

Они довольно близки к правильному ответу, а числа в скобках - это доверительные интервалы вокруг параметров. Помните, что каждый раз, когда вы генерируете новый набор баллов, вы получаете новый ответ для оценок.

Математически это максимальная вероятность для оценки как среднего, так и стандартного отклонения гауссианы. Вероятность означает (в данном случае) «вероятность данных заданных значений параметров». Максимальное правдоподобие означает «значения параметров, которые максимизируют вероятность генерации моих входных данных». Оценка максимального правдоподобия - это алгоритм для нахождения значений параметров, которые максимизируют вероятность генерации входных данных, и для некоторых распределений он может включать числовую оптимизацию алгоритмы. В R большая часть работы выполняется fitdistr , который в некоторых случаях будет вызывать optim .

Вы можете извлечь логарифмическую правдоподобие из ваших параметров следующим образом:

print( params$loglik )
[1] -139.5772

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

С такими вычислительными инструментами легко оценить параметры для любого распределения. Рассмотрим этот пример:

x = x[ x >= 0 ]

distributions = c("normal","exponential")

for ( dist in distributions ) {
    print( paste( "fitting parameters for ", dist ) )
    params = fitdistr( x, dist )
    print( params )
    print( summary( params ) )
    print( params$loglik )
}

Экспоненциальное распределение не генерирует отрицательные числа, поэтому я удалил их в первой строке. Вывод (который является стохастическим) выглядел так:

[1] "fitting parameters for  normal"
      mean          sd    
  0.72021836   0.54079027 
 (0.07647929) (0.05407903)
         Length Class  Mode   
estimate 2      -none- numeric
sd       2      -none- numeric
n        1      -none- numeric
loglik   1      -none- numeric
[1] -40.21074
[1] "fitting parameters for  exponential"
     rate  
  1.388468 
 (0.196359)
         Length Class  Mode   
estimate 1      -none- numeric
sd       1      -none- numeric
n        1      -none- numeric
loglik   1      -none- numeric
[1] -33.58996

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

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

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

10 голосов
/ 27 ноября 2010

Взгляните на fitdistrplus (http://cran.r -project.org / web / packages / fitdistrplus / index.html ).

Несколько быстрых замечаний:

  • Попробуйте использовать функцию descdist, которая предоставляет график асимметрии и эксцесса данных, а также показывает некоторые распространенные распределения.
  • fitdist позволяет вам соответствовать любым распределениям, которые вы можете определить в терминах плотности и cdf.
  • Затем вы можете использовать gofstat, который вычисляет статистику KS и AD, которые измеряют расстояниесоответствуют данным.
5 голосов
/ 27 ноября 2010

Это, вероятно, немного более общее, чем нужно, но может дать вам кое-что для продолжения.

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

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

M. Г. Кендалл и А. Стюарт, Продвинутая теория статистики, вып. 1, Чарльз Гриффин, 1963, был наиболее полным справочником, который я нашел для этого, с огромной целой страницей, посвященной этой теме; в большинстве других текстов было самое большее предложение или перечислялось расширение в терминах моментов вместо кумулянтов, что немного бесполезно. Удачи в поиске копии, тем не менее, я должен был отправить своего университетского библиотекаря в поездку в архив для этого ... но это было годы назад, так что, возможно, Интернет будет более полезным сегодня.

Самая общая форма вашего вопроса - это тема поля, известного как непараметрическая оценка плотности , где дано:

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

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

Лично мне немного не повезло в использовании непараметрической оценки плотности для чего-либо полезного, но если у вас есть постоянный запас здравомыслия, вы должны изучить это.

3 голосов
/ 27 ноября 2010

Вы, по сути, хотите сравнить данные реального мира с набором теоретических распределений.В базе R есть функция qqnorm(), которая сделает это для нормального распределения, но я предпочитаю функцию probplot в e1071, которая позволяет вам тестировать другие распределения.Вот фрагмент кода, который будет отображать ваши реальные данные по каждому из теоретических распределений, которые мы вставляем в список.Мы используем plyr для просмотра списка, но есть и несколько других способов просмотра списка.

library("plyr") 
library("e1071")

realData <- rnorm(1000) #Real data is normally distributed

distToTest <- list(qnorm = "qnorm", lognormal = "qlnorm", qexp =  "qexp")

#function to test real data against list of distributions above. Output is a jpeg for each distribution.
testDist <- function(x, data){
    jpeg(paste(x, ".jpeg", sep = ""))
    probplot(data, qdist = x)
    dev.off()
    }

l_ply(distToTest, function(x) testDist(x, realData))
2 голосов
/ 27 ноября 2010

Я не ученый, но если бы вы делали это с карандашом и бумагой, очевидным способом было бы сделать график, а затем сравнить его с одним из известных стандартных распределений.

* 1002Если идти дальше с этой мыслью, то «сравнение» означает, что кривые стандартного распределения и ваши совпадают.

Тригонометрия, касательные ... были бы моей последней мыслью.

IЯ не эксперт, просто еще один скромный веб-разработчик =)

0 голосов
/ 27 ноября 2010

Как бы то ни было, кажется, что вы, возможно, захотите взглянуть на распределение Пуассона.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...