실행 표준 편차를 효율적으로 계산하는 방법은 무엇입니까?
예를 들어 다음과 같은 숫자 목록이 있습니다.
[0] (0.01, 0.01, 0.02, 0.04, 0.03)
[1] (0.00, 0.02, 0.02, 0.03, 0.02)
[2] (0.01, 0.02, 0.02, 0.03, 0.02)
...
[n] (0.01, 0.00, 0.01, 0.05, 0.03)
내가하고 싶은 것은 모든 배열 요소에서 목록의 각 인덱스에서 평균과 표준 편차를 효율적으로 계산하는 것입니다.
평균을 내기 위해 배열을 반복하고 목록의 지정된 인덱스에서 값을 합산했습니다. 마지막으로 "평균 목록"의 각 값을 n(집단의 표본이 아닌 모집단으로 작업하고 있습니다)로 나눕니다 .
표준 편차를 수행하기 위해 평균을 계산 했으므로 다시 반복합니다.
나는 배열을 두 번, 평균에 대해 한 번, SD에 대해 한 번 (평균을 얻은 후) 피하고 싶습니다.
두 값을 모두 계산하는 효율적인 방법이 있습니까? 배열을 한 번만 통과합니까? 해석 된 언어 (예 : Perl 또는 Python) 또는 의사 코드로 된 코드는 괜찮습니다.
대답은 Welford의 알고리즘을 사용하는 것입니다. 이는 다음의 "순진한 방법"뒤에 매우 명확하게 정의되어 있습니다.
- Wikipedia : 분산 계산 알고리즘
다른 응답에서 제안 된 2- 패스 또는 온라인 단순 제곱합 수집기보다 수치 적으로 더 안정적입니다. 안정성 은 부동 소수점 문헌에서 " 치명적 취소 " 로 알려진 결과로 이어지기 때문에 서로 가까운 값이 많은 경우에만 실제로 중요합니다 .
또한 분산 계산 (제곱 편차)에서 샘플 수 (N)로 나누는 것과 N-1의 차이를 살펴볼 수 있습니다. N-1로 나누면 표본에서 편향되지 않은 분산 추정값이 나오는 반면, 평균 N으로 나누면 분산이 과소 평가됩니다 (표본 평균과 실제 평균 사이의 분산을 고려하지 않기 때문).
온라인에서 이전 값을 삭제하는 방법을 포함하여 자세한 내용을 다루는 주제에 대한 두 개의 블로그 항목을 작성했습니다.
내 Java 구현을 살펴볼 수도 있습니다. javadoc, 소스 및 단위 테스트는 모두 온라인입니다.
- Javadoc :
stats.OnlineNormalEstimator - 출처:
stats.OnlineNormalEstimator.java - JUnit 소스 :
test.unit.stats.OnlineNormalEstimatorTest.java - LingPipe 홈 페이지
기본적인 대답은 x ( 'sum_x1'이라고 부름)와 x 2 ( 'sum_x2'라고 부름)의 합계를 모으는 것입니다. 표준 편차의 값은 다음과 같습니다.
stdev = sqrt((sum_x2 / n) - (mean * mean))
어디
mean = sum_x / n
이것은 표본 표준 편차입니다. 제수로 'n-1'대신 'n'을 사용하여 모집단 표준 편차를 얻습니다.
큰 표본을 다루는 경우 두 큰 숫자의 차이를 가져 오는 수치 안정성에 대해 걱정해야 할 수도 있습니다. 자세한 내용은 다른 답변 (Wikipedia 등)의 외부 참조로 이동하십시오.
아마도 당신이 요청한 것은 아니지만 ... numpy 배열을 사용하면 효율적으로 작업을 수행합니다.
from numpy import array
nums = array(((0.01, 0.01, 0.02, 0.04, 0.03),
(0.00, 0.02, 0.02, 0.03, 0.02),
(0.01, 0.02, 0.02, 0.03, 0.02),
(0.01, 0.00, 0.01, 0.05, 0.03)))
print nums.std(axis=1)
# [ 0.0116619 0.00979796 0.00632456 0.01788854]
print nums.mean(axis=1)
# [ 0.022 0.018 0.02 0.02 ]
그건 그렇고,이 블로그 게시물에는 평균 및 분산 계산을위한 원 패스 방법에 대한 몇 가지 흥미로운 토론이 있습니다.
다음은 http://www.johndcook.com/standard_deviation.html 에서 Welford 알고리즘 구현의 문자 그대로 순수한 Python 번역입니다 .
https://github.com/liyanage/python-modules/blob/master/running_stats.py
class RunningStats:
def __init__(self):
self.n = 0
self.old_m = 0
self.new_m = 0
self.old_s = 0
self.new_s = 0
def clear(self):
self.n = 0
def push(self, x):
self.n += 1
if self.n == 1:
self.old_m = self.new_m = x
self.old_s = 0
else:
self.new_m = self.old_m + (x - self.old_m) / self.n
self.new_s = self.old_s + (x - self.old_m) * (x - self.new_m)
self.old_m = self.new_m
self.old_s = self.new_s
def mean(self):
return self.new_m if self.n else 0.0
def variance(self):
return self.new_s / (self.n - 1) if self.n > 1 else 0.0
def standard_deviation(self):
return math.sqrt(self.variance())
용법:
rs = RunningStats()
rs.push(17.0);
rs.push(19.0);
rs.push(24.0);
mean = rs.mean();
variance = rs.variance();
stdev = rs.standard_deviation();
파이썬 RUNSTATS 모듈 것은 단지 이런 종류입니다. PyPI에서 runstats 를 설치합니다 .
pip install runstats
Runstats 요약은 단일 데이터 패스에서 평균, 분산, 표준 편차, 왜도 및 첨도를 생성 할 수 있습니다. 이것을 사용하여 "실행중인"버전을 만들 수 있습니다.
from runstats import Statistics
stats = [Statistics() for num in range(len(data[0]))]
for row in data:
for index, val in enumerate(row):
stats[index].push(val)
for index, stat in enumerate(stats):
print 'Index', index, 'mean:', stat.mean()
print 'Index', index, 'standard deviation:', stat.stddev()
통계 요약은 Art of Computer Programming, Vol 2, p에 설명 된대로 표준 편차를 한 번에 계산하는 Knuth 및 Welford 방법을 기반으로합니다. 232, 3 판. 이것의 이점은 수치 적으로 안정적이고 정확한 결과입니다.
면책 조항 : 저는 Python runstats 모듈의 작성자입니다.
Statistics :: Descriptive 는 다음과 같은 유형의 계산을위한 매우 괜찮은 Perl 모듈입니다.
#!/usr/bin/perl
use strict; use warnings;
use Statistics::Descriptive qw( :all );
my $data = [
[ 0.01, 0.01, 0.02, 0.04, 0.03 ],
[ 0.00, 0.02, 0.02, 0.03, 0.02 ],
[ 0.01, 0.02, 0.02, 0.03, 0.02 ],
[ 0.01, 0.00, 0.01, 0.05, 0.03 ],
];
my $stat = Statistics::Descriptive::Full->new;
# You also have the option of using sparse data structures
for my $ref ( @$data ) {
$stat->add_data( @$ref );
printf "Running mean: %f\n", $stat->mean;
printf "Running stdev: %f\n", $stat->standard_deviation;
}
__END__
산출:
C:\Temp> g
Running mean: 0.022000
Running stdev: 0.013038
Running mean: 0.020000
Running stdev: 0.011547
Running mean: 0.020000
Running stdev: 0.010000
Running mean: 0.020000
Running stdev: 0.012566
PDL ( "piddle!"로 발음)을 살펴보십시오 .
이것은 고정밀 수학 및 과학 컴퓨팅을 위해 설계된 Perl 데이터 언어입니다.
다음은 수치를 사용한 예입니다 ....
use strict;
use warnings;
use PDL;
my $figs = pdl [
[0.01, 0.01, 0.02, 0.04, 0.03],
[0.00, 0.02, 0.02, 0.03, 0.02],
[0.01, 0.02, 0.02, 0.03, 0.02],
[0.01, 0.00, 0.01, 0.05, 0.03],
];
my ( $mean, $prms, $median, $min, $max, $adev, $rms ) = statsover( $figs );
say "Mean scores: ", $mean;
say "Std dev? (adev): ", $adev;
say "Std dev? (prms): ", $prms;
say "Std dev? (rms): ", $rms;
다음을 생성합니다.
Mean scores: [0.022 0.018 0.02 0.02]
Std dev? (adev): [0.0104 0.0072 0.004 0.016]
Std dev? (prms): [0.013038405 0.010954451 0.0070710678 0.02]
Std dev? (rms): [0.011661904 0.009797959 0.0063245553 0.017888544]
statsover 함수 에 대한 자세한 내용 은 PDL :: Primitive 를 참조하십시오 . 이것은 ADEV가 "표준 편차"임을 시사하는 것 같습니다.
그러나 PRMS (Sinan의 Statistics :: Descriptive 예제가 표시됨) 또는 RMS (ars의 NumPy 예제가 표시됨) 일 수 있습니다. 이 세 가지 중 하나가 맞을 것 같아요 ;-)
자세한 PDL 정보는 다음을 참조하십시오.
- pdl.perl.org (공식 PDL 페이지).
- PerlMonks에 대한 PDL 빠른 참조 가이드
- PDL에 대한 Dobb 박사의 기사
- PDL 위키
- PDL에 대한 Wikipedia 항목
- PDL 용 Sourceforge 프로젝트 페이지
어레이가 얼마나 큽니까? 길이가 무수히 많은 요소가 아니라면 두 번 반복하는 것에 대해 걱정하지 마십시오. 코드는 간단하고 쉽게 테스트됩니다.
내 선호는 numpy 배열 수학 확장을 사용하여 배열 배열을 numpy 2D 배열로 변환하고 표준 편차를 직접 얻는 것입니다.
>>> x = [ [ 1, 2, 4, 3, 4, 5 ], [ 3, 4, 5, 6, 7, 8 ] ] * 10
>>> import numpy
>>> a = numpy.array(x)
>>> a.std(axis=0)
array([ 1. , 1. , 0.5, 1.5, 1.5, 1.5])
>>> a.mean(axis=0)
array([ 2. , 3. , 4.5, 4.5, 5.5, 6.5])
이것이 옵션이 아니고 순수한 Python 솔루션이 필요한 경우 계속 읽으십시오 ...
어레이가
x = [
[ 1, 2, 4, 3, 4, 5 ],
[ 3, 4, 5, 6, 7, 8 ],
....
]
표준 편차는 다음과 같습니다.
d = len(x[0])
n = len(x)
sum_x = [ sum(v[i] for v in x) for i in range(d) ]
sum_x2 = [ sum(v[i]**2 for v in x) for i in range(d) ]
std_dev = [ sqrt((sx2 - sx**2)/N) for sx, sx2 in zip(sum_x, sum_x2) ]
배열을 한 번만 반복하기로 결정한 경우 누적 합계를 결합 할 수 있습니다.
sum_x = [ 0 ] * d
sum_x2 = [ 0 ] * d
for v in x:
for i, t in enumerate(v):
sum_x[i] += t
sum_x2[i] += t**2
위의 목록 이해 솔루션만큼 우아하지는 않습니다.
이 문제가 도움이 될 것 같습니다. 표준 편차
표준 편차 에 대한 Wikipedia 기사 , 특히 빠른 계산 방법에 대한 섹션을 볼 수 있습니다.
Python을 사용하는 기사도 있습니다. 많은 변경없이 코드를 사용할 수 있어야합니다. Subliminal Messages-Running Standard Deviations .
n=int(raw_input("Enter no. of terms:"))
L=[]
for i in range (1,n+1):
x=float(raw_input("Enter term:"))
L.append(x)
sum=0
for i in range(n):
sum=sum+L[i]
avg=sum/n
sumdev=0
for j in range(n):
sumdev=sumdev+(L[j]-avg)**2
dev=(sumdev/n)**0.5
print "Standard deviation is", dev
다음 답변에 설명되어 있습니다. pandas / scipy / numpy는 누적 표준 편차 함수를 제공합니까? Python Pandas 모듈에는 실행 또는 누적 표준 편차 를 계산하는 방법이 포함되어 있습니다 . 이를 위해 데이터를 pandas 데이터 프레임 (또는 1D 인 경우 시리즈)으로 변환해야하지만이를위한 함수가 있습니다.
다음은 함수형 프로그래밍 스타일로 여러 줄에 걸쳐있는 "한 줄짜리"입니다.
def variance(data, opt=0):
return (lambda (m2, i, _): m2 / (opt + i - 1))(
reduce(
lambda (m2, i, avg), x:
(
m2 + (x - avg) ** 2 * i / (i + 1),
i + 1,
avg + (x - avg) / (i + 1)
),
data,
(0, 0, 0)))
나는 다음과 같이 업데이트를 표현하고 싶습니다.
def running_update(x, N, mu, var):
'''
@arg x: the current data sample
@arg N : the number of previous samples
@arg mu: the mean of the previous samples
@arg var : the variance over the previous samples
@retval (N+1, mu', var') -- updated mean, variance and count
'''
N = N + 1
rho = 1.0/N
d = x - mu
mu += rho*d
var += rho*((1-rho)*d**2 - var)
return (N, mu, var)
1 회 통과 함수는 다음과 같습니다.
def one_pass(data):
N = 0
mu = 0.0
var = 0.0
for x in data:
N = N + 1
rho = 1.0/N
d = x - mu
mu += rho*d
var += rho*((1-rho)*d**2 - var)
# could yield here if you want partial results
return (N, mu, var)
note that this is calculating the sample variance (1/N), not the unbiased estimate of the population variance (which uses a 1/(N-1) normalzation factor). Unlike the other answers, the variable, var, that is tracking the running variance does not grow in proportion to the number of samples. At all times it is just the variance of the set of samples seen so far (there is no final "dividing by n" in getting the variance).
In a class it would look like this:
class RunningMeanVar(object):
def __init__(self):
self.N = 0
self.mu = 0.0
self.var = 0.0
def push(self, x):
self.N = self.N + 1
rho = 1.0/N
d = x-self.mu
self.mu += rho*d
self.var += + rho*((1-rho)*d**2-self.var)
# reset, accessors etc. can be setup as you see fit
This also works for weighted samples:
def running_update(w, x, N, mu, var):
'''
@arg w: the weight of the current sample
@arg x: the current data sample
@arg mu: the mean of the previous N sample
@arg var : the variance over the previous N samples
@arg N : the number of previous samples
@retval (N+w, mu', var') -- updated mean, variance and count
'''
N = N + w
rho = w/N
d = x - mu
mu += rho*d
var += rho*((1-rho)*d**2 - var)
return (N, mu, var)
'programing tip' 카테고리의 다른 글
| 지수 표기법에서 숫자 구문 분석 (0) | 2020.10.11 |
|---|---|
| Python 요청 라이브러리를 사용하여 게시 요청에서 쿠키를 보내는 방법은 무엇입니까? (0) | 2020.10.11 |
| AVAssetReader로 파형 그리기 (0) | 2020.10.11 |
| 회사 프록시에 의해 차단 된 Bower 통화 (0) | 2020.10.11 |
| "이 창을 닫으시겠습니까?"프롬프트를받지 않고 브라우저 창을 닫으려면 어떻게해야합니까? (0) | 2020.10.11 |