2010-11-27 1 views
26

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

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

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

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

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

+7

Вы не смогли описать наиболее важную часть вашего вопроса - что вы хотите сделать с моделью? – hadley 2010-11-27 09:25:30

+3

Этот вопрос лучше подходит для [stats.se] (http://stats.stackexchange.com/) – csgillespie 2010-11-27 10:43:17

+1

Эх, я согласен разрешить ему не учитывать, что он будет делать с параметрической моделью. Даже просто работа с синтетическими данными, полученными из адекватной параметрической модели, достаточно для того, чтобы задать такой вопрос. Бутстрап замечательный, но вы должны сохранить или отправить данные. – Iterator 2011-09-09 12:41:56

ответ

37

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

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

  1. Least Squares
  2. Maximum Likelihood

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

Вот конкретный пример того, как оценивать параметры в 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) 

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

Математически это использует максимальную вероятность для оценки среднего и стандартного отклонения гауссова. Вероятность означает (в данном случае) «вероятность данных при заданных значениях параметров». Максимальное правдоподобие означает «значения параметров, которые максимизируют вероятность генерации моих входных данных». Оценка максимального правдоподобия - это алгоритм для нахождения значений параметров, которые максимизируют вероятность генерации входных данных, а для некоторых распределений он может включать в себя алгоритмы numerical optimization. В 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). Из-за этого очень часто используется предварительное знание, чтобы ограничить ваши модели выбора подмножеством всех возможных моделей.

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

2

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

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

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

Я не эксперт, просто еще один скромный Web Developer =)

+4

Я ученый, и ваша идея построения графика ваших данных и сравнения его с известными дистрибутивами действительно хороша - это основа как максимального правдоподобия, так и наименьших квадратов. Разница между ними заключается в том, как они оценивают соответствие между вашими данными и дистрибутивами, но оба они основаны на вашей интуитивно привлекательной идее. :) – 2010-11-27 07:02:03

-4

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

3

Вы по существу хотите сравнить данные своего реального мира с набором теоретических распределений. В базе 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)) 
+0

Не могли бы вы рассказать мне, можно ли добавить в тестовый список распределение «Отрицательный биномный»? Я пытался, но не уверен, как поставить что-то с пространством, то есть ссылка на веб-сайте R говорит, что мне нужно поставить «Negative Binomial», но я не уверен, как добавить это в список. – Legend 2010-11-28 20:51:03

5

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

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

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

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

Наиболее общий вид вашего вопроса является темой поля, известного как оценки непараметрической плотности, где данный:

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

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

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

10

Посмотрите на fitdistrplus (http://cran.r-project.org/web/packages/fitdistrplus/index.html).

Несколько быстрых вещи, чтобы отметить:

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