p*********g
发帖数: 226 | 1 给了 n 个独立的 Bernolli 随机变量X_i，p(X_i = 1) = p_i. 现要采样（X_1, ...,
X_n) 使得它们的和为 m (m <=n). 如何高效地采样？
当然规定了和之后，它们就不再独立了。
简单的方法：独立采样所有 X_i, 如果和不为m，扔掉，重新采。当然效率太低了。 | a****m
发帖数: 693 | 2 p_1+p_2+...=m/n,
multiple normal distribution | p*********g
发帖数: 226 | 3 对不起，能具体点吗？
【在 a****m 的大作中提到】
: p_1+p_2+...=m/n,
: multiple normal distribution
| k*****u
发帖数: 1688 | 4 怎么觉得像sequential的问题？
美丑一次就检查一下，小雨m，继续抽样，直到到了m为止。 sequential test是不是就是这么做的？以前学的东西都忘了
你说不是ｍ就扔掉，这样太浪费了吧 | f***a
发帖数: 329 | 5 Here you are. coded in R. It is easy to read and get the idea directly from
the code. Have fun!
###############
n <- 10 # how many Bernolli variables
m <- 6 # how many "1"s you want
nout <- 100 # how many such samples you want
p <- runif(n) # pre-specified probability of Bernolli variables
pp <- p/sum(p)
out <- matrix(0,nout,n)
for(irun in 1:nout)
{
ind <- numeric(m)
iter <- 1
while(iter<=m)
{
res <- sum(runif(1) > cumsum(pp))+1
if(!(res %in% ind)) {ind[iter] <- res; iter <- iter+1}
}
out[irun,ind] <- 1
}
out # the sample! | p*********g
发帖数: 226 | 6 多谢。不过这里有个典型的陷阱，那就是这些随机变量在给定和之后就不再独立了。而
且每采完一个 X_i 之后，剩下的（条件）概率要随之调整。
举个反例吧。取 p1 = 1/2, p2 = 1/3, p3 = 1/4. 现要和 s = 2.
一方面
（*） p(x2 = 1 | s = 2) = p(x2 = 1, s = 2) / p(s = 2) = 2/3
（列举一下所有的8种可能，然后算这两个概率：
x1 x2 x3 prob
0 0 0 1/4
0 0 1 1/12
0 1 0 1/8
0 1 1 1/24
1 0 0 1/4
1 0 1 1/12
1 1 0 1/8
1 1 1 1/24
故 p(s = 2) = 1/24 + 1/12 + 1/8 = 1/4
p(x2 = 1, s = 2) = 1/24 + 1/8 = 1/6.)
另一方面，依照你的采样方法，列举一下采样到的概率：
先1 后2：6/13 * 4/7
先1 后3：6/13 * 3/7
先2 后1：4/13 * 2/3
先2 后3：4/13 * 1/3
先3 后1：3/13 * 3/5
先3 后2：3/13 * 2/5
故 p(x2 = 1 | s = 2)
= 6/13 * 4/7 + 4/13 * 2/3 + 4/13 * 1/3 + 3/13 * 2/5
= 302 / 455，
而非（*）里的 2/3。
from
【在 f***a 的大作中提到】
: Here you are. coded in R. It is easy to read and get the idea directly from
: the code. Have fun!
: ###############
: n <- 10 # how many Bernolli variables
: m <- 6 # how many "1"s you want
: nout <- 100 # how many such samples you want
: p <- runif(n) # pre-specified probability of Bernolli variables
: pp <- p/sum(p)
: out <- matrix(0,nout,n)
: for(irun in 1:nout)
| p*********g
发帖数: 226 | 7 俺先确认对你的算法理解正确：
S = empty
loop until |S| = m
sample i with prob pi / sum_{j=1}^n p_j
if i is not in S
add i to S
end
end
return S
n <- 10 # how many Bernolli variables
m <- 6 # how many "1"s you want
nout <- 100 # how many such samples you want
p <- runif(n) # pre-specified probability of Bernolli variables
pp <- p/sum(p)
out <- matrix(0,nout,n)
for(irun in 1:nout)
{
ind <- numeric(m)
iter <- 1
while(iter<=m)
{
res <- sum(runif(1) > cumsum(pp))+1
if(!(res %in% ind)) {ind[iter] <- res; iter <- iter+1}
}
out[irun,ind] <- 1
}
out # the sample!
【在 f***a 的大作中提到】
: Here you are. coded in R. It is easy to read and get the idea directly from
: the code. Have fun!
: ###############
: n <- 10 # how many Bernolli variables
: m <- 6 # how many "1"s you want
: nout <- 100 # how many such samples you want
: p <- runif(n) # pre-specified probability of Bernolli variables
: pp <- p/sum(p)
: out <- matrix(0,nout,n)
: for(irun in 1:nout)
| f***a
发帖数: 329 | 8 你理解是对的，不过我的算法是错的。。。。
譬如p_i=1的时候
有空再想，要走了要走了
【在 p*********g 的大作中提到】
: 俺先确认对你的算法理解正确：
: S = empty
: loop until |S| = m
: sample i with prob pi / sum_{j=1}^n p_j
: if i is not in S
: add i to S
: end
: end
: return S
:
| h***i
发帖数: 3844 | 9 去哪里？
【在 f***a 的大作中提到】
: 你理解是对的，不过我的算法是错的。。。。
: 譬如p_i=1的时候
: 有空再想，要走了要走了
| k*******a
发帖数: 772 | 10 没怎么看懂题目。。。可不可以这样，
第一步，X_i中找 和m/n 最近的那个p，抽样
第二步，如果第一个样品是1， 那么剩下的找和 m-1/n-1最近的p的抽样
反之，找和 m/n-1最近的p抽样
。。。。直到最后一个 | | | a****m
发帖数: 693 | 11 不好意思，可能自己的理解有偏差：
m/(m1*m2*m3...(m-m1-m2-m3-...))*(p1)^m1*p2^m2*p3^m3...(m/n-p1-p2-p3-...)
感觉好像是求能得到这样的概率，而不是怎么有效的抽样。
【在 p*********g 的大作中提到】
: 对不起，能具体点吗？
| p*********g
发帖数: 226 | 12 其实就是先定义一个(X1,...,Xn) 的联合分布，
p(X1,...,Xn) = \prod_i p(X_i)
p(X_i = 1) = pi.
然后按
p(X1, ..., Xn | \sum_i X_i = m)
抽样。
【在 k*******a 的大作中提到】
: 没怎么看懂题目。。。可不可以这样，
: 第一步，X_i中找 和m/n 最近的那个p，抽样
: 第二步，如果第一个样品是1， 那么剩下的找和 m-1/n-1最近的p的抽样
: 反之，找和 m/n-1最近的p抽样
: 。。。。直到最后一个
| p*********g
发帖数: 226 | 13 能给点 intuition 为什么这样能 work 吗？多谢。
【在 k*******a 的大作中提到】
: 没怎么看懂题目。。。可不可以这样，
: 第一步，X_i中找 和m/n 最近的那个p，抽样
: 第二步，如果第一个样品是1， 那么剩下的找和 m-1/n-1最近的p的抽样
: 反之，找和 m/n-1最近的p抽样
: 。。。。直到最后一个
| f***a
发帖数: 329 | 14 回来了回来了。
重新想了下，这个其实就是在一堆iid variables之间加了一个constraint。貌似
sample起来不难。
以最简单的n=2,m=1为例：
Without constraint, outputs space is {(0,0),(0,1),(1,0),(1,1)}.
The corresponding probability space is {(1-p1)*(1-p2), ..., p1*p2}.
With constraint, outputs space is O={(0,1),(1,0)}.
The corresponding probability space is P={(1-p1)*(p2), p1*(1-p2)}.
Under the constrain, standardize the probability space into
P.std={P1/(P1+P2),P2/(P1+P2)}.
Then under constrain, output (0,1) has the probability P1/(P1+P2) to be
sampled and (1,0) has the probability P2/(P1+P2) to be sampled.
Generalized this logic into arbitrary numbers for n and m. And the R code is
listed below.
##############
n <- 5 # how many Bernolli variables
m <- 3 # how many "1"s you want
set.seed(123)
p <- round(runif(n),2) # pre-specified probabilities for Bernolli variables
library(wle) #try: binary(15,dim=4)$binary
tt <- t(matrix(unlist(apply(data.frame(x=0:(2^n-1)),1,
function(t)binary(t,dim=n)$binary)),n,))
ind <- which(apply(tt,1,sum)==m)
sample.l <- tt[ind,]#this is the list of all possible outputs!
p.s <- matrix(rep(p,times=nrow(sample.l)),nrow(sample.l),,byrow=TRUE)
p.s1 <- p.s; p.s2 <- 1-p.s;
p.s1[sample.l==0] <- 0
p.s2[sample.l==1] <- 0
p.ss <- p.s1+p.s2
p.out <- apply(p.ss,1,prod)
p.out.std <- p.out/sum(p.out) #this is the list of probabilities
#w.r.t. all possbile outputs!
N <- 100 #how many such samples you want
ind.s <- sample(1:nrow(sample.l),N,replace =TRUE, prob =p.out.std )
res <- sample.l[ind.s,]# Here you are! | p*********g
发帖数: 226 | 15 给了 n 个独立的 Bernolli 随机变量X_i，p(X_i = 1) = p_i. 现要采样（X_1, ...,
X_n) 使得它们的和为 m (m <=n). 如何高效地采样？
当然规定了和之后，它们就不再独立了。
简单的方法：独立采样所有 X_i, 如果和不为m，扔掉，重新采。当然效率太低了。 | a****m
发帖数: 693 | 16 p_1+p_2+...=m/n,
multiple normal distribution | p*********g
发帖数: 226 | 17 对不起，能具体点吗？
【在 a****m 的大作中提到】
: p_1+p_2+...=m/n,
: multiple normal distribution
| k*****u
发帖数: 1688 | 18 怎么觉得像sequential的问题？
美丑一次就检查一下，小雨m，继续抽样，直到到了m为止。 sequential test是不是就是这么做的？以前学的东西都忘了
你说不是ｍ就扔掉，这样太浪费了吧 | f***a
发帖数: 329 | 19 Here you are. coded in R. It is easy to read and get the idea directly from
the code. Have fun!
###############
n <- 10 # how many Bernolli variables
m <- 6 # how many "1"s you want
nout <- 100 # how many such samples you want
p <- runif(n) # pre-specified probability of Bernolli variables
pp <- p/sum(p)
out <- matrix(0,nout,n)
for(irun in 1:nout)
{
ind <- numeric(m)
iter <- 1
while(iter<=m)
{
res <- sum(runif(1) > cumsum(pp))+1
if(!(res %in% ind)) {ind[iter] <- res; iter <- iter+1}
}
out[irun,ind] <- 1
}
out # the sample! | p*********g
发帖数: 226 | 20 多谢。不过这里有个典型的陷阱，那就是这些随机变量在给定和之后就不再独立了。而
且每采完一个 X_i 之后，剩下的（条件）概率要随之调整。
举个反例吧。取 p1 = 1/2, p2 = 1/3, p3 = 1/4. 现要和 s = 2.
一方面
（*） p(x2 = 1 | s = 2) = p(x2 = 1, s = 2) / p(s = 2) = 2/3
（列举一下所有的8种可能，然后算这两个概率：
x1 x2 x3 prob
0 0 0 1/4
0 0 1 1/12
0 1 0 1/8
0 1 1 1/24
1 0 0 1/4
1 0 1 1/12
1 1 0 1/8
1 1 1 1/24
故 p(s = 2) = 1/24 + 1/12 + 1/8 = 1/4
p(x2 = 1, s = 2) = 1/24 + 1/8 = 1/6.)
另一方面，依照你的采样方法，列举一下采样到的概率：
先1 后2：6/13 * 4/7
先1 后3：6/13 * 3/7
先2 后1：4/13 * 2/3
先2 后3：4/13 * 1/3
先3 后1：3/13 * 3/5
先3 后2：3/13 * 2/5
故 p(x2 = 1 | s = 2)
= 6/13 * 4/7 + 4/13 * 2/3 + 4/13 * 1/3 + 3/13 * 2/5
= 302 / 455，
而非（*）里的 2/3。
from
【在 f***a 的大作中提到】
: Here you are. coded in R. It is easy to read and get the idea directly from
: the code. Have fun!
: ###############
: n <- 10 # how many Bernolli variables
: m <- 6 # how many "1"s you want
: nout <- 100 # how many such samples you want
: p <- runif(n) # pre-specified probability of Bernolli variables
: pp <- p/sum(p)
: out <- matrix(0,nout,n)
: for(irun in 1:nout)
| | | p*********g
发帖数: 226 | 21 俺先确认对你的算法理解正确：
S = empty
loop until |S| = m
sample i with prob pi / sum_{j=1}^n p_j
if i is not in S
add i to S
end
end
return S
n <- 10 # how many Bernolli variables
m <- 6 # how many "1"s you want
nout <- 100 # how many such samples you want
p <- runif(n) # pre-specified probability of Bernolli variables
pp <- p/sum(p)
out <- matrix(0,nout,n)
for(irun in 1:nout)
{
ind <- numeric(m)
iter <- 1
while(iter<=m)
{
res <- sum(runif(1) > cumsum(pp))+1
if(!(res %in% ind)) {ind[iter] <- res; iter <- iter+1}
}
out[irun,ind] <- 1
}
out # the sample!
【在 f***a 的大作中提到】
: Here you are. coded in R. It is easy to read and get the idea directly from
: the code. Have fun!
: ###############
: n <- 10 # how many Bernolli variables
: m <- 6 # how many "1"s you want
: nout <- 100 # how many such samples you want
: p <- runif(n) # pre-specified probability of Bernolli variables
: pp <- p/sum(p)
: out <- matrix(0,nout,n)
: for(irun in 1:nout)
| f***a
发帖数: 329 | 22 你理解是对的，不过我的算法是错的。。。。
譬如p_i=1的时候
有空再想，要走了要走了
【在 p*********g 的大作中提到】
: 俺先确认对你的算法理解正确：
: S = empty
: loop until |S| = m
: sample i with prob pi / sum_{j=1}^n p_j
: if i is not in S
: add i to S
: end
: end
: return S
:
| h***i
发帖数: 3844 | 23 去哪里？
【在 f***a 的大作中提到】
: 你理解是对的，不过我的算法是错的。。。。
: 譬如p_i=1的时候
: 有空再想，要走了要走了
| k*******a
发帖数: 772 | 24 没怎么看懂题目。。。可不可以这样，
第一步，X_i中找 和m/n 最近的那个p，抽样
第二步，如果第一个样品是1， 那么剩下的找和 m-1/n-1最近的p的抽样
反之，找和 m/n-1最近的p抽样
。。。。直到最后一个 | a****m
发帖数: 693 | 25 不好意思，可能自己的理解有偏差：
m/(m1*m2*m3...(m-m1-m2-m3-...))*(p1)^m1*p2^m2*p3^m3...(m/n-p1-p2-p3-...)
感觉好像是求能得到这样的概率，而不是怎么有效的抽样。
【在 p*********g 的大作中提到】
: 对不起，能具体点吗？
| p*********g
发帖数: 226 | 26 其实就是先定义一个(X1,...,Xn) 的联合分布，
p(X1,...,Xn) = \prod_i p(X_i)
p(X_i = 1) = pi.
然后按
p(X1, ..., Xn | \sum_i X_i = m)
抽样。
【在 k*******a 的大作中提到】
: 没怎么看懂题目。。。可不可以这样，
: 第一步，X_i中找 和m/n 最近的那个p，抽样
: 第二步，如果第一个样品是1， 那么剩下的找和 m-1/n-1最近的p的抽样
: 反之，找和 m/n-1最近的p抽样
: 。。。。直到最后一个
| p*********g
发帖数: 226 | 27 能给点 intuition 为什么这样能 work 吗？多谢。
【在 k*******a 的大作中提到】
: 没怎么看懂题目。。。可不可以这样，
: 第一步，X_i中找 和m/n 最近的那个p，抽样
: 第二步，如果第一个样品是1， 那么剩下的找和 m-1/n-1最近的p的抽样
: 反之，找和 m/n-1最近的p抽样
: 。。。。直到最后一个
| f***a
发帖数: 329 | 28 回来了回来了。
重新想了下，这个其实就是在一堆iid variables之间加了一个constraint。貌似
sample起来不难。
以最简单的n=2,m=1为例：
Without constraint, outputs space is {(0,0),(0,1),(1,0),(1,1)}.
The corresponding probability space is {(1-p1)*(1-p2), ..., p1*p2}.
With constraint, outputs space is O={(0,1),(1,0)}.
The corresponding probability space is P={(1-p1)*(p2), p1*(1-p2)}.
Under the constrain, standardize the probability space into
P.std={P1/(P1+P2),P2/(P1+P2)}.
Then under constrain, output (0,1) has the probability P1/(P1+P2) to be
sampled and (1,0) has the probability P2/(P1+P2) to be sampled.
Generalized this logic into arbitrary numbers for n and m. And the R code is
listed below.
##############
n <- 5 # how many Bernolli variables
m <- 3 # how many "1"s you want
set.seed(123)
p <- round(runif(n),2) # pre-specified probabilities for Bernolli variables
library(wle) #try: binary(15,dim=4)$binary
tt <- t(matrix(unlist(apply(data.frame(x=0:(2^n-1)),1,
function(t)binary(t,dim=n)$binary)),n,))
ind <- which(apply(tt,1,sum)==m)
sample.l <- tt[ind,]#this is the list of all possible outputs!
p.s <- matrix(rep(p,times=nrow(sample.l)),nrow(sample.l),,byrow=TRUE)
p.s1 <- p.s; p.s2 <- 1-p.s;
p.s1[sample.l==0] <- 0
p.s2[sample.l==1] <- 0
p.ss <- p.s1+p.s2
p.out <- apply(p.ss,1,prod)
p.out.std <- p.out/sum(p.out) #this is the list of probabilities
#w.r.t. all possbile outputs!
N <- 100 #how many such samples you want
ind.s <- sample(1:nrow(sample.l),N,replace =TRUE, prob =p.out.std )
res <- sample.l[ind.s,]# Here you are! | I*****a
发帖数: 5425 | 29 sample(1:n, size = m, prob = c(p1, p2, ..., pn) )
.,
【在 p*********g 的大作中提到】
: 给了 n 个独立的 Bernolli 随机变量X_i，p(X_i = 1) = p_i. 现要采样（X_1, ...,
: X_n) 使得它们的和为 m (m <=n). 如何高效地采样？
: 当然规定了和之后，它们就不再独立了。
: 简单的方法：独立采样所有 X_i, 如果和不为m，扔掉，重新采。当然效率太低了。
| t*******2
发帖数: 384 | 30 Re
【在 I*****a 的大作中提到】
: sample(1:n, size = m, prob = c(p1, p2, ..., pn) )
:
: .,
|
|
