UFO ET IT

각 비트에 대해 주어진 확률이 0 또는 1 인 의사 랜덤 비트를 생성하는 빠른 방법

ufoet 2020. 12. 15. 19:35
반응형

각 비트에 대해 주어진 확률이 0 또는 1 인 의사 랜덤 비트를 생성하는 빠른 방법


일반적으로 난수 생성기는 각 위치에서 0 또는 1을 관찰 할 확률이 동일한 (즉, 50 %) 비트 스트림을 반환합니다. 이것을 편향되지 않은 PRNG라고합시다.

다음 속성을 사용하여 의사 랜덤 비트 문자열을 생성해야합니다. 각 위치에서 1을 볼 확률은 p입니다 (즉, 0을 볼 확률은 1-p입니다). 매개 변수 p는 0과 1 사이의 실수입니다. 내 문제에서 그것은 0.5 %의 해상도를 가지고 있습니다. 즉, 0 %, 0.5 %, 1 %, 1.5 %, ..., 99.5 %, 100 % 값을 취할 수 있습니다.

p는 정확한 분수가 아니라 확률입니다. n 비트 스트림에서 1로 설정된 실제 비트 수는 이항 분포 B (n, p)를 따라야합니다.

편향되지 않은 PRNG를 사용하여 각 비트 (의사 코드)의 값을 생성 할 수있는 순진한 방법이 있습니다.

generate_biased_stream(n, p):
  result = []
  for i in 1 to n:
    if random_uniform(0, 1) < p:
      result.append(1)
    else:
      result.append(0)
  return result

이러한 구현은 각 비트 당 한 번씩 난수 생성기 함수를 호출하기 때문에 편향되지 않은 스트림을 생성하는 것보다 훨씬 느립니다. 편향되지 않은 스트림 생성기는 단어 크기 당 한 번 호출합니다 (예 : 단일 호출로 32 또는 64 개의 임의 비트를 생성 할 수 있음).

나는 더 빠른 구현을 원하지만 무작위성을 약간 희생하더라도. 생각 나는 아이디어는 조회 테이블을 미리 계산하는 것입니다. p의 가능한 200 개 값 각각에 대해 느린 알고리즘을 사용하여 C 8 비트 값을 계산하고 테이블에 저장합니다. 그런 다음 고속 알고리즘은 이들 중 하나를 무작위로 선택하여 8 개의 왜곡 된 비트를 생성합니다.

얼마나 많은 메모리가 필요한지 알아보기위한 엔벨로프 계산의 뒷면 : C는 최소한 256 (가능한 8 비트 값의 수)이어야하며 샘플링 효과를 피하기 위해 더 많아야합니다. 1024라고 가정 해 봅시다. 아마도 숫자는 p에 따라 달라질 수 있지만 단순하게 유지하고 평균은 1024라고 가정하겠습니다. p의 값이 200 개이므로 => 총 메모리 사용량은 200KB입니다. 이것은 나쁘지 않으며 L2 캐시 (256KB)에 맞을 수 있습니다. 편향을 도입하는 샘플링 효과가 있는지 확인하기 위해 여전히 평가해야합니다.이 경우 C를 증가시켜야합니다.

이 솔루션의 결점은 많은 작업이있는 경우에도 한 번에 8 비트 만 생성 할 수있는 반면, 편향되지 않은 PRNG는 몇 개의 산술 명령만으로 한 번에 64 비트를 생성 할 수 있다는 것입니다.

조회 테이블 대신 비트 연산을 기반으로 더 빠른 방법이 있는지 알고 싶습니다. 예를 들어 난수 생성 코드를 직접 수정하여 각 비트에 대한 편향을 도입합니다. 이것은 편향되지 않은 PRNG와 동일한 성능을 달성합니다.


3 월 5 일 수정

여러분의 제안에 감사드립니다. 흥미로운 아이디어와 제안을 많이 받았습니다. 다음은 상위 항목입니다.

  • p가 1/200 대신 1/256의 해상도를 갖도록 문제 요구 사항을 변경하십시오. 이를 통해 비트를보다 효율적으로 사용할 수 있으며 최적화를위한 더 많은 기회를 얻을 수 있습니다. 이 변화를 만들 수있을 것 같습니다.
  • 산술 코딩을 사용하여 편향되지 않은 생성기에서 비트를 효율적으로 소비합니다. 위의 해상도 변경으로 훨씬 쉬워집니다.
  • 몇몇 사람들은 PRNG가 매우 빠르기 때문에 산술 코딩을 사용하면 도입 된 오버 헤드로 인해 실제로 코드가 느려질 수 있다고 제안했습니다. 대신 항상 최악의 경우 비트 수를 소비하고 해당 코드를 최적화해야합니다. 아래 벤치 마크를 참조하십시오.
  • @rici는 SIMD 사용을 제안했습니다. 이것은 우리가 항상 고정 된 수의 비트를 소비하는 경우에만 작동하는 좋은 아이디어입니다.

