Semiparametric Mixtures of Nonparametric Regressions Sijia Xiang, Weixin Yao %%%%%%%%% %% LEM %% %%%%%%%%% %% Step 1 function[est_p,est_mu,est_var]=backfitting_step1(x,y,x0,est_p,est_mu,est_var,h) [n,numc]=size(est_mu); m=length(x0); r=[];f=[];est_p_x0=[];est_mu_x0=[];est_var_x0=[]; rsd=1; likelihood=1; iter=1; acc=10^(-3); while(abs(rsd)>acc && iter<200) %E-step for j=1:numc f(:,j)=normpdf(y,est_mu(:,j),est_var(:,j).^0.5); end %M-step for j=1:numc r(:,j) = est_p(:,j).*f(:,j)./sum(f.*est_p,2); W=exp(-(repmat(x',m,1)-repmat(x0,1,n)).^2/2/h^2)/h/sqrt(2*pi); R=repmat(r(:,j),1,m)'; RY=repmat(r(:,j).*y,1,m)'; est_p_x0(:,j)=sum(W.*R,2)./sum(W,2); est_p(:,j)=interp1(x0,est_p_x0(:,j),x,'linear','extrap'); est_mu_x0(:,j)=sum(W.*RY,2)./sum(W.*R,2); est_mu(:,j)=interp1(x0,est_mu_x0(:,j),x,'linear','extrap'); RE=(repmat(y,1,m)'-repmat(est_mu_x0(:,j),1,n)).^2; %RE=repmat((y-est_mu(:,j)).^2.*r(:,j),1,m)'; est_var_x0(:,j)=sum(W.*RE.*R,2)./sum(W.*R,2); est_var(:,j)=interp1(x0,est_var_x0(:,j),x,'linear','extrap'); end rsd=likelihood-sum(log(sum(f.*est_p,2))); likelihood=sum(log(sum(f.*est_p,2))); iter=iter+1; end %% Step 2 function[est_p,est_var]=backfitting_step2(x,y,x0,est_p,est_mu,est_var) [n,numc]=size(est_mu); m=length(x0); r=[];f=[]; rsd=1; likelihood=1; iter=1; acc=10^(-3); while(abs(rsd)>acc && iter<200) %E-step for j=1:numc f(:,j)=normpdf(y,est_mu(:,j),est_var(:,j).^0.5); end %M-step for j=1:numc r(:,j) = est_p(:,j).*f(:,j)./sum(f.*est_p,2); est_p(:,j)=repmat(mean(r(:,j)),n,1); est_var(:,j)=repmat(sum((y-est_mu(:,j)).^2.*r(:,j))/sum(r(:,j)),n,1); end rsd=likelihood-sum(log(sum(f.*est_p,2))); likelihood=sum(log(sum(f.*est_p,2))); iter=iter+1; end %% Step 3 function[est_mu,est_mu_x0,likelihood]=backfitting_step3(x,y,x0,est_p,est_mu,est_var,h) [n,numc]=size(est_mu); m=length(x0); r=[];f=[];est_mu_x0=[]; rsd=1; likelihood=1; iter=1; acc=10^(-3); while(abs(rsd)>acc && iter<200) %E-step for j=1:numc f(:,j)=normpdf(y,est_mu(:,j),est_var(:,j).^0.5); end %M-step for j=1:numc r(:,j) = est_p(:,j).*f(:,j)./sum(f.*est_p,2); W=exp(-(repmat(x',m,1)-repmat(x0,1,n)).^2/2/h^2)/h/sqrt(2*pi); R=repmat(r(:,j),1,m)'; RY=repmat(r(:,j).*y,1,m)'; est_mu_x0(:,j)=sum(W.*RY,2)./sum(W.*R,2); est_mu(:,j)=interp1(x0,est_mu_x0(:,j),x,'linear','extrap'); end rsd=likelihood-sum(log(sum(f.*est_p,2))); likelihood=sum(log(sum(f.*est_p,2))); iter=iter+1; end %%%%%%%%% %% GEM %% %%%%%%%%% function[est_p,est_mu,est_var,est_mu_x0]=backfitting(x,y,x0,est_p,est_mu,est_var,h) [n,numc]=size(est_mu); m=length(x0); r=[];f=[];est_mu_x0=[]; rsd=1; likelihood=1; iter=1; acc=10^(-3); while(abs(rsd)>acc && iter<200) %E-step for j=1:numc f(:,j)=normpdf(y,est_mu(:,j),est_var(:,j).^0.5); end %M-step for j=1:numc r(:,j) = est_p(:,j).*f(:,j)./sum(f.*est_p,2); est_p(:,j)=repmat(sum(r(:,j))/n,n,1); est_var(:,j)=repmat(sum(r(:,j).*(y-est_mu(:,j)).^2)./sum(r(:,j)),n,1); W=exp(-(repmat(x',m,1)-repmat(x0,1,n)).^2/2/h^2)/h/sqrt(2*pi); R=repmat(r(:,j)',m,1); RY=repmat((r(:,j).*y)',m,1); est_mu_x0(:,j)=sum(W.*RY,2)./sum(W.*R,2); est_mu(:,j)=interp1(x0,est_mu_x0(:,j),x,'linear','extrap'); % est_mu(:,j)=interp1q(x0,est_mu_x0(:,j),x); end rsd=likelihood-sum(log(sum(f.*est_p,2))); likelihood=sum(log(sum(f.*est_p,2))); iter=iter+1; end