가우스를 통합하는 Python 함수를 작성하는 가장 좋은 방법은 무엇입니까?

StackOverflow https://stackoverflow.com/questions/509994

문제

가우스를 통합하기 위해 scipy의 쿼드 방법을 사용하려고 할 때(가우스라는 가우스 방법이 있다고 가정해 보겠습니다) 필요한 매개변수를 가우스에 전달하고 쿼드를 떠나 올바른 변수에 대한 통합을 수행하는 데 문제가 있었습니다.다차원 함수와 함께 쿼드를 사용하는 방법에 대한 좋은 예가 있는 사람이 있습니까?

그러나 이로 인해 일반적으로 가우스를 통합하는 가장 좋은 방법에 대한 더 큰 질문이 생겼습니다.놀랍게도 scipy에서 가우스 통합을 찾지 못했습니다.내 계획은 간단한 가우스 함수를 작성하여 이를 쿼드(또는 고정 폭 적분기)에 전달하는 것이었습니다.당신은 무엇을 하시겠습니까?

편집하다:고정 너비는 고정 dx를 사용하여 곡선 아래 영역을 계산하는 Trapz와 같은 것을 의미합니다.

지금까지 내가 찾아낸 것은 쿼드로 들어갈 수 있는 람다 함수를 반환하는 make___gauss 메서드입니다.이렇게 하면 통합하기 전에 필요한 평균과 분산으로 일반 함수를 만들 수 있습니다.

def make_gauss(N, sigma, mu):
    return (lambda x: N/(sigma * (2*numpy.pi)**.5) *
            numpy.e ** (-(x-mu)**2/(2 * sigma**2)))

quad(make_gauss(N=10, sigma=2, mu=0), -inf, inf)

일반 가우스 함수(x, N, mu 및 sigma로 호출해야 함)를 전달하고 다음과 같은 쿼드를 사용하여 일부 값을 채우려고 시도했을 때

quad(gen_gauss, -inf, inf, (10,2,0))

매개변수 10, 2, 0은 반드시 N=10, sigma=2, mu=0과 일치하지 않았으므로 더 확장된 정의가 필요했습니다.

scipy.special의 erf(z)를 사용하려면 처음에 t가 무엇인지 정확히 정의해야 하지만 그것이 거기에 있다는 것을 알 수 있어서 좋습니다.

도움이 되었습니까?

해결책

좋아요, 당신은 몇 가지 사항에 대해 상당히 혼란스러워하는 것 같습니다.처음부터 시작해 보겠습니다.당신은 "다차원 함수"에 대해 언급했지만 계속해서 일반적인 단일 변수 가우스 곡선에 대해 논의합니다.이것은 ~ 아니다 다차원 함수:이를 통합하면 하나의 변수(x)만 통합됩니다.있기 때문에 구별하는 것이 중요합니다. ~이다 진정한 다차원 함수인 "다변량 가우스 분포"라는 괴물은 통합할 경우 두 개 이상의 변수를 통합해야 합니다(이전에 언급한 값비싼 몬테 카를로 기술을 사용함).하지만 당신은 작업, 통합 등이 훨씬 쉬운 일반적인 단일 변수 가우스에 대해 이야기하고 있는 것 같습니다.

단일 변수 가우스 분포에는 두 개의 매개변수가 있습니다. sigma 그리고 mu, 이며 우리가 표시할 단일 변수의 함수입니다. x.또한 정규화 매개변수를 가지고 다니는 것 같습니다. n (이는 몇 가지 응용 프로그램에 유용합니다).정규화 매개변수는 일반적으로 ~ 아니다 마지막에 다시 추가할 수 있으므로 계산에 포함됩니다(적분은 선형 연산자라는 점을 기억하세요. int(n*f(x), x) = n*int(f(x), x) ).하지만 원하시면 가지고 다닐 수 있습니다.제가 정규 분포에 대해 좋아하는 표기법은 다음과 같습니다.

N(x | mu, sigma, n) := (n/(sigma*sqrt(2*pi))) * exp((-(x-mu)^2)/(2*sigma^2))

