TITLE: Bernoulli Uniform Bayes Estimator
QUESTION [1 upvotes]: My answer comes out as $(p|X)$~$BETA(x+1,-x+2)$, indicating that $p_{Bayes}=\frac{x+1}{3}$, but apparently the correct answer is $p_{Bayes}=\frac{\sum x_i+1}{n+2}$.
I don't understand where this answer comes from, however; can somebody here explain it?  Thanks.

REPLY [2 votes]: You are calculating the posterior distribution incorrectly. I'll use $\theta$ instead of $p$ to avoid confusing the notation, so that $\theta \sim U(0,1)$ and $p(\theta) = 1$:
\begin{align*}
p(\theta | X) &= p(X | \theta) p(\theta)\\
&= \prod_{i=1}^n p(X_i|\theta) p(\theta)\\
&= p(\theta) \prod_{i=1}^n \theta^{X_i}  (1- \theta)^{1-X_i}\\
&= 1 \times \theta^{n \bar{X}}(1- \theta)^{n-n\bar{X}}\\
&= \theta^{n \bar{X}}(1- \theta)^{n-n\bar{X}}
\end{align*}
which is the kernel of a beta distribution with parameters:
$$
B(\alpha = n \bar{X} + 1, \beta = n-n \bar{X} +1 )
$$
where $n \bar{X} = \sum_{i} X_i$.
Then, the bayes estimator of $\theta$ with squared error is simply the posterior mean (mean of a beta distribution with the above parameters):
\begin{align*}
\theta_{\text{bayes}} &= E(\theta|X) \\
&= \frac{\alpha}{ \alpha + \beta}\\
&= \frac{ n \bar{X} + 1}{ n \bar{X} + 1 + n-n \bar{X} +1}\\
&= \frac{n \bar{X} +1}{n+2}\\
& = \frac{\sum_i X_i +1}{n+2}\\
\end{align*}
To solve for the Bayes estimator of $\theta(1-\theta)$, apply what we just did. The Bayes estimator under SE loss is simply the posterior expectation of the thing we are trying to estimate.
$$
E(\theta(1-\theta)) = E(\theta) - E(\theta^2) = E(\theta) - (V(\theta) + E(\theta)^2)
$$
We already know what $E(\theta)$ is from above, and $V(\theta)$  is just:
$$
\frac{\alpha \beta}{ (\alpha + \beta)^2 ( \alpha + \beta + 1)}
$$


text:
\begin{document}
\maketitle

\section{Generalized Linear Model and EM}

Consider the generalized linear model (GLM) defined by
\begin{align}
\zbf &= \Abf\xbf,
\quad
\text{prior~} p(\xbf;\thetabf_x) = \prod_{j=1}^N p(x_j;\thetabf_x),
\quad
\text{and $\zbf$-likelihood~} p(\ybf|\zbf;\thetabf_z) = \prod_{i=1}^M p(y_i|z_i;\thetabf_z) ,
\label{eq:glm}
\end{align}
where
$\thetabf \defn [\thetabf_x,\thetabf_z]$ are unknown deterministic parameters,
and where $\zbf\in\R^M$ and $\xbf\in\R^N$.

The maximum-likelihood estimate of $\thetabf$ given $\ybf$ is
\begin{align}
\thetabfhat
&= \arg\min_\thetabf \big\{ -\ln p(\ybf;\thetabf) \big\} ,
\end{align}
where
\begin{align}
p(\ybf;\thetabf)
&= \int p(\ybf,\zbf,\xbf;\thetabf) \dif\zbf \dif\xbf
= \int p(\ybf|\zbf,\xbf;\thetabf) p(\zbf,\xbf;\thetabf)\dif\zbf \dif\xbf
= \int p(\ybf|\zbf;\thetabf_z) \delta(\zbf-\Abf\xbf) p(\xbf;\thetabf_x) \dif\zbf \dif\xbf
\label{eq:logy} .
\end{align}
Although the ``$\thetabf$-likelihood function'' $p(\ybf;\thetabf)$ is difficult to work with directly, the expectation-maximization (EM) algorithm \cite{NealHinton:98} offers an alternative.
There, the idea is to write
\begin{align}
-\ln p(\ybf;\thetabf)
&= J(b;\thetabf) - D(b \,\|\, p(\xbf|\ybf;\thetabf))
\end{align}
for an arbitrary belief $b(\xbf)$, where $D(\cdot\|\cdot)$ is KL divergence,
\begin{align}
J(b;\thetabf)
&\defn D(b \,\|\, p(\xbf;\thetabf_x)) + D(b \,\|\, p(\ybf|\xbf;\thetabf_z)) + H(b)
\end{align}
is known as the Gibbs free energy,
and $H(b)$ is the entropy of $b$.
Because $D(b \,\|\, p(\xbf|\ybf;\thetabf))\geq 0$ for any $b$, we have that $J(b;\thetabf)$ is an upper bound on $-\ln p(\ybf;\thetabf)$, the quantity that ML seeks to minimize.
Thus, if it is tractable to construct and minimize $J(b;\thetabf)$, it makes sense to iterate the following two steps (over $k=1,2,\dots$):
\begin{align}
\text{E step:~}& b^k(\xbf)
=p(\xbf|\ybf;\thetabfhat^k)
\label{eq:Estep}\\
\text{M step:~}& \thetabfhat^{\kp1}
=\arg\min_\thetabf J(b^k;\thetabf)
=\arg\min_\thetabf D(b^k \,\|\, p(\xbf;\thetabf_x)) + D(b^k \,\|\, p(\ybf|\xbf;\thetabf_z))
\label{eq:Mstep} ,
\end{align}
which together constitute the EM algorithm.
The ``E'' step creates an upper bound on $-\ln p(\ybf;\thetabf)$ that is tight at $\thetabf=\thetabfhat^k$, and the ``M'' step finds the estimate of $\thetabf$ that minimizes this bound.

Unfortunately, however, the posterior density required by the E-step \eqref{eq:Estep},
\begin{align}
p(\xbf|\ybf;\thetabf)
&= \frac{p(\xbf;\thetabf_x) p(\ybf|\xbf;\thetabf_z)}{p(\ybf;\thetabf)}
= \frac{p(\xbf;\thetabf_x) \int p(\ybf|\zbf;\thetabf_z) \delta(\zbf-\Abf\xbf) \dif\zbf}{\int p(\xbf;\thetabf_x) p(\ybf|\zbf;\thetabf_z) \delta(\zbf-\Abf\xbf) \dif\zbf\dif\xbf} ,
\end{align}
is difficult to compute due to the high-dimensional integration.
Thus we consider an approximation afforded by the GVAMP algorithm \cite{schniter2016vectorAsil}.
For this, we first reparameterize the GLM \eqref{eq:glm} as an standard linear model (SLM).

\section{An SLM Equivalent}

The GLM \eqref{eq:glm} can be written as an SLM using the following formulation:
\begin{align}
\ybfbar
&= \Abfbar\xbfbar + \wbfbar
~~~\text{with}~~~
\ybfbar\defn \zero,~
\Abfbar\defn \mat{\Abf & -\Ibf_M},~
\xbfbar\defn \mat{\xbf\\\zbf},~
\wbfbar\sim \mathcal{N}(\zero,\epsilon \Ibf_M)
~\text{s.t.}~
\epsilon\rightarrow 0
\label{eq:slm} .
\end{align}
Here, $\xbf$ is apriori independent of $\zbf$;
the dependence between $\xbf$ and $\zbf$ manifests only aposteriori, i.e., after the measurement $\ybfbar$ is observed.
For $\xbf$, we assign the prior $p(\xbf;\thetabf_x)$, and
for $\zbf$  we assign the \emph{improper} (i.e., unnormalized) prior $p(\ybf|\zbf;\thetabf_z)$.

The lack of normalization will not be an issue in GVAMP, because the ``prior'' $p(\ybf|\zbf;\thetabf_z)$ is used only to compute posteriors of the form
\begin{align}
  p(\zbf|\ybf;\pbfhat,\tau,\thetabf_z) \propto
  p(\ybf|\zbf;\thetabf_z) \mathcal{N}(\zbf;\pbfhat,\Ibf/\tau) ,
\end{align}
which are well defined because the right side is always integrable over $\zbf$.

Let us first consider direct ML estimation of $\thetabf$ in the above SLM.
The $\thetabf$-likelihood function is
\begin{align}
p(\ybfbar;\thetabf)
&= \int p(\ybfbar,\xbfbar;\thetabf) \dif\xbfbar
= \int p(\ybfbar|\xbfbar) p(\xbfbar;\thetabf) \dif\xbfbar
= \int \mathcal{N}(\ybfbar;\Abfbar\xbfbar,\epsilon \Ibf) p(\xbfbar;\thetabf) \dif\xbfbar \\
&= \int \underbrace{ \mathcal{N}(\zbf;\Abf\xbf,\epsilon \Ibf) }_{\displaystyle \rightarrow \delta(\zbf-\Abf\xbf)} p(\xbf;\thetabf_x) p(\ybf|\zbf;\thetabf_z)\dif\xbf \dif\zbf
\label{eq:logybar} ,
\end{align}
which is consistent with \eqref{eq:logy} as $\epsilon\rightarrow 0$.
Likewise, for any belief $b(\xbfbar)$, we can upper bound the negative log-likelihood by a Gibbs free energy $\bar{J}(b;\thetabf)$ of the form
\begin{align}
\bar{J}(b;\thetabf)
&\defn D(b \,\|\, p(\xbfbar;\thetabf)) + D(b \,\|\, p(\ybfbar|\xbfbar)) + H(b)
\label{eq:Jbar} ,
\end{align}
since
$-\ln p(\ybfbar;\thetabf)
= \bar{J}(b;\thetabf) - D(b \,\|\, p(\xbfbar|\ybfbar;\thetabf))$
with $D(b \,\|\, p(\xbfbar|\ybfbar;\thetabf))\geq 0$.
The corresponding EM algorithm is
\begin{align}
\text{E step:~}& b^k(\xbfbar)
=p(\xbfbar|\ybfbar;\thetabfhat^k)
\label{eq:Estepbar}\\
\text{M step:~}& \thetabfhat^{\kp1}
=\arg\min_\thetabf \bar{J}(b^k;\thetabf)
=\arg\min_\thetabf D(b^k\,\|\, p(\xbfbar;\thetabf))
\label{eq:Mstepbar} .
\end{align}
As before, the posterior density required by the E-step \eqref{eq:Estepbar}
\begin{align}
p(\xbfbar|\ybfbar;\thetabf)
&= \frac{p(\xbfbar;\thetabf) p(\ybfbar|\xbfbar)}{p(\ybfbar;\thetabf)}
= \frac{p(\xbf;\thetabf_x) p(\ybf|\zbf;\thetabf_z) \delta(\zbf-\Abf\xbf)}
       {\int p(\xbf;\thetabf_x) p(\ybf|\zbf;\thetabf_z) \delta(\zbf-\Abf\xbf) \dif\zbf\dif\xbf} ,
\end{align}
is difficult to compute due to the high-dimensional integral.
Thus we consider an approximation afforded by the GVAMP algorithm \cite{schniter2016vectorAsil}, as described in the next section.

