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

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

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

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,s2)) # function for probability distribution
  
}

g <- Vectorize(f,c("n1","s1","n2","s2")) 
set.seed(123)
res <- outer(df$mpg,df$sd,df$mpg,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 比另一辆车少的概率。但是概率的总和不会加到 1。鉴于可能会发生紧缩，我希望可以将其添加到 1 或更低。

例如，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),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,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,-1) # output

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

> subset(datalong_2,"car2"))
       p1   p2     value
1121 car1 car1 0.5000000
1122 car2 car1 0.3327904
1155 car1 car2 0.6672096
1156 car2 car2 0.5000000