<jupyter_start><jupyter_text>Define the experiment class<jupyter_code>class MLE:
    def __init__(self,dim,sigma):
        self.dim=dim
        self.T=lambda x: mvnorm.pdf(x=x,mean=np.zeros(dim))
        self.iP=lambda x: mvnorm.pdf(x=x,mean=np.zeros(dim),cov=sigma**2)
        self.iS=lambda size: mvnorm.rvs(size=size,mean=np.zeros(dim),cov=sigma**2)
        
    def __estimate(self,W,name,asym=True):
        Z=np.mean(W)
        err=np.abs(Z-1)
        if asym:
            aVar=np.var(W)
            aErr=np.sqrt(aVar/W.size)
            ESS=1/np.sum((W/np.sum(W))**2)
            print('{} est: {:.4f}; err: {:.4f}; a-var: {:.4f}; a-err: {:.4f}; ESS: {:.0f}/{}'\
                  .format(name,Z,err,aVar,aErr,ESS,W.size))
        else:
            print('{} est: {:.4f}; err: {:.4f}'.format(name,Z,err))
        
    def estimate_IS(self,size):
        S=self.iS(size)
        W=self.T(S)/self.iP(S)
        self.__estimate(W,'IS')
    
    def draw_TP(self,P,x,name,dim=0):
        X=np.zeros([x.size,self.dim])
        X[:,dim]=x
        ratio=np.reshape(self.T(np.zeros(self.dim))/P(np.zeros([1,self.dim])),1)[0]
        print('------------ pdf ratio at origin: {:.2f} ------------'.format(ratio))
        
        fig,ax=plt.subplots(figsize=(7,4))
        ax.plot(x,self.T(X))
        ax.plot(x,P(X))
        if name=='nonparametric':
            one=np.zeros(self.dim)
            one[dim]=1
            rW=np.array([self.h(one*loc,loc) for loc in self.rS])
            rW=rW/rW.max()*P(np.zeros([1,self.dim]))[0]
            rWmeans=np.ones_like(rW)*rW.mean()

            ax.plot(x,self.mP(X))
            ax.hist(self.rS[:,dim],bins=2*rW.size,weights=rWmeans)
            ax.hist(self.rS[:,dim],bins=2*rW.size,weights=rW)
            ax.legend(['target','nonparametric proposal','mixture proposal','centers','centers with weight'])
        elif name=='regression':
            G=self.G(X)
            rPO=self.regO.coef_.dot(G)+self.regO.intercept_*P(X)
            rPL=self.regL.coef_.dot(G)+self.regL.intercept_*P(X)
            mid=int(x.size/2)
            print('regression ratios: ordinary {:.4f}, lasso {:.4f}'\
                  .format(self.T(X[mid])/rPO[mid],self.T(X[mid])/rPL[mid]))
            
            ax.plot(x,rPO)
            ax.plot(x,rPL)
            ax.legend(['target','mixture proposal','ordinary regression','lasso regression'])
        else:
            ax.legend(['target','{} proposal'.format(name)])
            
        ax.set_title('{}-D target and {} proposal (cross-sectional view)'.format(self.dim,name))
        plt.show()
        
    def resample(self,size,ratio):
        S=self.iS(ratio*size)
        p=self.T(S)/self.iP(S)
        index=np.arange(S.shape[0])
        choice=np.random.choice(index,size,p=p/np.sum(p),replace=True)
        
        self.rS=S[choice]
        self.rSset=S[list(set(choice))]
        print('resampling rate: {}/{}'.format(self.rSset.shape[0],size))
        
    def estimate_NIS(self,size,rate,bdwth='scott',reg=0):
        if(type(bdwth)==str):
            method=bdwth
            tmp=sciKDE(self.rS.T,bw_method=method)
            bdwth=np.mean(np.sqrt(np.diag((tmp.covariance_factor()**2)*np.cov(self.rS.T))))
            print('bdwth: {:.4f} (based on {} ({:.4f}))'.format(bdwth,method,tmp.covariance_factor()))
        
        self.bdwth=bdwth
        self.kde=sklKDE(kernel='gaussian',bandwidth=bdwth).fit(self.rS)
        self.h=lambda x,loc: mvnorm.pdf(x=x,mean=loc,cov=self.bdwth**2)
        if reg==0:
            self.G=lambda x: np.array([self.h(x,loc) for loc in self.rSset])-self.iP(x)
        elif reg==1:
            self.G=lambda x: np.array([self.h(x,self.rSset[i])-self.h(x,self.rSset[0])\
                                       for i in range(1,self.rSset.shape[0])])
        elif reg==2:
            self.G=lambda x: np.array([self.h(x,self.rSset[i])-self.h(x,self.rSset[i-1])\
                                       for i in range(1,self.rSset.shape[0])])
        elif reg==3:
            def G(x):
                tmp=np.array([self.h(x,loc) for loc in self.rSset])
                return tmp[1:]-tmp.mean(axis=0)
                
            self.G=G
        
        self.nP=lambda x: np.exp(self.kde.score_samples(x))
        self.nS=lambda size: self.kde.sample(size)
        S=self.nS(size)
        W=self.T(S)/self.nP(S)
        self.__estimate(W,'NIS')
        
        self.mP=lambda x: (1-rate)*self.iP(x)+rate*self.nP(x)
        self.mS=lambda size: np.vstack([self.iS(size-round(rate*size)),self.nS(round(rate*size))])
        self.S=self.mS(size)
        W=self.T(self.S)/self.mP(self.S)
        self.__estimate(W,'MIS')
        
    def estimate_RIS(self,alpha,vif=False):
        X=(self.G(self.S)/self.mP(self.S)).T
        if vif:
            if X.shape[1]<=100:
                ind=np.arange(X.shape[1])
            else:
                ind=np.random.choice(np.arange(X.shape[1]),100,replace=False)

            tmp=np.hstack([X,np.ones([X.shape[0],1])])
            Vif=np.zeros(ind.size)
            for i in range(ind.size):
                Vif[i]=VIF(tmp,ind[i])

            print('VIF: (min {:.4f}, median {:.4f}, mean {:.4f}, max {:.4f}, [>5] {}/{})'\
                  .format(Vif.min(),np.median(Vif),Vif.mean(),Vif.max(),np.sum(Vif>5),Vif.size))
            
        Xn=X/np.linalg.norm(X,axis=0)
        lbd=np.linalg.eigvals(Xn.T.dot(Xn))
        tau=np.sqrt(lbd.max()/lbd)
        print('Condition index: (min {:.4f}, median {:.4f}, mean {:.4f}, max {:.4f}, [>30] {}/{})'\
              .format(tau.min(),np.median(tau),tau.mean(),tau.max(),np.sum(tau>30),tau.size))
        
        y=self.T(self.S)/self.mP(self.S)
        self.regO=Linear().fit(X,y)
        self.regL=Lasso(alpha).fit(X,y)
        print('Ordinary R2: {:.4f}; Lasso R2: {:.4f}'.format(self.regO.score(X,y),self.regL.score(X,y)))
        
        W=y-X.dot(self.regO.coef_)
        self.__estimate(W,'RIS(Ord)')
        W=y-X.dot(self.regL.coef_)
        self.__estimate(W,'RIS(Las)')
    
    def estimate_MLE(self,opt=False,init=0):
        mP=self.mP(self.S)
        G=self.G(self.S)
        target=lambda zeta: -np.mean(np.log(mP+zeta.dot(G)))
        gradient=lambda zeta: -np.mean(G/(mP+zeta.dot(G)),axis=1)
        hessian=lambda zeta: (G/(mP+zeta.dot(G))**2).dot(G.T)/G.shape[1]
        zeta0=np.zeros(G.shape[0])
        grad0=gradient(zeta0)
        print('Reference:')
        print('origin: value: {:.4f}; grad: (min {:.4f}, mean {:.4f}, max {:.4f}, std {:.4f})'\
              .format(target(zeta0),grad0.min(),grad0.mean(),grad0.max(),grad0.std()))
        
        print()
        print('Theoretical results:')
        X=(G/mP).T
        XX=X-X.mean(axis=0)
        zeta1=np.linalg.solve(XX.T.dot(XX),X.sum(axis=0))
        print('MLE(The) zeta: (min {:.4f}, mean {:.4f}, max {:.4f}, std {:.4f}, norm {:.4f})'\
              .format(zeta1.min(),zeta1.mean(),zeta1.max(),zeta1.std(),np.sqrt(np.sum(zeta1**2))))
        grad1=gradient(zeta1)
        print('theory: value: {:.4f}; grad: (min {:.4f}, mean {:.4f}, max {:.4f}, std {:.4f})'\
              .format(target(zeta1),grad1.min(),grad1.mean(),grad1.max(),grad1.std()))
        W=(self.T(self.S)/mP)*(1-XX.dot(zeta1))
        self.__estimate(W,'RIS(The)',asym=False)
        W=self.T(self.S)/(mP+zeta1.dot(G))
        self.__estimate(W,'MLE(The)',asym=False)
        
        if opt:
            zeta=zeta1 if init==1 else zeta0
            begin=dt.now()
            res=root(lambda zeta: (gradient(zeta),hessian(zeta)),zeta,method='lm',jac=True)
            end=dt.now()
            print()
            print('Optimization results (spent {} seconds):'.format((end-begin).seconds))
            if res['success']:
                zeta=res['x']
                print('MLE(Opt) zeta: (min {:.4f}, mean {:.4f}, max {:.4f}, std {:.4f}, norm {:.4f})'\
                      .format(zeta.min(),zeta.mean(),zeta.max(),zeta.std(),np.sqrt(np.sum(zeta**2))))
                print('Dist(zeta(Opt),zeta(The))={:.4f}'.format(np.sqrt(np.sum((zeta-zeta1)**2))))
                grad=gradient(zeta)
                print('optimal: value: {:.4f}; grad: (min {:.4f}, mean {:.4f}, max {:.4f}, std {:.4f})'\
                      .format(target(zeta),grad.min(),grad.mean(),grad.max(),grad.std()))
                W=self.T(self.S)/(mP+zeta.dot(G))
                self.__estimate(W,'MLE(Opt)',asym=False)
            else:
                print('MLE fail')<jupyter_output><empty_output><jupyter_text>**Limitations:**1. didn't consider self-normalized importance sampling2. only for symmetric normal target and symmetric normal initial proposal with only one mode3. only for normal KDE without weights and adaptive bandwidth4. only consider regression for proposal components based on mixture proposal  Run the experiments  Initial proposal and the curse of dimensionality<jupyter_code>mle=MLE(dim=2,sigma=2)