\section{GVAMP}

Recall that the exact posterior can (in principle) be found by solving the variational optimization problem
\begin{align}
p(\xbfbar|\ybfbar;\thetabf)
&= \arg\min_b D(b\,\|\, p(\xbfbar|\ybfbar;\thetabf)) \\
&= \arg\min_b \bar{J}(b;\thetabf) \label{eq:var1} \\
&= \arg\min_b D(b \,\|\, p(\xbfbar;\thetabf)) + D(b \,\|\, p(\ybfbar|\xbfbar)) + H(b)
\label{eq:var2} ,
\end{align}
where \eqref{eq:var1} follows from $\bar{J}(b;\thetabf)=D(b\,\|\, p(\xbfbar|\ybfbar;\thetabf))-\ln p(\ybfbar;\thetabf)$ and \eqref{eq:var2} follows from \eqref{eq:Jbar}.
But since the posterior computation problem is NP hard in general, \eqref{eq:var2} is no more tractable than any other approach.
The GVAMP algorithm computes a posterior approximation using the expectation-consistent (EC) method \cite{OppWin:05,fletcher2016expectation}.
In this application of EC, we first split $b(\xbfbar)$ into three copies, i.e.,
\begin{align}
p(\xbfbar|\ybfbar;\thetabf)
&= \arg\min_{b_1=b_2=q} D(b_1 \,\|\, p(\xbfbar;\thetabf)) + D(b_2 \,\|\, p(\ybfbar|\xbfbar)) + H(q) ,
\end{align}
and then relax the density-matching constraint $b_1=b_2=q$ to a moment-matching constraint:
\begin{align}
p(\xbfbar|\ybfbar;\thetabf)
&\approx \arg\min_{b_1,b_2,q} D(b_1 \,\|\, p(\xbfbar;\thetabf)) + D(b_2 \,\|\, p(\ybfbar|\xbfbar)) + H(q)
\label{eq:EC}
\\&\quad \nonumber
\text{~s.t.~} \Exp[\xbfbar|b_1]=\Exp[\xbfbar|b_2]=\Exp[\xbfbar|q]
\text{~~and~} \tr_2\{\Cov[\xbfbar|b_1]\}=\tr_2\{\Cov[\xbfbar|b_2]\}=\tr_2\{\Cov[\xbfbar|q]\} ,
\end{align}
where
$\Exp[\xbfbar|b_i]$ and $\Cov[\xbfbar|b_i]$ denote the expectation and covariance of $\xbfbar$ under $\xbfbar\sim b_i(\xbfbar)$, and where
\begin{align}
\tr_2\left(\mat{\Abf&\Bbf\\\Bbf\tran&\Cbf}\right)
&\defn \mat{\tr(\Abf)\\\tr(\Cbf)}
\text{~for~} \Abf\in\R^{N\times N}
\text{~and~} \Cbf\in\R^{M\times M} .
\end{align}
Essentially, $\tr_2\{\Cov[\xbfbar]\}$ separately computes the trace of the covariance of $\xbf$ and the trace of the covariance of $\zbf$.
The right side of \eqref{eq:EC} yields three different approximations of the posterior:
\begin{align}
b_1(\xbfbar;\thetabf)
&\propto  p(\xbfbar;\thetabf) \,\mathcal{N}\!\left(\xbfbar;\mat{\rbf_1\\\pbf_1},\mat{\Ibf_N/\gamma_1\\& \Ibf_M/\tau_1}\right) \label{eq:b1} \\
b_2(\xbfbar;\thetabf)
&\propto  p(\ybfbar|\xbfbar) \,\mathcal{N}\!\left(\xbfbar;\mat{\rbf_2\\\pbf_2},\mat{\Ibf_N/\gamma_2\\& \Ibf_M/\tau_2}\right) \label{eq:b2} \\
q(\xbfbar;\thetabf)
&\propto  \mathcal{N}\!\left(\xbfbar;\mat{\xbfhat\\\zbfhat},\mat{\Ibf_N/\eta\\& \Ibf_M/\zeta}\right) \label{eq:q} ,
\end{align}
where the form of \eqref{eq:b1}-\eqref{eq:q} can be deduced by analyzing the stationary points of the Lagrangian of \eqref{eq:EC}, as shown in \cite{OppWin:05}.

The GVAMP algorithm is an iterative approach to finding the values of $\rbf_1,\gamma_1,\pbf_1,\tau_1,\rbf_2,\gamma_2,\pbf_2,\gamma_2,\xbfhat,\eta,\zbfhat,\zeta$ under which the three beliefs in \eqref{eq:b1}-\eqref{eq:q} obey the moment constraints in \eqref{eq:EC}.
It is conjectured that the algorithm is guaranteed to converge for the GLM \eqref{eq:glm} when $\Abf$ is large and rotationally invariant, although this has not yet been rigorously proven.
Empirically, we find that the algorithm converges quickly in this scenario (e.g., on the order of $10$ iterations).

Note that the values of $\rbf_1,\gamma_1,\pbf_1,\tau_1,\rbf_2,\gamma_2,\pbf_2,\gamma_2,\xbfhat,\eta,\zbfhat,\zeta$ that satisfy the moment constraints are interdependent, and thus they all depend on the assumed value of $\thetabf$ through \eqref{eq:b1}.

\section{EM-GVAMP}

Recall that our current motivation for using GVAMP is to compute an approximation to the posterior $b^k(\xbfbar)=p(\xbfbar|\ybfbar;\thetabfhat^k)$ in the EM algorithm \eqref{eq:Estepbar}-\eqref{eq:Mstepbar}.
Of the three posterior approximations produced by GVAMP, the Gaussian approximation from \eqref{eq:q} is the simplest to use for this purpose.
Plugging the Gaussian approximation into \eqref{eq:Estepbar}-\eqref{eq:Mstepbar} yields
\begin{align}
\text{E step:~}& b^k(\xbfbar)
=\mathcal{N}\!\left(\xbfbar;\mat{\xbfhat^k\\\zbfhat^k},\mat{\Ibf_N/\eta^k\\& \Ibf_M/\zeta^k}\right)
\text{~found using GVAMP with $\thetabf=\thetabfhat^k$}
\label{eq:approxEstepbar}\\
\text{M step:~}& \thetabfhat^{\kp1}
=\arg\min_\thetabf D(b^k\,\|\, p(\xbfbar;\thetabf))
\label{eq:approxMstepbar} .
\end{align}
The difference between the EM algorithm
\eqref{eq:Estepbar}-\eqref{eq:Mstepbar}
and the EM algorithm
\eqref{eq:approxEstepbar}-\eqref{eq:approxMstepbar}
is that, in the former case, the bound is tight at each EM iteration $k$, whereas in the latter case the bound is only approximately tight.

Due to the form of $b^k(\xbfbar)$ in \eqref{eq:approxEstepbar}, the M-step is relatively easy to compute:
\begin{align}
\thetabfhat^{\kp1}
&=\arg\min_\thetabf D(b^k\,\|\, p(\xbfbar;\thetabf)) \\
&=\arg\min_\thetabf D\big(\mathcal{N}(\xbf;\xbfhat^k,\Ibf_N/\eta^k) \mathcal{N}(\zbf;\zbfhat^k,\Ibf_M/\zeta^k)\,\big\|\, p(\xbf;\thetabf_x) p(\ybf|\zbf;\thetabf_z)\big) \\
&=\arg\max_\thetabf \left\{
\int \mathcal{N}(\xbf;\xbfhat^k,\Ibf_N/\eta^k) \ln p(\xbf;\thetabf_x) \dif\xbf
+\int \mathcal{N}(\zbf;\zbfhat^k,\Ibf_N/\zeta^k) \ln p(\ybf|\zbf;\thetabf_z) \dif\zbf
\right\} \\
&=\arg\max_\thetabf \left\{
\sum_{j=1}^N \int \mathcal{N}(x_j;\xhat_j^k,1/\eta^k) \ln p(x_j;\thetabf_x) \dif x_j
+\sum_{i=1}^M \int \mathcal{N}(z_i;\zhat_i^k,1/\zeta^k) \ln p(y_i|z_i;\thetabf_z) \dif z_i
\right\}
\label{eq:approxMstepbar2} .
\end{align}
The resulting $(\thetabfhat_x^{\kp1},\thetabfhat_z^{\kp1})$ are necessarily values of $(\thetabf_x,\thetabf_z)$ that zero the gradient of the right side of \eqref{eq:approxMstepbar2} with respect to $\thetabf_x$ and to $\thetabf_z$.

\section{Application to Noise-Variance Estimation in Phase Retrieval}

Phase retrieval is a problem that can be formulated in the GLM setting \cite{Schniter:TSP:15}, allowing application of the GVAMP algorithm \cite{metzler2016prvamp}.
One way to model the $i$th measured intensity $y_i$ is via
\begin{align}
y_i &= \big|z_i+w_i\big|
\text{~for~i.i.d.~} w_i\sim \mathcal{N}(0,\nu_w)
\label{eq:yPR} ,
\end{align}
where $z_i,w_i\in\C$ and
$\mathcal{N}(w_i;\mu,\nu)=\frac{1}{\pi \nu}\exp(-|w_i-\mu|^2/\nu)$
represents a circular complex-Gaussian density with mean $\mu\in\C$ and variance $\nu>0$.
In this case, the measurement noise variance $\nu_w$ may be unknown in practice, and so we might try to estimate it using the methods described in this report.
In that case, the unknown $\zbf$-likelihood parameters ``$\thetabf_z$'' reduce to $\nu_w$.
In the sequel, we will use the notation $\nu_w$ instead of $\thetabf_z$.

It was shown \cite{Schniter:TSP:15} that, under \eqref{eq:yPR}, the $z_i$-likelihood function $p(y_i|z_i;\nu_w)$ takes the form
\begin{align}
p(y_i|z_i;\nu_w)
&= 1_{y_i\geq 0}~ y_i \int_0^{2\pi} \mathcal{N}(y_i e^{j\theta_i};z_i,\nu_w) \dif\theta_i \\
&= \frac{2y_i}{\nu_w} \exp\left(-\frac{y_i^2+|z_i|^2}{\nu_w}\right) I_0\left(\frac{2y_i|z_i|}{\nu_w}\right) 1_{y_i\geq 0} ,
\end{align}
where $I_0(\cdot)$ is the $0$th-order modified Bessel function of the first kind.
If we view $p(y_i|z_i;\nu_w)$ as a density on $y_i$, then $y_i$ is Rician (conditional on $z_i$).
Note that $\theta_i$ above denotes the (hidden) phase on $z_i+w_i$; it should not be confused with the statistical parameters $\thetabf$ described earlier in this report.

