2013-03-03 6 views
6

Jestem nowym użytkownikiem R i próbuję tworzyć wektory z liczbami losowo generowanymi na podstawie określonej dystrybucji (na przykład z komendą rnorm) z wektorami mającymi określoną sumę gęstości prawdopodobieństwa lub sumę rozkładów skumulowanych.Tworzenie wektorów RNG w R, które mają wstępnie zdefiniowaną sumę pdf lub sumy cdf

Na przykład, podczas generowania wektorów x1, x2 ... xn chcę ich słuchać albo

sum(pnorm(x1)) = sum(pnorm(x2)) = … sum(pnorm(xn))

lub

sum(pnorm(xi)) = ”fixed value”

lub zrobić to samo, ale z dnorm. Innymi słowy, czy istnieje możliwość ustawienia takich parametrów podczas korzystania z rnorm lub jakiegokolwiek innego RNG w R?

Porady i sugestie dotyczące strategii zamiast kompletnych rozwiązań również zostałyby docenione.

Z góry dziękujemy za poświęcony czas.

+0

Interesujące pytanie. Moim pierwszym pomysłem byłoby wykorzystanie pewnego rodzaju antytetycznego podejścia zmienności. Zobacz na przykład http://pl.wikipedia.org/wiki/Antitetyczne_wariaty –

+0

To wydaje się bardziej pytanie o statystyki/matematykę niż pytanie o programowanie. możesz chcieć przynieść to do stat.stackexchange.com lub math.sta ... a następnie, jeśli potrzebujesz pomocy przy implementacji, pytając tutaj o to –

+0

Czy zauważyłeś, że istnieje ograniczenie liczby wektorów tego typu? Na przykład, jeśli twoje x1, x2, .. są po prostu skalarami, istnieje tylko co najwyżej dwie wartości z identycznym prawdopodobieństwem w przypadku normalnej dystrybucji (x + mu i x-mu)? –

Odpowiedz

6

1. W przypadku rozkładu Gaussa, próbkowania od (X1,...,Xn) pod warunkiem, że X1+...+Xn=s jest tylko próbkowanie z conditional Gaussian distribution.

wektora (X1, X2, ... Xn X1 + ... + Xn) ma rozkład Gaussa, z zerową średnią i wariancją, matrycy

1 0 0 ... 0 1 
0 1 0 ... 0 1 
0 0 1 ... 0 1 
... 
0 0 0 ... 1 1 
1 1 1 ... 1 n. 

Możemy zatem próbka z nią następująco.

s <- 1 # Desired sum 
n <- 10 
mu1 <- rep(0,n) 
mu2 <- 0 
V11 <- diag(n) 
V12 <- as.matrix(rep(1,n)) 
V21 <- t(V12) 
V22 <- as.matrix(n) 
mu <- mu1 + V12 %*% solve(V22, s - mu2) 
V <- V11 - V12 %*% solve(V22,V21) 
library(mvtnorm) 
# Random vectors (in each row) 
x <- rmvnorm(100, mu, V) 
# Check the sum and the distribution 
apply(x, 1, sum) 
hist(x[,1]) 
qqnorm(x[,1]) 

W przypadku dystrybucji arbitralnej podejście to wymagałoby obliczenia rozkładu warunkowego, co może nie być łatwe.

2. Istnieje jeszcze jeden łatwy, specjalny przypadek: rozkład jednorodny.

Do próbki równomiernie n (pozytywnych) liczb, które sumują się do 1, można wziąć n-1 cyfry jednostajnie w [0,1], i sortować je: określają n przedziałów, których długości kolei sumują się do 1 i są równomiernie rozłożone.

Ponieważ te punkty tworzą proces Poissona, , można również wygenerować je z rozkładem wykładniczym.

x <- rexp(n) 
x <- x/sum(x) # Sums to 1, and each coordinate is uniform in [0,1] 

Ta idea jest wyjaśnione (z dużą ilością zdjęć) w następującym artykule: Portfolio Optimization for VaR, CVaR, Omega and Utility with General Return Distributions, (WT Shaw, 2011), strony 6 do 8.

3. (EDIT) Początkowo błędnie przeczytałem pytanie, które było pytaniem o sum(pnorm(x)), a nie sum(x). To okazuje się łatwiejsze.

Jeśli X ma rozkład Gaussa, to pnorm(X) ma jednorodny rozkład: problemem jest wówczas pobranie próbki z jednolitego rozkładu z określoną sumą.

n <- 10 
s <- 1 # Desired sum 
p <- rexp(n) 
p <- p/sum(p) * s # Uniform, sums to s 
x <- qnorm(p)  # Gaussian, the p-values sum to s 
+0

'apply (x, 1, sum)' powinno zwrócić tę samą wartość dla wszystkich wierszy? – agstudy

+1

Tak: jest to żądana suma, powinna wynosić (około) 1 dla każdego wiersza. –

+0

jest OK! Mam pomyłkę w niektórych sesjach zmiennych ..! +1! Przez pewien czas próbowałem tej niegrzecznej metody (generuję x2, aż stwierdzę warunek, x3 ...), ale jest to bardzo nieefektywne. – agstudy

Powiązane problemy