벤치 마크 (산술 디코딩 없음)

참고 : 많은 분들이 제안하신대로 해상도를 1/200에서 1/256으로 변경했습니다.

나는 단순히 8 개의 임의의 편향되지 않은 비트를 취하고 1 개의 편향된 비트를 생성하는 순진한 방법의 여러 구현작성 했습니다 .

  • SIMD없이
  • @rici 가 제안한대로 Agner Fog의 벡터 클래스 라이브러리를 사용하는 SIMD 사용
  • 내장 함수를 사용하는 SIMD 사용

나는 두 개의 편향되지 않은 의사 난수 생성기를 사용합니다.

또한 비교를 위해 편향되지 않은 PRNG의 속도를 측정합니다. 결과는 다음과 같습니다.


RNG: Ranvec1(Mersenne Twister for Graphics Processors + Multiply with Carry)

Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 16.081 16.125 16.093 [Gb/s]
Number of ones: 536,875,204 536,875,204 536,875,204
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency
Gbps/s: 0.778 0.783 0.812 [Gb/s]
Number of ones: 104,867,269 104,867,269 104,867,269
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 2.176 2.184 2.145 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 2.129 2.151 2.183 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical   : 104,857,600

SIMD는 스칼라 방법에 비해 성능을 3 배 향상시킵니다. 예상대로 편향되지 않은 생성기보다 8 배 느립니다.

가장 빠른 바이어스 생성기는 2.1Gb / s를 달성합니다.


RNG: xorshift128plus

Method: Unbiased with 1/1 efficiency (incorrect, baseline)
Gbps/s: 18.300 21.486 21.483 [Gb/s]
Number of ones: 536,867,655 536,867,655 536,867,655
Theoretical   : 104,857,600

Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 22.660 22.661 24.662 [Gb/s]
Number of ones: 536,867,655 536,867,655 536,867,655
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency
Gbps/s: 1.065 1.102 1.078 [Gb/s]
Number of ones: 104,868,930 104,868,930 104,868,930
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 4.972 4.971 4.970 [Gb/s]
Number of ones: 104,869,407 104,869,407 104,869,407
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 4.955 4.971 4.971 [Gb/s]
Number of ones: 104,869,407 104,869,407 104,869,407
Theoretical   : 104,857,600

xorshift의 경우 SIMD는 스칼라 방법에 비해 성능을 5 배 향상시킵니다. 편향되지 않은 생성기보다 4 배 느립니다. 이것은 xorshift의 스칼라 구현입니다.

가장 빠른 바이어스 생성기는 4.9Gb / s를 달성합니다.


RNG: xorshift128plus_avx2

Method: Unbiased with 1/1 efficiency (incorrect, baseline)
Gbps/s: 18.754 21.494 21.878 [Gb/s]
Number of ones: 536,867,655 536,867,655 536,867,655
Theoretical   : 104,857,600

Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 54.126 54.071 54.145 [Gb/s]
Number of ones: 536,874,540 536,880,718 536,891,316
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency
Gbps/s: 1.093 1.103 1.063 [Gb/s]
Number of ones: 104,868,930 104,868,930 104,868,930
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 19.567 19.578 19.555 [Gb/s]
Number of ones: 104,836,115 104,846,215 104,835,129
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 19.551 19.589 19.557 [Gb/s]
Number of ones: 104,831,396 104,837,429 104,851,100
Theoretical   : 104,857,600

이 구현은 AVX2를 사용하여 4 개의 편향되지 않은 xorshift 생성기를 병렬로 실행합니다.

가장 빠른 바이어스 발생기는 19.5Gb / s를 달성합니다.

산술 디코딩을위한 벤치 마크

간단한 테스트는 산술 디코딩 코드가 PRNG가 아니라 병목임을 보여줍니다. 그래서 저는 가장 비싼 PRNG만을 벤치마킹하고 있습니다.


RNG: Ranvec1(Mersenne Twister for Graphics Processors + Multiply with Carry)

Method: Arithmetic decoding (floating point)
Gbps/s: 0.068 0.068 0.069 [Gb/s]
Number of ones: 10,235,580 10,235,580 10,235,580
Theoretical   : 10,240,000

Method: Arithmetic decoding (fixed point)
Gbps/s: 0.263 0.263 0.263 [Gb/s]
Number of ones: 10,239,367 10,239,367 10,239,367
Theoretical   : 10,240,000

Method: Unbiased with 1/1 efficiency (incorrect, baseline)
Gbps/s: 12.687 12.686 12.684 [Gb/s]
Number of ones: 536,875,204 536,875,204 536,875,204
Theoretical   : 104,857,600

Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 14.536 14.536 14.536 [Gb/s]
Number of ones: 536,875,204 536,875,204 536,875,204
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency
Gbps/s: 0.754 0.754 0.754 [Gb/s]
Number of ones: 104,867,269 104,867,269 104,867,269
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 2.094 2.095 2.094 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical   : 104,857,600

Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 2.094 2.094 2.095 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical   : 104,857,600

