Find MLE and CI for $\psi = \mathbb{P}_\theta(X_1>0)$ when $X_1,...,X_n \sim \text{IID } \mathcal{N}(\theta,1)$

Question

Question: Let $X_1,\cdots,X_n \sim \text{IID }\mathcal{N}(\theta,1)$ where $\theta\in\mathbb{R}$ is unknown and let $\psi = \mathbb{P}_\theta(X_1>0)$. Find the maximum likelihood estimator $\hat{\psi}$ of $\psi$ and find a 95% confidence interval for $\psi$.

---

Comments: I think I understand the first part of the two part question. From the log-likelihood function, it can be shown that $\hat{\theta}_{MLE}=\bar{X}_n=:\sum_{i=1}^n X_i/n.$ Then, with the invariance property of the MLEs, we have $$\hat{\psi}_{MLE} = \mathbb{P}\left(\mathcal{N}(\bar{X}_n,1)>0\right).$$

My question is how do we construct the confidence interval? I am having difficulty calculating the any of the moments $\mathbb{E} \hat{\psi}_{MLE}$ and $\mathbb{E} \hat{\psi}_{MLE}^2$ and hence $\mathbb{V}(\hat{\psi}_{MLE}).$

Since $\bar{X}_n\sim \mathcal{N}(\theta,1/n)$, letting $f_{\bar{X}_n}(x)$ be the density function, I have tried using Tonelli's to switch the order of integration on the following expression $$\mathbb{E} \hat{\psi}_{MLE}=\int_{\mathbb{R}} \mathbb{P}\left(\mathcal{N}(x,1)>0\right) f_{\bar{X}_n}(x)dx=\int_{\mathbb{R}}\int_0^\infty \frac{1}{\sqrt{2\pi}}e^{-(y-x)^2/2}dy \frac{1}{\sqrt{2\pi n}}e^{-(x-\theta)^2/(2\sqrt{n})}dx,$$ but Tonelli's seems to not be the correct route. How do I obtain the variance for the confidence interval?

Edit: Further, I am interested in finding the (nondegenerate) asymptotic distribution of $\hat{\psi}_{MLE}$.

Answer 1

$\psi$ is a smooth function of $\theta$, so the delta-method says

$$\sqrt{n}(\psi(\hat\theta)-\psi(\theta))\stackrel{d}{\to} N(0, \sigma^2)$$ where $\sigma^2$ is the limiting value of $n\psi'(\theta)^2\mathrm{var}[\hat\theta]$

You know $\mathrm{var}[\hat\theta]=1/n$, so you just need $\psi'(\theta)$, which is $\phi(0-\theta)$ by the fundamental theorem of calculus (or something very close to it).

And a check

> theta<-1
>  n<-50
>  barx<-rnorm(100000,m=theta,s=1/sqrt(n))
>  psi<-function(theta) pnorm(0,theta,1)
> psi(1)
[1] 0.1586553
> mean(psi(barx))
[1] 0.1609813
> dnorm(0,1,1)^2
[1] 0.05854983
> var(psi(barx))*50
[1] 0.05881804