aSummary of statCompwritten by Haodong with MathThis is a brief summary for DATA130004, statistcal computing in the school of data science, Fudan University, fall 2021. The summary contains only the logic flow and the esssential parts of the course, and the "importance" is judged by how well I remember. Some details might be ignored. You are suppose to refer to this summary for previewing or reviewing for the course or other similar fields. You are not supposed to use it as text book or lecture note record word by word.The course focus on applying statistical methods on computer, especially for R. 1. generate random variables 1.1. inverse transform 1.2. acceptance-rejection method 2. Monte Carlo integration and variance reduction 2.1. Simple MC integration 2.2. Variance Reduction 2.2.1. antithetic variable 2.2.2. control variate 2.2.3. Antethetic variables as a control variate 2.2.4. Control variate and linear regression 2.2.5. importance sampling 2.2.6. stratified sampling 2.2.7. Stratified importance sampling 3. MC in statistical inference 3.1. point estimation 3.2. confidence interval 3.3. hypothesis testing 3.3.1. empirical type I error rate 3.3.2. Power of a test 4. Bootstrap 4.1. bootstrap estimate of distribution 4.2. point estimation 4.3. confidence interval 4.3.1. standard normal distribution 4.3.2. percentile CI 4.3.3. Basic bootstrap CI 4.3.4. Bootstrap t CI 5. Jacknife 5.1. Bias 5.2. Standard Error 6. Bayesian statistics and MCMC 6.1. Bayesian problem set-up 6.2. Markov chain Monte Carlo 6.2.1. Metropolis-Hasting sampler 6.2.2. Metropolis sampler 6.2.3. Random Walk sampler 6.2.4. independent sampler 6.2.5. Gibbs sampler 6.3. Monitoring the convergence (Gelman-Rubin method) 7. EM algorithm 8. Variational inference 8.1. KL divergence 8.2. Evidence lower bound (ELBO) 8.3. The mean-field variational family 1. generate random variablesIn this part, we talk about the methods for generating random numbers. For some certain distributions, like normal distribution, we need to think of how could we use computer to generate numbers from the distributuion we want. This is something you should bare in mind, one of the key idea for this course is to think from the view of a computer programmer.1.1. inverse transform You might learned from your probability course one interesting theorem, which you don't know what for in a long time, that is,

Theorem1Probability Integral TransformationIf X is a continuous random variables with cdf FX(x), then U=FX(X)∼Uniform(0,1).

Thus, given a uniform random number generator and a cdf, we can get the random number using the inverse of the cdf.1. Derive the inverse function F-1X(u).2. Generate u∼Unif(0,1).3. x=F-1(u)This method requires the F-1 easy to compute.1.2. acceptance-rejection methodIn this method, given a generator with pdf Y∼g(t), we can generate from our target X∼f(t). We propose the algorithm, then we prove it.First we assume that c satisfies f(t)

g(t)≤c for all t∈R, then 1. Generate y∼g(t), u∼Unif(0,1).2. If u<f(y)

cg(y) then accept y and set x=yNow we prove x is actually from target f.a

P(X=y|accept)	=P(Y=y|accept)
	=P(accept|Y=y)P(Y=y) P(accept)

whereP(accept|Y=y)=P(u<f(y)

cg(y)|Y=y)=f(y)

cg(y)P(accept)= ∑yP(accept|Y=y)P(Y=y)= ∑yf(y)

cg(y)⋅g(y)=1

cThus P(X=y|accept)=f(y)

cg(y)⋅g(y)

1

c=f(y)The continuous case is shown in the homework.2. Monte Carlo integration and variance reductionIn this part, we solve the integration problem 𝜃=∫g(x)dx. For x1,⋯,xmiid

∼

X, the empirical averge EM[g(x)]=1

m∑g(xi) estimates the population mean, i.e., the expectaion, Eg(x)=∫g(x)f(x)dx, X∼f. Thus we introduce the Mote Carlo integration. The Monte Carlo method is a widely used method. When saying MC, we are simulating many times to estimate the target. For example, in this part, for estimate the integral, we simulate many random numbers to approximate the result.2.1. Simple MC integrationTo estimate 𝜃=1∫0g(x)dx, X∼Unif(0,1).1. Generate x1,⋯,xmiid

∼

X2. 𝜃=1

mm∑i=1g(xi).Notice that the integration is limited to the range 0 to 1. We can generalize to range [a,b] by change of variable y=x-a

