Блог

Я хочу понять, почему я не получаю распределение вероятностей, когда использую моделирование из случайного нормального распределения:

 library(tidyverse)
df <- mtcars # data

df$sd <- sd(df$mpg) # standard deviation of the sample

set.seed(123)
f <- function(n1, s1, n2, s2){
  mean(rnorm(10000, n1, s1) < rnorm(10000, n2, s2)) # function for probability distribution
  
}

g <- Vectorize(f, c("n1", "s1", "n2", "s2")) 
set.seed(123)
res <- outer(df$mpg, df$sd, df$mpg, df$sd, FUN = g)
dimnames(res) <- list(row.names(df), row.names(df))
res <- data.frame(res)
res <- tibble::rownames_to_column(res, 'p1')

datalong_2 <- tidyr::gather(res, 'p2', 'value', 2:33) # output
 

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

Например, вероятность того, что Mazda Rx4 mpg будет ниже, чем Mazda Rx4 wag 0,5094, в то время как вероятность того, что Mazda Rx4 wag mpg будет ниже, чем Mazda Rx4 0,5029, сумма этой вероятности равна 1,0123. Как я могу изменить этот код, чтобы получить фактическое распределение вероятности того, что у одного автомобиля mpg ниже, чем у другого автомобиля?

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

Мы предполагаем, что X~N(u1,s1) и Y~N(u2,s2) где s1 и s2 есть отклонения.

Также мы знаем , что P(X<Y) = P(X-Y<0) , где X-Y ~ N(u1-u2,s1 s2) . Исходя из этого, мы можем точно рассчитать вероятности:

 df <- mtcars # data
df$sd <- sd(df$mpg) # standard deviation of the sample

f <- function(n1, n2){
  pnorm(0, mean = n1 - n2, sd = sqrt(2*df$sd^2))
}

res <- outer(X = df$mpg, Y = df$mpg, FUN = f)
dimnames(res) <- list(row.names(df), row.names(df))
res <- data.frame(res)
res <- tibble::rownames_to_column(res, 'p1')

datalong_2 <- tidyr::gather(res, 'p2', 'value', 2:33) # output

> datalong_2
                     p1                p2      value
1             Mazda RX4         Mazda.RX4 0.50000000
2         Mazda RX4 Wag         Mazda.RX4 0.50000000
3            Datsun 710         Mazda.RX4 0.41637203
4        Hornet 4 Drive         Mazda.RX4 0.48128464
5     Hornet Sportabout         Mazda.RX4 0.60636049
..                   ..                ..         ..
 

Кроме того, я думаю , что ваша главная проблема заключалась в функции outer() , для которой требуется 2 входа X и Y . Это сработало для меня, как только я изменил его.

Правки 2 и 3:

 df1 <- mtcars; df1$rownames = rownames(df1)
df2 <- mtcars; df2$rownames = rownames(df2)
df2$mpg = df2$mpg   rnorm(nrow(df2),0,3)
data = rbind(df1, df2)


df = ddply(data,~rownames,summarise,mean=mean(mpg),sd=sd(mpg))
df = rbind(df, c("car1",-1.02, 2.66))
df = rbind(df, c("car2",0.13, 0.06))
df$mean <- as.numeric(df$mean)
df$sd <- as.numeric(df$sd)

f <- function(x, y){
  n1 = df$mean[x]; n2 = df$mean[y]; sd1 = df$sd[x]; sd2 = df$sd[y]
  pnorm(0, mean = n1 - n2, sd = sqrt(sd1^2   sd2^2))
}

res <- outer(X = 1:nrow(df), Y = 1:nrow(df), f)
dimnames(res) <- list(df$rownames, df$rownames)
res <- data.frame(res)
res <- tibble::rownames_to_column(res, 'p1')

datalong_2 <- tidyr::gather(res, 'p2', 'value', -1) # output

subset(datalong_2, p1 %in% c("car1","car2") amp; p2 %in% c("car1","car2"))

> subset(datalong_2, p1 %in% c("car1","car2") amp; p2 %in% c("car1","car2"))
       p1   p2     value
1121 car1 car1 0.5000000
1122 car2 car1 0.3327904
1155 car1 car2 0.6672096
1156 car2 car2 0.5000000