This repository has been archived by the owner on Oct 26, 2022. It is now read-only.
/
snr.py
209 lines (160 loc) · 9.08 KB
/
snr.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
#! /usr/bin/env python2
# encoding:UTF-8
# Copyright (c) 2012 Victor Terron. All rights reserved.
# Institute of Astrophysics of Andalusia, IAA-CSIC
#
# This file is part of LEMON.
#
# LEMON is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
from __future__ import division
import math
import operator
import numpy
def snr_to_error(snr):
"""Signal-to-noise ratio to error in magnitudes conversion.
The method computes the error in magnitudes to which a determined SNR is
equivalent. This is achieved by comparing the mean number of counts, c, to
the maximum and minimum values, both of which are returned, induced by the
noise, i.e, Δm = -2.5 log10((c +- c/SNR)/c) = -2.5 log10(1 ± 1/SNR). The
two values derived from this formula are returned as a two-element tuple:
the first one will correspond to the solution which the addition and the
second one to the subtraction of the 'c/SNR' part of the equation.
Note that this means that a S/N of 100 implies an observational error or
0.01. This formula is taken from 'Astronomical Photometry', by Arne A.
Henden and Ronald H. Kaitchuc [1982, page 79]
The domain of the formula is the set of all non-negative numbers, as the
base-10 logarithm of zero or a negative value cannot be calculated. Thus,
ValueError is raised if the signal-to-noise ratio is not above one.
"""
if not numpy.any(snr > 1):
raise ValueError("SNR cannot be less than or equal to one")
operators = (operator.add, operator.sub)
return tuple([-2.5 * numpy.log10(f(1, 1 / snr)) for f in operators])
def error_to_snr(error):
"""Error in magnitudes to signal-to-noise ratio conversion.
If we find the value of Δm for a given SNR, we can also obtain the
signal-to-noise ration to which a determined error in magnitudes is
equivalent. To do this we need to isolate in Henden's equation the
value of SNR:
1. Δm = -2.5 log10(1 ± 1/SNR)
2. Δm / -2.5 = log10(1 ± 1/SNR)
3. 10 ** (Δm / -2.5) = 10 ** (log10(1 ± 1/SNR))
4. 10 ** (Δm / -2.5) = 1 ± 1/SNR
5. 10 ** (Δm / -2.5) -1 = ± 1/SNR
6. ± 1 / (10 ** (Δm / -2.5) -1) = SNR
Note that the equation has a plus-minus sign, while the function only
receives one error in magnitudes. How to handle this? Easily. Examining
Henden's equation we can notice that when plus sign is used the noise is
always negative, while a minus makes the returned value to be positive.
Therefore, which sign the equation must use can be straight-forwardly
determined: a plus if the error is negative and a minus if it turns out
to be negative.
"""
return (1 if error < 0 else -1) / (math.pow(10, error / -2.5) - 1)
def difference_error(*errors):
"""Return the absolute error of the difference of a series of errors.
The error of a combination (addition / subtraction) of *independent* errors
is given by their addition in quadrature. A pessimistic approach could
consist in simply adding them, obtaining the maximum and minimum probable
errors, but since they are independent quantities the first error is just
as likely to (partially) cancel the second error as it is to increase it.
The correct way of combining independent errors, therefore, is to add them
in quadrature; in other words, the absolute error in a combination of sums
and differences is obtained by taking the root of the sum of the squares of
the absolute errors in the original quantities. For example, if we had two
measures whose errors were 0.10 and 0.05, then the error of the addition or
difference of these measures would be be sqrt(0.1 ** 2 + 0.05 ** 2) = 0.11.
Note: the name of this error is so because, in our code, it is expected to
be used only when calculating the error of the difference in magnitude
between two or more stars. However, it may be perfectly used for additions
too, and of course also for a combination of additions and subtractions.
"""
return math.sqrt(sum(e ** 2 for e in errors))
def difference_snr(*snrs):
"""Return the SNR of the difference of a series of SNRs.
The method returns the signal-to-noise ratio of the difference of a series
of signal-to-noise ratios. This is internally done by converting the SNRs
to errors in magnitudes, computing the absolute error and converting the
resulting value back to its equivalent signal-to-noise ratio.
As it is the case with difference_error, this method, despite its name, may
also be used to compute the resulting signal-to-noise ratio of the addition
of different SNRs, as well as the combination of additions and subtractions.
"""
# We take the second of the two errors (i.e., the positive one) returned by
# the conversion from signal-to-noise ratio. Using the negative value (the
# minimum noise) would not work, as the addition is quadrature would make
# the errors become positive, and as such they would be considered when
# converting back to SNR.
errors = [snr_to_error(s)[1] for s in snrs]
assert all(e >= 0 for e in errors)
error = difference_error(*errors)
return error_to_snr(error)
def mean_error(errors, weights=None):
"""Return the absolute error of the arithmetic mean of errors.
If we have a series of independent values (for example, instrumental
magnitudes) and their corresponding uncertainty (the error in magnitudes),
what is the uncertainty (error) if we calculate their arithmetic mean? By
default all the values are assumed to have contributed equally to the mean,
but the method can also calculate the error resulting from a weighted one.
The arithmetic mean is just a scaled version of the sum, so the error
scales as the quantity itself under scaling. Therefore, the error in the
arithmetic mean of a series on n independent values is given by the
equation \sqtr{e_{1}^{2} + e_{2}^{2} + ... + e_{n}^{2}} / n
As for the weighted mean of n independent values with their corresponding
weights, c_{1}, c_{2}... c{n}, and where \sum{i=0}^n{c_{i}} = 1, the error
is \sqrt{c_{1}^2 e_{1}^2 + c_{2}^2 e_{2}^2 + ... + c_{n}^2 e_{n}^2). This
is the equation implemented by the method, indeed, where the coefficients
default to 1/n if no weights are given.
Thanks so much to the people at Math Stack Exchange for their help:
http://math.stackexchange.com/q/123276/
Keyword arguments:
weights - the coefficients of the weighted mean. The i-th weight is
interpreted to correspond to the i-th error received by the
method. The coefficients do not have to sum up to one; they will
internally normalized to guarantee this. Therefore, the weights
[0.5, 0.5], [1.0, 1.0] and [2.6, 2.6], e.g., are equivalent.
"""
if weights is None:
# All the values contribute equally (and weights sum up to one)
weights = [1 / len(errors)] * len(errors)
elif len(weights) != len(errors):
raise ValueError("number of weights must equal that of errors")
else:
# Normalize the values so that they sum up to one
weights = [w / math.fsum(weights) for w in weights]
return math.sqrt(math.fsum(((w ** 2) * (e ** 2) for e, w in zip(errors, weights))))
def mean_snr(snrs, weights=None):
"""Return the SNR of the arithmetic mean of a series of SNSRs.
The method returns the signal-to-noise ratio of the arithmetic or weighted
mean of a series of signal-to-noise ratios. This is internally done by
converting the SNRs to errors in magnitudes, computing the absolute error
and converting the resulting value back to its equivalent SNR.
Keyword arguments:
weights - the coefficients of the weighted mean. The i-th weight is
interpreted to correspond to the i-th error received by the
method. The coefficients do not have to sum up to one; they will
internally normalized to guarantee this. Therefore, the weights
[0.5, 0.5], [1.0, 1.0] and [2.6, 2.6], e.g., are equivalent.
"""
# We take the second of the two errors (i.e., the positive one) returned by
# the conversion from signal-to-noise ratio. Using the negative value (the
# minimum noise) would not work, as the addition is quadrature would make
# the errors become positive, and as such they would be considered when
# converting back to SNR.
errors = [snr_to_error(s)[1] for s in snrs]
assert all(e >= 0 for e in errors)
error = mean_error(errors, weights=weights)
return error_to_snr(error)
if __name__ == "__main__":
pass