전력 법칙 분포를 생성하는 임의 번호 생성기?
문제
C ++ Command Linux 앱에 대한 테스트를 작성하고 있습니다. 전력 법/긴 꼬리 분포가있는 수많은 정수를 생성하고 싶습니다. 즉, 나는 일부 숫자를 매우 자주 얻지 만 대부분은 비교적 드물게됩니다.
이상적으로는 rand () 또는 stdlib 랜덤 함수 중 하나와 함께 사용할 수있는 마법 방정식이 있습니다. 그렇지 않다면 사용하기 쉬운 C/C ++ 덩어리가 좋을 것입니다.
감사!
해결책
이것 Wolfram Mathworld의 페이지 균일 한 분포 (대부분의 임의의 숫자 생성기가 제공하는 것)에서 전력법 분포를 얻는 방법에 대해 논의합니다.
짧은 답변 (위의 링크에서 파생) :
x = [(x1^(n+1) - x0^(n+1))*y + x0^(n+1)]^(1/(n+1))
어디 와이 균일 한 변수입니다. N 분포 전력입니다. x0 그리고 x1 분포 범위를 정의하고 엑스 전력 법률 분산 변수입니다.
다른 팁
원하는 분포 (PDF (Propability Distribution Function)라고 함)를 알고 있고 올바르게 정규화되면 누적 분포 함수 (CDF)를 얻기 위해 통합 한 다음 CDF (가능하면)를 뒤집어 변환을 얻을 수 있습니다. 유니폼이 필요합니다 [0,1]
원하는 분포.
따라서 원하는 분포를 정의하는 것으로 시작합니다.
P = F(x)
([0,1]의 X의 경우) 그런 다음 통합되어
C(y) = \int_0^y F(x) dx
이것이 반전 될 수 있다면 당신은 얻을 수 있습니다
y = F^{-1}(C)
그러니 전화 rand()
결과를 AS에 연결합니다 C
마지막 줄에서 y를 사용하십시오.
이 결과를 샘플링의 기본 정리라고합니다. 이것은 정규화 요구 사항과 기능을 분석적으로 반전시켜야하기 때문에 번거 로움입니다.
또는 거부 기술을 사용할 수 있습니다. 원하는 범위에 균일하게 숫자를 던지고 다른 숫자를 던지고 첫 번째 던지기가 빚진 위치에서 PDF와 비교하십시오. 두 번째 던지기가 PDF를 초과하는지 거부하십시오. 긴 꼬리가있는 사람들과 같이 확률 영역이 낮은 PDF에 비효율적 인 경향이 있습니다 ...
중간 접근 방식은 Brute Force에 의해 CDF를 반전시키는 것과 관련이 있습니다. CDF를 조회 테이블로 저장하고 결과를 얻기 위해 리버스 조회를 수행합니다.
여기의 진짜 악취는 간단합니다 x^-n
분포는 범위에서 정규화 할 수 없습니다 [0,1]
, 따라서 샘플링 정리를 사용할 수 없습니다. try (x+1)^-n 대신 ...
전력법 분포를 생성하는 데 필요한 수학에 대해서는 언급 할 수 없습니다 (다른 게시물에는 제안이 있습니다). <random>
. 이것들은보다 더 많은 기능을 제공합니다 std::rand
그리고 std::srand
. 새로운 시스템은 발전기, 엔진 및 분배를위한 모듈 식 API를 지정하고 많은 사전 설정을 공급합니다.
포함 된 분포 사전 설정은 다음과 같습니다.
uniform_int
bernoulli_distribution
geometric_distribution
poisson_distribution
binomial_distribution
uniform_real
exponential_distribution
normal_distribution
gamma_distribution
전력법 분배를 정의 할 때는 기존 발전기 및 엔진과 함께 연결할 수 있어야합니다. 그 책 C ++ 표준 라이브러리 확장 Pete Becker는 훌륭한 장을 가지고 있습니다 <random>
.
여기에 기사가 있습니다 다른 배포판을 만드는 방법에 대해 (Cauchy, Chi-Squared, Student T 및 Snedecor f의 예제 포함)
(정당하게) 수락 된 답변에 대한 보완으로 실제 시뮬레이션을 수행하고 싶었습니다. r에서는 코드가 (의사) -pseudo 코드가 될 정도로 간단합니다.
하나의 작은 차이 Wolfram Mathworld 공식 허용 된 답변 및 기타, 아마도 더 일반적인 방정식에서 전력 법률 지수 n
(일반적으로 알파로 표시됨)는 명시 적 부정적인 부호를 가지고 있지 않습니다. 따라서 선택된 알파 값은 음수이어야하며 일반적으로 2와 3 사이입니다.
x0
그리고 x1
분포의 하한 및 상한을 나타냅니다.
그래서 여기에 있습니다 :
x1 = 5 # Maximum value
x0 = 0.1 # It can't be zero; otherwise X^0^(neg) is 1/0.
alpha = -2.5 # It has to be negative.
y = runif(1e5) # Number of samples
x = ((x1^(alpha+1) - x0^(alpha+1))*y + x0^(alpha+1))^(1/(alpha+1))
hist(x, prob = T, breaks=40, ylim=c(0,10), xlim=c(0,1.2), border=F,
col="yellowgreen", main="Power law density")
lines(density(x), col="chocolate", lwd=1)
lines(density(x, adjust=2), lty="dotted", col="darkblue", lwd=2)
또는 로그 척도로 플롯 팅 :
h = hist(x, prob=T, breaks=40, plot=F)
plot(h$count, log="xy", type='l', lwd=1, lend=2,
xlab="", ylab="", main="Density in logarithmic scale")
데이터 요약은 다음과 같습니다.
> summary(x)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.1000 0.1208 0.1584 0.2590 0.2511 4.9388