size=1000000
mle.estimate_IS(size)
x=np.linspace(-4,4,101)
mle.draw_TP(mle.iP,x,'initial')

print('=======================================================')
print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++')
print('=======================================================')

mle=MLE(dim=4,sigma=2)
mle.estimate_IS(size)
mle.draw_TP(mle.iP,x,'initial')

print('=======================================================')
print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++')
print('=======================================================')

mle=MLE(dim=8,sigma=2)
mle.estimate_IS(size)
mle.draw_TP(mle.iP,x,'initial')<jupyter_output>IS est: 1.0007; err: 0.0007; a-var: 1.2873; a-err: 0.0011; ESS: 437529/1000000
------------ pdf ratio at origin: 4.00 ------------<jupyter_text>**Summary:*** The cumulative phenomenon of scale difference of all dimension. **Future:*** What about the comparison between multi-normal and multi-t?   KDE with different initial proposal and bandwidth in 6-dimension<jupyter_code>np.random.seed(1234)
mle=MLE(dim=6,sigma=1)
size=50000
mle.estimate_IS(size)
x=np.linspace(-4,4,101)
mle.draw_TP(mle.iP,x,'initial')
print('=======================================================')
mle.resample(1000,20)
mle.estimate_NIS(size,0.9)
mle.draw_TP(mle.nP,x,'nonparametric')