단순 고정 소수점 방법은 0.25Gb / s를 달성하는 반면 순진 스칼라 방법은 3 배 빠르며 순진 SIMD 방법은 8 배 빠릅니다. 산술 디코딩 방법을 더 최적화 및 / 또는 병렬화하는 방법이있을 수 있지만 그 복잡성으로 인해 여기서 멈추고 순진한 SIMD 구현을 선택하기로 결정했습니다.

도움을 주셔서 감사합니다.


p256 개의 가능한 값을 기준으로 근사화 할 준비가 되어 있고 개별 비트가 서로 독립적 인 균일 한 값을 생성 할 수있는 PRNG가있는 경우 벡터화 된 비교를 사용하여 단일 난수에서 여러 편향된 비트를 생성 할 수 있습니다. .

이는 (1) 난수 품질에 대해 걱정하고 (2) 동일한 바이어스를 가진 많은 수의 비트가 필요할 경우에만 가치가 있습니다. 두 번째 요구 사항은 다음과 같이 제안 된 솔루션을 비판하는 원래 질문에 의해 암시 된 것 같습니다. "이 솔루션의 결함은 한 번에 8 비트 만 생성 할 수 있다는 것입니다. 작업이 많은 경우에도 편향되지 않은 PRNG 몇 가지 산술 명령으로 한 번에 64 개를 생성 할 수 있습니다. " 여기서 의미 는 단일 호출에서 바이어스 된 비트의 큰 블록을 생성하는 것이 유용 하다는 것입니다 .

난수 품질은 어려운 주제입니다. 측정이 불가능하지는 않더라도 어렵 기 때문에 다른 사람들은 "무작위성"의 다른 측면을 강조하거나 평가 절하하는 다른 측정 항목을 제안 할 것입니다. 일반적으로 낮은 "품질"을 위해 난수 생성 속도를 절충하는 것이 가능합니다. 이것이 가치가 있는지 여부는 정확한 응용 프로그램에 달려 있습니다.

난수 품질에 대한 가장 간단한 테스트는 개별 값의 분포와 생성기의주기 길이를 포함합니다. C 라이브러리 rand및 Posix random함수 의 표준 구현 은 일반적으로 배포 테스트를 통과하지만주기 길이는 장기 실행 애플리케이션에 적합하지 않습니다.

그러나 이러한 생성기는 일반적으로 매우 빠릅니다. glibc 구현 random에는 몇 개의 사이클 만 필요한 반면, 고전적인 선형 합동 생성기 (LCG)는 곱셈과 덧셈이 필요합니다. (또는 glibc 구현의 경우, 위의 3 개가 31 비트를 생성합니다.) 품질 요구 사항에 충분하다면, 특히 편향 확률이 자주 변경되는 경우 최적화를 시도 할 필요가 없습니다.

주기 길이는 예상되는 샘플 수보다 훨씬 길어야합니다. 이상적으로는이 숫자의 제곱보다 커야하므로 기가 바이트의 임의 데이터를 생성하려는 경우 주기 길이가 2 31 인 선형 합동 생성기 (LCG) 는 적합하지 않습니다. 주기 길이가 약 2 35 라고 주장되는 Gnu 3 항 비선형 가산 피드백 생성기조차도 수백만 개의 샘플이 필요한 애플리케이션에 사용해서는 안됩니다.

테스트하기 훨씬 더 어려운 또 다른 품질 문제는 연속 샘플의 독립성과 관련이 있습니다. 반복이 시작되면 생성 된 난수가 과거 값과 정확하게 연관되기 때문에 짧은주기 길이는이 메트릭에서 완전히 실패합니다. Gnu 삼항 알고리즘은주기가 더 길지만 생성 된 i 번째 난수 r i 는 항상 두 값 r i −3 + r i −31 또는 r 중 하나이기 때문에 명확한 상관 관계가 있습니다. 나는 −3 + r 나는 −31 +1. 이것은 놀랍 거나 적어도 혼란 스러울 수 있습니다. 특히 베르누이 실험의 결과.

다음은 Agner Fog의 유용한 벡터 클래스 라이브러리를 사용하여 구현 한 것 입니다. SSE 내장 함수의 성가신 세부 정보를 추상화하고 빠른 벡터화 된 난수 생성기 ( 아카이브 special.zip내부에 있음) 와 함께 제공되어 vectorclass.zip256 비트를 생성 할 수 있습니다. 256 비트 PRNG에 대한 8 번의 호출. 메르 센 트위스터조차도 품질 문제가있는 이유에 대한 Dr. Fog의 설명과 제안 된 솔루션을 읽을 수 있습니다. 나는 실제로 논평 할 자격이 없지만 적어도 내가 시도한 Bernoulli 실험에서 예상되는 결과를 제공하는 것처럼 보입니다.

