\documentclass[notitlepage]{article}
\usepackage{Math}
%\usepackage{../../Math598}
\usepackage{listings}
\usepackage{numprint}
\usepackage{enumerate}
\usepackage{bm}
\usepackage{bbm}
\usepackage{verbatim}
\usepackage{mathtools}
\usepackage{dsfont}
\usepackage{scalerel}

\usepackage{amsfonts,amsmath,amssymb}

\usepackage{tikz}
\usetikzlibrary{shapes,decorations,arrows,calc,arrows.meta,fit,positioning}


\usepackage{pgfplots}
\pgfplotsset{compat=1.17}

\newcommand*\circled[1]{\tikz[baseline=(char.base)]{
            \node[shape=circle,draw,inner sep=2pt] (char) {#1};}}

\usepackage{xcolor}

\newtheorem{alg}{Algorithm}
\newtheorem{ex}{Example}
\newtheorem{note}{Note}
\newenvironment{proof}{\par\noindent{\normalfont{\bf{Proof
}}}}
{\hfill \small$\blacksquare$\\[2mm]}


\setlength{\parindent}{0in}

\def\E{\mathbbmss{E}}
\def\Var{\text{Var}}
\def\Cov{\text{Cov}}
\def\Corr{\text{Corr}}

\def\bW{\mathbf{W}}
\def\bH{\mathbf{H}}

\newcommand{\bs}[1]{\bm{#1}}
\newcommand{\Rho}{\mathrm{P}}
\newcommand{\cif}{\text{ if }}
\newcommand{\Ra}{\Longrightarrow}

\def\bphi{\boldsymbol \phi}


\lstloadlanguages{R}

\definecolor{keywordcolor}{rgb}{0,0.6,0.6}
\definecolor{delimcolor}{rgb}{0.461,0.039,0.102}
\definecolor{Rcommentcolor}{rgb}{0.101,0.043,0.432}

\lstdefinestyle{Rsettings}{
  basicstyle=\ttfamily,
  breaklines=true,
  showstringspaces=false,
  keywords={if, else, function, theFunction, tmp}, % Write as many keywords
  otherkeywords={},
  commentstyle=\itshape\color{Rcommentcolor},
  keywordstyle=\color{keywordcolor},
  moredelim=[s][\color{delimcolor}]{"}{"},
}

\lstset{basicstyle=\ttfamily, numbers=none, literate={~} {$\sim$}{2}}

\def\Xb{\mathbf{X}_{\beta}}
\def\Xp{\mathbf{X}_{\psi}}
\def\beps{\boldsymbol{\epsilon}}
\def\betahat{\widehat \beta}
\def\psihat{\widehat \psi}

\begin{document}

<<setup, cache=FALSE, echo=FALSE>>=
library(knitr)
# global chunk options
opts_chunk$set(cache=TRUE, autodep=TRUE)
options(scipen=6)
options(repos=c(CRAN="https://cloud.r-project.org/"))
knit_hooks$set(document = function(x) {sub('\\usepackage[]{color}', '\\usepackage{xcolor}', x, fixed = TRUE)})
@
\begin{center}
{\textsclarge{Project 1: Solutions}}
\end{center}

\tikzset{
    -Latex,auto,node distance =1 cm and 1 cm,semithick,
    state/.style ={circle, draw, minimum width = 0.7 cm},
    box/.style ={rectangle, draw, minimum width = 0.7 cm, fill=lightgray},
    point/.style = {circle, draw, inner sep=0.08cm,fill,node contents={}},
    bidirected/.style={Latex-Latex,dashed},
    el/.style = {inner sep=3pt, align=left, sloped}

}

We use the DAG 
\begin{figure}[ht]
\centering
\begin{tikzpicture}[scale=1.5]
    % x node set with absolute coordinates
    \node[state] (x1) at (0,2) {${X_1}$};
    \node[state] (x2) at (-1,1){${X_2}$};
    \node[state] (x3) at (1,1) {${X_3}$};
    \node[state] (x4) at (0,0) {${X_4}$};
    \node[state] (x5) at (0,-1) {${X_5}$};

    \path (x1) edge (x2);
    \path (x1) edge (x3);
    \path (x2) edge (x4);
    \path (x3) edge (x4);
    \path (x4) edge (x5);

 \end{tikzpicture}
\end{figure}

which implies the factorization
\[
f_{X_1}(x_1) f_{X_2|X_1}(x_2|x_1) f_{X_3|X_1}(x_3|x_1) f_{X_4|X_2,X_3}(x_4|x_2,x_3) f_{X_5|X_4}(x_5|x_4)
\]
with
\begin{itemize}
  \item $X_1 \sim Normal(1,1)$
  \item $X_2 |X_1 = x_1 \sim Normal(x_1,1)$ and $X_3 |X_1 = x_1 \sim Normal(2 x_1,1)$ independent
  \item $X_4 | X_2=x_2, X_3 = x_3 \sim Normal(x_3 - x_2, 1)$
  \item $X_5 |X_4 = x_4 \sim Normal(x_4,1)$
\end{itemize}
yielding a jointly multivariate Normal distribution.


<<C1,comment='+', fig.width=6, fig.height=6, fig.align='center'>>=
#Calculation for large N
library(mvnfast)
set.seed(22087)
n<-10000
X1<-rnorm(n,1,1)
X2<-rnorm(n,X1,1)
X3<-rnorm(n,2*X1,1)
X4<-rnorm(n,X3-X2,1)
X5<-rnorm(n,X4,1)
X<-cbind(X1,X2,X3,X4,X5)
cor(X)                    #Correlation
par(mar=c(2,2,2,2))
pairs(X,pch=19,cex=0.6)
Xa<-X
@

\begin{enumerate}[(a)]

\item \textbf{Blocking and unblocking via conditioning:}  By construction, $X_2$ and $X_3$ are conditionally independent given $X_1$, but are marginally dependent .  We can demonstrate this by regression, or by indpecting the correlation or partial correlation between the variables: first, for dependence:
<<C2,comment='+', fig.width=6, fig.height=6, fig.align='center'>>=
cor(X2,X3)                  #Correlation not zero
coef(summary(lm(X3~X2)))    #Regression coefficient not zero
coef(summary(lm(X2~X3)))    #Regression coefficient not zero
@
Secondly, for conditional independence, we can inspect the partial correlation of $(X_2,X_3)$ given $X_1$ by looking at the $(2,3)$ element of the inverse correlation matrix formed by considering the three variables. 
<<C3,comment='+', fig.width=6, fig.height=6, fig.align='center'>>=
solve(cor(cbind(X1,X2,X3))) #almost zero in the (2,3) element
coef(summary(lm(X3~X2+X1))) #conditioning on X1 makes the coefficient for X2 almost zero.
@

\item \textbf{d-separation and conditional independence:}  From the DAG, we can see that
\[
X_1 \independent X_4 \given (X_2,X_3)
\]
so that $X_1$ and $X_4$ are d-separated by $S = (X_2,X_3)$.
<<C4,comment='+', fig.width=6, fig.height=6, fig.align='center'>>=
solve(cor(cbind(X1,X2,X3,X4))) #almost zero in the (1,4) element
coef(summary(lm(X1~X4+X2+X3))) #conditioning on X2,X3 makes the coefficient for X4 zero.
coef(summary(lm(X4~X1+X2+X3))) #conditioning on X2,X3 makes the coefficient for X1 zero.
@
\item \textbf{the impact of intervention:} for the intervention case we consider intervening to set $X_3$ deterministically to its marginal mean value, that is, $X_3 = 2$ with probability 1.
<<C5,comment='+', fig.width=6, fig.height=6, fig.align='center'>>=
#Calculation for large N
library(mvnfast)
set.seed(22087)
n<-10000
X1<-rnorm(n,1,1)
X2<-rnorm(n,X1,1)
X3<-rep(2,n)
X4<-rnorm(n,X3-X2,1)
X5<-rnorm(n,X4,1)
X<-cbind(X1,X2,X4,X5)
cor(X)                    #Correlation
par(mar=c(2,2,2,2))
pairs(X,pch=19,cex=0.6)
Xb<-X
@
Note that under the intervention, the joint distribution of $(X_2,X_4,X_5)$ is now different; in the original DAG, $(X_2,X_4)$ and $(X_2,X_5)$ are marginally uncorrelated

<<C6,comment='+', fig.width=6, fig.height=6, fig.align='center'>>=
cor(Xa[,c(2,4)])
cor(Xa[,c(2,5)])
@
whereas under the intervention, they are not.
<<C7,comment='+', fig.width=6, fig.height=6, fig.align='center'>>=
cor(Xb[,c(2,3)])
cor(Xb[,c(2,4)])
@
This is due to the specific settings of the parameters used in the simulation; in general the correlation will change, but the lack of correlation in the original simulation does not hold in general, as the next simulation indicates: for example if we change the model to
\begin{itemize}
  \item $X_1 \sim Normal(1,1)$
  \item $X_2 |X_1 = x_1 \sim Normal(x_1,1)$ and $X_3 |X_1 = x_1 \sim Normal(x_1,1)$ independent
  \item $X_4 | X_2=x_2, X_3 = x_3 \sim Normal(x_3 - x_2, 1)$
  \item $X_5 |X_4 = x_4 \sim Normal(x_4,1)$
\end{itemize}
then the correlation pattern changes.

<<C8,comment='+', fig.width=6, fig.height=6, fig.align='center'>>=
#Calculation for large N
library(mvnfast)
set.seed(22087)
n<-10000
X1<-rnorm(n,1,1)
X2<-rnorm(n,X1,1)
X3<-rnorm(n,X1,1)
X4<-rnorm(n,X3-X2,1)
X5<-rnorm(n,X4,1)
X<-cbind(X1,X2,X3,X4,X5)
cor(X)                    #Correlation
par(mar=c(2,2,2,2))
pairs(X,pch=19,cex=0.6)
@
\end{enumerate}

\end{document}