function [index_ini,index,theta,theta0,Test,runlength] = thom4(data,trueA0,trueA,m,sigma,niter,w,sigma9,lambda0,controlimit,sigma0,v) [n,p] = size(data); [~,q0] = size(trueA0); [~,q] = size(trueA); index = zeros(m,n); B0=zeros(m,q0,n); B1=zeros(m,q,n); X0=zeros(m,n); XB=zeros(q,n); BB=zeros(q,n); BBjk=zeros(q,q,n); theta=zeros(q,n); theta0=zeros(q0,n); Test = zeros(1,n); runlength = n; %% for the first n=1 index_ini =1:p; B0_ini = trueA0(index_ini,:); B1_ini = trueA(index_ini,:); BB(:,1)=diag(B1_ini'*B1_ini); BBjk(:,:,1)=B1_ini'*B1_ini; s2=1./(BB(:,1)/sigma^2+1./sigma9.^2); a=0.5*ones(q,1); u=zeros(q,1); X0_ini = data(1,index_ini)'; mu_tmp = zeros(q,1); inv_B0B0 = inv(B0_ini'*B0_ini+1e-6.*eye(q0)); theta0(:,1) = inv_B0B0*(B0_ini'*(X0_ini-B1_ini*mu_tmp)); X1_ini = X0_ini - B0_ini*theta0(:,1); for h=1:niter XB(:,1) = (X1_ini'*B1_ini)'; ss=zeros(q,1); for j=1:q u(j)=(XB(j,1)-(BBjk(j,:,1)*(a.*u)-BBjk(j,j,1)*a(j)*u(j)))*s2(j)/sigma^2; %ss(j)=log(w(j)/(1-w(j)))+1/2-(s2(j)+u(j)^2)/(2*sigma9(j)^2)+log(sqrt(s2(j))/sigma9(j))+u(j)^2/s2(j)-BB(j,1)*(u(j)^2+s2(j))/(2*sigma^2); ss(j)=log(w(j)/(1-w(j)))+u(j)^2/(2*sigma9(j)^2)+BB(j,1)*(u(j)^2-s2(j)+v*s2(j))/(2*sigma^2); a(j)=exp(ss(j))/(1+exp(ss(j))); if isnan(a(j))==1 a(j)=1; end end mu_tmp = u.*a; %theta0(:,1) = inv_B0B0*(B0_ini'*(X0_ini-B1_ini*mu_tmp)); sigma_b(:,:,1)=(B0_ini'*(sigma^(-2).*eye(p))*B0_ini+sigma0^(-1).*eye(q0))^(-1); mu_theta0(:,1)=sigma_b(:,:,1)*B0_ini'*(sigma^(-2).*eye(p))*(X0_ini-B1_ini*mu_tmp); X1_ini = X0_ini - B0_ini*mu_theta0(:,1); end r=zeros(q,1); for j=1:q r(j)=random('Binomial',1,a(j)); if r(j)==1 theta(j,1)=normrnd(u(j),s2(j)^0.5); else %theta(j,1)=0; theta(j,1)=normrnd(0,(v*s2(j))^0.5); end end for j=1:q for k=1:q A(j,k)=a(j)*a(k); end A(j,j)=A(j,j)^(1/2); end % X=normrnd(trueA*theta(:,1),sigma); % thom=2*trueA*theta(:,1).*X-(trueA*theta(:,1)).^2; % shuffle = randsample(1:p,p,0); % [~,ind] = sort(thom(shuffle),'descend'); % index(:,2)=shuffle(ind(1:m)); theta_a = u.*a; allset=nchoosek(1:p,m); [num,~]=size(allset); thom=zeros(num,1); for k=1:num B1_tmp=trueA(allset(k,:),:); B0_tmp=trueA0(allset(k,:),:); X_tmp=normrnd(B1_tmp*theta(:,1),sigma); H_tmp=B0_tmp*inv(B0_tmp'*B0_tmp)*B0_tmp'; %thom(k)=theta(:,1)'*B1_tmp'*H*B1_tmp*theta(:,1); %thom(k)=X'*H*X; thom(k)=2*X_tmp'*(eye(m)-H_tmp)*B1_tmp*theta_a-u'*(B1_tmp'*(eye(m)-H_tmp)*B1_tmp.*A)*u+theta_a'*B1_tmp'*H_tmp*B1_tmp*theta_a; end matmp=allset(thom==max(thom),:); myrow=size(allset(thom==max(thom),:),1); index(:,2)=matmp(randsample(1:myrow,1),:); H=B0_ini*inv(B0_ini'*B0_ini)*B0_ini'; Test(1)=2*X1_ini'*squeeze(B1_ini)*theta_a-u'*(squeeze(B1_ini)'*squeeze(B1_ini).*A)*u... -2*X1_ini'*H*squeeze(B1_ini)*theta_a+theta_a'*squeeze(B1_ini)'*H*squeeze(B1_ini)*theta_a; %% for later n = 2, 3,.... for i = 2:n B0(:,:,i)= trueA0(index(:,i),:); B1(:,:,i)= trueA(index(:,i),:); % X0(:,i)= data(i,index(:,i))'; BB(:,i)=BB(:,i-1)*lambda0+diag(B1(:,:,i)'*B1(:,:,i)); BBjk(:,:,i)=BBjk(:,:,i-1)*lambda0+B1(:,:,i)'*B1(:,:,i); s2=1./(BB(:,i)/sigma^2+1./sigma9.^2); a=0.5*ones(q,1); u=zeros(q,1); X0(:,i)= data(i,index(:,i))'; mu_tmp = zeros(q,1); inv_B0B0 = inv(B0(:,:,i)'*B0(:,:,i)+1e-6.*eye(q0)); theta0(:,i) = inv_B0B0*(B0(:,:,i)'*(X0(:,i)-B1(:,:,i)*mu_tmp)); X1(:,i)= X0(:,i) - B0(:,:,i)*theta0(:,i); for h=1:niter if i ==1 XB(:,i) = (X1(:,i)'*B1(:,:,i))'; else XB(:,i)=XB(:,i-1)*lambda0+(X1(:,i)'*B1(:,:,i))'; end ss=zeros(q,1); for j=1:q u(j)=(XB(j,i)-(BBjk(j,:,i)*(a.*u)-BBjk(j,j,i)*a(j)*u(j)))*s2(j)/sigma^2; % fid=fopen('test3.txt','at'); % fprintf(fid,'%f ',u(j)); % fclose(fid); %ss(j)=log(w(j)/(1-w(j)))+1/2-(s2(j)+u(j)^2)/(2*sigma9(j)^2)+log(sqrt(s2(j))/sigma9(j))+u(j)^2/s2(j)-BB(j,i)*(u(j)^2+s2(j))/(2*sigma^2); ss(j)=log(w(j)/(1-w(j)))+u(j)^2/(2*sigma9(j)^2)+BB(j,i)*(u(j)^2-s2(j)+v*s2(j))/(2*sigma^2); a(j)=exp(ss(j))/(1+exp(ss(j))); if isnan(a(j))==1 a(j)=1; end end mu_tmp = u.*a; %theta0(:,i) = inv_B0B0*(B0(:,:,i)'*(X0(:,i)-B1(:,:,i)*mu_tmp)); sigma_b(:,:,i)=(B0(:,:,i)'*(sigma^(-2).*eye(m))*B0(:,:,i)+sigma0^(-1).*eye(q0))^(-1); mu_theta0(:,i)=sigma_b(:,:,i)*B0(:,:,i)'*(sigma^(-2).*eye(m))*(X0(:,i)-B1(:,:,i)*mu_tmp); X1(:,i) = X0(:,i) - B0(:,:,i)*mu_theta0(:,i); % fid=fopen('test3.txt','at'); % fprintf(fid,'\n'); % fprintf(fid,'\n'); % fclose(fid); end r=zeros(q,1); for j=1:q r(j)=random('Binomial',1,a(j)); if r(j)==1 theta(j,i)=normrnd(u(j),s2(j)^0.5); else %theta(j,i)=0; theta(j,i)=normrnd(0,(v*s2(j))^0.5); end end for j=1:q for k=1:q A(j,k)=a(j)*a(k); end A(j,j)=A(j,j)^(1/2); end % X=normrnd(trueA*theta(:,i),sigma); % thom=2*trueA*theta(:,i).*X-(trueA*theta(:,i)).^2; % shuffle = randsample(1:p,p,0); % [~,ind] = sort(thom(shuffle),'descend'); % index(:,i+1)=shuffle(ind(1:m)); theta_a = u.*a; allset=nchoosek(1:p,m); [num,~]=size(allset); thom=zeros(num,1); for k=1:num B1_tmp=trueA(allset(k,:),:); B0_tmp=trueA0(allset(k,:),:); X_tmp=normrnd(B1_tmp*theta(:,i),sigma); %thom(k)=2*theta(:,i)'*B1_tmp'*X-theta(:,i)'*B1_tmp'*B1_tmp*theta(:,i); H_tmp=B0_tmp*inv(B0_tmp'*B0_tmp)*B0_tmp'; thom(k)=2*X_tmp'*(eye(m)-H_tmp)*B1_tmp*theta_a-u'*(B1_tmp'*B1_tmp.*A)*u+theta_a'*B1_tmp'*H_tmp*B1_tmp*theta_a; end matmp=allset(thom==max(thom),:); myrow=size(allset(thom==max(thom),:),1); index(:,i+1)=matmp(randsample(1:myrow,1),:); H=B0(:,:,i)*inv(B0(:,:,i)'*B0(:,:,i))*B0(:,:,i)'; Test(i)=2*X1(:,i)'*squeeze(B1(:,:,i))*theta_a-u'*(squeeze(B1(:,:,i))'*squeeze(B1(:,:,i)).*A)*u... -2*X1(:,i)'*H*squeeze(B1(:,:,i))*theta_a+theta_a'*squeeze(B1(:,:,i))'*H*squeeze(B1(:,:,i))*theta_a; if Test(i)>controlimit runlength = i; break; end end end