(이것을 "정규 분포"라고 읽으십시오. x 주어진 sigma, mu, 그리고 n ...에 의해 제공됩니다.") 지금까지는 매우 좋습니다.이것은 당신이 가지고 있는 기능과 일치합니다.유일한 실제 변수 여기는 x:나머지 세 가지 매개변수는 다음과 같습니다. 결정된 특정 가우스에 대해.

이제 수학적 사실을 살펴보겠습니다.모든 가우스 곡선이 동일한 모양을 갖고 있다는 것은 증명된 사실입니다. 단지 약간만 이동했을 뿐입니다.그래서 우리는 함께 일할 수 있습니다 N(x|0,1,1), "표준 정규 분포"라고 하며 결과를 다시 일반 가우스 곡선으로 변환합니다.그래서 만약 당신이 적분을 가지고 있다면 N(x|0,1,1), 가우시안의 적분을 간단하게 계산할 수 있습니다.이 적분은 너무 자주 나타나므로 특별한 이름이 있습니다.그만큼 오류 기능 erf.일부 오래된 관례 때문에 그렇지 않습니다. 정확히 erf;몇 가지 덧셈과 곱셈 요소도 함께 사용됩니다.

만약에 Phi(z) = integral(N(x|0,1,1), -inf, z);그건, Phi(z) 는 마이너스 무한대에서 최대 표준 정규 분포의 적분입니다. z, 그러면 오류 함수의 정의에 따르면 이는 참입니다.

Phi(z) = 0.5 + 0.5 * erf(z / sqrt(2)).

마찬가지로 만약에 Phi(z | mu, sigma, n) = integral( N(x|sigma, mu, n), -inf, z);그건, Phi(z | mu, sigma, n) 주어진 매개변수에 정규 분포를 적분한 것입니다. mu, sigma, 그리고 n 마이너스 무한대부터 z, 그러면 오류 함수의 정의에 따르면 이는 참입니다.

Phi(z | mu, sigma, n) = (n/2) * (1 + erf((x - mu) / (sigma * sqrt(2)))).

보세요 일반 CDF의 Wikipedia 기사 이 사실에 대한 더 자세한 내용이나 증거를 원한다면.

좋아요, 배경 설명은 이것으로 충분하겠습니다.(수정된) 게시물로 돌아갑니다."scipy.special의 erf(z)를 사용하려면 처음에 t가 무엇인지 정확히 정의해야 합니다"라고 말합니다.나는 이것이 무엇을 의미하는지 전혀 모릅니다.어디 있니? t (시간?) 이것에 전혀 들어가나요?위의 설명을 통해 오류 함수에 대한 이해가 조금 풀렸기를 바라며 이제 오류 함수가 해당 작업에 적합한 함수인 이유가 더 명확해졌기를 바랍니다.

귀하의 Python 코드는 괜찮지만 저는 람다보다 클로저를 선호합니다.

def make_gauss(N, sigma, mu):
    k = N / (sigma * math.sqrt(2*math.pi))
    s = -1.0 / (2 * sigma * sigma)
    def f(x):
        return k * math.exp(s * (x - mu)*(x - mu))
    return f

클로저를 사용하면 상수를 미리 계산할 수 있습니다. k 그리고 s, 따라서 반환된 함수는 호출될 때마다 더 적은 작업을 수행해야 합니다. 이는 함수를 통합하는 경우 중요할 수 있으며 이는 여러 번 호출된다는 의미입니다.또한 지수 연산자의 사용을 피했습니다. **, 이는 단순히 제곱을 작성하는 것보다 느리고 내부 루프에서 나누기를 끌어올려 곱셈으로 대체했습니다.Python에서의 구현을 전혀 살펴보지 않았지만 원시 x87 어셈블리를 사용하여 순수한 속도를 위해 내부 루프를 조정한 마지막 시간부터 더하기, 빼기 또는 곱하기에는 각각 약 4 CPU 사이클이 소요되고 나누는 것을 기억하는 것 같습니다. 36이고 지수는 약 200입니다.그것은 몇 년 전의 일이므로 그 숫자를 소금 한 알로 받아들이십시오.그래도 상대적인 복잡성을 보여줍니다.계산도 그렇고 exp(x) 무차별 방식은 매우 나쁜 생각입니다.좋은 구현을 작성할 때 취할 수 있는 요령이 있습니다. exp(x) 일반보다 훨씬 빠르고 정확합니다. a**b 스타일 지수화.

나는 상수 pi와 e의 numpy 버전을 사용한 적이 없습니다.나는 항상 평범하고 오래된 수학 모듈 버전을 고수해 왔습니다.왜 둘 중 하나를 선호하는지 모르겠습니다.

당신이 그걸로 무엇을 하려는지 잘 모르겠습니다 quad() 부르다. quad(gen_gauss, -inf, inf, (10,2,0)) 마이너스 무한대에서 플러스 무한대까지 재정규화된 가우스를 적분해야 하며, 가우스가 실제 선에 대해 1로 적분되므로 항상 10(정규화 인자)을 내야 합니다.10에서 멀리 떨어진 대답은 (나는 기대하지 않습니다 정확히 10 이후 quad() 결국 근사치일 뿐입니다) 어딘가에 문제가 있다는 의미입니다...실제 반환 값과 내부 작동 방식을 알지 못하면 무엇이 문제인지 말하기 어렵습니다. quad().

이것이 혼란을 어느 정도 해소하고 오류 함수가 문제에 대한 정답인 이유와 궁금한 경우 직접 수행하는 방법을 설명했기를 바랍니다.설명이 명확하지 않은 경우 먼저 Wikipedia를 잠깐 살펴보는 것이 좋습니다.여전히 궁금한 점이 있으면 주저하지 말고 문의하세요.

다른 팁

Scipy는 "오류 함수", AKA Gaussian Integral과 함께 배송합니다.

import scipy.special
help(scipy.special.erf)

나는 당신이 다변량 가우시안을 처리하고 있다고 가정합니다. 그렇다면 Scipy는 이미 원하는 기능을 가지고 있습니다. MVNDIST ( "다변량 정규 분포)라고합니다. SCIPY 문서는 그 어느 때보다도 끔찍하므로 기능이 묻힌 위치조차 찾을 수는 없습니다. 어딘가에 있습니다. 이 문서는 Scipy의 최악의 부분이며 과거에는 끝나지 않았습니다.

단일 변수 가우시안은 많은 구현을 사용할 수있는 좋은 오래된 오류 기능을 사용합니다.

James Thompson이 언급 한 것처럼 일반적으로 문제를 공격하는 경우, 자신의 가우스 분포 기능을 작성하여 Quad ()에 공급하고 싶습니다. 그러나 일반화 된 통합을 피할 수 있다면 그렇게하는 것이 좋습니다. 특정 기능을위한 특수 통합 기술 (MVNDIST 사용과 같은)은 표준 Monte Carlo 다중 차원 통합보다 훨씬 빠릅니다. 높은 정확도를 위해.

가우스 분포를 정규 분포라고도합니다. Scipy Norm 모듈의 CDF 기능은 원하는 작업을 수행합니다.

from scipy.stats import norm
print norm.cdf(0.0)
>>>0.5

http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html#scipy.stats.norm

왜 항상 -infinity에서 +infinity로의 통합을하지 않으므로 항상 답을 알 수 있습니까? (농담!)

내 생각에 Scipy에 통조림 가우시안 기능이 아직 없다는 유일한 이유는 그것이 쓰는 것이 사소한 기능이기 때문입니다. 자신의 기능을 작성하고 쿼드로 전달하여 탁월한 소리를냅니다. 이 작업을 수행하기 위해 허용되는 Scipy 도구를 사용합니다. 최소한의 코드 노력이며 Scipy를 본 적이 없어도 다른 사람들에게는 매우 읽을 수 있습니다.

고정형 통합기가 정확히 무엇을 의미합니까? Quadpack이 사용하는 것과 다른 알고리즘을 사용한다는 것을 의미합니까?

편집 : 완전성을 위해 다음은 평균 0과 표준 편차가 0에서 +무한대의 가우스를 위해 시도한 것과 같은 것입니다.

from scipy.integrate import quad
from math import pi, exp
mean = 0
sd   = 1
quad(lambda x: 1 / ( sd * ( 2 * pi ) ** 0.5 ) * exp( x ** 2 / (-2 * sd ** 2) ), 0, inf )

가우스 기능이 약간 길지만 여전히 글을 쓰는 것이 매우 사소하기 때문에 약간 추악합니다.

라이센스 : CC-BY-SA ~와 함께 속성
제휴하지 않습니다 StackOverflow
scroll top