我有一个生物化学问题,可以简化为两辊骰子实验(我认为...)。
假定骰子的面数不均匀,有10个面,即各个面的概率不为1/10。我们想知道这些概率。
但是,我们拥有的给定数据集是将(相同)骰子滚动两次的总和的直方图。因此,观察到的垃圾箱范围是2-20(2 = 1 + 1; 3 = 1 + 2、2 + 1、4 = 2 + 2、1 + 3、3 + 1等)。
相加的面孔的概率是各个概率的乘积(s:相加的面孔的观测概率; p:各个面的概率),可以写成如下:
s2 ~ p1^2
s3 ~ 2*p1*p2
s4 ~ 2*p1*p3 + p2^2
s5 ~ 2*p1*p4 + 2*p2*p3
s6 ~ 2*p1*p5 + 2*p2*p4 + p3^2
s7 ~ 2*p1*p6 + 2*p2*p5 + 2*p3*p4
s8 ~ 2*p1*p7 + 2*p2*p6 + 2*p3*p5 + p4^2
s9 ~ 2*p1*p8 + 2*p2*p7 + 2*p3*p6 + 2*p4*p5
s10 ~ 2*p1*p9 + 2*p2*p8 + 2*p3*p7 + 2*p4*p6 + p5^2
s11 ~ 2*p1*p10 + 2*p2*p9 + 2*p3*p8 + 2*p4*p7 + 2*p5*p6
s12 ~ 2*p2*p10 + 2*p3*p9 + 2*p4*p8 + 2*p5*p7 + p6^2
s13 ~ 2*p3*p10 + 2*p4*p9 + 2*p5*p8 + 2*p6*p7
s14 ~ 2*p4*p10 + 2*p5*p9 + 2*p6*p8 + p7^2
s15 ~ 2*p5*p10 + 2*p6*p9 + 2*p7*p8
s16 ~ 2*p6*p10 + 2*p7*p9 + p8^2
s17 ~ 2*p7*p10 + 2*p8*p9
s18 ~ 2*p8*p10 + p9^2
s19 ~ 2*p9*p10
s20 ~ p10^2
在这种情况下,有20-1 = 19个已知变量和10个未知变量,因此系统是超定的。使用代数手工求解也很容易。据我所记得:二次项将为每张面孔带来2种可能的解决方案。概率总是正的,因此实际上应该有一个解决方案。对吧?
有没有办法解决R中的这个系统?我对R中的线性逆问题很熟悉,但是我不知道如何在R中处理这个(二次方)系统。
以下是一些模拟问题的代码:
options(stringsAsFactors = FALSE) library(gtools) library(dplyr) dice <- data.frame(face = 1:10) ### functions split_dice_faces <- function(summed_face){ face_face <- strsplit(x = as.character(summed_face),split = "[/_\\|]")[[1]] names(face_face) <- c("face1","face2") as.numeric(face_face) } sum_dice_faces <- function(face_face){ sapply(face_face, function(face_face_i){ face1 <- split_dice_faces(face_face_i)[1] face2 <- split_dice_faces(face_face_i)[2] sum(c(face1[1], face2[1])) }) } simulate_2_rolls <- function(dice_pool){ dice_perm <- data.frame(permutations(n = dim(dice_pool)[1], r = 2, v = as.character(dice_pool$face), repeats.allowed = T )) dice_perm$face_face <- paste(dice_perm[[1]],"|",dice_perm[[2]], sep = "") dice_perm$prob <- dice_pool$prob[match(dice_perm[[1]], dice_pool$face)]*dice_pool$prob[match(dice_perm[[2]], dice_pool$face)] dice_perm$summed_face <- sum_dice_faces(dice_perm$face_face) dice_perm <- dice_perm %>% arrange(summed_face) %>% select(one_of(c("face_face", "summed_face","prob"))) dice_perm } summarise_2_rolls_experiment <- function(simulate_2_rolls_df){ simulate_2_rolls_df %>% group_by(summed_face) %>% summarise(prob = sum(prob)) } from_face_probs_to_summed_observations <- function(face_probs){ face_probs %>% data.frame(face = dice$face, prob = .) %>% simulate_2_rolls() %>% summarise_2_rolls_experiment() %>% pull(prob) } generate_formulas <- function() { output <- dice_sum_probs %>% group_by(summed_face) %>% group_split() %>% sapply(function(i){ left_hand <- paste("s",i$summed_face[1],sep="") right_hand <- sapply(strsplit(i$face_face, "\\|") , function(row){ row_i <- as.numeric(row) row_i <- row_i[order(row_i)] row_i <- paste("p",row_i,sep = "") if(row_i[1] == row_i[2]){ paste(row_i[1],"^2",sep="") } else { paste(row_i,collapse="*") } }) right_hand <- paste(sapply(unique(right_hand), function(right_hand_i){ fact <- sum(right_hand == right_hand_i) if(fact > 1){fact <- paste(fact,"*",sep = "")} else {fact <- ""} paste(fact,right_hand_i,sep = "") }), collapse = " + ") paste(left_hand, "~", right_hand) }) return(output) }
模拟数据集:
### random individual probabilites dice_probs <- data.frame(face = dice$face, prob = runif(n = dim(dice)[1]) %>% (function(x){x / sum(x)})) dice_probs ### simulate infinite amount of trials, observations expressed as probabilities dice_sum_probs <- simulate_2_rolls(dice_probs) dice_sum_probs ### sum experiment outcomes with the same sum dice_sum_probs_summary <- dice_sum_probs %>% summarise_2_rolls_experiment() ### plot, this is the given dataset with(data = dice_sum_probs_summary, barplot(prob, names.arg = summed_face)) ### how to calculate / approach p1, p2, ..., p10?
谢谢!
我有一个生物化学问题,可以简化为两辊骰子实验(我认为...)。假设有10张面的骰子参差不齐,即单个面的概率不为1/10。我们...
[如果我们创建概率outer(p, p)
的乘法表,然后使用outer(1:10, 1:10, "+")
将超过tapply
的恒定值的那些求和,则会出现以下非线性回归问题: