--- title: "Generic Simplified Log Likelihood" output: rmarkdown::html_vignette: fig_caption: yes keep_md: true vignette: > %\VignetteIndexEntry{Generic Simplified Log Likelihood} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} bibliography: references.bib --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.dim = c(6, 4) ) library(knitr) local({ formatEng <- function( x, digits=4) { # from https://stat.ethz.ch/pipermail/r-help/2006-July/108808.html # Hans-Joerg Bibiko # Max Planck Institute for Evolutionary Anthropology # Department of Linguistics # Deutscher Platz 6 phone: +49 (0) 341 3550 341 # D-04103 Leipzig fax: +49 (0) 341 3550 333 # Germany e-mail: bibiko@eva.mpg.de # modified by Phil Shea phil1shea@duck.go s <- as.numeric( strsplit( format( x, scientific=TRUE, digits=digits), "e")[[ 1]]) mant <- s[1] * 10^(s[2] %% 3) expo <- as.integer( s[2] - (s[2] %% 3)) return( ifelse( expo==0, paste( mant), paste( mant, expo, sep="e"))) } inline_hook <- function (x) { if (is.numeric(x)) { # ifelse does a vectorized comparison # If integer, print with commas; otherwise print two places res <- ifelse( is.integer( x) || ((x == round(x)) && (x < .Machine$integer.max)), noquote( prettyNum( x, big.mark=",")), formatEng( x)) paste(res, collapse = ", ") } else paste( x) } knit_hooks$set(inline = inline_hook) }) ``` ```{r } 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: ```{r data} QPSKdB.B2B <- B2BConvert( QPSKdB) s <- 0:20 O2 <- 3 B2 <- Inf N <- 1000000 (r2 <- rbinom( length( s), N, QPSKdB.B2B( s, B2, O2))) mle6 <- mleB2B( Errors=r2, N=N, f=QPSKdB.B2B, fparms=list( x=s), start=c(b2b=20, offset=0)) summary( mle6) ``` This estimated a B2B of `r coef( mle6)["b2b"]` 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 $\mathrm{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 $\gamma$ is the linear SNR. We also note that R does not have the $\mathrm{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: $$ \begin{aligned} Off_{dB} &= -dB( a) \\ B2B_{dB} &= -dB( b) \\ s_{dB} &= dB( \gamma) \end{aligned} $$ 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`). ```{r} 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) ``` This 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 $\alpha$ and $\beta$, but we will want $a = \alpha^2$ and $b = \beta^2$. ```{r Q_ab} 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) ``` Note that the small `beta` value of `r coef(mle8)["beta"]` is equivalent to a `r -2*dB( abs(coef(mle8)["beta"]))` dB B2BQ. Below is a plot of that result. ```{r plotQab} mle8coef <- coef( mle8) plot( s, y=r2/N, log='y', type='p', panel.first = grid()) 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 $b^2$ is greater than $b^2$, so it is reasonable to let $b = 0$. Now let's try it with a detectable B2B. ```{r} O4 <- 3 B4 <- 15 (r4 <- rbinom( length( s), N, QPSKdB.B2B( s, B4, O4))) mle9 <- mleB2B( N=N, Errors=r4, f=Q_ab, start=c( alpha=1, beta=0), fparms=list(gamma=gamma)) summary( mle9) mle9coef <- coef( mle9) (mle9sd <- sqrt( diag( vcov( mle9)))) plot( s, r4/N, log='y',panel.first = grid()) 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: ```{r} -2*dB(mle9coef) ``` [^1]: see Abramowitz and Stegun 29.2.29