# Is there mistake in derivation/simulation of Maximum A Posterori (MAP) estimation?

I'm trying to estimate $$X$$ using MAP Estimator. If $$Y = aX + N$$, where $$a$$ is a constant, $$X \sim \mathcal{N}(0,\sigma_{x}^{2})$$ and $$N \sim \mathcal{N}(0,\sigma_{n}^{2})$$.

The likelihood function for $$Y = aX + N$$ can be evaluated as $$$$\label{llf} f_{Y|X}(y|x)=\frac{1}{\sqrt{2\pi}\sigma^2}\exp\Bigg\{-\frac{1}{2}\frac{(y-\mu)^2}{\sigma^2}\Bigg\}$$$$ where $$\mu$$ and $$\sigma^2$$ are mean and variance of $$f_{Y|X}(y|x)$$, which can be evaluated as \begin{align}\label{mean} \mu&=\mathbb{E}[Y\vert X] \notag\\ &= \mathbb{E}[aX + N\vert X] \notag\\ &= a\mathbb{E}[X\vert X]+\mathbb{E}[N\vert X] \notag\\ &=ax \end{align} and the variance can be evaluated as \begin{align}\label{var} \sigma^2 &= \mathbb{E}[(Y - \mu)^2\vert X]\notag\\ &= \mathbb{E}[Y^2\vert X] - \mu^2 \notag\\ &= \mathbb{E}[(aX + N)^2\vert X] - a^2x^2\notag\\ &= a^2\mathbb{E}[X^2\vert X]+\mathbb{E}[N^2\vert X] + 2a\mathbb{E}[X N\vert X] - a^2x^2 \notag\\ &= a^2x^2 + \sigma_{n}^2 + 2a\mathbb{E}[X\vert X]\mathbb{E}[N\vert X]-a^2x^2\notag\\ &= \sigma_{n}^2 \end{align} And a priori density function is $$$$\label{ppdf} f_{X}(x)=\frac{1}{\sqrt{2\pi}\sigma_{x}^2}\exp\Bigg\{-\frac{x^2}{2\sigma_{x}^2}\Bigg\}$$$$ Therefore the Bayesian log-likelihood function will be \begin{align} \mathcal{L}(X,Y) &= \ln f_{Y\vert X}(y\vert x)f_{X}(x)\notag\\ &=\ln\frac{1}{2\pi\sigma_{n}^2\sigma_{x}^2}\exp\Bigg\{-\frac{1}{2}\Bigg[\frac{(y-ax)^2}{\sigma_{n}^2} + \frac{x^2}{\sigma_{x}^2}\Bigg]\Bigg\}\notag\\ &= \ln\frac{1}{2\pi\sigma_{n}^2\sigma_{x}^2} -\frac{1}{2}\Bigg[\frac{(y-ax)^2}{\sigma_{n}^2} + \frac{x^2}{\sigma_{x}^2}\Bigg] \end{align} On differentiating above equation with respect to $$x$$, we get \begin{align} \frac{\partial}{\partial x}\mathcal{L}(X,Y)=\frac{\partial}{\partial x}\ln\frac{1}{2\pi\sigma_{n}^2\sigma_{x}^2} -\frac{1}{2}\frac{\partial}{\partial x}\Bigg[\frac{(y-ax)^2}{\sigma_{n}^2} + \frac{x^2}{\sigma_{x}^2}\Bigg] \end{align} And now equating the above equation for $$x$$ equal to zero. \begin{align} \frac{1}{2\sigma_{n}^2}\frac{\partial}{\partial x}(y-ax)^2 - \frac{1}{2\sigma_{x}^2}\frac{\partial}{\partial x}x^2 &= 0\notag\\ -\frac{(a^2x - ay)}{\sigma_{n}^2} - \frac{x}{\sigma_{x}^2}&=0\notag\\ x\Bigg(\frac{a^2\sigma_{x}^2 + \sigma_{n}^2}{\sigma_{x}^2\sigma_{n}^2}\Bigg) &= \frac{ay}{\sigma_{n}^2}\notag\\ x &= \frac{a\sigma_{x}^2y}{a^2\sigma_{x}^2 + \sigma_{n}^2} \end{align} Hence $$$$\hat{X}_{MAP} = \frac{a\sigma_{x}^2y}{a^2\sigma_{x}^2 + \sigma_{n}^2}$$$$

Below is the MATLAB code which I'm using to simulate and validate the theoretical result obtained.

%% Release memory and clear screen
clear;clc;

%% Theoretical
varX = 3;                   % Variance of X
varN = 1;                   % Variance of N
a = 2;                      % Constant
L = 1e3;                    % Number of points
x = linspace(-10,10,L);

%% Generate Gaussianly Distributed Random Number
X = varX*randn;
N = varN*randn;
Y = a*X + N;
% MAP Estimator of X
Xhatsim = a*varX*Y/(a*a*varX + varN);
mu = a*X;
fx = normpdf(x,0,sqrt(varX));     % Density of X
fn = normpdf(x,0,sqrt(varN));     % Density of N
fyx = normpdf(x,mu,sqrt(varN));   % Density of Y|X
fxy = fyx.*fx;                    % Density of X|Y
[Xhattheo,loc] = max(fxy);

%% Plot Results
plot(x,fx,'r',x,fn,'b',x,fyx,'k',x,fxy,'m',...
x(loc),Xhattheo,'o',x(loc),Xhatsim,'*');
axis([min(x),max(x),...
min(min(min(min(fx,fn),fyx),fxy)),...
max(max(max(max(fx,fn),fyx),fxy))]);
str = {'$$f_X(x)$$','$$f_N(n)$$','$$f_{Y\vert X}(y\vert x)$$',...
'$$f_{X\vert Y}(x\vert y)$$','Theoretical','Simulated'};
legend(str,'Interpreter', 'latex')

But as it can be seen from the result shown in the figure doesn't match. Where am I doing a mistake? Is it in the derivation or in Simulation code?