mẫu thống nhất một cách ngẫu nhiên từ một đơn vị n-chiều simplex

Question

Lấy mẫu thống nhất một cách ngẫu nhiên từ một đơn vị n-chiều simplex là cách ưa thích để nói rằng bạn muốn số n ngẫu nhiên mà

• họ là tất cả phi phủ định,
• chúng được tính thành một và
• mọi vectơ có thể có của n số không âm có tổng là một số có khả năng giống nhau.

Trong trường hợp n = 2, bạn muốn lấy mẫu thống nhất từ ​​phân đoạn của dòng x + y = 1 (tức là y = 1-x) ở góc phần tư dương. Trong n = 3 trường hợp bạn đang lấy mẫu từ phần hình tam giác của mặt phẳng x + y + z = 1 có nghĩa là trong chòm sao bát nhân tích cực của R3:

(Hình ảnh từ http://en.wikipedia.org/wiki/Simplex.)

Lưu ý rằng chọn n số ngẫu nhiên đồng nhất và sau đó bình thường hóa chúng để chúng tổng hợp thành một không hoạt động. Bạn kết thúc với một thiên vị đối với số ít cực đoan hơn.

Tương tự, chọn số ngẫu nhiên đồng nhất n-1 và sau đó lấy số thứ n là một trừ đi tổng số tiền cũng giới thiệu sai lệch.

Wikipedia cung cấp hai thuật toán để thực hiện điều này một cách chính xác: http://en.wikipedia.org/wiki/Simplex#Random_sampling (Mặc dù câu hỏi thứ hai hiện chỉ đúng trong thực tế, chứ không phải theo lý thuyết. Tôi hy vọng làm sạch hoặc làm rõ nó khi tôi hiểu điều này tốt hơn. Ban đầu tôi bị mắc kẹt trong một "CẢNH BÁO: giấy như vậy và tuyên bố sau đây là sai" trên trang Wikipedia đó và một người nào khác biến nó thành "chỉ hoạt động trong thực tế" báo trước.)

Cuối cùng, câu hỏi: Bạn xem xét việc thực hiện tốt nhất việc lấy mẫu simplex trong Mathematica (tốt nhất là với xác nhận thực nghiệm là nó đúng)?

câu hỏi liên quan

• Generating a probability distribution
• java random percentages

Answer 1

mã này có thể làm việc:

samples[n_] := Differences[Join[{0}, Sort[RandomReal[Range[0, 1], n - 1]], {1}]] 

Về cơ bản, bạn chỉ cần chọn n - 1 các địa điểm trong khoảng thời gian [0,1] để chia nhỏ nó rồi lấy kích thước của từng miếng bằng cách sử dụng Differences.

Chạy nhanh Timing cho thấy điều này nhanh hơn một chút so với câu trả lời đầu tiên của Janus.

Answer 2

Sau một chút đào bới xung quanh, tôi thấy this page mà đưa ra một thực hiện tốt đẹp của các phân phối Dirichlet. Từ đó nó có vẻ như nó sẽ được khá đơn giản để làm theo phương pháp của Wikipedia 1. Điều này có vẻ như cách tốt nhất để làm điều đó.

Là một thử nghiệm:

In[14]:= RandomReal[DirichletDistribution[{1,1}],WorkingPrecision->25] 
Out[14]= {0.8428995243540368880268079,0.1571004756459631119731921} 
In[15]:= Total[%] 
Out[15]= 1.000000000000000000000000 

Một lô 100 mẫu:

alt text http://www.public.iastate.edu/~zdavkeos/simplex-sample.png

Answer 3

Tôi với zdav: phân phối Dirichlet có vẻ là cách dễ nhất phía trước, và các thuật toán để lấy mẫu phân phối Dirichlet mà zdav đề cập đến cũng được trình bày trên trang Wikipedia trên Dirichlet distribution.

Thực hiện, đó là một chút chi phí để thực hiện phân phối Dirichlet đầy đủ trước tiên, vì tất cả những gì bạn thực sự cần là n ngẫu nhiên Gamma[1,1] mẫu. Hãy so sánh dưới đây
đơn giản thực hiện

SimplexSample[n_, opts:OptionsPattern[RandomReal]] := 
    (#/Total[#])& @ RandomReal[GammaDistribution[1,1],n,opts] 

Dirichlet thực hiện đầy đủ

DirichletDistribution/:Random`DistributionVector[ 
DirichletDistribution[alpha_?(VectorQ[#,Positive]&)],n_Integer,prec_?Positive]:= 
    Block[{gammas}, gammas = 
     Map[RandomReal[GammaDistribution[#,1],n,WorkingPrecision->prec]&,alpha]; 
     Transpose[gammas]/Total[gammas]] 

SimplexSample2[n_, opts:OptionsPattern[RandomReal]] := 
    (#/Total[#])& @ RandomReal[DirichletDistribution[ConstantArray[1,{n}]],opts] 

Timing

Timing[Table[SimplexSample[10,WorkingPrecision-> 20],{10000}];] 
Timing[Table[SimplexSample2[10,WorkingPrecision-> 20],{10000}];] 
Out[159]= {1.30249,Null} 
Out[160]= {3.52216,Null} 

Vì vậy, toàn bộ Dirichlet là một yếu tố của 3 chậm hơn. Nếu bạn cần m> 1 điểm mẫu tại một thời điểm, bạn có thể có thể giành chiến thắng thêm bằng cách thực hiện (#/Total[#]&)/@RandomReal[GammaDistribution[1,1],{m,n}].

Answer 4

Dưới đây là một việc thực hiện tốt chính xác của thuật toán thứ hai từ Wikipedia:

SimplexSample[n_] := Rest@# - Most@# &[Sort@Join[{0,1}, RandomReal[{0,1}, n-1]]] 

Đó là chuyển thể từ đây: http://www.mofeel.net/1164-comp-soft-sys-math-mathematica/14968.aspx (Ban đầu nó có Liên minh thay vì Sắp xếp @ Tham gia - sau này là nhanh hơn một chút.)

(Xem ý kiến ​​đối với một số bằng chứng cho thấy điều này là đúng!)

Answer 5

Tôi đã tạo một thuật toán để tạo ngẫu nhiên đồng đều trên một đơn vị. Bạn có thể tìm thấy các chi tiết trong bài báo trong các liên kết sau đây: http://www.tandfonline.com/doi/abs/10.1080/03610918.2010.551012#.U5q7inJdVNY

nói ngắn gọn, bạn có thể sử dụng sau đây công thức đệ quy để tìm các điểm ngẫu nhiên trên simplex n chiều:

x = 1- R^1/n-1

x _k = (1- Σ _{i = 1} ^k x _i) (1- R _k^1/nk), k = 2, ..., N-1

x _n = 1- Σ _{i = 1}^n-1 x _i

đâu R_i của là số ngẫu nhiên giữa 0 và 1.

Bây giờ tôi đang cố gắng tạo một thuật toán để tạo các mẫu đồng nhất ngẫu nhiên từ ràng buộc simplex.that là giao nhau giữa một đơn giản và một cơ thể lồi.