b-a, or just generate X∼Unif(a,b), and 𝜃=(b-a)1

mm∑i=1g(xi).2.2. Variance ReductionAlthough the simple MC integration is unbiased eastimator of 𝜃, we can use better estimators with smaller variance.2.2.1. antithetic variableIn simple MC, we use independent random variables. However, with usage of dependent variables, we might reduce the variables.Suppose Y and Z have same distribution with X but dependent. Var(1

2(Y+Z))=1

4aVar(Y)+Var(Z)+2Cov(Y,Z)a, if Cov(Y,Z)<0, the variance can be reduced.Known that if U∼Unif(0,1), then 1-U∼Unif(0,1) and U and 1-U are negatively correlated. And we can expect that

Corollary1 If g(X)=g(X1,⋯,Xn) is monotone, then Y=g(F-1X(u1),⋯,F-1X(un)) and Y=g(F-1X(1-u1),⋯,F-1X(1-un)) are negatively correlated.

Proof is ignored here.Then, instead of generating m Unif(0,1) random varables, we needs m

2 generations, and for j=1,⋯,m

2, we define Yj=g(F-1X(u(j)1),⋯,F-1X(u(j)n)) and Yj'=g(F-1X(1-u(j)1),⋯,F-1X(1-u(j)n)), then 𝜃=1