#include "vectorclass/vectorclass.h"
#include "vectorclass/ranvec1.h"

class BiasedBits {
  public:
    // Default constructor, seeded with fixed values
    BiasedBits() : BiasedBits(1)  {}
    // Seed with a single seed; other possibilities exist.
    BiasedBits(int seed) : rng(3) { rng.init(seed); }

    // Generate 256 random bits, each with probability `p/256` of being 1.
    Vec8ui random256(unsigned p) {
      if (p >= 256) return Vec8ui{ 0xFFFFFFFF };
      Vec32c output{ 0 };
      Vec32c threshold{ 127 - p };
      for (int i = 0; i < 8; ++i) {
        output += output;
        output -= Vec32c(Vec32c(rng.uniform256()) > threshold);
      }
      return Vec8ui(output);
    }

  private:
    Ranvec1 rng;
};

내 테스트에서 260ms에서 268435456 비트 또는 나노 초당 1 비트를 생성하고 계산했습니다. 테스트 머신은 i5이므로 AVX2가 없습니다. YMMV.

실제 사용 사례에서에 대해 201 개의 가능한 값을 사용 p하면 8 비트 임계 값 계산이 성가 시게 정확하지 않습니다. 이러한 부정확성이 바람직하지 않은 경우, 두 배의 난수를 생성하는 대신 16 비트 임계 값을 사용하도록 위의 내용을 조정할 수 있습니다.

또는 10 비트 임계 값을 기반으로 벡터화를 핸드 롤링 할 수 있습니다. 그러면 10 비트마다 차용을 확인하여 벡터화 된 임계 값 비교를 수행하는 표준 비트 조작 해킹을 사용하여 0.5 % 증분까지 매우 좋은 근사치를 얻을 수 있습니다. 값 벡터와 반복 임계 값을 뺀 값입니다. 예를 들어,와 결합 std::mt19937_64하면 64 비트 난수마다 평균 6 비트가 제공됩니다.


한 가지 할 수있는 것은 기본 무 편향 생성기에서 여러 번 샘플링하여 32 비트 또는 64 비트 단어를 여러 번 얻은 다음 비트 부울 산술을 수행하는 것입니다. 예를 들어 4 단어의 b1,b2,b3,b4경우 다음 분포를 얻을 수 있습니다.

    표현 | p (비트는 1)
    ----------------------- + -------------
    b1 및 b2 및 b3 및 b4 | 6.25 %
    b1 및 b2 및 b3 | 12.50 %
    b1 및 b2 및 (b3 | b4) | 18.75 %
    b1 및 b2 | 25.00 %
    b1 | (b2 & (b3 | b4)) | 31.25 %
    b1 & (b2 | b3) | 37.50 %
    b1 & (b2 | b3 | b4)) | 43.75 %
    b1 | 50.00 %

더 미세한 해상도를 위해 유사한 구조를 만들 수 있습니다. 약간 지루하고 여전히 더 많은 생성기 호출이 필요하지만 적어도 비트 당 하나는 필요하지 않습니다. 이것은 a3f의 대답과 비슷하지만 구현하기가 더 쉽고 단어를 스캔하는 것보다 빠릅니다 0xF.

원하는 0.5 % 해상도의 경우 하나의 편향된 단어에 대해 8 개의 편향되지 않은 단어가 필요하며, 이는 (0.5 ^ 8) = 0.390625 %의 해상도를 제공합니다.


정보 이론적 관점에서, 편향된 비트 스트림 (가있는 p != 0.5)은 편향되지 않은 스트림보다 정보 으므로 이론적 으로는 편향되지 않은 입력의 1 비트 미만 을 (평균적으로) 단일 비트를 생성해야합니다. 편향된 출력 스트림의. 예를 들어, 베르누이 랜덤 변수 엔트로피p = 0.1-0.1 * log2(0.1) - 0.9 * log2(0.9)비트 주변에있는 0.469비트입니다. 이는 경우에 p = 0.1대해 편향되지 않은 입력 비트 당 출력 스트림의 2 비트를 약간 넘게 생성 할 수 있어야 함을 시사합니다 .

아래에서는 편향된 비트를 생성하는 두 가지 방법을 제공합니다. 둘 다 가능한 한 입력 무 편향 비트를 적게 요구한다는 의미에서 최적의 효율성에 가깝습니다.

방법 1 : 산술 (디코딩) 코딩

실용적인 방법은 이미 alexis답변에 설명 된 것처럼 산술 (디) 코딩을 사용하여 편향되지 않은 입력 스트림을 디코딩 하는 것입니다 . 이 간단한 경우에는 코드를 작성하는 것이 어렵지 않습니다. 다음은이를 수행하는 최적화되지 않은 의사 코드 ( cough, Python )입니다.

