如何从汽车对汽车的模拟中获得正态概率分布？

Question

我想了解为什么当我使用随机正态分布的模拟时我没有得到概率分布：

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

我做了这个模拟，但出于某种原因，我没有得到实际的概率分布，我的目标是评估一辆汽车的 mpg 低于另一辆汽车的概率。但是概率之和并不等于一。我希望这可以增加到一个或更低，因为可能会发生紧张。

例如，Mazda Rx4 的 mpg 低于 Mazda Rx4 wag 的概率为 0.5094，而 Mazda Rx4 wag 的 mpg 低于 Mazda Rx4 的概率为 0.5029，这个概率的总和是 1.0123。我如何更改此代码以获得一辆汽车的实际概率分布低于另一辆汽车的 mpg？

Answer 1

除非你绝对需要运行模拟，否则你可以使用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") & p2 %in% c("car1","car2"))

> subset(datalong_2, p1 %in% c("car1","car2") & 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

如何从汽车对汽车的模拟中获得正态概率分布？

How to get normal probability distribution from a simulation on car agains car?

simulation

r

normal-distribution

probability