Generic Simplified Log Likelihood
library( binfunest)
library( stats4) # need this in the search path to work with mle results.
set.seed( 31394)If you don’t have substantial samples of the high-SNR tail of a bit error rate function, the offset and back-to-back (B2B) can get confused. Also, the form requires an infinite estimate for B2B if the system is performing more-or-less theoretically. Using functions in Decibels yields the following:
QPSKdB.B2B <- B2BConvert( QPSKdB)
s <- 0:20
O2 <- 3
B2 <- Inf
N <- 1000000
(r2 <- rbinom( length( s), N, QPSKdB.B2B( s, B2, O2)))
#> [1] 158193 130640 103416 78594 56209 37571 22730 12447 6019 2451
#> [11] 700 211 23 3 0 0 0 0 0 0
#> [21] 0
mle6 <- mleB2B( Errors=r2, N=N, f=QPSKdB.B2B, fparms=list( x=s),
start=c(b2b=20, offset=0))
summary( mle6)
#> Maximum likelihood estimation
#>
#> Call:
#> stats4::mle(minuslogl = function (b2b = 20, offset = 0)
#> -sum(sort(stats::dbinom(Errors, N, (function (x, B2B = Inf, offset = 0)
#> QPSKdB(-dB((undB(-x) + undB(-B2B))) - offset))(x = 0:20, b2b,
#> offset), log = TRUE))), start = c(b2b = 20, offset = 0),
#> method = "Nelder-Mead", nobs = 21L)
#>
#> Coefficients:
#> Estimate Std. Error
#> b2b 46.114366 35.950957636
#> offset 2.996244 0.004852486
#>
#> -2 log L: 163.4187This estimated a B2B of 46.11 dB.We need a method that is numerically better behaved in more conditions. We can cast the BER in terms of the Q function: $Q( s) = (1/2)\mathrm{Erfc}(s/\sqrt{2})$ where Erfc(x) is the complementary error function 1. Using this we can create a generic BER function as:
$$ BER( \gamma, a, b) = Q\left( \sqrt{ \frac{a \gamma} {1 + b \gamma}} \right) $$
where a and b are the linear versions of the
offset and B2B and γ is the
linear SNR. We also note that R does not have the Erfc, however, the following is true:
Erfc( x) = 2 * pnorm( x * sqrt( 2), lower=FALSE), and the
package includes the function Q_(x). With this formulation,
we have the following conversions:
$$
We can use this directly, but we have to limit b to be greater than zero or
sqrt might be passed a negative number. See documentation
on optim for how this is done (. mleB2B
parameters in the ... are passed to mle, which
are passed to optim).
gamma <- undB( s)
Q_ab <- function( gamma, a, b)
Q_( sqrt( a * gamma / (1 + b * gamma)))
mle7 <- mleB2B( N=N, Errors=r2, f=Q_ab, start=c( a=1, b=0),
fparms=list(gamma=gamma), method="L-BFGS-B", lower=list(b=0))
summary( mle7)
#> Maximum likelihood estimation
#>
#> Call:
#> stats4::mle(minuslogl = function (a = 1, b = 0)
#> -sum(sort(stats::dbinom(Errors, N, (function (gamma, a, b)
#> Q_(sqrt(a * gamma/(1 + b * gamma))))(gamma = c(1, 1.25892541179417,
#> 1.58489319246111, 1.99526231496888, 2.51188643150958, 3.16227766016838,
#> 3.98107170553497, 5.01187233627272, 6.30957344480193, 7.94328234724282,
#> 10, 12.5892541179417, 15.8489319246111, 19.9526231496888, 25.1188643150958,
#> 31.6227766016838, 39.8107170553497, 50.1187233627273, 63.0957344480193,
#> 79.4328234724282, 100), a, b), log = TRUE))), start = c(a = 1,
#> b = 0), method = "L-BFGS-B", nobs = 21L, lower = list(b = 0))
#>
#> Coefficients:
#> Estimate Std. Error
#> a 1.00322 0.0015879675
#> b 0.00000 0.0004025817
#>
#> -2 log L: 163.4075This doesn’t always work though, so sometimes you might want to allow a negative b. In order to allow the parameters to be negative without getting a negative square root, square the parameters in the objective function. Since our data was created with b = 0 (i.e., the B2B Q is infinite) we want the gradient search to allow negative a and b without failing.
$$ BER_2( \gamma, \alpha, \beta) = Q\left( \sqrt{ \frac{\alpha^2 \gamma} {1 + \beta^2 \gamma}} \right) $$
We generate some observations of QPSK data with a 3 dB offset and an
infinite B2B Q. Since QPSK follows the Q function plus 3 dB, the fit should
find a equal to one. The mleB2B function will
estimate α and β, but we will want a = α2 and b = β2.
gamma <- undB( s)
Q_ab2 <- function( gamma, alpha, beta)
Q_( sqrt( alpha^2 * gamma / (1 + beta^2 * gamma)))
mle8 <- mleB2B( N=N, Errors=r2, f=Q_ab2, start=c( alpha=1, beta=0),
fparms=list(gamma=gamma))
summary( mle8)
#> Maximum likelihood estimation
#>
#> Call:
#> stats4::mle(minuslogl = function (alpha = 1, beta = 0)
#> -sum(sort(stats::dbinom(Errors, N, (function (gamma, alpha, beta)
#> Q_(sqrt(alpha^2 * gamma/(1 + beta^2 * gamma))))(gamma = c(1,
#> 1.25892541179417, 1.58489319246111, 1.99526231496888, 2.51188643150958,
#> 3.16227766016838, 3.98107170553497, 5.01187233627272, 6.30957344480193,
#> 7.94328234724282, 10, 12.5892541179417, 15.8489319246111, 19.9526231496888,
#> 25.1188643150958, 31.6227766016838, 39.8107170553497, 50.1187233627273,
#> 63.0957344480193, 79.4328234724282, 100), alpha, beta), log = TRUE))),
#> start = c(alpha = 1, beta = 0), method = "Nelder-Mead", nobs = 21L)
#>
#> Coefficients:
#> Estimate Std. Error
#> alpha 1.001609e+00 0.0004542291
#> beta -2.414827e-05 0.0675007113
#>
#> -2 log L: 163.4075Note that the small beta value of -24.15e-6 is
equivalent to a 92.34 dB B2BQ. Below is a plot of that result.
mle8coef <- coef( mle8)
plot( s, y=r2/N, log='y', type='p', panel.first = grid())
#> Warning in xy.coords(x, y, xlabel, ylabel, log): 7 y values <= 0 omitted from
#> logarithmic plot
lines( s, QPSKdB( s))
lines( s, Q_(sqrt( gamma)), col='red')
lines( s, y=Q_ab( gamma, mle8coef[1], mle8coef[2]), col="green")
legend( "bottomleft",
legend=c( "Data", "QPSK Theory", "Q w/ 0 dB Offset",
"Estimated"),
lty=c( NA, 1, 1, 5), col=c( 'black', 'black', 'red', 'green'),
pch=c( 1, NA, NA, NA))
That generated very nice estimates. Note that the standard deviation of b2 is greater than b2, so it is reasonable to let b = 0.
Now let’s try it with a detectable B2B.
O4 <- 3
B4 <- 15
(r4 <- rbinom( length( s), N, QPSKdB.B2B( s, B4, O4)))
#> [1] 161655 134889 109241 85045 62946 44893 29886 18680 11005 5826
#> [11] 2874 1340 597 227 100 30 14 2 1 0
#> [21] 0
mle9 <- mleB2B( N=N, Errors=r4, f=Q_ab, start=c( alpha=1, beta=0),
fparms=list(gamma=gamma))
#> Warning in sqrt(a * gamma/(1 + b * gamma)): NaNs produced
summary( mle9)
#> Maximum likelihood estimation
#>
#> Call:
#> stats4::mle(minuslogl = function (alpha = 1, beta = 0)
#> -sum(sort(stats::dbinom(Errors, N, (function (gamma, a, b)
#> Q_(sqrt(a * gamma/(1 + b * gamma))))(gamma = c(1, 1.25892541179417,
#> 1.58489319246111, 1.99526231496888, 2.51188643150958, 3.16227766016838,
#> 3.98107170553497, 5.01187233627272, 6.30957344480193, 7.94328234724282,
#> 10, 12.5892541179417, 15.8489319246111, 19.9526231496888, 25.1188643150958,
#> 31.6227766016838, 39.8107170553497, 50.1187233627273, 63.0957344480193,
#> 79.4328234724282, 100), alpha, beta), log = TRUE))), start = c(alpha = 1,
#> beta = 0), method = "Nelder-Mead", nobs = 21L)
#>
#> Coefficients:
#> Estimate Std. Error
#> alpha 1.00463983 0.0014578468
#> beta 0.03211088 0.0003437894
#>
#> -2 log L: 195.4187
mle9coef <- coef( mle9)
(mle9sd <- sqrt( diag( vcov( mle9))))
#> alpha beta
#> 0.0014578468 0.0003437894
plot( s, r4/N, log='y',panel.first = grid())
#> Warning in xy.coords(x, y, xlabel, ylabel, log): 2 y values <= 0 omitted from
#> logarithmic plot
lines( s, QPSKdB( s))
lines( s, y=QPSKdB.B2B( s, Inf, O4), col="blue")
lines( s, y=QPSKdB.B2B( s, B4, O4), col="red")
lines( s, y=Q_ab( gamma, mle9coef[1], mle9coef[2]), col="green", lty=5)
legend( "bottomleft",
legend=c( "Data", "Theory", "Theory + 3 dB", "3 dB Offset + 15 dB B2B",
"Estimated"),
lty=c( NA, 1, 1, 1, 5),
col=c( 'black', 'black', 'blue', 'red', 'green'),
pch=c( 1, NA, NA, NA, NA))
We can get the coefficients in Decibels from the following expression:
see Abramowitz and Stegun 29.2.29↩︎