import random

def random_bits():
    """
    Infinite generator generating a stream of random bits,
    with 0 and 1 having equal probability.
    """
    global bit_count  # keep track of how many bits were produced
    while True:
        bit_count += 1
        yield random.choice([0, 1])

def bernoulli(p):
    """
    Infinite generator generating 1-bits with probability p
    and 0-bits with probability 1 - p.
    """
    bits = random_bits()

    low, high = 0.0, 1.0
    while True:
        if high <= p:
            # Generate 1, rescale to map [0, p) to [0, 1)
            yield 1
            low, high = low / p, high / p
        elif low >= p:
            # Generate 0, rescale to map [p, 1) to [0, 1)
            yield 0
            low, high = (low - p) / (1 - p), (high - p) / (1 - p)
        else:
            # Use the next random bit to halve the current interval.
            mid = 0.5 * (low + high)
            if next(bits):
                low = mid
            else:
                high = mid

다음은 사용 예입니다.

import itertools
bit_count = 0

# Generate a million deviates.
results = list(itertools.islice(bernoulli(0.1), 10**6))

print("First 50:", ''.join(map(str, results[:50])))
print("Biased bits generated:", len(results))
print("Unbiased bits used:", bit_count)
print("mean:", sum(results) / len(results))

위는 다음과 같은 샘플 출력을 제공합니다.

First 50: 00000000000001000000000110010000001000000100010000
Biased bits generated: 1000000
Unbiased bits used: 469036
mean: 0.100012

약속 한대로 소스 비 편향 스트림에서 오십 만 개 미만을 사용하여 출력 편향 스트림의 백만 비트를 생성했습니다.

최적화를 위해이를 C / C ++로 변환 할 때 부동 소수점 대신 정수 기반 고정 소수점 산술을 사용하여이를 코딩하는 것이 합리적 일 수 있습니다.

방법 2 : 정수 기반 알고리즘

정수를 직접 사용하도록 산술 디코딩 방법을 변환하는 대신 여기에 더 간단한 방법이 있습니다. 더 이상 산술 디코딩은 아니지만 완전히 관련이 없으며 위의 부동 소수점 버전과 동일한 출력 편향 비트 / 입력 편향 비트 비율에 가깝습니다. 모든 수량이 부호없는 32 비트 정수에 맞도록 구성되어 있으므로 C / C ++로 쉽게 변환 할 수 있습니다. 이 코드 p는가의 정확한 배수 인 경우에 특화되어 1/200있지만이 방법은 p분모가 합리적으로 작은 유리수로 표현할 수있는 모든 경우에 적용됩니다.

def bernoulli_int(p):
    """
    Infinite generator generating 1-bits with probability p
    and 0-bits with probability 1 - p.

    p should be an integer multiple of 1/200.
    """
    bits = random_bits()
    # Assuming that p has a resolution of 0.05, find p / 0.05.
    p_int = int(round(200*p))

    value, high = 0, 1
    while True:
        if high < 2**31:
            high = 2 * high
            value = 2 * value + next(bits)
        else:
            # Throw out everything beyond the last multiple of 200, to
            # avoid introducing a bias.
            discard = high - high % 200
            split = high // 200 * p_int
            if value >= discard:  # rarer than 1 time in 10 million
                value -= discard
                high -= discard
            elif value >= split:
                yield 0
                value -= split
                high = discard - split
            else:
                yield 1
                high = split

핵심 관찰은 while루프 의 시작에 도달 할 때마다의 value모든 정수에 균일하게 분포 [0, high)되고 이전에 출력 된 모든 비트와 독립적이라는 것입니다. 당신은 완벽한 정확성보다 더 빠른 속도에 대해 걱정하는 경우, 당신은 제거 할 수 discardvalue >= discard지점 : 바로 거기가 우리의 출력을 보장 0하고 1딱 맞는 확률로합니다. 그 복잡함은 그대로두면 대신 거의 올바른 확률을 얻을 수 있습니다. 또한 해상도를 보다 p같게 만들면 잠재적으로 시간이 많이 걸리는 분할 및 모듈로 연산이 비트 연산으로 대체 될 수 있습니다.1/2561/200

이전과 동일한 테스트 코드를 사용하지만 bernoulli_int대신 사용 bernoulli하면 다음과 같은 결과가 나타납니다 p=0.1.

First 50: 00000010000000000100000000000000000000000110000100
Biased bits generated: 1000000
Unbiased bits used: 467997
mean: 0.099675

1이 나타날 확률이 6,25 % (1/16)라고 가정 해 보겠습니다. 4 비트 번호에 대해 16 개의 가능한 비트 패턴이 있습니다 0000,0001, ..., 1110,1111.

이제 예전과 같은 난수를 생성하고 1111니블 경계에서 1모든 것을 a 로 바꾸고 나머지는 모두 0.

