Р: Как получить значение р (для получения х <2) для популяции (моделируется комбинацией из трех нормальных распределений)? - PullRequest
2 голосов
/ 27 июня 2011

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

https://stats.stackexchange.com/questions/10062/which-r-package-to-use-to-calculate-component-parameters-for-a-mixture-model

После подгонки моегоданные, параметры для трех нормальных распределений следующие:

      pi      mu sigma
1 0.5552 -0.4868 2.044
2 0.2739  8.3846 1.399
3 0.1709 12.5317 1.036

Чтобы проверить соответствие между моими данными (x) и модельным распределением (ee), я сделал следующие шаги:

e1 <- rnorm(5552, mean=-0.4868, sd=2.044)
e2 <- rnorm(2739, mean=8.3846, sd=1.399)
e3 <- rnorm(1709, mean=12.5317, sd=1.036)
ee <- c(e1,e2,e3)
qqplot(x, ee)

Я получил qqplot следующим образом: (http://i.stack.imgur.com/3favy.png)

Вроде бы неплохо, поэтому я хочу вычислить p-значение для получения значения, равного или меньшего 2,0 для этой моделиНаселение. Не могли бы вы научить меня, как рассчитать это p-значение, используя R?

График плотности модели населения ee прилагается (http://i.stack.imgur.com/pExhF.png).

Ответы [ 2 ]

1 голос
/ 27 июня 2011

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

Отдельные вклады "слева от 2":

> probs <- with(param, pi*pnorm(2, mu, sigma) )
> probs
[1] 4.930888e-01 6.883473e-07 2.409615e-25

Итого:

> prob <- with(param, sum(pi*pnorm(2, mu, sigma)) )
> prob
[1] 0.4930895

Просто взглянув на вывод, вы догадаетесь, что почти весь вклад будет от первого компонента, поскольку средства двух других находятся далеко направо.В значении первого вклада преобладает оценка pi (пропорции), поскольку большая часть его массы находится «слева от 2».

0 голосов
/ 27 июня 2011

Я думаю, вы просто хотите вычислить вероятность области хвоста с помощью чего-то вроде.

> sum(ee <= 2) / length(ee)
> [1] 0.4936
...