From \eqref{eq:approxMstepbar2}, we see that the EM estimate $\hat{\nu}_w^{\kp1}$ of $\nu_w$ must obey
\begin{align}
0
&=\frac{\partial}{\partial \nu_w}
\sum_{i=1}^M \int_\C \mathcal{N}(z_i;\hat{z}_i^k,1/\eta^k) \ln p(y_i|z_i;\hat{\nu}_w^{\kp1}) \dif z_i \\
&= \sum_{i=1}^M \int_\C \mathcal{N}(z_i;\hat{z}_i^k,1/\eta^k)
\frac{\partial}{\partial \nu_w}
\ln \int_0^{2\pi} \mathcal{N}(y_i e^{j\theta_i};z_i,\hat{\nu}_w^{\kp1}) \dif\theta_i \dif z_i \\
&= \sum_{i=1}^M \int_\C \mathcal{N}(z_i;\hat{z}_i^k,1/\eta^k)
\frac{
\int_0^{2\pi} \frac{\partial}{\partial \nu_w}
\mathcal{N}(y_i e^{j\theta_i};z_i,\hat{\nu}_w^{\kp1}) \dif\theta_i
}{\int_0^{2\pi} \mathcal{N}(y_i e^{j\theta_i'};z_i,\hat{\nu}_w^{\kp1}) \dif\theta'_i}
\dif z_i
\label{eq:zero} .
\end{align}

Plugging the derivative expression (see \cite{Vila:TSP:13})
\begin{align}
\frac{\partial}{\partial \nu_w}
        \mathcal{N}(y_i e^{j\theta_i};z_i,\hat{\nu}_w^{\kp1})
&= \frac{\mathcal{N}(y_i e^{j\theta_i};z_i,\hat{\nu}_w^{\kp1})}{(\hat{\nu}_w^{\kp1})^2}
        \big(
        |y_i e^{j\theta_i}-z_i|^2-\hat{\nu}_w^{\kp1}
        \big) ,
\end{align}
into \eqref{eq:zero} and multiplying both sides by $(\hat{\nu}_w^{\kp1})^2$, we find
\begin{align}
\hat{\nu}_w^{\kp1}
&= \frac{1}{M}\sum_{i=1}^M \int_\C \mathcal{N}(z_i;\widehat{z}_i^k,1/\eta^k)
        \frac{
        \int_0^{2\pi}
        |y_i e^{j\theta_i}-z_i|^2
        \mathcal{N}(y_i e^{j\theta_i};z_i,\hat{\nu}_w^{\kp1})
        \dif\theta_i}
        {\int_0^{2\pi}
        \mathcal{N}(y_i e^{j\theta_i'};z_i,\hat{\nu}_w^{\kp1})
        \dif\theta_i'}
        \dif z_i \\
&= \frac{1}{M}\sum_{i=1}^M
        \int_\C
        \mathcal{N}(z_i;\hat{z}_i^k,1/\eta^k)
        \int_0^{2\pi}
        |y_i e^{j\theta_i}-z_i|^2
        p(\theta_i;z_i,\hat{\nu}_w^{\kp1})
        \dif\theta_i \dif z_i
                \label{eq:EM1}
\end{align}
with the newly defined pdf
\begin{align}
p(\theta_i;z_i,\hat{\nu}_w^{\kp1})
&\triangleq \frac{ \mathcal{N}(y_i e^{j\theta_i};z_i,\hat{\nu}_w^{\kp1}) }
        {\int_0^{2\pi} \mathcal{N}(y_i e^{j\theta_i'};z_i,\hat{\nu}_w^{\kp1})
        \dif\theta_i'}
\propto
 \exp\Big(-\frac{|z_i-y_i e^{j\theta_i}|^2}{\hat{\nu}_w^{\kp1}}\Big)\\
&\propto
 \exp\big(\kappa_i \cos(\theta_i-\phi_i)\big)
 \text{~for~}\kappa_i \triangleq \frac{2|z_i|y_i}{\hat{\nu}_w^{\kp1}},
        \label{eq:ptheta}
\end{align}
where $\phi_i$ denotes the phase of $z_i$.
The proportionality \eqref{eq:ptheta} identifies this pdf as a von Mises distribution \cite{Mardia:Book:00}, which can be stated in normalized form as
\begin{align}
p(\theta_i;z_i,\widehat{\nu}_w^{\kp1})
&= \frac{\exp(\kappa_i \cos(\theta_i-\phi_i))}{2\pi I_0(\kappa_i)} .
        \label{eq:ptheta2}
\end{align}
Expanding the quadratic in \eqref{eq:EM1} and plugging in \eqref{eq:ptheta2}, we get
\begin{align}
\widehat{\nu}_w^{\kp1}
&= \frac{1}{M}\sum_{i=1}^M
        \int_{\C}
        \mathcal{N}(z_i;\widehat{z}_i^k,1/\eta^k)
        \bigg( y_i^2 + |z_i|^2

        - 2y_i |z_i|
        \int_0^{2\pi}
        \cos(\theta_i-\phi_i)\,
        \frac{\exp(\kappa_i \cos(\theta_i-\phi_i))}{2\pi I_0(\kappa_i)}
        \dif\theta_i
        \bigg)
        \dif z_i \\
&= \frac{1}{M}\sum_{i=1}^M
        \int_{\C}
        \mathcal{N}(z_i;\widehat{z}_i^k,1/\eta^k)
        \bigg( y_i^2 + |z_i|^2
        - 2y_i |z_i|
        R_0\bigg(\frac{2|z_i|y_i}{\widehat{\nu}_w^{\kp1}}\bigg)
        \bigg)
        \dif z_i
\label{eq:EM2} ,
\end{align}
where $R_0(\cdot)$ is the modified Bessel function ratio
$R_0(\kappa_i)\triangleq I_1(\kappa_i)/I_0(\kappa_i)$
and \eqref{eq:EM2} follows from \cite[9.6.19]{Abramowitz:Book:64}.

Simplifying approximations of \eqref{eq:EM2} could be taken as needed.
For example, in the high-SNR case, the expansion
$R_0(\kappa) = 1 - \frac{1}{2\kappa} - \frac{1}{8\kappa^2} - \frac{1}{8\kappa^3}+ o(\kappa^{-3})$
from \cite[Lemma~5]{Robert:SaPL:90}
could be used to justify
\begin{align}
R_0(\kappa)
&\approx 1 - \frac{1}{2\kappa} ,
\end{align}
which, when applied to \eqref{eq:EM2}, yields
\begin{align}
\widehat{\nu}_w^{\kp1}
&\approx \frac{2}{M}\sum_{i=1}^M
        \int_{\mathbb{C}}
        \big( y_i - |z_i| \big)^2
        \mathcal{N}(z_i;\widehat{z}_i^k,1/\eta^k)
        \dif z_i .
        \label{eq:EM3}
\end{align}
Approximation \eqref{eq:EM3} can be reduced to an expression that involves the mean of a Rician distribution.
In particular, using $z_i=\rho_i e^{j\phi_i}$, the integral in \eqref{eq:EM3} can be converted to polar coordinates as follows:
\begin{align}
\int_0^\infty \big( y_i - \rho_i \big)^2
        \hspace{-11mm}
        \underbrace{
        \int_0^{2\pi} \mathcal{N}(\rho_i e^{j\phi_i};\widehat{z}_i^k,1/\eta^k)
        \dif \phi_i \rho_i
        }_{\displaystyle
\frac{2\rho_i}{1/\eta^k}
\exp\left(-\frac{\rho_i^2+|\hat{z}_i^k|^2}{1/\eta^k}\right)
I_0\left(\frac{2\rho_i|\hat{z}_i^k|}{1/\eta^k}\right)
1_{\rho_i\geq 0}
        }
        \hspace{-11mm}
\dif\rho_i
&= y_i^2 - 2y_i \Exp[\rho_i] + \Exp[\rho_i^2] ,
\end{align}
where, for the expectations, $\rho_i$ has the Ricean density under the brace.
For this density, it is known that
\begin{align}
\Exp[\rho_i]
&= \sqrt{\frac{\pi}{4\eta^k}} \,L_{1/2}\!\left(-\eta^k|\hat{z}_i^k|^2\right) \\
\Exp[\rho_i^2]
&= 1/\eta^k + |\hat{z}_i^k|^2 ,
\end{align}
where the Laguerre polynomial $L_{1/2}(x)$ can be computed as
\begin{align}
L_{1/2}(x)
&= \exp\left(\frac{x}{2}\right)\left[(1-x)I_0\left(-\frac{x}{2}\right)-xI_1\left(-\frac{x}{2}\right)\right] .
\end{align}
Note that, for reasons of numerical precision, $\exp(x/2)I_d(-x/2)$ should be computed using ``\texttt{besseli($d$,$-x/2$,1)}'' in Matlab, not using ``\texttt{exp($x/2$).*besseli($d,-x/2$)}.''

\bibliographystyle{ieeetr}
\bibliography{../bibl.bib}
\end{document}

config: arxiv


TITLE: Graph Realization Problem
QUESTION [3 upvotes]: Given degrees of nodes, is it possible to construct a graph with those degrees, and if yes devise the algorithm?
I know of Handshaking Lemma that describes my problem, and another algorithm for it as follows:
The following algorithm decides if a simple graph can be constructed with given node degrees:

Sort the degrees in descending order.
If the first degree is $0$ (i.e. all degrees are $0$), then obviously such a graph can be formed (no edges) and you are done.
If the first degree has value $d$, then the following $d$ degrees must be greater than $0$. If not, you are done: no such graph can be formed.
Take away the first degree (value $d$) and reduce the following $d$ degrees by one (i.e. draw the requested number of edges from the node with highest degree to the nodes with highest degrees among the remaining ones), then continue with step 1 (with now one node less).

This algorithm takes $O(n^2 \log n)$ time where $n$ is number of nodes. Hence computationally very expensive!!  Is there an algorithm which solves this problem in $O(n\log n)$ or $O(n)$??
NOTE:Parallel edges and loops are not allowed

REPLY [1 votes]: The degree sequence problem asks essentially this, does there exist a graph that achieves the specified degrees of nodes?  Also called the graph realization problem.
The Havel-Hakimi algorithm constructs a solution, if one exists, recursively.  It has complexity $\Theta(n^2)$, but provides at most a single solution.
The Erdős–Gallai theorem can be recast, according to recent results by Tripathi, Venugopalanb, and West, as a construction of $O(n^3)$ worst case running time: "A short constructive proof of the Erdős-Gallai characterization of graphic lists", Discrete Math. 310, (4) (2010) 833–834.
More recently the Erdős-Gallai Linear algorithm was proposed by Iványi, Lucz, Móri, and Sótér (2011) claimed to have $\Theta(n)$ worst case running time.


text:
\begin{document}
\title{\bf A subexponential upper bound for the van der Waerden numbers $W(3,k)$}

\author{ By\\  \\{\sc Tomasz Schoen}}

\date{}

\maketitle

\begin{abstract} We show an improved upper estimate for the van der Waerden number $W(3,k):$ there is  an absolute  constant $c>0$ such that if  $\{1,\dots,N\}=X\cup Y$ is a partition such that $X$  does not contain any  arithmetic progression of length $3$ and $Y$ does not contain any  arithmetic progression of length $k$ then
$$N\le \exp(O(k^{1-c}))\,.$$
\end{abstract}

\section{Introduction}\label{s:intro}

Let $k$ and $l$ be positive integers. The van der Waerden number $W(k,l)$ is the smallest positive integer $N$ such that in any partition
$\{1,\dots, N\}=X\cup Y$ there is an arithmetic progression of length $k$ in $X$ or an arithmetic progression of length $l$ in $Y$. The existence of such numbers was established by van der Waerden \cite{vdW}, however the order of magnitude of $W(k,l)$ is unknown for  $k,l\ge 3.$ Clearly,  $W(k,l)$ is related to Szemer\'edi's theorem on arithmetic progressions \cite{szemeredi} and any effective estimate in this theorem leads to an upper bound on the van der Waerden numbers. Currently best known bounds in the most important
 diagonal case are
$$(1-o(1))\frac{2^{k-1}}{ek}\le W(k,k)\le {2^{2^{2^{2^{2^{k+9}}}}}}\,.$$
The upper bound follows from the famous work of Gowers \cite{gowers} and the lower bound was proved by Szab\'o \cite{szabo} using probabilistic argument.
Furthermore, Berlekamp \cite{berlekamp} showed that if $k-1$ is a prime number then
$$W(k,k)\ge (k-1)2^{k-1}\,.$$

 Another very intriguing instance  the problem are numbers $W(3,k)$ as they are related to Roth's theorem \cite{roth} that provides  more efficient estimates for sets avoiding three-terms arithmetic progressions. Let us denote by $r(N)$ the size of the largest progression-free subset of  $\{1,\dots, N\}$. We know that
\begin{equation}\label{roth}
r(N)\ll \frac{N}{{\log N}^{1-o(1)}}\,,
\end{equation}
see \cite{bloom, bloom-sisask,sanders,schoen}.
However this bound is not strong enough to imply a subexponential estimate for  $W(3,k)$.

Green \cite{green} proposed a  very clever argument based on  arithmetic properties of sumsets to bound $W(3,k).$ Building on this method and applying results from \cite{crs} it was showed in \cite{cwalina-schoen} that
$$W(3,k)\le \exp(O(k\log k))\,.$$
The best known lower bound was obtained by Li and Shu \cite{ls} (see also \cite{blr}), who showed that
$$W(3,k)\gg \Big(\frac{k}{\log k}\Big )^2.$$
The purpose of this paper is to prove a subexponential bound on $W(3,k)$.

\begin{theorem}\label{main} There are absolute constants $C,c>0$ such that for every $k$ we have
$$W(3,k)\le \exp(Ck^{1-c})\,.$$
\end{theorem}

Our argument is based on the method of \cite{schoen}, which explores in details the structure of a large spectrum.  This method can be partly applied (see Lemma \ref{increment}) in our approach and it deals only with a progression-free  partition class. The second part of the proof exploits the structure of both partition classes  and in this case   the argument of  \cite{schoen} has to be significantly modified.

\section{Notation} \label{s:notation}

The Fourier coefficients of a function $f:\Z/N\Z\to\C$ are defined
by
$$
\h f(r)=\sum_{x=0}^{N-1}f(x)e^{-2\pi ixr/N},
$$
where $r\in\Z/N\Z$. The inversion formula states that
$$f(x)=\frac1N\sum_{x=0}^{N-1}\h f(r)e^{2\pi ixr/N}.$$ We denote by $1_A(x)$  the indicator function of  set $A$.
Thus using the inversion formula and the fact that $\h{(1_A*1_B)}(r)=\hA(r)\hB(r)$ one can express the number of three--term arithmetic progressions (including trivial ones) by
$$\frac1{N}\sum_{r=0}^{N-1}\hA(r)^2\hA(-2r)=|A|\,.$$

 Parseval's identity asserts in particular  that
$$\sum_{r=0}^{N-1}|\hA(r)|^{2}=|A|N\,.$$

Let  $\th\ge 0$ be a real number,  the $\th-$spectrum of  $A$ is defined by
$$
\D_{\th}(A)=\big \{r\in\Z/N\Z:|\hA(r)|\ge\th|A|\big \}.
$$
If $A$ is specified then we  write $\D_{\th}$ instead of $\D_{\th}(A).$

By the span of a finite set $S$ we mean
$$\Span(S)=\Big\{\sum_{s\in S}\e_ss: \e_s\in \{-1,0,1\} \text{ for all } s\in S\Big\}$$
and the dimension of $A$ is defined by
$$\dim(A)=\min\big\{|S|: A\sbeq \Span (S) \big\}\,.$$
Chang's Spectral Lemma provides an upper bound
for the dimension of a spectrum.

\begin{lemma} {\rm\cite{chang}}\label{chang}
Let $A\subseteq \Z/N\Z$ be a set of size $|A|=\d N$ and let $\th>0.$ Then
$$\dim(\D_\th(A))\ll \th^{-2}\log(1/\d)\,.$$
\end{lemma}
We are going to use  Bohr sets \cite{bourgain-1/2} to prove the main result.

Let $\G\sbeq\h G$ and   $\g\in(0,\frac{1}{2}]$ then the Bohr set generated by $\G$  with radius $\g$  is
$$
B(\G,\g)
  =  \big\{ x\in \Z/N\Z: \ \|{tx}/{N}\|\le\g \text { for all } t\in\G\big\}\,,$$
where $\left\Vert x\right\Vert =\min_{y\in\Z}|x-y|$.
 The rank of $B$ is the size of $\G$ and we denote it by $\rk(B).$
Given $\eta>0$ and a
Bohr set $B=B(\G,\g)$, by $B_{\eta}$ we mean the Bohr set $B(\G,\eta\g).$
We will use two basic properties of Bohr sets concerning its size and regularity, see \cite{tao-vu}.

\begin{lemma} {\rm\cite{bourgain-1/2}}\label{bohr-size}
 For every $\g\in(0,\frac{1}{2}]$ we have
$$
\g^{|\G|}N\le|B(\G,\g)|\le 8^{|\G|+1}|B(\G,\g/2)|\,.
$$
\end{lemma}
We call a Bohr set $B(\G,\g)$ {regular} if for every $\eta$,
where $|\eta|\le 1/(100|\G|)$ we have
$$
(1-100|\G||\eta|)|B|\le|B_{1+\eta}|\le(1+100|\G||\eta|)|B|.
$$
Bourgain \cite{bourgain-1/2} showed that regular Bohr sets are ubiquitous.
\begin{lemma} {\rm\cite{bourgain-1/2}}\label{bohr-regular} For every Bohr set $B(\G,\g),$ there exists
$\g'$ such that $\frac{1}{2}\g\le\g'\le\g$ and $B(\G,\g')$ is
regular.
\end{lemma}

\section{Proof of Theorem \ref{main}}
\label{s:large-coeff}

Our main tool is the next lemma, which can be extracted from \cite{schoen} (see Lemmas 7, 9, 12 and 13). Its proof  makes use of the deep result by Bateman and Katz in \cite{bateman-katz, bateman-katz-nonsmooth} describing  the structure of the large
spectrum.

\begin{lemma} {\rm\cite{schoen}}\label{increment}
There exists an absolute constant $c>0$ such that   the following holds.
 Let $A\subseteq \Z/N\Z,\, |A|=\d N$ be a set without any non-trivial arithmetic progressions of length three and such that
\begin{equation}\label{mid}
\sum_{r:\, \d^{1+c} |A|\le |\hA(r)|\le \d^{1/10} |A|}|\hA(r)|^3\ge \frac1{10}\d^{c/5}|A|^3\,.
\end{equation}
Then  there is a regular Bohr set $B$ with
$\rk( B)\ll \d^{-1+c}$  and  radius $\Omega( \d^{1-c})$ such that     for some $t$
    $$|(A+t)\cap B|\gg \d^{1-c}|B|.$$
\end{lemma}

Furthermore, we  apply Bloom's iterative lemma, that provides a density increment by  a constant factor greater than $1$ for progression-free sets and 	Sanders' lemma on a containment of long arithmetic progressions in dense subsets of regular Bohr sets.

\begin{lemma}\label{bloom} {\rm\cite{bloom}} There exists an absolute constant $c_1 > 0$ such that the following
holds. Let $B\sbeq \Z/N\Z$ be a regular Bohr set of rank $d$. Let
$A_1 \sbeq B$ and $A_2 \sbeq B_\e,$ each with relative densities
$\a_i$. Let $\a = \min(c_1, \a_1, \a_2)$ and assume that   $d \le
\exp(c_1(\log^2(1/\a)).$ Suppose that   $B_\e$ is also regular and
$c_1\a/(4d) \le \e \le c_1\a/d.$ Then either
\begin{itemize}
\item[(i)] there is a regular Bohr set $B'$ of rank $\rk(B') \le d +
O(\a^{-1}\log(1/\a))$ and size
$$ |B'| \ge \exp\big(-O(\log^2(1/\a)(d + \a^{-1}\log(1/\a)))\big)|B|$$
such that
$$|(A_1+t)\cap B'|\gg (1+c_1)\a_1|B'|$$
for some $t\in \zn$;
\item[(ii)] or there are $\Omega(\a_1^2\a_2 |B| |B_\e|)$ three-term arithmetic progressions $x+y=2z$ with $x,y\in A_1, z\in A_2$;
\end{itemize}
\end{lemma}

\begin{lemma} {\rm\cite{sanders}}\label{sanders}
Let $B(\G,\g)\sbeq \zn$  be a regular Bohr set of rank $d$   and let $\e$ be a positive number satisfying  $\e^{-1}\ll \g d^{-1}N^{1/d}$. Suppose
that  $A\sbeq B$ contains at least a proportion $1-\e$ of $B(\G, \g).$ Then
  $A$ contains an arithmetic
progression of length at least $1/(4\e).$
\end{lemma}

\no {\bf Proof of Theorem \ref{main}.}  Put $M=W(3,k)-1$ and let  $\{1,\dots,M\}=X\cup Y$ be a partition such that $X$ and $Y$ avoid $3$ and $k$-term arithmetic progressions respectively. Clearly, we may assume that $M\ge 100k$ hence
$$|Y|\le M-\lf M/k\rf\le M-M/(2k),$$
as no block of $k$ consecutive numbers is contained in $Y$ and therefore $|X|\ge M/(2k).$
Let $N$ be any prime number satisfying $2M<N\le 4M.$ We embed $\{1,\dots,M\}=X\cup Y$ in $\zn$ in a natural way and observe that $X$ and $Y$ possess the same properties. Put $|X|=\d N$.

First let us assume that
\begin{equation}\label{mid}
\sum_{r:\, \d^{1+\m} |X|\le |\hX(r)|\le \d^{1/10} |X|}|\hX(r)|^3\ge \d^{\m/5}|X|^3\,,
\end{equation}
Then by Lemma \ref{increment}   there is  $t\in \zn$ and a regular Bohr set $B^0$ with
$\rk( B^0)\ll \d^{-1+c}$  and  radius $\Omega( \d^{1-c})$ such that
    $$|(X+t)\cap B^0|\gg \d^{1-c}|B|$$
		for some absolute constant $c>0.$
 Writing  $A_0=(X+t)\cap B_0$ we have
$$|X_0\cap B^0|\gg \a|B^0|\,,$$
where
$$\a\gg \d^{1-c}\,.$$
By Lemma \ref{bohr-size} we have
$$|B^0|\ge \exp\big (-O(\d^{-1+c}\log(1/\d))\big)N\,.$$
Next, we
iteratively  apply  Lemma \ref{bloom}. Since after each step the density increases by factor $1+c_1$ it follows that  after  $l\ll \log (1/\a)$ steps  case $(ii)$ of Lemma \ref{bloom} holds.   Let $B^i$ be Bohr sets obtained in the iterative procedure and observe that $\rk(B^i)\ll \a^{-1}\log^2(1/\a)$ for every $i\le l$.  Therefore, there are
$$\Omega( \a^3 |B^k||B^l_\e|)$$
three-term arithmetic progressions in $X,$ where $\e\ge c_1\a/(4\rk(B^l))\gg \a^2\log^2(1/\a).$
By Lemma \ref{bloom} and  Lemma \ref{bohr-size}
we have
$$|B^l|\ge \exp\big(-O(\a^{-1}\log^4(1/\a))\big)N\ge  \exp\big(-O(\d^{-1+c}\log^4(1/\d))\big)N\,,$$
and
\begin{eqnarray*}
|B^l_\e|&\ge& \exp\big(-O(\a^{-1}\log^3(1/\a))\big)\exp\big(-O(\a^{-1}\log^4(1/\a))\big)N\\
&\ge&  \exp\big(-O(\d^{-1+c}\log^4(1/\d))\big)N\,.
\end{eqnarray*}
Thus, $X$ contains
 $$\Omega (\d^{3-3c}
\exp\big(-O(\d^{-1+c}\log^4(1/\d))\big)N^2)$$
arithmetic progressions of length three.
Since there are only $|X|$ trivial progressions in $X$ it follows that
$$|X|\gg \d^{3-3c} \exp\big(-O(\d^{-1+c}\log^4(1/\d))\big)N^2\,,$$
so
$$W(3,k)\ll N\ll \exp\big(O(\d^{-1+c}\log^4(1/\d)\big)\le \exp\big(O(k^{1-c}\log^4k)\big)\,.$$

Next let us assume that \eqref{mid} does not hold.
By Chang's lemma
$$d:=\dim(\D_{\d^{1/10}}(X))\ll \d^{-1/5}\log (1/\d)$$
 hence there is a set $\La$ such that $|\La|=d$ and
$\D_{\d^{1/10}}\sbeq \Span(\La)$. By Lemma \ref{bohr-regular} there is   a regular Bohr set $B=B(\Lambda,\g)$   with radius $\g\gg \d^3$. Let $\b=\frac1{|B|}1_B$ then for every $r\in \D_{\d^{1/10}}$ we have
\begin{equation}\label{beta}
\big|\h \b(r)-1\big|\le \frac1{|B|}\sum_{b\in B} |e^{-2\pi i\lambda b/N}-1|\le \frac{2\pi}{|B|}\sum_{b\in B} \sum_{\lambda\in \La}\|rb/N\|\le 2\pi \d^2\,,
\end{equation}
and similarly $|\h \b(2r)-1| \ll \d^2.$
For $t\in \zn$ put
$$f(t)=\b*1_X(t)\,$$
and note that if for some $t\in [4\g M,(1-4\g)M]$ we have   $f(t)=\frac1{|B|}|X\cap (B+t)|\le \d^{1+c'},$ where  $c'=c/20,$ then since $B+t\sbeq [1,M]$ it follows that
$$|Y\cap (B+t)|\ge 1-\d^{1+c'}\,.$$
Therefore, by Lemma \ref{sanders} either $\d^{-1-c'}\gg \g d^{-1}N^{1/d}$ or $Y$ contains an arithmetic progression of length
$\frac14\d^{-1-c'}$. The former inequality implies that
$$k^{1+c'}\gg\g d^{-1}N^{1/d}\gg \d^4N^{O(\d^{1/5}\log^{-1} (1/\d))}\gg k^{-4} N^{O(k^{-1/5}\log^{-1} k)}\,,$$
so
$$W(3,k)\ll N\ll \exp\big(O(k^{1/5}\log^2 k)\big)\,.$$ If the second alternative holds then
$$\frac14\d^{-1-c'}<k$$
hence by \eqref{roth}
$$k^{-1/(1+c')}\ll \d\ll (\log N)^{-1+o(1)}$$
so
$$W(3,k)\le N\ll \exp\big(O(k^{\frac1{1+c'}+o(1)})\big)\,.$$
Finally we can assume that for every $t\in [4\g M,(1-4\g)M]$ we have   $f(t)\ge \d^{1+c'}.$
Let $T(X)$ denote the number of three-term arithmetic progression in $X$ and let
$$T(f)=\sum_{x+y=2z}f(x)f(y)f(z)\,.$$
Then clearly
\begin{equation}\label{f}
T(f)\gg \d^{3+3c'}M^2\gg \d^{3+c/6}N^2
\end{equation}
and we will show that  $T(X)$  does not differ  much from $T(f)$
\begin{eqnarray}\label{t}
\big|T(X)-T(f)\big|&=&\frac1{N}\big|\sum_{r=0}^{N-1} \h {1_X}(r)^2\hX(-2r)-\sum_{r=0}^{N-1} \h f(r)^2\h f(-2r)\big|\nonumber\\
&\le&\frac1{N}\sum_{r=0}^{N-1} |\hX(r)^2\hX(-2r)(1-\h \b(r)^2\h \b (-2r))|\\
&=& S_1+S_2+S_3,\nonumber
\end{eqnarray}
where $S_1,S_2$ and $S_3$ summations of \eqref{t} respectively over $\D_{\d^{1/10}},\D_{\d^{1+c}}\setminus \D_{\d^{1/10}}$ and $\zn\setminus \D_{\d^{1+c}}.$
By \eqref{beta}, \eqref{mid}, Parseval's formula and  H\"older's inequality  we have
$$S_1\ll \d^2 \frac1{N}\sum_{r\in \D_{\d^{1/10}}}|\hX(r)|^3\le \d^3\sum_{r=0}^{N-1} |\h {1_X}(r)|^2= \d^2|X|^2\,,$$
$$S_2\le \d^{1+c/5}|X|^2$$
and
$$S_3\le \frac2{N} \d^{1+c}|X| \sum_{r=0}^{N-1} |\h {1_X}(r)|^2=\d^{1+c}|X|^2\,.$$
Thus,
$$\big|T(X)-T(f)\big|\ll \d^{3+c/5}N^2,$$
so by \eqref{f} and the fact that $X$ avoids non-trivial three-term arithmetic progression we have
$$|X|=T(X)\gg \d^{3+c/6}N^2\,,$$
hence
$$W(3,k)\le N\ll \d^{-2-c/6}\ll k^3$$
which concludes the proof.$\hfill\Box$


text:
\begin{document}

\title[Classification of barrier options]{
Classification of barrier options}
\author[J.C. Ndogmo]{J.C. Ndogmo}

\address{P O Box 2446\\
Bellville 7535\\
South Africa.}

\email{ndogmoj@yahoo.com}

\keywords{ }

\begin{abstract}
For a given level of accuracy in option prices, the paper considers
the problem of deciding when exactly, as one or more of the pricing
parameters change, a barrier option degenerates into a simpler type
of option. This problem is meaningful in the real world where option
prices are always determined within a certain level of accuracy. The
problem is reduced to finding certain critical values of the initial
stock price, and this is achieved through a probability-based
approach.
\end{abstract}

\keywords{Critical stock price,  Barrier options classification,
non-dimensional approximate parameter}

\subjclass[2000]{91b24, 65c30}

\maketitle

\section{Introduction}
\label{s:intro}
One of the topics that financial mathematics has been concerned with
over the last decades is that of option pricing, and barrier option
pricing features prominently among these topics, in particular
because barriers can be added to any existing option. In this
regard, Merton \cite{merton} obtained in $1973$ the first pricing
formula for a European down-and-out call option by solving
explicitly the associated differential equation with the appropriate
boundary conditions. Latter on in 1985, Cox and Rubistein~
\cite{cox} obtained a similar formula for the up-and-out barrier.
Using a probabilistic approach based on certain theoretical results
derived earlier on by Levy~ \cite{levy} and by Anderson~
\cite{ander}, Kunitomo and Ikeda~ \cite{kuni} obtained more general
pricing formulas for European double barrier options with  curved
barriers and for a variety of path-dependent options and corporate
securities, some of which had already been independently obtained
for restricted cases in \cite{coxb,inger, connell}.\par
 Attempts to extend these pricing formulas to the more complex case of
American-style options have been undertaken in papers by Broadie and
Detemple~ \cite{broadie}, Gao {\em et al}~ \cite{gao}, and Haug~
\cite{haug2}. All of these extensions to American-style options have
however only lead to close-form approximations. In parallel with the
determination of these pricing formulas, numerical methods have been
used for pricing barrier options, especially in those cases where
analytical pricing solutions remain unavailable, such as for
discrete barrier options \cite{MC, TT, MCh, fusai}, or for
American-style options that, in some cases, comprise transaction
costs or stochastic volatilities \cite{barraq, gao, clarke,
ikonen}.\par

An important problem about barrier options that doesn't seem to have
been considered is that of deciding when exactly, for a given level
of accuracy in option prices, a double barrier option degenerates
into a much elementary type of barrier option, as one or more of the
pricing parameters change with time. Although this problem is
meaningless in terms of the Brownian motion of the of the underlying
asset price, it becomes fully meaningful in the real word where all
prices are known only to within a finite number of precision digits.
It is also a practical problem in the sense that it can lead to a
more efficient pricing of various types of barrier options, and to a
better insight into the dynamics of options, and more specifically
barrier options.\par

In this paper, we consider the problem of determining when exactly,
as the initial stock price or other pricing parameters evolve, and
for a fixed level of accuracy in option prices, a barrier option
with a single or a double barrier degenerates into a simpler type of
option, which may be, for instance in the case of a double barrier
option, a down-and-out or an up-and-out barrier option, or just a
vanilla. We consider in this analysis more general types of barriers
with time-dependent absorbing boundaries. The problem is reduced to
finding certain critical values of the initial asset price, and we
use a probabilistic approach to derive close-form approximate
solutions for these values. Numerical experiments based on
analytical option pricing results such as those from \cite{kuni} are
then implemented, to check the validity of the formulas
obtained.\par

In the next section, we derive formulas for the classifying critical
asset prices which in  Section $3,$ are applied to the problem of
classification of barrier options. We present the numerical test of
these formulas in  Section $4,$ and give some concluding remarks in
Section~ \ref{s:conclusion}. Due to the usual parity considerations,
we focus our analysis solely on knock-out barrier options.

\section{Classifying critical asset prices}\label{s:critical}

We make the usual assumption that the stock price $S_t$ follows a
geometric Brownian motion, or equivalently, that it is lognormally
distributed. In this case, the stochastic differential equation
satisfied by $S_t= S(t)$ is given in terms of the standard Wienner
process $W$ by

\begin{equation} \label{eq:gbm} d S_t= S_t \mu dt + S_t \sigma d W, \end{equation}

where $\mu$ and $\sigma$ are the constant drift and volatility
parameters, respectively. Consider two continuous curves $S= B_l(t)$
and $S=B_u(t)$ over the time interval $I=[0,T],$ such that $B_l <
S_0 < B_u,$ where $B_l=B_l(0)$,  $B_u= B_u(0),$  $S_0=S(0)$ and $T$
is a given time. Assume furthermore that the two curves are
non-overlapping, so that we have $B_l(t)< B_u (t)$ for all $t \in
I.$  A sample path of stock price $S_t$ together with the two
adjacent curves are depicted in Figure~ \ref{fg:path1}.\par

\begin{figure}[ht]
\begin{center}
\caption{\label{fg:path1} \protect \footnotesize Sample asset price
path together with lower and upper barriers. \vspace{2mm}}
\includegraphics[width= \textwidth]{dblepath1.eps}
\end{center}
\end{figure}

 Denote by $\theta= \theta(f)$ the accuracy
in a given option price $f.$ For simplicity, we may write
$\theta(f)$ in the form $\theta(f)= 10^{- \mathsf{m}},$ where
$\mathsf{m}$ is the number of significant digits to the right of the
decimal point in $f.$ Suppose now that $f$ is the price of a double
knock-out option contingent on $S_t,$ with curved barriers $B_u(t)$
and $B_l(t).$ For any option under consideration, unless otherwise
stated we denote by $S_0$ the initial price of the underlying asset,
by $K$ its strike price and by $T$ its expiry date. Let $\mathsf{v}$
be a variable symbol that may take on the value $u$ to mean "upper"
or $l$ to mean "lower". Let $E_\mathsf{v}$ be the event that the
$\mathsf{v}$ barrier is breached before the other barrier by time
$T,$ and let ${\rm R}_\mathsf{v}$ be the corresponding rebate paid
at the expiry date of the option. If we denote the risk-free
interest rate by $r,$ by $\prob$ the probability operator, and by
$\Esp$ the risk neutral expectation operator, then the price $f$ of
the double knock-out option is given by

\begin{equation}\label{eq:gndbl}
f= e^{-r T}\Esp \left[ {\rm R}_u \prob (E_u)+ {\rm R}_l \prob(E_l) +
(1-\prob(E_u)-\prob(E_l)) \prob (S_T > K) (S_T -K)\right]
\end{equation}

   In terms of the geometric Brownian motion of the underlying, the
probability that the stock price breaches a barrier before any given
time $T$ is never zero. Note that the operator $\Esp$ is additive
with respect to each term in its argument in \eqref{eq:gndbl}. By
the same properties of the geometric Brownian motion, by choosing
$S_0$ sufficiently far away from an initial barrier value
$B_\mathsf{v},$ we can make $\prob (E_\mathsf{v})$ so small that if
we let $\Delta f$ be the difference between the corresponding value
of $f$ and the value of $f$ obtained by setting $\prob
(E_\mathsf{v})=0,$ then $\Delta f < 5\times 10^{-1} \theta (f).$
This in terms of the accuracy of $f$ means that $\Delta (f)=0,$ and
thus the barrier $B_\mathsf{v}(t)$ has therefore no effect on the
value of $f$, and the double knock-out option degenerates into a
simpler type of option, which may be a down-and-out barrier option
if $\mathsf{v}=u,$ or an up-and-out barrier option if
$\mathsf{v}=l.$ In the real world markets where all prices are
determined only in terms of a finite number of precision digits, we
may thus assume that the probability of a stock price breaching a
barrier before a given time $T$ can take on the value zero for
appropriately chosen values of the initial price $S_0$.\par

It is therefore meaningful to consider the following two
complementary problems, for a given accuracy $\theta$ of option
prices.
\begin{description}
\itemsep = 1mm
\topsep = 2.5mm

\item[{\small Problem 1}] Find the minimum value $S_{ml}$ of $S_0$ which ensures that
the stock price will not breach the lower curve before time $T.$
\item[{\small Problem 2}] Find the maximum value $S_{mu}$ of $S_0$ which ensures
that the stock price will not breach the upper curve before time
$T.$
\end{description}

\begin{defn}
We say a curve $S=g(t)$ is worthless with respect to the stock price
movement over a given time interval and for a given accuracy
$\theta$  in stock prices if, with probability $1,$ the path
followed by $S_t$ will not breach the curve (either from above or
from below) over that time interval.
\end{defn}

In other words, the problem we are interested in is that of finding
the minimum (maximum) value $S_{ml}\; (S_{mu})$ of the initial stock
price that will ensure that the lower (upper) curve  is worthless
over $[0,T].$  Let $P^{\, l}_t= \prob (S_t > B_l(t))$ be the
probability that $S_t > B_l(t).$ Then by Eq.~ \eqref{eq:gbm} we
clearly have

\begin{align}
\label{eq:pl0}  P^{\, l}_t &= \Phi(\omega^{\,l}_t) \notag \\[-4mm]
\intertext{ where }
                \omega^{\,l}_t &= \frac{(\mu- \sigma^2/2)t + \ln \left(S_0/ B_l(t)\right)}{\sigma
                \sqrt{t}},
\end{align}
and where $\Phi $ is the standard normal distribution function. It
is well known that $\Phi$ is a positive definite increasing function
satisfying $\lim_{x \ra \infty} \Phi(x)=1.$  However, this
convergence of $\Phi$ to $1$ is very fast and if we denote by $\pi=
\pi (\Phi)$ the accuracy of $\Phi,$ then it is easy to see that we
have for instance $\Phi (4.753424)=1,$ for $\pi= 10^{-6}.$ For a
given value of $\pi,$ let $\nu$ be the smallest number for which the
approximation $\Phi(\nu)=1$ holds, that is, such that $\phi(\nu)
\geq 1- \pi.$ We therefore have $P_t^l=1$ if and only if $\omega_t^l
\geq \nu,$ and solving this inequality for $S_0$ shows that

\begin{subequations}\label{eq:pl1}
\begin{align}
 P^{\, l}_t &=1 \Leftrightarrow   S_0 \geq S_l(t)\label{eq:pl1.1} \\[-4mm]
 \intertext{ where }
S_l (t)&= B_l (t)  e^{(\nu \sigma \sqrt{t}- \mu_1 t)} \label{eq:slt}
\\[-4mm]
\intertext{ and where }

\mu_1&= \mu - \sigma^2 /2 \label{eq:mu1} \\[-2mm] \notag
\end{align}
\end{subequations}
 Now, let $S_m$ be the maximum value of the continuous function $S_l(t)$ achieved
over the interval $I.$ Then by compactness we have $S_m= S_l(t_m)$
for some $t_m \in I.$ For $S_0 \geq S_m,$ Eq.~ \eqref{eq:pl1.1}
shows that $P^{\, l}_t=1$ for all $t \in I,$ that is the lower curve
is worthless. On the other hand, if $S_0 < S_m,$ then
\eqref{eq:pl1.1} shows again that $P^{\, l}_{t_m}\neq 1,$ and thus
the lower barrier isn't worthless. Consequently, $S_m$ is the
minimum value of the initial stock price above which the lower curve
becomes worthless, that is, $S_m= S_ml$ is the solution to Problem
1.\par

Similarly, let $P^u_t = \prob (S_t < B_u (t))$ be the probability
that $S_t < B_u (t).$ Then we have
\begin{align}
\label{eq:pu0}  P^u_t &= \Phi(\omega^u_t), \notag \\
\intertext{ where }
                \omega^u_t &= \frac{\ln\left(B^u(t)/S_0\right) - \mu_1 t}{\sigma
                \sqrt{t}},
\end{align}
and using again the fact that $P^u_t=1$ if and only if $\omega^u_t
\geq \nu,$ it readily follows that

\begin{subequations}\label{eq:pu1}
\begin{align}
 P^u_t &=1 \Leftrightarrow   S_0 \leq S_u(t),\label{eq:pu1.1} \\[-4mm]
 \intertext{ where }
S_u (t)&= B_u (t)  e^{-(\nu \sigma \sqrt{t}+ \mu_1 t)}.
\label{eq:sut}
\end{align}
\end{subequations}

Consequently, if we let $S_{\sc M}$ be the minimum value of the
continuous function $S_u(t)$ over the interval $I,$ then it can be
shown in a similar manner as in the preceding case for the lower
curve that $S_{\sc M}$ is the maximum value of the initial stock
price $S_0$ below which the upper curve is worthless. That is,
$S_{\sc M}= S_{mu}$ is the solution to Problem $2.$ We have thus
obtained the following result.

\begin{thm}\label{th:critic}
Suppose that the stock price $S_t$ follows a geometric Brownian
motion of the form \eqref{eq:gbm}. Suppose also that the two
functions $S_l (t)$ and $S_u(t)$  given by \eqref{eq:slt} and
\eqref{eq:sut}, respectively, are continuous and considered over the
time interval $I=[0,T].$ Let $S_{ml}$ be the maximum value of
$S_l(t)$ and $S_{mu}$ the minimum value of  $S_u(t).$ Then we have

\begin{enumerate}
\item[(a)] $S_{ml}$ is the minimum value of the initial stock price  $S_0$ above
which the curve $S= B_l(t)$ becomes worthless.
\item[(b)] $S_{mu}$ is the maximum value of the initial stock price
$S_0$ below which the curve $S= B_u (t)$ becomes worthless.
\end{enumerate}
\end{thm}

\section{Applications to barrier options}
\label{s:applications} We now consider more explicitly the problem
of determining when exactly a barrier option contingent on a stock
price $S_t$ will degenerate into a simpler type of option, as a
result of a change in the option's pricing parameters.  All barrier
options will be assumed to have, in general curved barriers, and the
underlying asset price $S_t$ is supposed to follow a geometric
Brownian motion of the form \eqref{eq:gbm}. We shall also denote by
$B_l(t)$ and by $B_u (t)$ the lower barrier and the upper barrier,
respectively, and if any, for a given barrier option.

\begin{defn}
We say that a barrier option of a given type is a typical barrier of
that type if none of its barriers can be considered worthless for
the corresponding parameter set. That is, if it cannot be priced as
a barrier option of a different type with the same parameter set.
\end{defn}

The equations \eqref{eq:slt} and \eqref{eq:sut} show that in
addition to the parameters of the curved barriers, the critical
values $S_{ml}$ and $S_{mu}$ of Theorem~\ref{th:critic} depend only
on the drift and volatility parameters $\mu$ and $\sigma$ of the
stock price $S_t,$ and on the expiry time $T.$ In a risk neutral
world, we have $\mu= r - \rho,$ where $r$ is the interest rate and
$\rho$ is the continuous dividend yield, if any, paid by the
underlying asset. However, the value of $\rho,$ is of no importance
for our analysis and we may assume without loss of generality that
$\rho=0.$ By the fixed parameter set for a given option, we shall
therefore refer to the other parameters $T, \sigma,$ and $r,$ as
well as the initial asset price $S_0$, the pay-off function
parameters such the strike price $K,$ and the curved barriers'
parameters.

\begin{thm} \label{th:dao}
Consider a down-and-out barrier option with  {\rm (}lower {\rm )}
barrier $S~=~B_l(t),$ and a given parameter set.
\begin{enumerate}
\item[(a)] If $S_0 \geq S_{ml},$ the barrier is worthless and the
option is equivalent to a vanilla option with the same parameter
set.
\item[(b)] The option is a typical down-and-out barrier option if and only
if \\$B_l(0)< S_0 < S_{ml}.$
\end{enumerate}
\end{thm}

\begin{proof}
Part (a) of the theorem readily follows from the definition of
$S_{ml},$ and the fact that a down-and-out option with a worthless
 barrier is naturally equivalent to a vanilla option with the same parameter
set. On the other hand, if $B_l(0)< S_0 < S_{ml},$ there is no
certainty that the barrier will not be breached before the expiry
date $T,$ and thus the option is necessarily a typical down-and-out
barrier option.
\end{proof}

\begin{thm} \label{th:uao}
Consider an up-and-out barrier option with {\rm (}upper{\rm )}
barrier $S~=~B_u(t),$ and a given parameter set.
\begin{enumerate}
\item[(a)] If $S_0 \leq S_{mu},$ the barrier is worthless and the
option is equivalent to a vanilla option with the same parameter
set.
\item[(b)] The option is a typical up-and-out barrier option if and only
if \\$S_{mu}< S_0 < B_u(0).$
\end{enumerate}
\end{thm}

\begin{proof}
The proof also follows from the definition of $S_{mu},$ and is
similar to that of Theorem~\ref{th:dao}.
\end{proof}

The following result is also straightforward from the definitions of
$S_{ml}$ and $S_{mu}$ given by Theorem~ \ref{th:critic}.

\begin{thm} \label{th:db}
Consider a double barrier option with lower barrier $S~=~B_l(t)$ and
upper barrier $S~=~B_u(t),$ and a fixed parameter set.
\begin{enumerate}
\topsep = 2.5mm
\itemsep = 1mm
\item[(a)] If $S_0 < S_{ml}, S_{mu},$ the option degenerates into a
down-and-out barrier option.
\item[(b)] If $S_{ml}, S_{mu} < S_0,$ the option degenerates into an
up-and-out barrier option.
\item[(c)] The option is equivalent to the vanilla option with the
           same parameter set if and only if  $S_{ml} \leq S_0 \leq S_{mu}.$
\item[(d)] The option is  a typical double barrier
option if and only if $S_{mu} < S_0 < S_{ml}.$
\end{enumerate}
\end{thm}

It should be noted that the results of the  theorems \ref{th:dao} -
\ref{th:db} do not depend on the pay-off function of the options
considered, and consequently they may be applied to all options
endowed with barriers, and in particular they apply to both barrier
call and barrier put options, European as well as American options,
and options including a rebate payment. The importance of these
theorems also stem from the fact that barriers can be added to
virtually all sorts of existing options, and this is one of the
reasons why barrier options have become so popular. An important
particular case is that of constant barriers. In such case, the
lower and upper curves, if any, for a given option can be denoted by
$S=B_l$ and $S=B_u,$ respectively, where $B_l$ and $B_u$ are
constants, and we can therefore give  more explicit expressions for
the critical values $S_{ml}$ and $S_{mu}.$ It is easy to see that
the  functions $S_l(t)$ and $S_u(t)$ given by the equations
\eqref{eq:slt} and \eqref{eq:sut} have in this particular case a
common expression for the turning point $t_p,$ when it exists,
given by

\begin{equation} \label{eq:tp}  t_p= \left(  \frac{\nu \sigma}{2 \mu_1 }\right)^2.
\end{equation}
Consequently we have the following explicit expressions for the
critical asset prices $S_{ml}$ and $S_{mu}.$

\begin{subequations}\label{eq:fcritic}
\begin{align}
 S_{ml} &= \begin{cases} B_{l \,} e^{(\nu \sigma
\sqrt{T} - \mu_1 T)},& \text{ if  $\mu_1 \leq 0$ or $ t_p \geq T$ }\\
B_{l \,}  e^{(\nu \sigma \sqrt{t_{p}} - \mu_1 t_{p})}, &\text{
otherwise}
\end{cases} \label{eq:smlf}\\
\intertext{ and }

S_{mu} &= \begin{cases} B_{u \,}e^{-(\nu \sigma \sqrt{T} + \mu_1 T
)}, & \text{ if  $ \mu_1 \geq 0 $
or $ t_{p} \geq T$ }\\
B_{u \,}e^{-(\nu \sigma \sqrt{t_{p}} + \mu_1 t_{p} )}, & \text{
otherwise.}
\end{cases} \label{eq:smuf}
\end{align}
\end{subequations}
Although to the best of our knowledge the critical value $S_{ml}$
and $S_{mu}$ are considered for curved barriers only for the first
time in this paper, it should be noted that equations of the form
\eqref{eq:fcritic} have been used in \cite{ndogbar1} for an optimal
determination of the solution domain in the numerical pricing of
barrier options with flat barriers. It's however, also the first
time that these critical values are applied to the classification of
barrier options.

\section{Numerical test}
\label{s:test}  The choice of the value of $\nu$ in the equations
\eqref{eq:slt} and \eqref{eq:sut} is clearly crucial. Indeed, these
equations show that the larger $\nu$ is, the farther both $S_{ml}$
and $S_{mu}$ will be with respect to their corresponding barriers,
and vice versa. This means that, as expected, larger values of $\nu$
yield critical values $S_{ml}$ and $S_{mu}$ which ensure with a
better certainty that the barriers are worthless, but which are
however not optimal. Inversely, smaller values of $\nu$ will yield
critical asset prices which are wrong by being too close to the
barrier, and this is just a consequence of the fact that the
equality $\Phi (\nu)=1$ is wrong for such values of $\nu$. Ideally,
$\nu$ should be as indicated the smallest number for which the
equality $\Phi (\nu)=1$ holds. However, for any numerical
determination, the constant $\nu,$ and consequently both $S_{ml}$
and $S_{mu},$ are very sensitive to the level of accuracy required
in option prices.\par

 A numerical value for $S_{ml}$ and $S_{mu}$ can be found by trial and
error on computing systems such as {\sc mathematica}. For example in
the case of a down-and-out barrier option, $S_{ml}$ is the smallest
value of the initial price $S_0$ for which the option price is the
same as that of the vanilla option with the same parameters. We can
then use such a numerical value to find the value of $\nu$ that
yields the same value of $S_{ml}$ given by Theorem~ \ref{th:critic}.
Next, using this choice of $\nu,$ and the same accuracy in option
prices, we can check whether  numerically computed critical asset
prices match those given by Theorem~ \ref{th:critic}, for different
sets of parameters.\par

\begin{table}[htb]
\parbox[]{4.1in}{\caption{\label{tb:test1} \protect \footnotesize Numerical test
of $S_{ml}$ for a down-and-out barrier option. The fixed parameter
set is $K=100, B_l=70,$ and  $r=0.10.$}}

\begin{tabular*}{3.1in}{ c c c c c c}
\hline \\[0.3mm]
T && $ \sigma$ & \parbox[c]{1 in}{ $\quad {\rm An} (S_{ml}$) \\
with
$\nu = 4.9 $}& $ {\rm N} (S_{ml})$ & $\nu$ \\[4mm] \hline \\[0.5mm]
0.25 &&  0.15   & 98.87      & 94.852       &4.347       \\
0.25 &&  0.30   &144.00      &156.744       &5.465       \\
0.50 &&  0.15   &112.60      &112.000       &4.850       \\
0.50 &&  0.30   &192.00      &229.999       &6.129       \\ \hline
\\[1.0mm]

\end{tabular*}
\parbox[t]{3.1in}{  \protect \footnotesize For $x=S_{ml},$  ${\rm An}(x)$
is the analytical value of $x$ given by Theorem \ref{th:critic}, and
$N(x)$ is the corresponding numerical value computed by trial and
error.}
\end{table}

First, we show by an example the effect of accuracy in  option
prices on the critical asset prices and on $\nu,$ by considering the
case of a down-and-out barrier option. We assume that the parameter
set is $T=1/4, K=100, B_l= 70, \sigma = 0.30$ and $r=0.10.$ Setting
$x= S_{ml},$ we see that when accuracy in option prices is to within
$2$ digits, we have $(x, \nu)=(77.182, 0.76),$ to within $6$ digits,
$(x, \nu)= (97.000, 2.267),$ and to within $17$ digits we have $(x,
\nu)=(156.744, 5.465).$ Intuitively, this clearly shows that  higher
accuracy yields  better values for both $\nu$ and the critical
values.\par

Table \ref{tb:test1} reports the test of $S_{ml}$ for a down-and-out
barrier option with a flat barrier and with fixed parameter set
$K=100, B_l =70, r=0.10$ and for variable values of the expiry date
$T$ and volatility parameters $\sigma.$ In the table,  ${\rm An}
(x)$ is the analytical value of $x=S_{ml}$ as given by Theorem
\ref{th:critic}, while $N(x)$ is the corresponding numerical value
computed by trial and error. The third column of the table contains
analytical values of $S_{ml}$ computed with a choice of $\nu=4.9$
according to the result of Theorem~ \ref{th:critic}, and more
specifically by using Eq.~ \eqref{eq:smlf}. Column $4$ gives the
actual numerical value of $S_{ml}$ for each parameter set, while the
last column gives the value of $\nu$ for which the analytical value
of  $S_{ml}$ given by Eq.~  \eqref{eq:smlf} would agree with the
numerical value given in the fourth column.\par
  The table shows that the value of $\nu=4.9,$ which is closely the
smallest number  $y$ for which {\sc mathematica} gives $\Phi (y)=1,$
is an appropriate value for Eq.~  \eqref{eq:slt} for relatively
small volatilities such as $\sigma=0.15.$ It also shows however that
$\nu$ should be larger for larger volatilities. Unfortunately, we
don't know how exactly the choice of $\nu$ changes with the
volatility, and with other pricing parameters. Nevertheless, Theorem
 \ref{th:critic} can be used efficiently when only approximate
 values of $S_{ml}$ and $S_{mu}$ are requested~ \cite{ndogbar1}. Thanks to
 the non dimensional nature of the parameter $\nu,$ the equations
 ~\eqref{eq:slt} and \eqref{eq:sut} together with Theorem~
 \ref{th:critic} provide all the properties of the critical assets
 prices and show how they depend in particular on the parameters $T,
 \sigma,$ and $r.$

\section{Concluding Remarks}\label{s:conclusion}
In this paper, we've found with very simple arguments in Theorem~
\ref{th:critic} critical stock prices that determine when exactly a
given stock price becomes worthless with respect to some given
curves, and we've given some applications of these critical values
to barrier options. This simple determination gives fully explicit
expressions for  the critical stock prices $S_{ml}$ and $S_{mu},$
although it remains approximate due to the presence of the
approximate parameter $\nu$ in these expressions. However, the
explicit nature of this determination, as expressed for example in
the particular case of flat barriers in  \eqref{eq:fcritic} ,
provide us with simple means to understand the properties of this
critical asset prices. Indeed, Table~ \ref{tb:test1} intuitively
confirms all the effects of the expiry date $T$ and volatility
$\sigma$ on these critical stock prices.\par
   Although the comparison in Table \ref{tb:test1} between actual
numerical values and corresponding semi-analytical values of the
critical prices displays some discrepancies, such values can be
safely and efficiently used in certain cases such as the optimal
truncation of the solution domain in a numerical scheme, and
especially in the case of the difficult problem of discrete barrier
options pricing. Indeed, in the case of constant barriers, an
application of the critical asset prices given by \eqref{eq:fcritic}
to the pricing of discrete barrier options  has given  in
\cite{ndogbar1}  results closest to the only semi-analytical one
available \cite{fusai}, when compared with similar results obtained
with five other numerical methods.\par

There are ways of finding exact numerical values for $S_{ml}$ and
$S_{mu}.$ For instance, the probability $Q(S, t, T)$ that the stock
price will breach a given barrier before time $T,$ given that it has
value $S$ at time $t$ can be formulated as a first time exit
problem, and it is well-known \cite{wilmott1} that $Q(S, t, T)$
satisfies the Kolmogorov Backward equation

\begin{equation}\label{eq:kolmo} \frac{\pd Q}{\pd t} +
B(S, t) \frac{\pd Q}{\pd S} + \frac{1}{2} A(S, t) \frac{\pd ^2
Q}{\pd S^2}=0, \end{equation}

with  $A(S,t)= \sigma S$ and $B(S,t)= \mu S,$ and with some boundary
conditions that are easy to determine. For instance we have $Q(S, T,
T)=0,$ and $Q(S, t, T)=1$ on the boundary curves.  The corresponding
boundary value problem can be solved explicitly, at least by
reducing Eq.~  \eqref{eq:kolmo} to a simpler equation. But after all
such lengthy calculations the corresponding expressions for $S_{ml}$
and $S_{mu}$ are not likely to be found explicitly, although their
exact numerical values can be found in this way. Equations of the
form \eqref{eq:fcritic} therefore have the advantage of providing in
a much quicker and simpler manner, in terms of the non-dimensional
parameter $\nu,$ all the information needed for the properties of
the critical values $S_{ml}$ and $S_{mu}.$\par

By replacing the two curves in Theorem \ref{th:critic} by the
straight line determined by the strike price $K$ of a given European
option, Theorem \ref{th:critic} can also be used to determine the
initial asset price range beyond which the option becomes worthless.
There are certainly several other applications of this theorem in
option pricing and finance.


TITLE: Minimum value of $a\cdot b+2a\cdot c+b\cdot c$
QUESTION [0 upvotes]: find the minimum value of $$S=\vec{a}\cdot\vec{b}+2 \vec{a}\cdot\vec{c}+\vec{b}\cdot\vec{c}$$ if $|a|=1,|b|=2,|c|=3$

The solution is pretty easy: $$|2\vec{a} +2\vec{c}+\vec{b}|^2\ge 0$$ $$\iff 4+36+4+4S\ge 0$$ $$\iff S\ge -11$$ However if i want to solve this like this it becomes difficult
Let $p,q,r$ be the angles btw the vectors $(a,b),(b,c),(c,a)$ respectively.Then we have to minimise $$S=2\cos p+6\cos q+6\cos r$$ We cannot take $p=q=r=-\pi$ as that would mean all vectors make an angle of $-\pi$ with each other which is insane!
How would i go about solving from this method??
UpdateAs Calvin Lin pointed out My first meethod is wrong

REPLY [0 votes]: You need the additional constraint that $ p + q + r \leq 2 \pi $ (where we take the angle in the range $ [ 0, \pi ] $).
This is because the angles at the vertex/origin sum up to at most $2 \pi$.
Furthermore, because the coefficients are positive, we will always want to increase $p/q/r$, so the minimum will occur when $ p + q + r = 2 \pi$.
At this point, we can use Lagrange Multipliers to calculate the minimum.


TITLE: physical meaning of heat equation
QUESTION [1 upvotes]: consider the heat equation
$u_t=a(t)u_{xx}+f(x,t)$, $0<x<L$, $0<t<T$
subject to the initial condition
$u(x,0)=g(x)$
and boundary conditions
$u(1,t)=0,$ $u_x(0,t)+hu(0,t)=0$
where the functions $f,g$ and the constant$h$ is known.
My question is about meaning of this model in diverse scientific field. Such as in physics, chemistry, biology etc. can you give me some example about what it models in these area ? What does the boundary conditions model? and Is there any relationship between $a(t)$ and $h$?
Please help me about this problem. I really wonder the answers of these questions.
Thanks in advance!

REPLY [1 votes]: For example, it can model a process of heating a rod with time-depending heat conductivity coefficient (it can either rise or decay), for example, there is a chemical reaction in a rod, and a material (consequently, conductivity coefficient) changes with time. The boundary condition tells us that on one side of rod the temperature is constant, and on another a heat flow is constant (you might have made a mistake in a text, because in this moment your problem is ill-posed), which says us, for example, that there is a constant heat exchange with something. Here $a$ and $h$ can be both dependent and independent, because $a $ is a characteristic of rod, and $h$ is a characteristic of media, in which this rod is situated


TITLE: ODE of the form $y'=\exp(-(\cos(2\pi y))$
QUESTION [2 upvotes]: Consider the function $h:[0,1]\to \mathbb{R}$
$$h(\theta):=\sum_{k\geq 1}\frac{a_{k}}{\sqrt{k}}\cos(2\pi k \theta)+\frac{b_{k}}{\sqrt{k}}\sin(2\pi k \theta),$$
where $a_{k},b_{k}\in\mathbb{R}$. For simplicity let's take $\big(\frac{a_{k}}{\sqrt{k}}\big),\big(\frac{b_{k}}{\sqrt{k}}\big)\in \ell^{1}$, so that $h\in H^{1}(S^{1})$.
The ODE we have is
$$y'=\exp(-h(y))\;\Leftrightarrow \;y(x)=\int_{0}^{x}\exp(-h(y(t)))dt.$$
Does this fit anywhere into the ODE theory (e.g. fixed points of integral equations)? Maybe the simpler problem with a truncated series
$$y'=\exp(-h^{N}(y)),$$
where $h^{N}(\theta)=\sum_{k=1}^{N}\frac{a_{k}}{\sqrt{k}}\cos(2\pi k \theta)+\frac{b_{k}}{\sqrt{k}}\sin(2\pi k \theta)$. Or even simpler with a single term:
$$y'=\exp\Big(-\big(\frac{a_{k}}{\sqrt{k}}\cos(2\pi k y)+\frac{b_{k}}{\sqrt{k}}\sin(2\pi k y)\big)\Big).$$
Or even more simply
$$y'=\exp(-(\cos(2\pi  y)).$$
Our main question is to obtain a (series) approximation solution to $y$ (for any of the above ODEs)?

REPLY [1 votes]: It can fall under " nonlinear Volterra integral equations":
$$u(x)=\int_{0}^{x} K(x,t)F(u(t))dt,$$
and so many approaches are possible:

successive approximation
series approximation as suggested from the other answer.
Adomian Decomposition Method

See that and many other methods in the wonderful reference: "Linear and nonlinear integral equations"


TITLE: Calculating the resulting phase from its complex $x$ and $y$ field components
QUESTION [0 upvotes]: Suppose I have multiple light sources in a 2 dimension plane (X and Y), for which I only know the real and imaginary component of the electric field in the X and Y directions within the 2D plane.
This image is an example of such a case. There is 3 sources S1,s2,s3 and I only know the electric fields $\vec{Ex}=A+Bi$ and $\vec{Ey}=C+Di$ at points within the plane where A,B,C and D are constants. Suppose I wish to find out what the phase of the resulting waveform is at (x1,y1), is this possible with the available information? Calculating the phase of individual components I believe would be tan inverse of the imaginary/real component. Would the components of the wave($\vec{Ex}$ and $\vec{Ey}$) have a different phase and what would they mean? What would this phase be in relation to phase of the two components?

REPLY [0 votes]: I am going to assume that your three light sources emit coherent waves which combine to produce oscillating E field components at each point in the region, and the phase of these oscillations (relative to some chosen reference time) is represented by vectors rotating in the complex number plane.  The A,B,C, and D represent instantaneous values for the components of these vectors.  These imaginary vectors are similar to those used to represent voltages in an AC circuit. The magnitude of the oscillations of the electric field components at each point would be given by the square root of the real part of the phase vector squared plus the imaginary part squared, and the phase by the arc-tangent as you suggested. The resultant electric field vector would rotate  the xy plane. If the phase difference between the two electric field components is $90^o$, the resultant vector will sweep out an ellipse. If the two are in phase, the result is a straight line. Other phase differences give intermediate shapes.


TITLE: What is the decoding error percentage of (7,4) Hamming code?
QUESTION [3 upvotes]: The decoding error occurs where there are more than 1 bits that are flipped. Hence, if the probability that a bit is flip is $10\%$, then shouldn't the decoding error percentage be $1-0.1\cdot0.9^6\cdot7-0.9^7 = 0.149$? Why my textbook (by MacKay) claims that it is $7\%$?

REPLY [2 votes]: Let $E$ denote the event that the output of the decoder is a codeword that
is not the same as the transmitted codeword, and let $E_i$ denote the event
that the $i$-th data bit in the output codeword differs from the $i$-th data bit of the transmitted codeword. Clearly we have that
$$E = E_1 \cup E_2 \cup \cdots \cup E_k.\tag{1}$$ since if none of the $E_i$
occur, then the data bits in the output codeword are all correct and so
the output codeword must be the same as the transmitted codeword.
Now, in general, the events $E_i$ are neither independent events nor equiprobable events.  The Bit Error Rate (BER) for the decoded data bits, denoted by $P_b$ here, is defined as the arithmetic average of the $k$ $P(E_i)$ values. It is generally quite time-consuming to determine the values of the $P(E_i)$ or their average $P_b$ except for codes of very short blocklengths, whereas $P(E)$, the probability that the decoder output is not the same as the transmitted codeword, is easier to determine (or overbound for bounded-distance decoders that are prone to
decoder failure). It is straightforward to show that
$$\frac 1k P(E) \leq P_b \leq P(E)\tag{2}$$
which in the instance that is puzzling the OP tells us that the BER
must be between $3.725\%$ and $14.9\%$.  More can be said in this
particular case, however, the codewords of a $[7,4]$ Hamming code
are very closely related to what is called a biorthogonal signal set
in the communications literature (the codewords of the
$[8,4]$ extended Hamming code are exactly a biorthogonal signal set).
For a biorthogonal signal set, the BER is very nearly $\frac 12P(E)$, and so the $7\%$ BER claimed in the Mackay text can be explained
as being consistent with this result, with some differences because the $[7,4]$ Hamming code is not quite a biorthogonal signal set, and because $P_b$
is not exactly $\frac 12 P(E)$ for a biorthogonal signal set.
Masochists who really want to delve into the details can find information
about orthogonal and biorthogonal  signal sets and their relationship to first-order Reed-Muller codes (of which class the $[8,4]$ extended Hamming code is an example) can find some basic information on pp. 161-180 of this ancient Lecture Note of mine.