다른 확률에 맞게 조정하십시오.


이론적으로 최적의 동작을 얻을 수 있습니다. 즉 , 난수 생성기를 최소한으로 사용 하고 산술 코딩을 사용하여 이에 접근하면 확률 p를 정확하게 모델링 할 수 있습니다 .

산술 코딩은 메시지를 숫자 범위의 하위 간격으로 나타내는 데이터 압축 형식입니다. 이론적으로 최적의 인코딩을 제공하며 각 입력 기호에 대해 소수 비트 수를 사용할 수 있습니다.

아이디어는 다음과 같습니다 . 확률이 p 인 1 인 무작위 비트 시퀀스가 ​​있다고 가정 해보십시오 . 편의상 대신 비트가 0이 될 확률을 사용 합니다. ( q = 1-p ). 산술 코딩은 숫자 범위의 각 비트 부분에 할당됩니다. 첫 번째 비트의 경우 입력이 0이면 간격 [0, q)를 할당하고 입력이 1이면 간격 [q, 1)을 할당합니다. 이후 비트는 현재 범위의 비례 하위 간격을 할당합니다. 예를 들어, q = 1/3 이라고 가정합니다 . 입력 1 0 0 은 다음과 같이 인코딩됩니다.q

Initially       [0, 1),             range = 1
After 1         [0.333, 1),         range = 0.6666        
After 0         [0.333, 0.5555),    range = 0.2222   
After 0         [0.333, 0.407407),  range = 0.074074

첫 번째 숫자 인 1 은 범위의 상위 2/3 (1-q)를 선택합니다. 두 번째 숫자 인 0 아래 1/3을 선택하는 식 입니다. 첫 번째 단계와 두 번째 단계 후에 간격은 중간 지점에 걸쳐 있습니다. 그러나 세 번째 단계 이후에는 완전히 중간 점 아래에 있으므로 첫 번째 압축 된 숫자가 출력 될 수 있습니다 0.. 프로세스가 계속되고 특수 EOF기호가 종료 자로 추가됩니다.

이것이 당신의 문제와 무슨 관련이 있습니까? 압축 된 출력은 동일한 확률로 임의의 0과 1을 갖습니다. 따라서 확률이 p비트를 얻으 려면 RNG의 출력이 위와 같이 산술 코딩의 결과라고 가정 하고 디코더 프로세스를 적용하십시오. 즉, 라인 간격을 더 작은 조각으로 세분화하는 것처럼 비트를 읽습니다. 예를 들어 01RNG에서 읽은 후에는 [0.25, 0.5) 범위에있게됩니다 . 충분한 출력이 "디코딩"될 때까지 비트를 계속 읽으십시오. 압축 해제를 모방하고 있기 때문에 입력 한 것보다 더 많은 무작위 비트를 얻을 수 있습니다.산술 코딩은 이론적으로 최적이기 때문에 임의성을 희생하지 않고 RNG 출력을 더 편향된 비트로 전환 할 수있는 방법은 없습니다. 진정한 최대 값을 얻습니다.

몇 줄의 코드로는이 작업을 수행 할 수 없으며, 제가 가리킬 수있는 라이브러리를 알지 못합니다 (사용할 수있는 라이브러리가 있어야 함). 그래도 꽤 간단합니다. 위의 문서는 C. 그것은 꽤 간단에서 범용 인코더 및 디코더에 대한 코드를 제공하며, 임의의 확률로 다수의 입력 심볼을 지원합니다; 귀하의 경우 에는 확률 모델이 사소하기 때문에 훨씬 더 간단한 구현이 가능합니다 ( Mark Dickinson의 대답이 지금 보여 주듯이). 확장 된 사용의 경우 각 비트에 대해 많은 부동 소수점 계산을 수행하지 않는 강력한 구현을 생성하려면 조금 더 많은 작업이 필요합니다.

또한 Wikipedia 에는 작업을 보는 또 다른 방법 인 기수 변경으로 간주되는 산술 인코딩에 대한 흥미로운 논의가 있습니다.


어, 의사 난수 생성기는 일반적으로 매우 빠릅니다. 이것이 어떤 언어인지 (아마도 Python) 확실하지 않지만 "result.append"(거의 확실히 메모리 할당을 포함 함)는 "random_uniform"(약간의 수학 만 수행함)보다 느립니다.

이 코드의 성능을 최적화하려면 :

  1. 문제인지 확인하십시오. 최적화는 약간의 작업이며 코드를 유지하기 어렵게 만듭니다. 필요한 경우가 아니면하지 마십시오.
  2. 프로파일 링하십시오. 코드의 어떤 부분이 실제로 가장 느린 지 확인하려면 몇 가지 테스트를 실행하십시오. 속도를 높이는 데 필요한 부분입니다.
  3. 변경하고 실제로 더 빠른지 확인하십시오. 컴파일러는 꽤 똑똑합니다. 종종 명확한 코드는 더 빨리 보이는 것보다 복잡한 코드로 컴파일됩니다.