mm/2∑i=1(Yj+Yj').2.2.2. control variate In this part, we still try to use the benefits of correlation. Suppose there is a f with 𝜇=Ef(x) known and f is correlated with g.Then define 𝜃c=g(x)+c(f(x)-𝜇). 𝜃c is still an unbiased estimator of 𝜃 and Var(𝜃c)=Var(g(x))+2c⋅Cov(g(x),f(x))+c2Var(f(x))Let c*=-Cov(g(x),f(x))

Var(f(x)), we minimize the Var(𝜃c*)=Var(g(x))-Cov2(g(x),f(x))

Var(f(x)).The percent of reduction is Var(g(x))-Var(𝜃c*)

Var(g(x))=Cor2(g(x),f(x))×100%.2.2.3. Antethetic variables as a control variateWe combine the two methods, formulate control variate as linear combination of two unbiased estimator 𝜃c=c𝜃1+(1-c)𝜃2. Suppose 𝜃1 and 𝜃2 have identical distributions and r=Cor(𝜃1, 𝜃2)<0. Then Var(𝜃c)=c2Var(𝜃1)+2c(1-c)Cov(𝜃1, 𝜃2)+(1-c)2Var(𝜃2)=Var(𝜃1)ac2+2c(1-c)r+(1-c)2a=Var(𝜃1)a(2-2r)c2-(2-2r)c+1ac*=1

2.2.2.4. Control variate and linear regressionIn control variate method, suppose we have n samples (f(x1),g(x1)),⋯,(f(xn),g(xn)). When applying the linear regression g(x)=𝛼+𝛽f(x)+𝜀, we have the following four important properties.1. the OLS estimatorsa

𝛼=⏨⏨g(x)-𝛽⏨⏨f(x)
𝛽=Cov(f(x),g(x)) Var(f(x))=-c*

2. The predicted value at 𝜇=Ef(x) is the control variate estimator of the target integration𝛼+𝛽𝜇=⏨⏨g(x)-𝛽(⏨⏨f(x)-𝜇)=⏨⏨g(x)+c*(⏨⏨f(x)-𝜇)=𝜃c*3. The variance of the control variate estimator is the residual mean squared error (MSE).Var(⏨⏨g(x)+c(⏨⏨f(x)-𝜇))=1

nVar(g(x)+c(f(x)-𝜇))=1

nVar(g(x)-𝛽f(x)-𝛼)=𝜎2𝜀

n4. The percetage of improvement is {Cor(g(x),f(x))}×100% is the coefficient of determintion.2.2.5. importance samplingSimple MC integration, b-a

mm∑i=1g(Xi), weight the interval [a,b] uniformly. The replicates X1,⋯,Xm are uniformly distributed on [a,b]. Then we consider other weight functions.

Algorithm1Importance sampling1. decide a "envolope" f(x)2. For i=1,⋯,m(a) generate xi∼f(b) record g(xi)

f(xi)3. 𝜃IS=1

mm∑i=1g(xi)

f(xi)

Let X be a r.v. with density f such that f(x)>0 on the support of g. Set Y=g(x)

f(x) then 𝜃=∫g(x)dx=∫g(x)

f(x)f(x)dx=Ef[g(x)

f(x)]Thus we can use 𝜃IS=1

mm∑i=1g(xi)

f(xi) to estimate 𝜃.Now we analysis the variance.Var(𝜃IS)=Var(1

mm∑i=1Yi)=1

mVar(Y)=1

mVar(g(x)

f(x))=1

m{E(g(x)

f(x))2-(Eg(x)

f(x))2}where Eg(x)

f(x)=∫g(x)

f(x)f(x)dx=𝜃 and E(g(x)

f(x))2=∫g2(x)

f(x)dx.Var(𝜃IS)=1

m{∫g2(x)

f(x)dx-𝜃2}where ∫g2(x)

f(x)dx={∫g2(x)

f(x)dx}{∫f(x)dx}≥(∫g(x)dx)2. The equation holds iff f(x)∝|g(x)|.2.2.6. stratified sampling

Algorithm2Stratified Sampling1. divide [0,1] into k strata where the j-th strata is Ij=(j-1

k,j

k)2. On each strata, for i=1,⋯mj(a) generate x(j)i∼Unif(Ij) by density fj(x)=k⋅1(x∈Ij)(b) 𝜃j=1

mjmj∑i=1g(x(j)i)3. 𝜃S=1

kk∑j=1𝜃j

Note that E𝜃j=Eg(x(j))=∫g(x)k⋅1(x∈Ij)dx=k∫Ijg(x)dx. That's why we need 1

k.Now we show that the variance of importance sampling is smaller than simple MC integration.Denote 𝜃M be the simple MC estimation. For simplicity, we suppose in importance sampling each strata has equal numer, m, of replicates, and total number M=mk. Denote 𝜃j=E{g(u)|u∈Ij} and 𝜎2j=Var{g(u)|u∈Ij}.Var(𝜃S)=Var(1

kk∑j=1𝜃j)=1

k2k∑j=1Var(𝜃j)=1

k2k∑j=1𝜎2j

m=1

Mkk∑j=1𝜎2jSuppose a two-step experiment, J is discrete uniform on {1,⋯,k}. For i=1,⋯,M,1. draw J2. generate from IJVar(𝜃M)=1

MVar(g(u))=1

MaVaraE(g(u)aJ)a+EaVar(g(u)aJ)aa=1

MaVar(𝜃J)+1

kk∑j=1𝜎2ja=1

Mkk∑j=1𝜎2j+1

MVar(𝜃J)≥Var(𝜃S)The equation holds iff Var(𝜃J)=0, i.e., 𝜃1=⋯=𝜃k.2.2.7. Stratified importance sampling

Algorithm3Stratified importance sampling1. choose an importance function f.2. divide the real line into k strata where the j-th strata is Ij=(aj-1,aj), where a0=-∞, aj=F-1(j

k), ak=+∞.3. On strata j define gj(x)=g(x)1(x∈Ij) and fj(x)=fx|Ij(x|Ij)=f(x)1(x∈Ij)

∫f(x)1(x∈Ij)dx=kf(x)1(x∈Ij), for i=1,⋯m(a) generate x(j)i∼fj(x)(b) 𝜃j=1

mm∑i=1g(x(j)i)4. 𝜃SIS=k∑j=1𝜃j

Note that 𝜃j=∫gj(x)dx=∫Ijg(x)dx, thus 𝜃=∑𝜃j.Now we prove that Var(𝜃SIS)≤Var(𝜃IS). In Ij, denote 𝜃j=∫Ijg(x)dx=E{gj(x)

fj(x)} and 𝜎2j=Var{gj(x)

fj(x)}, where x∼fj.Var(𝜃SIS)=Var(∑𝜃j)=∑Var(𝜃j)=∑𝜎2j

m=k

M∑𝜎2jVar(𝜃IS)=1

MVar(Y)=𝜎2

M where 𝜎2=Var(g(x)

f(x))Next we show that 𝜎2-k∑𝜎2j≥0. We consider the two-stage experiment, For J=j, we geberate x* from fj and set Y*=gj(x*)

fj(x*)=g(x*)

kf(x*). x*d

=

x and kY*d

=

Y. Var(Y*)=E{Var(Y*|J)}+Var{E(Y*|J)}where E{Var(Y*|J)}=E(𝜎2J)=1

k∑𝜎2jVar{E(Y*|J)}=Var(𝜃J)Then 𝜎2=Var(Y)=k2Var(Y*)=k2Var(𝜃J)+k∑𝜎2j≥k∑𝜎2j3. MC in statistical inference 3.1. point estimation𝜃=1

m∑𝜃(j)se(⏨x)=a

1 n{∑(xi-⏨x)2}1 2
1 n{1 n-1∑(xi-⏨x)2}1 2

mse=1

m∑(𝜃(j)-𝜃)23.2. confidence interval

Algorithm4Monte Carlo Confindence interval1. For each replicate j=1,⋯,m(a) generate the j-th random sample X(j)1,⋯,X(j)n(b) compute the confidence interval Cj for the j-th sample(c) Compute yj=1(𝜃∈Cj) for the j-the sample2. Compute the empirical confidence level ⏨y=1

m∑yj

3.3. hypothesis testing3.3.1. empirical type I error rate

Algorithm5MC type I error rate1. For each replicate, indexed by j=1,⋯,m(a) Generate the j-th random sample x(j)1,⋯,x(j)n from the null distribution. (b) Compute the test statistic Tj from the j-th sample.(c) Record the test decision Ij=1 if H0 is rejected at significance level 𝛼 and otherwise Ij=0.2. Compute the proportion of significant tests 1

m∑Ij. This proportion is the observed Type I error rate.

3.3.2. Power of a test

Algorithm6MC power of a test1. select a particular 𝜃1∈𝛩12. For each replicate, indexed by j=1,⋯,m(a) Generate the j-th random sample x(j)1,⋯,x(j)n under 𝜃1. (b) Compute the test statistic Tj from the j-th sample.(c) Record the test decision Ij=1 if H0 is rejected at significance level 𝛼 and otherwise Ij=0.3. Compute the proportion of significant tests 𝜋(𝜃1)=1

m∑Ij.

4. Bootstrap4.1. bootstrap estimate of distribution

Algorithm7bootstrap estimate of distribution1. Foe each Bootstrap replicate b=1,⋯,B(a) Generate samle x*(b)=(x*1(b),⋯,x*m(b)) by sample with replacement from the observation x1, ⋯,xn(b) Compute the b-th replicate 𝜃(b) using x*(b).2. The bootstrap estimate of F𝜃 is the empirical distribution of 𝜃(1),⋯,𝜃(B)

4.2. point estimation1. se of 𝜃. seB≜se(𝜃*)=1

B-1B∑b=1(𝜃(b)-⏨𝜃*)2 where ⏨𝜃*=1

BB∑b=1𝜃(b).2. bias of 𝜃. bias(𝜃)=1

BB∑b=1𝜃(b)-𝜃.4.3. confidence intervalNow we use Bootstrap to estimate confidence intervals.4.3.1. standard normal distributionUse approximately normal poperty, [𝜃±1.96seB(𝜃)]4.3.2. percentile CIUse the sample quantile [𝜃*𝛼/2,𝜃*1-𝛼/2].4.3.3. Basic bootstrap CISuppose (L,U) is the confidence interval, i.e., P(L≥𝜃)=P(U≤𝜃)=𝛼

2. Then 𝛼

2=P(L≥𝜃)=P(L-𝜃≥𝜃-𝜃)=P(𝜃-𝜃≥𝜃-L)Thus L-𝜃 is the 1-𝛼

2 percentile of 𝜃-𝜃. We can estimate the 1-𝛼

2 quantiles of 𝜃 using bootstrap replicate 𝜃*1-𝛼/2, then 𝜃*1-𝛼/2-𝜃 is approximately equal to the 1-𝛼

2 quantile of 𝜃-𝜃. Set 𝜃-L=𝜃*1-𝛼/2-𝜃, we have L=2𝜃-𝜃*1-𝛼/2.Similarly, we have U=2𝜃-𝜃*𝛼/2. Thus the CI is [2𝜃-𝜃*1-𝛼/2,2𝜃-𝜃*𝛼/2].4.3.4. Bootstrap t CI

Algorithm8Bootstrap t confidence interval1. Compute 𝜃 from the observed data.2. For each bootstrap replicate b=1,⋯,B(a) Sample with replacement x(b)=(x(b)1,⋯,x(b)n)(b) compute 𝜃(b) from x(b)(c) estimate se(𝜃(b)). (another layer of Bootstrap, resample from x(b), not from x)(d) compute the b-th replicate of the t* distribution t(b)=𝜃(b)-𝜃

se(𝜃(b)).3. find sample quantiles t*𝛼

2 and t*1-𝛼

2 from {t(b)}Bb=1.4. compute se(𝜃) from {𝜃(b)}Bb=1.5. confidence interval [𝜃-t*1-𝛼

2se(𝜃),𝜃-t*𝛼

2se(𝜃)]

5. Jacknife Jacknife is similar to the leave-one-out method. In each sample, the Jacknife except one obvervation, i.e., the i-th Jacknife sample is x(i)=(x1,⋯,xi-1,xi+1,⋯,xn). The i-th Jacknife estimator is 𝜃(i)=f(x(i)).5.1. BiasThe jacknife bias is biasjack=(n-1)(⏨𝜃(⋅)-𝜃)where ⏨𝜃(⋅)=1

nn∑i=1𝜃(i).why n-1? For example, 𝜃=1

nn∑i=1(xi-⏨x)2, bias(𝜃)=E(𝜃)-𝜎2=n-1

n𝜎2-𝜎2=-1

n𝜎2. whileE(𝜃(i)-𝜃)=E(𝜃(i)-𝜃)-E(𝜃-𝜃)=-1

n-1𝜎2+1

n𝜎2=-𝜎2

n(n-1)=bias(𝜃)

n-15.2. Standard ErrorThe jacknife standard error is sejack=n-1

nn∑i=1(𝜃(i)-⏨𝜃(⋅))2Why n-1

n? For example, if 𝜃=⏨x, then 𝜃(i)=n⏨x-xi

n-1, ⏨𝜃(⋅)=⏨x.n∑i=1(𝜃(i)-⏨𝜃(⋅))2=n∑i=1(n⏨x-xi

n-1-⏨x)2=1

(n-1)2n∑i=1(⏨x-xi)2=1

n-1s26. Bayesian statistics and MCMCBayesian statistics looks problems in a differen way from frenquetists. In frequentist, the parameter is fixed but unknown, and he experiment is repeatable. However, in Bayesian, the parameter is random and unknown, and the experiment is fixed, or to say, the target of Bayesian is the posteriorp(𝜃|data), rather than the likelihood p(data|𝜃). 6.1. Bayesian problem set-up

parameter	𝜃
data	X
prior distribution of 𝜃	𝜋(𝜃)
sampling model for X	f(x|𝜃)
posterior distribution for 𝜃|X	p(𝜃|X)=𝜋(𝜃)f(X|𝜃) ∫𝛩𝜋(𝜃)f(X|𝜃)d𝜃∝prior×likelihood

6.2. Markov chain Monte CarloConstruct a Markov chain {Xt:t=0,1,⋯} whose stationary distribution is target distribution.6.2.1. Metropolis-Hasting samplerNow supppose our target distribution is f, to move from one state to another, we have a proposal Y∼g(y|Xt).

Algorithm9MH sampler1. Choose a proper proposal g(⋅|Xt)2. Initialize X0 and repeat until the chain converges. At time t,(a) Generate Y from g(⋅|Xt)(b) Compute the acceptance rate r(Xt,Y)=f(Y)g(Xt|Y)

f(Xt)g(Y|Xt)(c) Let 𝛼(Xt, Y)=min{r(Xt,Y),1} be the acceptance ratio. Set Xt+1=a

Y		with probability 𝛼(Xt,Y)
Xt		with probability 1-𝛼(Xt,Y)

Now we show how the MH sampler works.Let s, t be two different state. Without loss of generality, we assume that f(s)g(r|s)≥f(r)g(s|r). Then 𝛼(r,s)=1 and 𝛼(s,r)=f(r)g(s|r)

f(s)g(r|s), and a

P(Xt=s,Xt+1=r)	=P(Xt=s)P(Xt+1=r|Xt=s)
	=f(s)g(r|s)𝛼(s,r)
	=f(r)g(s|r)
P(Xt=r,Xt+1=s)	=P(Xt=r)P(Xt+1=s|Xt=r)
	=f(r)g(s|r)𝛼(r,s)
	=f(r)g(s|r)

⟹ P(Xt=s,Xt+1=r)=P(Xt=r,Xt+1=s)Thus a

P(Xt=r)	= ∑sP(Xt=r,Xt+1=s)
	= ∑sP(Xt=s,Xt+1=r)
	=P(Xt+1=r)

The Markov chain is stationary at f.□We can also prove it in a kernel point of view. Define K(r,s)=g(s|r)𝛼(r,s) be the transition kernel from r to s. We call f is the stationary distriibution if the balance condition holds ∫f(r)K(r,s)dr=f(s)Under MH sampler, a

f(r)K(r,s)	=f(r)g(s|r)𝛼(r,s)
	=f(r)g(s|r)min{1,f(s)g(r|s) f(r)g(s|r)}
	=min{f(r)g(s|r), f(s)g(r|s)}
	=f(s)g(r|s)𝛼(s,r)
	=f(s)K(s,r)

We have the detailed balanced consition f(r)K(r,s)=f(s)K(s,r). Then ∫f(r)K(r,s)dr=∫f(s)K(s,r)dr=f(s)∫K(s,r)dr=f(s)6.2.2. Metropolis samplerWith symmetric proposal g(X|Y)=g(Y|X). r(Xt, Y)=f(Y)

f(Xt).6.2.3. Random Walk sampler Proposal g(Y|Xt)=g(|Y-Xt|). For example Y=Xt+N(0,𝜎2). r(Xt, Y)=f(Y)

f(Xt).6.2.4. independent samplerProposal g(Y|Xt)=g(Y). r(Xt, Y)=f(Y)g(Xt)

f(Xt)g(Y).6.2.5. Gibbs samplerNow suppose we generate samples from a multivariate f(x) where x∈𝜒⊂Rd. We partition the d-dimension vector x into K disjoint blocks, denoted by x=(x1,⋯,xK)T where K≤d. Let fk(xk|x-k)=fk(xk|x1,⋯,xk-1,xk+1,⋯,xK), k=1,⋯,K

Algorithm10Gibbs sampler1. Starting with an arbitrary point x(0)∈𝜒 with f(x(0))>02. At time t, (1) generate x(t)1∼f1(x1|x(t-1)2,⋯,x(t-1)K) ⋮ (k) generate x(t)k∼fk(xk|x(t)1,⋯,x(t)k-1,x(t-1)k+1,⋯,x(t-1)K) ⋮ (K) generate x(t)K∼fK(xK|x(t-1)2,⋯,x(t-1)K-1) 3. Set x(t)=(x(t)1,⋯,x(t)K)T

Gibbs sampler is a special case ofMH sampler, with accept ratio 1. Proof can be found in homework 7.6.3. Monitoring the convergence (Gelman-Rubin method)In MCMC, the chain may be traped in some local mode, we can run multiple chains from different initial point to check the convergence. Recall ANOVA in DATA130046 statistics II. Suppose there are J chains, and n is the number in each chain after discarding the burn-in period. 𝜓 is a function of data. Write𝜓(j)i=𝜓(X(j)1,⋯,X(j)i), i=1,⋯,n;j=1,⋯,J⏨𝜓(j)=1

nn∑i=1𝜓(j)i⏨𝜓=1

nn∑i=1⏨𝜓(j)iThen the Between sequence variance is B=n

J-1J∑j=1(⏨𝜓(j)-⏨𝜓)2 and the Within sequence variance is W=1

JJ∑j=1s2j where s2j=1

n-1n∑i=1(𝜓(j)i-⏨𝜓(j))2. The Gelman-Rubin statistic is R=n-1

nW+1

nB

Wwhich should decrease and converge to 1 if chain converges well. A recommended threshold is 1.1.7. EM algorithm The expectation-maximization algorithm is used to overcome the difficulties in maximizing likelihoods. Suppose except observed data Y, there are unobserved data U. f(Y|𝜃) is not easily evaluate, but f(Y,U|𝜃) is easy to calculate. Then we write the total log-likelihood aslogf(Y,U|𝜃)=logf(Y|𝜃)+logf(U|Y,𝜃)Take expactation with respect to u, and given Y=y, 𝜃=𝜃*, we have Eualogf(Y,U|𝜃)|Y=y,𝜃*a=ℓ(𝜃)+Eualogf(U|Y,𝜃)|Y=y,𝜃*aWrite it as(1)Q(𝜃,𝜃*)=ℓ(𝜃)+C(𝜃,𝜃*)Let 𝜃=𝜃*, we have (2)Q(𝜃*,𝜃*)=ℓ(𝜃*)+C(𝜃*,𝜃*)(1)-(2) we get ℓ(𝜃)-ℓ(𝜃*)≥Q(𝜃,𝜃*)-Q(𝜃*,𝜃*)-{C(𝜃,𝜃*)-C(𝜃*,𝜃*)}.The right part is a

C(𝜃,𝜃*)-C(𝜃*,𝜃*)	=E{logf(U|Y,𝜃) f(U|Y,𝜃*)|Y=y,𝜃*}
	≤logE{f(U|Y,𝜃) f(U|Y,𝜃*)|Y=y,𝜃*}
	=log{∫f(u|y,𝜃) f(u|y,𝜃*)f(u|y,𝜃*)du}
	=0

Thus ℓ(𝜃)-ℓ(𝜃*)=Q(𝜃,𝜃*)-Q(𝜃*,𝜃*). Given 𝜃*, if we find 𝜃 that Q(𝜃,𝜃*)≥Q(𝜃*,𝜃*), then ℓ(𝜃)≥ℓ(𝜃*).

Algorithm11The EM algorithm1. E-step: compute Q(𝜃,𝜃(t))=EUalogf(Y,U|𝜃)aY=y,𝜃(t)a2. M-step: 𝜃(t+1)=argmax𝜃Q(𝜃,𝜃(t))

8. Variational inferenceTo approximate the posterior p(z|x)=𝜋(z)f(x|z)

f(x) in bayesian statistics, rather than generating random samples using MCMC, we can solve the problem by optimization, i.e., find the closest density to p(x|z). To judge the distance between our density q and the target p, we can use the Kullback-Leibler divergence. 8.1. KL divergenceKL(f‖g)=Ef{logf(x)

g(x)}=+∞∫-∞logf(x)

g(x)f(x)dxProperties:1. KL divergence is non-negative and KL=0 iff f=g2. KL is not a distance, i.e., KL(f||g)≠KL(g||f)Proof of Property 1:-KL(f||g)=Ef{logg(x)

f(x)}Jensen's inequality

≤

log{Efg(x)

f(x)}=log+∞∫-∞g(x)

f(x)f(x)dx=08.2. Evidence lower bound (ELBO)Suppose we have a family of densities over latent variables, then we wantq*(z)=argminq∈QKL(q(z) || p(x|z))a

KL(q(z) || p(x|z))	=Eq{logq(z) p(z|x)}
	=Eq{log q(z)}-Eq{log p(z|x)}
	=Eq{logq(z)}-Eq{logf(x,z)}+Eq{logf(x)}
	=logf(x)-ELBO(q)

Define evidence lower bound (ELBO) as ELBO(q)≜Eq{logf(x,z)}-Eq{logq(z)}Then logf(x)=KL(q(z)||p(z|x))=KL(q(z) || p(x|z))+ELBO(q)≥ELBO(q)To minimize KL(q(z) || p(x|z)), we maximize ELBO(q).a

ELBO(q)	=Eq{logf(x,z)}-Eq{logq(z)}
	=Eq{logf(x|z)}+Eq{log𝜋(z)}-Eq{logq(z)}
	=Eq{logf(x|z)}-Eq{logq(z) 𝜋(z)}
	=Eq{logf(x|z)}-KL(q(z)||𝜋(z))

To maximize ELBO(q), on one hand, we maximize Eq{logf(x|z)}, i.e., find q that fits the data. On the other hand, we minimize KL(q(z)||𝜋(z)), i.e., find q that fits the prior.8.3. The mean-field variational familyNow we consider a special case that the family of the q is Q≜aq: q(z)=m∏j=1qj(zj)a. It can be treated as the elements of the q is "indepdent".a

ELBO(q)	=Eq{logf(x,z)}-Eq{logq(z)}
	=∫m∏j=1qjalogf(x,z)-logqjadz

We use a method called coordinate ascend to maximize ELBO(q), i.e., fix all others, climb up on qk. Then ELBO(q)=∫qk{∫logf(x,z) ∏j≠kqjdzj}dzk-∫qklogqkdzk+c1, where c1 is the constant have nothing with qk.Define logp(x,zk)=∫logf(x,z) ∏j≠kqjdzj+c2, thus ELBO(q)=∫qklogp(x,zk)dzk-∫qklogqkdzk+c1=-KL(qk||p(x,zk))+c3q*k=p∝expaE-klogf(x,z)a