print('=======================================================')
print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++')
print('=======================================================')

np.random.seed(1234)
mle=MLE(dim=6,sigma=2)
mle.estimate_IS(size)
mle.draw_TP(mle.iP,x,'initial')
print('=======================================================')
mle.resample(1000,20)
np.random.seed(123456)
mle.estimate_NIS(size,0.9)
mle.draw_TP(mle.nP,x,'nonparametric')
np.random.seed(123456)
mle.estimate_NIS(size,0.9,0.6)
mle.draw_TP(mle.nP,x,'nonparametric')

print('=======================================================')
print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++')
print('=======================================================')

np.random.seed(123456)
mle.estimate_NIS(size,0.9,0.5)
mle.draw_TP(mle.nP,x,'nonparametric')

print('=======================================================')
print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++')
print('=======================================================')

np.random.seed(123456)
mle.estimate_NIS(size,0.9,0.4)
mle.draw_TP(mle.nP,x,'nonparametric')

print('=======================================================')
print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++')
print('=======================================================')

np.random.seed(123456)
mle.estimate_NIS(size,0.8,0.4)
mle.draw_TP(mle.nP,x,'nonparametric')<jupyter_output>NIS est: 0.9986; err: 0.0014; a-var: 0.3562; a-err: 0.0027; ESS: 36840/50000
MIS est: 1.0006; err: 0.0006; a-var: 0.4574; a-err: 0.0030; ESS: 34322/50000
------------ pdf ratio at origin: 2.70 ------------<jupyter_text>**Summary:*** Resampling ratio 20:1 can give us good 'samples' from the target. (We can increase it for better resample rate in higher dimension. )* Only conservative bandwidth can work well in higher dimension, because little bandwidth can raise serious tail problem. **Future*** What about KDE with weights? * What about the adaptive bandwidth?   Regression performance visualization<jupyter_code>np.random.seed(1234)
mle=MLE(dim=8,sigma=2)
size=100000
mle.estimate_IS(size)
x=np.linspace(-4,4,101)
mle.draw_TP(mle.iP,x,'initial')
print('=======================================================')
mle.resample(1000,20)
mle.estimate_NIS(size,0.9)
mle.draw_TP(mle.nP,x,'nonparametric')
print('=======================================================')
mle.estimate_RIS(0.1)
mle.draw_TP(mle.mP,x,'regression')
np.random.seed(1234)
mle=MLE(dim=10,sigma=2)
mle.estimate_IS(size)
mle.draw_TP(mle.iP,x,'initial')
print('=======================================================')
np.random.seed(123456)
mle.resample(1000,20)
mle.estimate_NIS(size,0.9)
mle.draw_TP(mle.nP,x,'nonparametric')
print('=======================================================')
mle.estimate_RIS(0.1)
mle.draw_TP(mle.mP,x,'regression')
np.random.seed(123456)
mle.resample(1000,100)
mle.estimate_NIS(size,0.9)
mle.draw_TP(mle.nP,x,'nonparametric')
print('=======================================================')
mle.estimate_RIS(0.1)
mle.draw_TP(mle.mP,x,'regression')<jupyter_output>resampling rate: 787/1000
bdwth: 0.7891 (based on scott (0.6105))
NIS est: 0.9955; err: 0.0045; a-var: 1.5623; a-err: 0.0040; ESS: 38811/100000
MIS est: 1.0036; err: 0.0036; a-var: 1.8405; a-err: 0.0043; ESS: 35369/100000
------------ pdf ratio at origin: 12.95 ------------<jupyter_text>**Summary:*** Regression improve the estimation performance by mimic the shape of the target. * In higher dimension, improve the quality of kernels by increasing resampling ratio can slightly improve KDE but more dramatically improve regression. **Future*** The problem of singularity? (What about multi-t initial? )* What about regression for kernels not used in KDE? * What about regression with KDE with weights? * What about regression for kernels with adaptive bandwidth? (Will that improve the kernels' quality? )  MLE method investigation<jupyter_code>np.random.seed(1234)
mle=MLE(dim=10,sigma=2)
size=100000
mle.estimate_IS(size)
x=np.linspace(-4,4,101)
mle.draw_TP(mle.iP,x,'initial')
print('=======================================================')
mle.resample(1000,100)
mle.estimate_NIS(size,0.9)
mle.draw_TP(mle.nP,x,'nonparametric')
print('=======================================================')
mle.estimate_RIS(0.1)
mle.draw_TP(mle.mP,x,'regression')
print('=======================================================')
mle.estimate_MLE(True)
np.random.seed(1234)
mle=MLE(dim=10,sigma=1)
mle.estimate_IS(size)
mle.draw_TP(mle.iP,x,'initial')
print('=======================================================')
mle.resample(1000,100)
mle.estimate_NIS(size,0.9)
mle.draw_TP(mle.nP,x,'nonparametric')
print('=======================================================')
mle.estimate_RIS(0.1)
mle.draw_TP(mle.mP,x,'regression')
print('=======================================================')
mle.estimate_MLE(True)
np.random.seed(1234)
mle=MLE(dim=5,sigma=2)
size=20000
mle.estimate_IS(size)
x=np.linspace(-4,4,101)
mle.draw_TP(mle.iP,x,'initial')
print('=======================================================')
mle.resample(1000,100)
mle.estimate_NIS(size,0.9)
mle.draw_TP(mle.nP,x,'nonparametric')
print('=======================================================')
mle.estimate_RIS(0.1)
mle.draw_TP(mle.mP,x,'regression')
print('=======================================================')
mle.estimate_MLE(True,0)

print('=======================================================')
print('+++++++++++++++++++++++++++++++++++++++++++++++++++++++')
print('=======================================================')

mle.estimate_MLE(True,1)<jupyter_output>Condition index: (min 1.0000, median 32.3609, mean 85.9046, max 1758.5111, [>30] 504/963)
Ordinary R2: 0.9973; Lasso R2: 0.9411
RIS(Ord) est: 1.0005; err: 0.0005; a-var: 0.0010; a-err: 0.0002; ESS: 19979/20000
RIS(Las) est: 1.0003; err: 0.0003; a-var: 0.0228; a-err: 0.0011; ESS: 19555/20000
------------ pdf ratio at origin: 2.57 ------------
regression ratios: ordinary 0.9980, lasso 0.9457<jupyter_text>**Summary:*** Theoretical zeta is very close to optimal zeta, so we can use it as the initial guess for optimization. * Root finding for the gradient is far more efficient than maximization. * Theoretical RIS give exact same result as regression. * Theoretical MLE can achieve similar results as MLE based on optimization, so we can use it as an alternative method. * Theoretical MLE can't give exact result when target can be expressed by proposals, because its gradients are not small enough. * Sometimes (low-dimension, many kernels and small sample size? ) the theoretical zeta can be out of the feasible region, which indicates that the optimal zeta is very close to the boundary. * A bad theoretical zeta still is a good initial guess for optimization, which may result in infeasible but usable optimization result. **Future*** Try more detailed experiments for MLE.   Regression and MLE only with kernels<jupyter_code>np.random.seed(1234)
mle=MLE(dim=8,sigma=2)
size=20000
mle.estimate_IS(size)
x=np.linspace(-4,4,101)
mle.draw_TP(mle.iP,x,'initial')
print('=======================================================')
mle.resample(500,50)<jupyter_output>IS est: 0.9448; err: 0.0552; a-var: 23.0217; a-err: 0.0339; ESS: 747/20000
------------ pdf ratio at origin: 256.00 ------------<jupyter_text>Control variables: $K_i(x)-q_0(x)$, $i\ge 1$<jupyter_code>np.random.seed(1997)
mle.estimate_NIS(size,0.9,reg=0)
mle.draw_TP(mle.nP,x,'nonparametric')
print('=======================================================')
mle.estimate_RIS(0.1)
mle.draw_TP(mle.mP,x,'regression')
print('=======================================================')
mle.estimate_MLE(True)<jupyter_output>bdwth: 0.7982 (based on scott (0.5958))
NIS est: 1.0001; err: 0.0001; a-var: 1.3666; a-err: 0.0083; ESS: 8452/20000
MIS est: 1.0054; err: 0.0054; a-var: 1.5742; a-err: 0.0089; ESS: 7821/20000
------------ pdf ratio at origin: 9.55 ------------<jupyter_text>Control variables: $K_i(x)-K_1(x)$, $i\ge 2$<jupyter_code>np.random.seed(1997)
mle.estimate_NIS(size,0.9,reg=1)
mle.draw_TP(mle.nP,x,'nonparametric')
print('=======================================================')
mle.estimate_RIS(0.1)
mle.draw_TP(mle.mP,x,'regression')
print('=======================================================')
mle.estimate_MLE(True)<jupyter_output>bdwth: 0.7982 (based on scott (0.5958))
NIS est: 1.0001; err: 0.0001; a-var: 1.3666; a-err: 0.0083; ESS: 8452/20000
MIS est: 1.0054; err: 0.0054; a-var: 1.5742; a-err: 0.0089; ESS: 7821/20000
------------ pdf ratio at origin: 9.55 ------------<jupyter_text>Control variables: $K_i(x)-K_{i-1}(x)$, $i\ge 2$<jupyter_code>np.random.seed(1997)
mle.estimate_NIS(size,0.9,reg=2)
mle.draw_TP(mle.nP,x,'nonparametric')
print('=======================================================')
mle.estimate_RIS(0.1)
mle.draw_TP(mle.mP,x,'regression')
print('=======================================================')
mle.estimate_MLE(True)<jupyter_output><empty_output><jupyter_text>Control variables: $K_i(x)-\overline{K}(x)$, $i\ge 2$<jupyter_code>np.random.seed(1997)
mle.estimate_NIS(size,0.9,reg=3)
mle.draw_TP(mle.nP,x,'nonparametric')
print('=======================================================')
mle.estimate_RIS(0.1)
mle.draw_TP(mle.mP,x,'regression')
print('=======================================================')
mle.estimate_MLE(True)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.300331379046648, 'fraction_numerical': 0.047616619933724194, 'mean_word_length': 3.908908908908909, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 39, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '282386', 'license': 'MIT', 'n_tokens_mistral': 7835, 'n_tokens_neox': 7033, 'n_words': 1113, 'path': 'MLE and REG/.ipynb_checkpoints/MLE-1-checkpoint.ipynb', 'repo_name': 'IanFla/Two-Stage', 'size': '19615'}