컴파일 된 언어 (JIT 컴파일 포함)로 작업하는 경우 모든 제어 전송 (if, while, 함수 호출 등)에 대해 성능 저하가 발생합니다. 당신이 할 수있는 것을 제거하십시오. 메모리 할당도 (보통) 상당히 비쌉니다.

통역 언어로 작업하는 경우 모든 베팅이 해제됩니다. 가장 간단한 코드가 최고 일 가능성이 높습니다. 인터프리터의 오버 헤드는 당신이하고있는 모든 일을 축소 시키므로 가능한 한 작업을 줄이십시오.

성능 문제가 어디에 있는지 추측 할 수 있습니다.

  1. 메모리 할당. 배열을 전체 크기로 미리 할당하고 나중에 항목을 채 웁니다. 이렇게하면 항목을 추가하는 동안 메모리를 재 할당 할 필요가 없습니다.
  2. 가지. 결과 또는 유사한 것을 캐스팅하여 "if"를 피할 수 있습니다. 이것은 컴파일러에 따라 많이 달라집니다. 어셈블리 (또는 프로파일)를 확인하여 원하는 작업을 수행하는지 확인하십시오.
  3. 숫자 유형. 난수 생성기가 기본적으로 사용하는 유형을 찾고 해당 유형으로 산술을 수행하십시오. 예를 들어 생성기가 자연스럽게 32 비트 부호없는 정수를 반환하는 경우 먼저 "p"를 해당 범위로 조정 한 다음 비교에 사용합니다.

그건 그렇고, 가능한 한 최소한의 무작위성을 사용하고 싶다면 "산술 코딩"을 사용하여 임의 스트림을 디코딩하십시오. 빠르지 않을 것입니다.


정확한 결과를 제공하는 한 가지 방법은 먼저 이항 분포 다음의 1 비트 수를 k 비트 블록에 대해 무작위로 생성 한 다음 여기 에있는 방법 중 하나를 사용하여 정확히 그 수의 비트로 k 비트 워드를 생성하는 입니다. 예를 들어 mic006에 의한 방법은 k-bit 난수에 대한 로그 k 비트 만 필요하고 광산은 하나만 필요합니다.


p가 0에 가까우면 n 번째 비트가 1 인 첫 번째 비트 일 확률을 계산할 수 있습니다. 그런 다음 0과 1 사이의 난수를 계산하고 그에 따라 n을 선택합니다. 예를 들어 p = 0.005 (0.5 %)이고 난수가 0.638128이면 n = 321을 계산할 수 있으므로 321 0 비트와 1 비트 세트로 채 웁니다.

p가 1에 가까우면 p 대신 1-p를 사용하고 1 비트와 1 개의 0 비트를 설정합니다.

p가 1 또는 0에 가까우 지 않으면 8 비트로 구성된 256 개 시퀀스의 표를 만들고 누적 확률을 계산 한 다음 난수를 구하고 누적 확률 배열에서 이진 검색을 수행하고 8 비트를 설정할 수 있습니다. .


임의 비트 생성기에 액세스 할 수 있다고 가정하면 p비트별로 비교할 값을 생성하고 생성 된 값이보다 작거나 크거나 같음을 증명할 수있는 즉시 중단 할 수 있습니다 p.

주어진 확률로 스트림에 하나의 항목을 생성하려면 다음과 같이 진행하십시오 p.

  1. 0.바이너리로 시작
  2. 임의의 비트를 추가하십시오. a 1가 그려 졌다고 가정 하면0.1
  3. 결과 (이진 표기법)가보다 작음이 입증 p되면 다음을 출력합니다.1
  4. 결과가보다 크거나 같으면 다음을 p출력합니다.0
  5. 그렇지 않으면 (둘 다 배제 할 수없는 경우) 2 단계로 진행하십시오.

p이진 표기법에서 다음과 같다고 가정 해 봅시다 0.1001101.... 이 공정은 임의의 발생하는 경우 0.0, 0.1000, 0.10010, ..., 값이 커지거나 같 수 p이상; 어느 한 경우에 0.11, 0.101, 0.100111, ... 생성 된 값보다 작아 질 수 없다 p.

나에게이 방법은 예상대로 약 두 개의 임의 비트를 사용하는 것처럼 보입니다. 산술 코딩 (Mark Dickinson의 답변 에서 볼 수 있음 )은 고정 된 경우 편향된 비트 (평균) 당 최대 1 개의 임의 비트를 소비합니다 p. 수정 비용 p은 명확하지 않습니다.


그것이하는 일

이 구현은 "/ dev / urandom"특수 문자 파일의 인터페이스를 통해 임의의 장치 커널 모듈을 단일 호출하여 주어진 해상도에서 모든 값을 나타내는 데 필요한 임의의 데이터 수를 가져옵니다. 가능한 최대 해상도는 1 / 256 ^ 2이므로 0.005는 다음과 같이 나타낼 수 있습니다.

328 / 256 ^ 2,

즉 :

해결 : 256 * 256

x : 328

오류 0.000004883.

어떻게 하는가

구현 bits_per_byte은 주어진 해상도를 처리하는 데 필요한 균일하게 분산 된 비트 수인 비트 수를 계산 @resolution합니다. 즉, 모든 값을 나타냅니다 . 그런 다음 임의 화 장치에 대한 단일 호출을 수행합니다 ( URANDOM_DEVICE정의 된 경우 "/ dev / urandom" , 그렇지 않으면 비트에 엔트로피가 충분하지 않은 경우 차단할 수있는 "/ dev / random"호출을 통해 장치 드라이버의 추가 노이즈를 사용합니다). 균일하게 분산 된 필요한 수의 바이트를 얻고 rnd_bytes임의의 바이트 배열 채 웁니다 . 마지막으로 rnd_bytes 배열의 각 bytes_per_byte 바이트에서 각 Bernoulli 샘플 당 필요한 비트 수를 읽고이 비트의 정수 값을에서 주어진 단일 Bernoulli 결과의 성공 확률과 비교합니다 x/resolution. 값이 맞으면, 즉x/resolution 길이를 임의로 [0, x / 해상도) 세그먼트로 선택한 다음 성공을 확인하고 결과 배열에 1을 삽입합니다.


임의의 장치에서 읽기 :

/* if defined use /dev/urandom (will not block),
 * if not defined use /dev/random (may block)*/
#define URANDOM_DEVICE 1

/*
 * @brief   Read @outlen bytes from random device
 *          to array @out.
 */
int
get_random_samples(char *out, size_t outlen)
{
    ssize_t res;
#ifdef URANDOM_DEVICE
    int fd = open("/dev/urandom", O_RDONLY);
    if (fd == -1) return -1;
    res = read(fd, out, outlen);
    if (res < 0) {
        close(fd);
        return -2;
    }
#else
    size_t read_n;
    int fd = open("/dev/random", O_RDONLY);
    if (fd == -1) return -1;
    read_n = 0;
    while (read_n < outlen) {
       res = read(fd, out + read_n, outlen - read_n);
       if (res < 0) {
           close(fd);
           return -3;
       }
       read_n += res;
    }
#endif /* URANDOM_DEVICE */
    close(fd);
    return 0;
}

Bernoulli 샘플로 구성된 벡터 채우기 :

/*
 * @brief   Draw vector of Bernoulli samples.
 * @details @x and @resolution determines probability
 *          of success in Bernoulli distribution
 *          and accuracy of results: p = x/resolution.
 * @param   resolution: number of segments per sample of output array 
 *          as power of 2: max resolution supported is 2^24=16777216
 * @param   x: determines used probability, x = [0, resolution - 1]
 * @param   n: number of samples in result vector
 */
int
get_bernoulli_samples(char *out, uint32_t n, uint32_t resolution, uint32_t x)
{
    int res;
    size_t i, j;
    uint32_t bytes_per_byte, word;
    unsigned char *rnd_bytes;
    uint32_t uniform_byte;
    uint8_t bits_per_byte;

    if (out == NULL || n == 0 || resolution == 0 || x > (resolution - 1))
        return -1;

    bits_per_byte = log_int(resolution);
    bytes_per_byte = bits_per_byte / BITS_PER_BYTE + 
                        (bits_per_byte % BITS_PER_BYTE ? 1 : 0);
    rnd_bytes = malloc(n * bytes_per_byte);
    if (rnd_bytes == NULL)
        return -2;
    res = get_random_samples(rnd_bytes, n * bytes_per_byte);
    if (res < 0)
    {
        free(rnd_bytes);
        return -3;
    }

    i = 0;
    while (i < n)
    {
        /* get Bernoulli sample */
        /* read byte */
        j = 0;
        word = 0;
        while (j < bytes_per_byte)
        {
            word |= (rnd_bytes[i * bytes_per_byte + j] << (BITS_PER_BYTE * j));
            ++j;
        }
        uniform_byte = word & ((1u << bits_per_byte) - 1);
        /* decision */
        if (uniform_byte < x)
            out[i] = 1;
        else
            out[i] = 0;
        ++i;
    }

    free(rnd_bytes);    
    return 0;
}

용법:

int
main(void)
{
    int res;
    char c[256];

    res = get_bernoulli_samples(c, sizeof(c), 256*256, 328); /* 328/(256^2) = 0.0050 */
    if (res < 0) return -1;

    return 0;
}

완전한 코드 , 결과 .

참조 URL : https://stackoverflow.com/questions/35795110/fast-way-to-generate-pseudo-random-bits-with-a-given-probability-of-0-or-1-for-e

반응형