diff --git a/BNP_covreg.m b/matlab/BNP_covreg.m similarity index 100% rename from BNP_covreg.m rename to matlab/BNP_covreg.m diff --git a/BNP_covreg_varinds.m b/matlab/BNP_covreg_varinds.m similarity index 100% rename from BNP_covreg_varinds.m rename to matlab/BNP_covreg_varinds.m diff --git a/flu_US.mat b/matlab/flu_US.mat similarity index 100% rename from flu_US.mat rename to matlab/flu_US.mat diff --git a/runstuff_BNPcovreg.m b/matlab/runstuff_BNPcovreg.m similarity index 99% rename from runstuff_BNPcovreg.m rename to matlab/runstuff_BNPcovreg.m index b1cb711..4642c7a 100644 --- a/runstuff_BNPcovreg.m +++ b/matlab/runstuff_BNPcovreg.m @@ -1,7 +1,7 @@ %% SCRIPT FOR RUNNING BASIC SAMPLER ON SIMULATED DATA EXAMPLE %% % Load simulated data: -load simData.mat +load simData % y is a p x N matrix specifying the N observations of the p-dimensional data [p N] = size(y); diff --git a/runstuff_varinds_flu.m b/matlab/runstuff_varinds_flu.m similarity index 100% rename from runstuff_varinds_flu.m rename to matlab/runstuff_varinds_flu.m diff --git a/simData.mat b/matlab/simData.mat similarity index 100% rename from simData.mat rename to matlab/simData.mat diff --git a/utilities/SIMplots.m b/matlab/utilities/SIMplots.m similarity index 100% rename from utilities/SIMplots.m rename to matlab/utilities/SIMplots.m diff --git a/utilities/calculate_hpd.m b/matlab/utilities/calculate_hpd.m similarity index 100% rename from utilities/calculate_hpd.m rename to matlab/utilities/calculate_hpd.m diff --git a/python/BNP_covreg.py b/python/BNP_covreg.py new file mode 100644 index 0000000..fb19a1c --- /dev/null +++ b/python/BNP_covreg.py @@ -0,0 +1,438 @@ +# Autogenerated with SMOP version +# /usr/local/bin/smop BNP_covreg.m + +from __future__ import division +try: + from runtime import * +except ImportError: + from smop.runtime import * + +def BNP_covreg_(y=None,prior_params=None,settings=None,restart=None,true_params=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 5-[y,prior_params,settings,restart,true_params].count(None)+len(args) + + T0=1 + Tf=1 + n_adapt=1000 + temp=copy_(T0) + empBayes=1 + p,N=size_(y,nargout=2) + if exist_(char('true_params'),char('var')): + cov_true_diag=zeros_(p,N) + for tt in arange_(1,N).reshape(-1): + cov_true_diag[:,tt]=diag_(true_params.cov_true(arange_(),arange_(),tt)) + if not exist_(char('restart'),char('var')): + restart=0 + y=init_y_(y,settings,true_params) + sample_K_flag=settings.sample_K_flag + latent_mean=settings.latent_mean + k=settings.k + L=settings.L + Niter=settings.Niter + trial=settings.trial + if not restart: + Stats[1:settings.saveEvery / settings.storeEvery]=struct_(char('zeta'),zeros_(L,k,N),char('psi'),zeros_(k,N),char('invSig_vec'),zeros_(1,p),char('theta'),zeros_(p,L),char('eta'),zeros_(k,N),char('phi'),zeros_(p,L),char('tau'),zeros_(1,L),char('K_ind'),0,char('y_heldout'),zeros_(1,sum_(inds2impute_(arange_())))) + store_counter=1 + delta=zeros_(1,L) + delta[1]=gamrnd_(prior_params.hypers.a1,1) + delta[2:L]=gamrnd_(prior_params.hypers.a2 * ones_(1,L - 1),1) + tau=exp_(cumsum_(log_(delta))) + phi=gamrnd_(prior_params.hypers.a_phi * ones_(p,L),1) / prior_params.hypers.b_phi + theta=zeros_(p,L) + for pp in arange_(1,p).reshape(-1): + theta[pp,:]=chol_(diag_(1.0 / (phi[pp,:].dot(tau)))).T * randn_(L,1) + invSig_vec=gamrnd_(prior_params.sig.a_sig * ones_(1,p),1) / prior_params.sig.b_sig + psi=zeros_(k,N) + if empBayes: + xi=sample_xi_init_(y,invSig_vec,psi,temp) + else: + xi=randn_(k,N) + eta=psi + xi + if sample_K_flag == 1 or sample_K_flag == 2: + Pk=cumsum_(prior_params.K.c_prior) + K_ind=1 + sum_(Pk[end()] * rand_(1) > Pk) + else: + K_ind=1 + if empBayes: + zeta=zeros_(L,k,N) + for ii in arange_(1,10).reshape(-1): + zeta,Sig_est=initialize_zeta_(zeta,y,theta,invSig_vec,nargout=2) + zeta=sample_zeta_(y,theta,eta,invSig_vec,zeta,prior_params.K.invK(arange_(),arange_(),K_ind),temp) + theta=sample_theta_(y,eta,invSig_vec,zeta,phi,tau,temp) + invSig_vec=sample_sig_(y,theta,eta,zeta,prior_params.sig,temp) + else: + zeta=zeros_(L,k,N) + cholK=chol_(prior_params.K.K) + for ll in arange_(1,L).reshape(-1): + for kk in arange_(1,k).reshape(-1): + zeta[ll,kk,:]=cholK.T * randn_(N,1) + zeta=sample_zeta_(y,theta,eta,invSig_vec,zeta,prior_params.K.invK(arange_(),arange_(),K_ind),temp) + if latent_mean: + psi=sample_psi_margxi_(y,theta,invSig_vec,zeta,psi,prior_params.K.invK(arange_(),arange_(),K_ind),temp) + else: + psi=zeros_(size_(xi)) + xi=sample_xi_(y,theta,invSig_vec,zeta,psi,temp) + eta=psi + xi + if not isempty_(inds2impute): + y=sample_y_(y,theta,invSig_vec,zeta,psi,inds2impute) + if not exist_(settings.saveDir,char('file')): + mkdir_(settings.saveDir) + if isfield_(settings,char('filename')): + settings_filename=strcat_(settings.saveDir,char('/'),settings.filename,char('_info4trial'),num2str_(trial)) + init_stats_filename=strcat_(settings.saveDir,char('/'),settings.filename,char('initialStats_trial'),num2str_(trial)) + else: + settings_filename=strcat_(settings.saveDir,char('/info4trial'),num2str_(trial)) + init_stats_filename=strcat_(settings.saveDir,char('/initialStats_trial'),num2str_(trial)) + if nargin > 3: + save_(settings_filename,char('y'),char('settings'),char('prior_params'),char('true_params')) + else: + save_(settings_filename,char('y'),char('settings'),char('prior_params')) + save_(init_stats_filename,char('zeta'),char('psi'),char('eta'),char('theta'),char('invSig_vec'),char('phi'),char('tau'),char('K_ind')) + num_iters=1 + nstart=1 + else: + load_([settings.saveDir,char('/BNP_covreg_stats'),char('iter'),num2str_(settings.lastIter),char('trial'),num2str_(settings.trial)]) + theta=Stats_(end()).theta + eta=Stats_(end()).eta + zeta=Stats_(end()).zeta + phi=Stats_(end()).phi + tau=Stats_(end()).tau + psi=Stats_(end()).psi + nstart=settings.lastIter + 1 + num_iters=1 + store_counter=1 + for nn in arange_(nstart,Niter).reshape(-1): + invSig_vec=sample_sig_(y,theta,eta,zeta,prior_params.sig,temp) + phi,tau=sample_hypers_(theta,phi,tau,prior_params.hypers,nargout=2) + theta=sample_theta_(y,eta,invSig_vec,zeta,phi,tau,temp) + if sample_K_flag == 1: + K_ind=sample_K_marg_zeta_(y,theta,eta,invSig_vec,prior_params.K,K_ind) + else: + if sample_K_flag == 2: + K_ind=sample_K_cond_zeta_(zeta,prior_params.K) + else: + K_ind=1 + if latent_mean: + psi=sample_psi_margxi_(y,theta,invSig_vec,zeta,psi,prior_params.K.invK(arange_(),arange_(),K_ind),temp) + else: + psi=zeros_(size_(xi)) + xi=sample_xi_(y,theta,invSig_vec,zeta,psi,temp) + eta=psi + xi + for ii in arange_(1,num_iters).reshape(-1): + zeta=sample_zeta_(y,theta,eta,invSig_vec,zeta,prior_params.K.invK(arange_(),arange_(),K_ind),temp) + if not isempty_(inds2impute): + y=sample_y_(y,theta,invSig_vec,zeta,psi,inds2impute) + num_iters=1 + if nn <= n_adapt: + temp=T0 * ((Tf / T0) ** (nn / n_adapt)) + if rem_(nn,settings.storeEvery) == 0 and nn >= settings.saveMin: + Stats[store_counter].zeta=zeta + Stats[store_counter].eta=eta + Stats[store_counter].theta=theta + Stats[store_counter].invSig_vec=invSig_vec + Stats[store_counter].psi=psi + Stats[store_counter].phi=phi + Stats[store_counter].tau=tau + Stats[store_counter].K_ind=K_ind + y_tmp=y[inds2impute] + Stats[store_counter].y_heldout=y_tmp[:] + store_counter=store_counter + 1 + if rem_(nn,settings.saveEvery) == 0: + if isfield_(settings,char('filename')): + filename=strcat_(settings.saveDir,char('/'),settings.filename,char('iter'),num2str_(nn),char('trial'),num2str_(settings.trial)) + else: + filename=strcat_(settings.saveDir,char('/BNP_covreg_stats'),char('iter'),num2str_(nn),char('trial'),num2str_(settings.trial)) + save_(filename,char('Stats')) + store_counter=1 + if not rem_(nn,100): + display_([char('Iter: '),num2str_(nn),char(', K_ind: '),num2str_(K_ind)]) + if exist_(char('true_params'),char('var')): + cov_est=zeros_(p,N) + for tt in arange_(1,N).reshape(-1): + cov_est[:,tt]=diag_(theta * zeta[:,:,tt] * zeta[:,:,tt].T * theta.T + diag_(1.0 / invSig_vec)) + plot_(cov_true_diag.T,char('LineWidth'),2) + hold_(char('on')) + plot_(cov_est.T,char('--'),char('LineWidth'),2) + hold_(char('off')) + drawnow + return +def sample_zeta_(y=None,theta=None,eta=None,invSig_vec=None,zeta=None,invK=None,temp=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 7-[y,theta,eta,invSig_vec,zeta,invK,temp].count(None)+len(args) + + p,L=size_(theta,nargout=2) + k,N=size_(eta,nargout=2) + mu_tot=zeros_(p,N) + for nn in arange_(1,N).reshape(-1): + mu_tot[:,nn]=theta * zeta[:,:,nn] * eta[:,nn] + for ll in arange_(1,L).reshape(-1): + theta_ll=theta[:,ll] + for kk in randperm_(k).reshape(-1): + eta_kk=eta[kk,:] + zeta_ll_kk=squeeze_(zeta[ll,kk,:]).T + mu_tot=mu_tot - theta_ll[:,ones_(1,N)].dot(eta_kk[ones_(p,1),:]).dot(zeta_ll_kk[ones_(p,1),:]) + A_lk_invSig_A_lk=diag_((eta[kk,:] ** 2) * ((theta[:,ll] ** 2).T * invSig_vec.T)) + theta_tmp=theta[:,ll].T.dot(invSig_vec) + ytilde=y - mu_tot + theta_lk=eta[kk,:].T.dot((theta_tmp * ytilde).T) + cholSig_lk_trans=numpy.linalg.solve(chol_(invK + A_lk_invSig_A_lk),eye_(N)) + m_lk=cholSig_lk_trans * (cholSig_lk_trans.T * theta_lk) + zeta[ll,kk,:]=m_lk + temp * cholSig_lk_trans * randn_(N,1) + zeta_ll_kk=squeeze_(zeta[ll,kk,:]).T + mu_tot=mu_tot + theta_ll[:,ones_(1,N)].dot(eta_kk[ones_(p,1),:]).dot(zeta_ll_kk[ones_(p,1),:]) + return zeta +def sample_psi_margxi_(y=None,theta=None,invSig_vec=None,zeta=None,psi=None,invK=None,temp=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 7-[y,theta,invSig_vec,zeta,psi,invK,temp].count(None)+len(args) + + p,L=size_(theta,nargout=2) + k,N=size_(psi,nargout=2) + Sigma_0=diag_(1.0 / invSig_vec) + mu_tot=zeros_(p,N) + Omega=zeros_(p,k,N) + invOmegaOmegaSigma0=zeros_(p,p,N) + for nn in arange_(1,N).reshape(-1): + Omega[:,:,nn]=theta * zeta[:,:,nn] + temp=numpy.linalg.solve((Omega[:,:,nn] * Omega[:,:,nn].T + Sigma_0),eye_(N)) + OmegaInvOmegaOmegaSigma0[:,:,nn]=Omega[:,:,nn].T * temp + mu_tot[:,nn]=Omega[:,:,nn] * psi[:,nn] + if sum_(sum_(mu_tot)) == 0: + numTotIters=100 + else: + numTotIters=5 + for numIter in arange_(1,numTotIters).reshape(-1): + for kk in randperm_(k).reshape(-1): + Omega_kk=squeeze_(Omega[:,kk,:]) + psi_kk=psi[kk,:] + mu_tot=mu_tot - Omega_kk.dot(psi_kk[ones_(p,1),:]) + theta_k=diag_(squeeze_(OmegaInvOmegaOmegaSigma0[kk,:,:]).T * (y - mu_tot)) + Ak_invSig_Ak=diag_(squeeze_(OmegaInvOmegaOmegaSigma0[kk,:,:]).T * Omega_kk) + cholSig_k_trans=numpy.linalg.solve(chol_(invK + diag_(Ak_invSig_Ak)),eye_(N)) + m_k=cholSig_k_trans * (cholSig_k_trans.T * theta_k) + psi[kk,:]=m_k + temp * cholSig_k_trans * randn_(N,1) + psi_kk=psi[kk,:] + mu_tot=mu_tot + Omega_kk.dot(psi_kk[ones_(p,1),:]) + return psi +def sample_xi_(y=None,theta=None,invSig_vec=None,zeta=None,psi=None,temp=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 6-[y,theta,invSig_vec,zeta,psi,temp].count(None)+len(args) + + p,N=size_(y,nargout=2) + L,k=size_(zeta[:,:,1],nargout=2) + invSigMat=invSig_vec[ones_(k,1),:] + xi=zeros_(k,N) + for nn in arange_(1,N).reshape(-1): + theta_zeta_n=theta * zeta[:,:,nn] + y_tilde_n=y[:,nn] - theta_zeta_n * psi[:,nn] + zeta_theta_invSig=theta_zeta_n.T.dot(invSigMat) + cholSig_xi_n_trans=numpy.linalg.solve(chol_(eye_(k) + zeta_theta_invSig * theta_zeta_n),eye_(k)) + m_xi_n=cholSig_xi_n_trans * (cholSig_xi_n_trans.T * (zeta_theta_invSig * y_tilde_n)) + xi[:,nn]=m_xi_n + temp * cholSig_xi_n_trans * randn_(k,1) + return xi +def sample_theta_(y=None,eta=None,invSig_vec=None,zeta=None,phi=None,tau=None,temp=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 7-[y,eta,invSig_vec,zeta,phi,tau,temp].count(None)+len(args) + + p,N=size_(y,nargout=2) + L=size_(zeta,1) + theta=zeros_(p,L) + eta_tilde=zeros_(L,N) + for nn in arange_(1,N).reshape(-1): + eta_tilde[:,nn]=zeta[:,:,nn] * eta[:,nn] + eta_tilde=eta_tilde.T + for pp in arange_(1,p).reshape(-1): + chol_Sig_theta_p_trans=numpy.linalg.solve(chol_(diag_(phi[pp,:].dot(tau)) + invSig_vec[pp] * (eta_tilde.T * eta_tilde)),eye_(L)) + m_theta_p=invSig_vec[pp] * (chol_Sig_theta_p_trans * chol_Sig_theta_p_trans.T) * (eta_tilde.T * y[pp,:].T) + theta[pp,:]=m_theta_p + temp * chol_Sig_theta_p_trans * randn_(L,1) + return theta +def sample_sig_(y=None,theta=None,eta=None,zeta=None,prior_params=None,temp=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 6-[y,theta,eta,zeta,prior_params,temp].count(None)+len(args) + + p,N=size_(y,nargout=2) + a_sig=prior_params.a_sig + b_sig=prior_params.b_sig + invSig_vec=zeros_(1,p) + for pp in arange_(1,p).reshape(-1): + sq_err=0 + for nn in arange_(1,N).reshape(-1): + sq_err=sq_err + (y[pp,nn] - theta[pp,:] * zeta[:,:,nn] * eta[:,nn]) ** 2 + a_temp=1 + (a_sig + 0.5 * N - 1) / temp + b_temp=(b_sig + 0.5 * sq_err) / temp + invSig_vec[pp]=gamrnd_(a_temp,1) / b_temp + return invSig_vec +def sample_hypers_(theta=None,phi=None,tau=None,prior_params=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 4-[theta,phi,tau,prior_params].count(None)+len(args) + + p,L=size_(theta,nargout=2) + a1=prior_params.a1 + a2=prior_params.a2 + a_phi=prior_params.a_phi + b_phi=prior_params.b_phi + a=matlabarray([a1,a2 * ones_(1,L - 1)]) + delta=exp_([log_(tau[1]),diff_(log_(tau))]) + for numIter in arange_(1,50).reshape(-1): + phi=gamrnd_((a_phi + 0.5) * ones_(p,L),1) / (b_phi + 0.5 * tau[ones_(1,p),:].dot((theta ** 2))) + sum_phi_theta=sum_(phi.dot((theta ** 2)),1) + for hh in arange_(1,L).reshape(-1): + tau_hh=exp_(cumsum_(log_(delta))).dot([zeros_(1,hh - 1),ones_(1,L - hh + 1) / delta[hh]]) + delta[hh]=gamrnd_(a[hh] + 0.5 * p * (L - hh + 1),1) / (1 + 0.5 * sum_(tau_hh.dot(sum_phi_theta))) + tau=exp_(cumsum_(log_(delta))) + return phi,tau +def sample_K_marg_zeta_(y=None,theta=None,eta=None,invSig_vec=None,prior_params=None,K_ind=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 6-[y,theta,eta,invSig_vec,prior_params,K_ind].count(None)+len(args) + + k,N=size_(eta,nargout=2) + p,L=size_(theta,nargout=2) + c_prior=prior_params.c_prior + K=prior_params.K + grid_size=length_(c_prior) + Pk=- Inf * ones_(1,grid_size) + Pk_tmp=- Inf * ones_(1,grid_size) + nbhd=copy_(Inf) + nbhd_vec=matlabarray([arange_(max_(K_ind - nbhd,1),min_(K_ind + nbhd,grid_size))]) + tmp_mat_init=diag_(repmat_(1.0 / invSig_vec,1,N)) + if p < sqrt_(N): + for cc in nbhd_vec.reshape(-1): + Kc=K[:,:,cc] + tmp_mat=copy_(tmp_mat_init) + for ll in arange_(1,L).reshape(-1): + theta_ll=theta[:,ll] + for kk in arange_(1,k).reshape(-1): + eta_kk=eta[kk,:] + tmp_mat=tmp_mat + kron_((eta_kk * eta_kk.T).dot(Kc),theta_ll * theta_ll.T) + Pk[cc]=normpdfln_(y[:],zeros_(p * N,1),tmp_mat) + else: + for cc in nbhd_vec.reshape(-1): + Kc=K[:,:,cc] + Kc=sparse_(Kc) + tmp_mat=copy_(tmp_mat_init) + for ll in arange_(1,L).reshape(-1): + theta_ll=theta[:,ll] + for kk in arange_(1,k).reshape(-1): + eta_kk=eta[kk,:] + eta_mat_kk=diag_(eta_kk) + eta_theta_ll_kk=kron_(eta_mat_kk,theta_ll) + tmp_mat=tmp_mat + eta_theta_ll_kk * Kc * eta_theta_ll_kk.T + Pk[cc]=normpdfln_(y[:],zeros_(p * N,1),tmp_mat) + nbhd_mask=zeros_(1,length_(prior_params.c_prior)) + nbhd_mask[nbhd_vec]=1 + Pk=Pk + log_(prior_params.c_prior).dot(nbhd_mask) + Pk=cumsum_(exp_(Pk - max_(Pk))) + K_ind=1 + sum_(Pk[end()] * rand_(1) > Pk) + return K_ind +def sample_K_cond_zeta_(zeta=None,prior_params=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 2-[zeta,prior_params].count(None)+len(args) + + L,k,N=size_(zeta,nargout=3) + c_prior=prior_params.c_prior + invK=prior_params.invK + zeta_tmp=reshape_(zeta,[L * k,N]).T + Pk=zeros_(1,length_(c_prior)) + for ii in arange_(1,length_(c_prior)).reshape(-1): + Pk[ii]=- sum_(diag_(zeta_tmp.T * (invK[:,:,ii] * zeta_tmp))) + Pk=Pk - 0.5 * (L * k) * prior_params.logdetK + log_(prior_params.c_prior) + Pk=cumsum_(exp_(Pk - max_(Pk))) + K_ind=1 + sum_(Pk[end()] * rand_(1) > Pk) + return K_ind +def sample_y_(y=None,theta=None,invSig_vec=None,zeta=None,psi=None,inds2impute=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 6-[y,theta,invSig_vec,zeta,psi,inds2impute].count(None)+len(args) + + p,N=size_(y,nargout=2) + inds_vec=matlabarray([arange_(1,p)]) + for nn in arange_(1,N).reshape(-1): + Sigma_nn=theta * zeta[:,:,nn] * zeta[:,:,nn].T * theta.T + diag_(1.0 / invSig_vec) + mu_nn=theta * zeta[:,:,nn] * psi[:,nn] + inds2impute_nn=inds_vec[inds2impute[:,nn]] + J=length_(inds2impute_nn) + reg_inds_nn=setdiff_(inds_vec,inds2impute_nn) + ordered_inds=matlabarray([inds2impute_nn,reg_inds_nn]) + ordered_cov_rows=zeros_(p) + ordered_mean=zeros_(p,1) + for jj in arange_(1,p).reshape(-1): + ordered_cov_rows[jj,:]=Sigma_nn[ordered_inds[jj],:] + ordered_mean[jj]=mu_nn[ordered_inds[jj]] + ordered_cov=zeros_(size_(ordered_cov_rows)) + for jj in arange_(1,p).reshape(-1): + ordered_cov[:,jj]=ordered_cov_rows[:,ordered_inds[jj]] + Sig11=ordered_cov[1:J,1:J] + Sig12=ordered_cov[1:J,J + 1:p] + Sig12invSig22=Sig12 / ordered_cov[J + 1:p,J + 1:p] + mu1=ordered_mean[1:J] + mu2=ordered_mean[J + 1:p] + reg_mean=mu1 + Sig12invSig22 * (y[reg_inds_nn,nn] - mu2) + reg_cov=Sig11 - Sig12invSig22 * Sig12.T + y[inds2impute_nn,nn]=reg_mean + chol_(reg_cov).T * randn_(J,1) + return y +def init_y_(y=None,settings=None,true_params=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 3-[y,settings,true_params].count(None)+len(args) + + inds2impute=settings.inds2impute + p,N=size_(y,nargout=2) + y[inds2impute]=0 + x=linspace_(- 1,1,5) + mu=0 + sig=1 + tmp=normpdf_(x,mu,sig) + y_conved=copy_(y) + for pp in arange_(1,p).reshape(-1): + y_conved[pp,:]=conv_(y[pp,:],tmp,char('same')) + y[inds2impute]=y_conved[inds2impute] + return y +def initialize_zeta_(zeta=None,y=None,theta=None,invSig_vec=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 4-[zeta,y,theta,invSig_vec].count(None)+len(args) + + p,N=size_(y,nargout=2) + L,k,tmp=size_(zeta,nargout=3) + Sig_mat=diag_(1.0 / invSig_vec) + Sig_est=zeros_(p,N) + x=floor_(linspace_(1,N,20)) + zeta_knots=zeros_(L,k,length_(x)) + for ii in arange_(1,length_(x)).reshape(-1): + inds_ii=matlabarray([arange_(max_(1,x[ii] - N / 10),min_(N,x[ii] + N / 10))]) + cov_est_ii=cov_(y[:,inds_ii].T) + 0 * ones_(p) + C=chol_(cov_est_ii).T + C_ii=C[:,1:k] + zeta_knots[:,:,ii]=numpy.linalg.solve(theta,C_ii) + Sig_est[:,ii]=diag_(C_ii * C_ii.T) + zeta_spline=spline_(x,zeta_knots) + zeta=ppval_([arange_(1,N)],zeta_spline) + return zeta,Sig_est +def sample_xi_init_(y=None,invSig_vec=None,psi=None,temp=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 4-[y,invSig_vec,psi,temp].count(None)+len(args) + + p,N=size_(y,nargout=2) + k,N=size_(psi,nargout=2) + invSigMat=invSig_vec[ones_(k,1),:] + xi=zeros_(k,N) + for nn in arange_(1,N).reshape(-1): + inds_nn=matlabarray([arange_(max_(1,nn - N / 10),min_(N,nn + N / 10))]) + cov_est_nn=cov_(y[:,inds_nn].T) + 0 * ones_(p) + C=chol_(cov_est_nn).T + C_nn=C[:,1:k] + theta_zeta_n=copy_(C_nn) + y_tilde_n=y[:,nn] - theta_zeta_n * psi[:,nn] + zeta_theta_invSig=theta_zeta_n.T.dot(invSigMat) + cholSig_xi_n_trans=numpy.linalg.solve(chol_(eye_(k) + zeta_theta_invSig * theta_zeta_n),eye_(k)) + m_xi_n=cholSig_xi_n_trans * (cholSig_xi_n_trans.T * (zeta_theta_invSig * y_tilde_n)) + xi[:,nn]=m_xi_n + temp * cholSig_xi_n_trans * randn_(k,1) + return xi +def normpdfln_(x=None,mu=None,Sigma=None,iSigma=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 4-[x,mu,Sigma,iSigma].count(None)+len(args) + + if nargin < 4: + q=size_(Sigma,1) + iSigma=numpy.linalg.solve(Sigma,eye_(q)) + logdet_iSigma=- 2 * sum_(log_(diag_(chol_(Sigma)))) + else: + logdet_iSigma=2 * sum_(log_(diag_(chol_(iSigma)))) + logp=0.5 * logdet_iSigma - 0.5 * (x - mu).T * (iSigma * (x - mu)) + return logp diff --git a/python/BNP_covreg_varinds.py b/python/BNP_covreg_varinds.py new file mode 100644 index 0000000..e90b37c --- /dev/null +++ b/python/BNP_covreg_varinds.py @@ -0,0 +1,324 @@ +# Autogenerated with SMOP version +# /usr/local/bin/smop BNP_covreg_varinds.m -o ../python/BNP_covreg_varinds.py + +from __future__ import division +try: + from runtime import * +except ImportError: + from smop.runtime import * + +def BNP_covreg_varinds_(y=None,prior_params=None,settings=None,restart=None,true_params=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 5-[y,prior_params,settings,restart,true_params].count(None)+len(args) + + p,N=size_(y,nargout=2) + if exist_(char('true_params'),char('var')): + cov_true_diag=zeros_(p,N) + for tt in arange_(1,N).reshape(-1): + cov_true_diag[:,tt]=diag_(true_params.cov_true(arange_(),arange_(),tt)) + if not exist_(char('restart'),char('var')): + restart=0 + sample_K_flag=settings.sample_K_flag + latent_mean=settings.latent_mean + inds_y=settings.inds_y + y[not inds_y]=0 + k=settings.k + L=settings.L + Niter=settings.Niter + trial=settings.trial + if not restart: + Stats[1:settings.saveEvery / settings.storeEvery]=struct_(char('zeta'),zeros_(L,k,N),char('psi'),zeros_(k,N),char('invSig_vec'),zeros_(1,p),char('theta'),zeros_(p,L),char('eta'),zeros_(k,N),char('phi'),zeros_(p,L),char('tau'),zeros_(1,L),char('K_ind'),0) + store_counter=1 + delta=zeros_(1,L) + delta[1]=gamrnd_(prior_params.hypers.a1,1) + delta[2:L]=gamrnd_(prior_params.hypers.a2 * ones_(1,L - 1),1) + tau=exp_(cumsum_(log_(delta))) + phi=gamrnd_(prior_params.hypers.a_phi * ones_(p,L),1) / prior_params.hypers.b_phi + theta=zeros_(p,L) + for pp in arange_(1,p).reshape(-1): + theta[pp,:]=chol_(diag_(1.0 / (phi[pp,:].dot(tau)))).T * randn_(L,1) + xi=randn_(k,N) + psi=zeros_(k,N) + eta=psi + xi + invSig_vec=gamrnd_(prior_params.sig.a_sig * ones_(1,p),1) / prior_params.sig.b_sig + if sample_K_flag == 1 or sample_K_flag == 2: + Pk=cumsum_(prior_params.K.c_prior) + K_ind=1 + sum_(Pk[end()] * rand_(1) > Pk) + else: + K_ind=1 + zeta=zeros_(L,k,N) + zeta=sample_zeta_(y,theta,eta,invSig_vec,zeta,prior_params.K.invK(arange_(),arange_(),K_ind),inds_y) + if not exist_(settings.saveDir,char('file')): + mkdir_(settings.saveDir) + if isfield_(settings,char('filename')): + settings_filename=strcat_(settings.saveDir,char('/'),settings.filename,char('_info4trial'),num2str_(trial)) + init_stats_filename=strcat_(settings.saveDir,char('/'),settings.filename,char('initialStats_trial'),num2str_(trial)) + else: + settings_filename=strcat_(settings.saveDir,char('/info4trial'),num2str_(trial)) + init_stats_filename=strcat_(settings.saveDir,char('/initialStats_trial'),num2str_(trial)) + if nargin > 4: + save_(settings_filename,char('y'),char('settings'),char('prior_params'),char('true_params')) + else: + save_(settings_filename,char('y'),char('settings'),char('prior_params')) + save_(init_stats_filename,char('zeta'),char('psi'),char('eta'),char('theta'),char('invSig_vec'),char('phi'),char('tau'),char('K_ind')) + num_iters=1 + nstart=1 + else: + load_([settings.saveDir,char('/BNP_covreg_stats'),char('iter'),num2str_(settings.lastIter),char('trial'),num2str_(settings.trial)]) + theta=Stats_(end()).theta + eta=Stats_(end()).eta + zeta=Stats_(end()).zeta + phi=Stats_(end()).phi + tau=Stats_(end()).tau + psi=Stats_(end()).psi + nstart=settings.lastIter + 1 + num_iters=1 + store_counter=1 + for nn in arange_(nstart,Niter).reshape(-1): + invSig_vec=sample_sig_(y,theta,eta,zeta,prior_params.sig,inds_y) + phi,tau=sample_hypers_(theta,phi,tau,prior_params.hypers,nargout=2) + theta=sample_theta_(y,eta,invSig_vec,zeta,phi,tau,inds_y) + if sample_K_flag == 1: + error_(char('need to code up K sampling for inds_y')) + else: + if sample_K_flag == 2: + error_(char('need to code up K sampling for inds_y')) + else: + K_ind=1 + if latent_mean: + psi=sample_psi_margxi_(y,theta,invSig_vec,zeta,psi,prior_params.K.invK(arange_(),arange_(),K_ind),inds_y) + else: + psi=zeros_(size_(xi)) + xi=sample_xi_(y,theta,invSig_vec,zeta,psi,inds_y) + eta=psi + xi + for ii in arange_(1,num_iters).reshape(-1): + zeta=sample_zeta_(y,theta,eta,invSig_vec,zeta,prior_params.K.invK(arange_(),arange_(),K_ind),inds_y) + num_iters=1 + if rem_(nn,settings.storeEvery) == 0 and nn >= settings.saveMin: + Stats[store_counter].zeta=zeta + Stats[store_counter].eta=eta + Stats[store_counter].theta=theta + Stats[store_counter].invSig_vec=invSig_vec + Stats[store_counter].psi=psi + Stats[store_counter].phi=phi + Stats[store_counter].tau=tau + Stats[store_counter].K_ind=K_ind + store_counter=store_counter + 1 + if rem_(nn,settings.saveEvery) == 0: + if isfield_(settings,char('filename')): + filename=strcat_(settings.saveDir,char('/'),settings.filename,char('iter'),num2str_(nn),char('trial'),num2str_(settings.trial)) + else: + filename=strcat_(settings.saveDir,char('/BNP_covreg_stats'),char('iter'),num2str_(nn),char('trial'),num2str_(settings.trial)) + save_(filename,char('Stats')) + store_counter=1 + if not rem_(nn,100): + display_([char('Iter: '),num2str_(nn),char(', K_ind: '),num2str_(K_ind)]) + if exist_(char('true_params'),char('var')): + cov_est=zeros_(p,N) + for tt in arange_(1,N).reshape(-1): + cov_est[:,tt]=diag_(theta * zeta[:,:,tt] * zeta[:,:,tt].T * theta.T + diag_(1.0 / invSig_vec)) + plot_(cov_true_diag.T,char('LineWidth'),2) + hold_(char('on')) + plot_(cov_est.T,char('--'),char('LineWidth'),2) + hold_(char('off')) + drawnow + return +def sample_zeta_(y=None,theta=None,eta=None,invSig_vec=None,zeta=None,invK=None,inds_y=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 7-[y,theta,eta,invSig_vec,zeta,invK,inds_y].count(None)+len(args) + + p,L=size_(theta,nargout=2) + k,N=size_(eta,nargout=2) + mu_tot=zeros_(p,N) + error_nn=zeros_(L,N) + for nn in arange_(1,N).reshape(-1): + mu_tot[:,nn]=theta * zeta[:,:,nn] * eta[:,nn] + error_nn[:,nn]=(theta ** 2).T * (invSig_vec.T.dot((1 - inds_y[:,nn]))) + numiter=1 + for nn in arange_(1,numiter).reshape(-1): + for ll in arange_(1,L).reshape(-1): + theta_ll=theta[:,ll] + for kk in randperm_(k).reshape(-1): + eta_kk=eta[kk,:] + zeta_ll_kk=squeeze_(zeta[ll,kk,:]).T + mu_tot=mu_tot - theta_ll[:,ones_(1,N)].dot(eta_kk[ones_(p,1),:]).dot(zeta_ll_kk[ones_(p,1),:]) + A_lk_invSig_A_lk=(eta[kk,:] ** 2).dot(((theta[:,ll] ** 2).T * invSig_vec.T - error_nn[ll,:])) + theta_tmp=theta[:,ll].T.dot(invSig_vec) + ytilde=(y - mu_tot).dot(inds_y) + theta_lk=eta[kk,:].T.dot((theta_tmp * ytilde).T) + cholSig_lk_trans=numpy.linalg.solve(chol_(invK + diag_(A_lk_invSig_A_lk)),eye_(N)) + m_lk=cholSig_lk_trans * (cholSig_lk_trans.T * theta_lk) + zeta[ll,kk,:]=m_lk + cholSig_lk_trans * randn_(N,1) + zeta_ll_kk=squeeze_(zeta[ll,kk,:]).T + mu_tot=mu_tot + theta_ll[:,ones_(1,N)].dot(eta_kk[ones_(p,1),:]).dot(zeta_ll_kk[ones_(p,1),:]) + return zeta +def sample_psi_margxi_(y=None,theta=None,invSig_vec=None,zeta=None,psi=None,invK=None,inds_y=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 7-[y,theta,invSig_vec,zeta,psi,invK,inds_y].count(None)+len(args) + + p,L=size_(theta,nargout=2) + k,N=size_(psi,nargout=2) + Sigma_0=diag_(1.0 / invSig_vec) + mu_tot=zeros_(p,N) + Omega=zeros_(p,k,N) + OmegaInvOmegaOmegaSigma0=zeros_(k,p,N) + for nn in arange_(1,N).reshape(-1): + Omega[inds_y[:,nn],:,nn]=theta[inds_y[:,nn],:] * zeta[:,:,nn] + temp=numpy.linalg.solve((Omega[inds_y[:,nn],:,nn] * Omega[inds_y[:,nn],:,nn].T + Sigma_0[inds_y[:,nn],inds_y[:,nn]]),eye_(sum_(inds_y[:,nn]))) + OmegaInvOmegaOmegaSigma0[:,inds_y[:,nn],nn]=Omega[inds_y[:,nn],:,nn].T * temp + mu_tot[:,nn]=Omega[:,:,nn] * psi[:,nn] + if sum_(sum_(mu_tot)) == 0: + numTotIters=50 + else: + numTotIters=5 + for numIter in arange_(1,numTotIters).reshape(-1): + for kk in randperm_(k).reshape(-1): + Omega_kk=squeeze_(Omega[:,kk,:]) + psi_kk=psi[kk,:] + mu_tot=mu_tot - Omega_kk.dot(psi_kk[ones_(p,1),:]) + theta_k=diag_(squeeze_(OmegaInvOmegaOmegaSigma0[kk,:,:]).T * (y - mu_tot)) + Ak_invSig_Ak=diag_(squeeze_(OmegaInvOmegaOmegaSigma0[kk,:,:]).T * Omega_kk) + cholSig_k_trans=numpy.linalg.solve(chol_(invK + diag_(Ak_invSig_Ak)),eye_(N)) + m_k=cholSig_k_trans * (cholSig_k_trans.T * theta_k) + psi[kk,:]=m_k + cholSig_k_trans * randn_(N,1) + psi_kk=psi[kk,:] + mu_tot=mu_tot + Omega_kk.dot(psi_kk[ones_(p,1),:]) + return psi +def sample_xi_(y=None,theta=None,invSig_vec=None,zeta=None,psi=None,inds_y=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 6-[y,theta,invSig_vec,zeta,psi,inds_y].count(None)+len(args) + + p,N=size_(y,nargout=2) + L,k=size_(zeta[:,:,1],nargout=2) + xi=zeros_(k,N) + for nn in arange_(1,N).reshape(-1): + theta_zeta_n=theta * zeta[:,:,nn] + y_tilde_n=y[:,nn] - theta_zeta_n * psi[:,nn] + invSigMat=invSig_vec.dot(inds_y[:,nn].T) + invSigMat=invSigMat[ones_(k,1),:] + zeta_theta_invSig=theta_zeta_n.T.dot(invSigMat) + cholSig_xi_n_trans=numpy.linalg.solve(chol_(eye_(k) + zeta_theta_invSig * theta_zeta_n),eye_(k)) + m_xi_n=cholSig_xi_n_trans * (cholSig_xi_n_trans.T * (zeta_theta_invSig * y_tilde_n)) + xi[:,nn]=m_xi_n + cholSig_xi_n_trans * randn_(k,1) + return xi +def sample_theta_(y=None,eta=None,invSig_vec=None,zeta=None,phi=None,tau=None,inds_y=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 7-[y,eta,invSig_vec,zeta,phi,tau,inds_y].count(None)+len(args) + + p,N=size_(y,nargout=2) + L=size_(zeta,1) + theta=zeros_(p,L) + eta_tilde=zeros_(L,N) + for nn in arange_(1,N).reshape(-1): + eta_tilde[:,nn]=zeta[:,:,nn] * eta[:,nn] + eta_tilde=eta_tilde.T + for pp in arange_(1,p).reshape(-1): + inds_y_p=inds_y[pp,:].T + eta_tilde_p=eta_tilde.dot(inds_y_p[:,ones_(1,L)]) + chol_Sig_theta_p_trans=numpy.linalg.solve(chol_(diag_(phi[pp,:].dot(tau)) + invSig_vec[pp] * (eta_tilde_p.T * eta_tilde_p)),eye_(L)) + m_theta_p=invSig_vec[pp] * (chol_Sig_theta_p_trans * chol_Sig_theta_p_trans.T) * (eta_tilde_p.T * y[pp,:].T) + theta[pp,:]=m_theta_p + chol_Sig_theta_p_trans * randn_(L,1) + return theta +def sample_sig_(y=None,theta=None,eta=None,zeta=None,prior_params=None,inds_y=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 6-[y,theta,eta,zeta,prior_params,inds_y].count(None)+len(args) + + p,N=size_(y,nargout=2) + a_sig=prior_params.a_sig + b_sig=prior_params.b_sig + inds_vec=matlabarray([arange_(1,N)]) + invSig_vec=zeros_(1,p) + for pp in arange_(1,p).reshape(-1): + sq_err=0 + for nn in inds_vec[inds_y[pp,:]].reshape(-1): + sq_err=sq_err + (y[pp,nn] - theta[pp,:] * zeta[:,:,nn] * eta[:,nn]) ** 2 + invSig_vec[pp]=gamrnd_(a_sig + 0.5 * sum_(inds_y[pp,:]),1) / (b_sig + 0.5 * sq_err) + return invSig_vec +def sample_hypers_(theta=None,phi=None,tau=None,prior_params=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 4-[theta,phi,tau,prior_params].count(None)+len(args) + + p,L=size_(theta,nargout=2) + a1=prior_params.a1 + a2=prior_params.a2 + a_phi=prior_params.a_phi + b_phi=prior_params.b_phi + a=matlabarray([a1,a2 * ones_(1,L - 1)]) + delta=exp_([log_(tau[1]),diff_(log_(tau))]) + for numIter in arange_(1,50).reshape(-1): + phi=gamrnd_((a_phi + 0.5) * ones_(p,L),1) / (b_phi + 0.5 * tau[ones_(1,p),:].dot((theta ** 2))) + sum_phi_theta=sum_(phi.dot((theta ** 2)),1) + for hh in arange_(1,L).reshape(-1): + tau_hh=exp_(cumsum_(log_(delta))).dot([zeros_(1,hh - 1),ones_(1,L - hh + 1) / delta[hh]]) + delta[hh]=gamrnd_(a[hh] + 0.5 * p * (L - hh + 1),1) / (1 + 0.5 * sum_(tau_hh.dot(sum_phi_theta))) + tau=exp_(cumsum_(log_(delta))) + return phi,tau +def sample_K_marg_zeta_(y=None,theta=None,eta=None,invSig_vec=None,prior_params=None,K_ind=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 6-[y,theta,eta,invSig_vec,prior_params,K_ind].count(None)+len(args) + + k,N=size_(eta,nargout=2) + p,L=size_(theta,nargout=2) + c_prior=prior_params.c_prior + K=prior_params.K + grid_size=length_(c_prior) + Pk=- Inf * ones_(1,grid_size) + nbhd=copy_(Inf) + nbhd_vec=matlabarray([arange_(max_(K_ind - nbhd,1),min_(K_ind + nbhd,grid_size))]) + tmp_mat_init=diag_(repmat_(1.0 / invSig_vec,1,N)) + if p < sqrt_(N): + for cc in nbhd_vec.reshape(-1): + Kc=K[:,:,cc] + tmp_mat=copy_(tmp_mat_init) + for ll in arange_(1,L).reshape(-1): + theta_ll=theta[:,ll] + for kk in arange_(1,k).reshape(-1): + eta_kk=eta[kk,:] + tmp_mat=tmp_mat + kron_((eta_kk * eta_kk.T).dot(Kc),theta_ll * theta_ll.T) + Pk[cc]=normpdfln_(y[:],zeros_(p * N,1),tmp_mat) + else: + for cc in nbhd_vec.reshape(-1): + Kc=K[:,:,cc] + Kc=sparse_(Kc) + tmp_mat=copy_(tmp_mat_init) + for ll in arange_(1,L).reshape(-1): + theta_ll=theta[:,ll] + for kk in arange_(1,k).reshape(-1): + eta_kk=eta[kk,:] + eta_mat_kk=diag_(eta_kk) + eta_theta_ll_kk=kron_(eta_mat_kk,theta_ll) + tmp_mat=tmp_mat + eta_theta_ll_kk * Kc * eta_theta_ll_kk.T + Pk[cc]=normpdfln_(y[:],zeros_(p * N,1),tmp_mat) + nbhd_mask=zeros_(1,length_(prior_params.c_prior)) + nbhd_mask[nbhd_vec]=1 + Pk=Pk + log_(prior_params.c_prior).dot(nbhd_mask) + Pk=cumsum_(exp_(Pk - max_(Pk))) + K_ind=1 + sum_(Pk[end()] * rand_(1) > Pk) + return K_ind +def sample_K_cond_zeta_(zeta=None,prior_params=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 2-[zeta,prior_params].count(None)+len(args) + + L,k,N=size_(zeta,nargout=3) + c_prior=prior_params.c_prior + invK=prior_params.invK + zeta_tmp=reshape_(zeta,[L * k,N]).T + Pk=zeros_(1,length_(c_prior)) + for ii in arange_(1,length_(c_prior)).reshape(-1): + Pk[ii]=- sum_(diag_(zeta_tmp.T * (invK[:,:,ii] * zeta_tmp))) + Pk=Pk - 0.5 * (L * k) * prior_params.logdetK + log_(prior_params.c_prior) + Pk=cumsum_(exp_(Pk - max_(Pk))) + K_ind=1 + sum_(Pk[end()] * rand_(1) > Pk) + return K_ind +def normpdfln_(x=None,mu=None,Sigma=None,iSigma=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 4-[x,mu,Sigma,iSigma].count(None)+len(args) + + if nargin < 4: + q=size_(Sigma,1) + iSigma=numpy.linalg.solve(Sigma,eye_(q)) + logdet_iSigma=- 2 * sum_(log_(diag_(chol_(Sigma)))) + else: + logdet_iSigma=2 * sum_(log_(diag_(chol_(iSigma)))) + logp=0.5 * logdet_iSigma - 0.5 * (x - mu).T * (iSigma * (x - mu)) + return logp diff --git a/python/__init__.py b/python/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/python/flu_US.mat b/python/flu_US.mat new file mode 100644 index 0000000..3ade014 Binary files /dev/null and b/python/flu_US.mat differ diff --git a/python/runstuff_BNPcovreg.py b/python/runstuff_BNPcovreg.py new file mode 100644 index 0000000..f512e9d --- /dev/null +++ b/python/runstuff_BNPcovreg.py @@ -0,0 +1,74 @@ +# Autogenerated with SMOP version +# /usr/local/bin/smop runstuff_BNPcovreg.m -o ../python/runstuff_BNPcovreg.py + +from __future__ import division +try: + from runtime import * +except ImportError: + from smop.runtime import * + +load(char('simData')) +p,N=size(y,nargout=2) +inds_y=ones(size(y)) +inds_y=inds_y > 0 +x=[arange_(1,N)] / N +c=10 +d=1 +r=1e-05 +K=zeros(N) +for ii in arange_(1,N).reshape(-1): + for jj in arange_(1,N).reshape(-1): + dist_ii_jj=_abs(x(ii) - x(jj)) + K[ii,jj]=d * exp(- c * (dist_ii_jj ** 2)) +K=K + diag(r * ones(1,N)) +invK=inv(K) +logdetK=2 * _sum(log(diag(chol(K)))) +prior_params.K.c_prior=1 +prior_params.K.invK=invK +prior_params.K.K=K +prior_params.K.logdetK=logdetK +prior_params.sig.a_sig=1 +prior_params.sig.b_sig=0.1 +prior_params.hypers.a_phi=1.5 +prior_params.hypers.b_phi=1.5 +prior_params.hypers.a1=2 +prior_params.hypers.a2=2 +settings.L=5 +settings.k=4 +settings.Niter=10000 +settings.saveEvery=100 +settings.storeEvery=10 +settings.saveMin=1 +settings.saveDir=char('test') +settings.trial=1 +settings.sample_K_flag=3 +settings.latent_mean=1 +settings.inds_y=inds_y +restart=0 +if restart: + settings.lastIter=10000 +BNP_covreg_varinds(y,prior_params,settings,restart,true_params) +c=10 +d=1 +r=1e-05 +c_vec=copy_(c) +Nc=length(c_vec) +invK=zeros(N,N,Nc) +K=zeros(N,N,Nc) +logdetK=zeros(1,Nc) +for kk in arange_(1,Nc).reshape(-1): + K_tmp=zeros(N,N) + for ii in arange_(1,N).reshape(-1): + for jj in arange_(1,N).reshape(-1): + dist_ii_jj=_abs(x(ii) - x(jj)) + K_tmp[ii,jj]=d * exp(- c_vec(kk) * (dist_ii_jj ** 2)) + K_tmp=K_tmp + diag(r * ones(1,N)) + K[arange_(),arange_(),kk]=K_tmp + invK[arange_(),arange_(),kk]=inv(K_tmp) + logdetK[kk]=2 * _sum(log(diag(chol(K_tmp)))) +prior_params.K.c_prior=ones(1,Nc) / Nc +prior_params.K.invK=invK +prior_params.K.K=K +prior_params.K.logdetK=logdetK +settings.sample_K_flag=1 +settings.sample_K_flag=2 diff --git a/python/runstuff_varinds_flu.py b/python/runstuff_varinds_flu.py new file mode 100644 index 0000000..f3771de --- /dev/null +++ b/python/runstuff_varinds_flu.py @@ -0,0 +1,92 @@ +# Autogenerated with SMOP version +# /usr/local/bin/smop runstuff_varinds_flu.m -o ../python/runstuff_varinds_flu.py + +from __future__ import division +try: + from runtime import * +except ImportError: + from smop.runtime import * + +import os +import numpy as np +from numpy import diagonal +from numpy.linalg import inv, cholesky as chol +from .BNP_covreg_varinds import BNP_covreg_varinds_ + +FILE_DIR = os.path.dirname( os.path.abspath(__file__) ) +data = load_(open(os.path.join( FILE_DIR, 'flu_US.mat' ))) + +month_names=[ 'January','February','March', + 'April','May','June','July','August','September','October','November','December' ] + +flu = data['data'].T + +for ii in xrange(size_(flu,1)): + start_date_ii = 0 + for tt in xrange(size_(flu,2)): + if np.isnan(flu[ii,tt]): + start_date_ii += 1 + flu[ii,0:start_date_ii] = 0 + +q,T=size_(flu,nargout=2) + +# normalize the flu data +flu_norms = np.sqrt( np.max(flu, axis=1) ) +flu /= flu_norms[:, np.newaxis] +y = flu * 1.75 # why? + +# Start only after date where there is data for all y vecs +tmp = np.cumsum(np.sum(y,axis=-1), axis=-1) +tmp = np.where(tmp == 0)[0] +if np.any(tmp): + start_time = tmp[-1] + 1 + y = y[:][ start_time:-1 ] + +# Boolean vector containing True for all +# y vectors that have a value greater than zero +inds_y = np.ones(y.shape) +inds_y[np.where(y == 0)[0]] = 0 +inds_y = inds_y > 0 + +p,N = y.shape +x = np.arange(1, N+1) / N +c=100 +d=1 +r=1e-05 +K=np.zeros((N,N)) + +for ii in xrange(N): + for jj in xrange(N): + dist_ii_jj = np.abs(x[ii] - x[jj]) + K[ii,jj] = d * np.exp(- c * (dist_ii_jj ** 2)) + +K = K + diagonal(r * np.ones((1,N))) +print diagonal(K) +invK = inv(K) +logdetK=2 * np.sum(np.log(diagonal( chol(K)))) + +prior_params.K.c_prior=1 +prior_params.K.invK=invK +prior_params.K.K=K +prior_params.K.logdetK=logdetK +prior_params.sig.a_sig=1 +prior_params.sig.b_sig=0.1 +prior_params.hypers.a_phi=1.5 +prior_params.hypers.b_phi=1.5 +prior_params.hypers.a1=10 +prior_params.hypers.a2=10 + +settings.L=10 +settings.k=20 +settings.Niter=10000 +settings.saveEvery=100 +settings.storeEvery=10 +settings.saveMin=1 +settings.saveDir=char('flu') +settings.trial=1 +settings.init2truth=0 +settings.sample_K_flag=3 +settings.latent_mean=1 +settings.inds_y=inds_y + +BNP_covreg_varinds_(y,prior_params,settings,0) diff --git a/python/simData.mat b/python/simData.mat new file mode 100644 index 0000000..6871d46 Binary files /dev/null and b/python/simData.mat differ diff --git a/python/utilities/SIMplots.py b/python/utilities/SIMplots.py new file mode 100644 index 0000000..8c8537f --- /dev/null +++ b/python/utilities/SIMplots.py @@ -0,0 +1,143 @@ +# Autogenerated with SMOP version +# /usr/local/bin/smop SIMplots.m -o ../../python/utilities/SIMplots.py + +from __future__ import division +try: + from runtime import * +except ImportError: + from smop.runtime import * + +latent_mean=settings.latent_mean +inds2impute=settings.inds2impute +sampleEvery=copy_(storeEvery) +var_mean=zeros(p,p,N) +var_var=zeros(p,p,N) +var_u=zeros(p,p,N) +var_l=zeros(p,p,N) +mu_mean=zeros(p,N) +mu_var=zeros(p,N) +mu_u=zeros(p,N) +mu_l=zeros(p,N) +cov_true=true_params.cov_true +mu_true=true_params.mu +for tt in arange_(1,N).reshape(-1): + theta_zeta_tt=zeros(p,k,(Niter - Nburn) / sampleEvery) + var_tt=zeros(p,p,(Niter - Nburn) / sampleEvery) + mu_tt=zeros(p,(Niter - Nburn) / sampleEvery) + m=1 + for nn in arange_(Nburn + 1,Niter,sampleEvery).reshape(-1): + n=nn + saveEvery - 1 + if rem(n,saveEvery) == 0 and n <= Niter: + filename=[saveDir,char('/BNP_covreg_statsiter'),num2str(n),char('trial'),num2str(trial),char('.mat')] + load(filename) + store_count=1 + theta_zeta_tt[arange_(),arange_(),m]=Stats(store_count).theta * Stats(store_count).zeta(arange_(),arange_(),tt) + var_tt[arange_(),arange_(),m]=Stats(store_count).theta * Stats(store_count).zeta(arange_(),arange_(),tt) * Stats(store_count).zeta(arange_(),arange_(),tt).T * Stats(store_count).theta.T + diag(1.0 / Stats(store_count).invSig_vec) + mu_tt[arange_(),m]=Stats(store_count).theta * Stats(store_count).zeta(arange_(),arange_(),tt) * Stats(store_count).psi(arange_(),tt) + m=m + 1 + store_count=store_count + 1 + var_mean[arange_(),arange_(),tt]=mean(var_tt,3) + var_var[arange_(),arange_(),tt]=var(var_tt,0,3) + mu_mean[arange_(),tt]=mean(mu_tt,2) + mu_var[arange_(),tt]=var(mu_tt,0,2) + for pp in arange_(1,p).reshape(-1): + for jj in arange_(pp,p).reshape(-1): + var_u(pp,jj,tt),var_l(pp,jj,tt)=calculate_hpd(var_tt(pp,jj,arange_()),0.95,nargout=2) + if latent_mean: + mu_u(pp,tt),mu_l(pp,tt)=calculate_hpd(mu_tt(pp,arange_()),0.95,nargout=2) + if not rem(tt,100): + display(num2str(tt)) +LineWidth=1.5 +fs=20 +figure +for pp in arange_(1,p).reshape(-1): + for jj in arange_(pp,p).reshape(-1): + plot(squeeze(var_mean(pp,jj,arange_())),char('LineWidth'),LineWidth) + hold(char('all')) +xlim([0,N]) +ylabel(char('Variance/Covariance'),char('FontSize'),fs) +xlabel(char('Time'),char('FontSize'),fs) +_set(gca,char('FontSize'),16) +ylim([- 2,2]) +title(char('BNP Covariance Regression'),char('Fontsize'),20) +figure +for pp in arange_(1,p).reshape(-1): + for jj in arange_(pp,p).reshape(-1): + plot(squeeze(true_params.cov_true(pp,jj,arange_())),char('LineWidth'),LineWidth) + hold(char('all')) +xlim([0,N]) +ylabel(char('Variance/Covariance'),char('FontSize'),fs) +xlabel(char('Time'),char('FontSize'),fs) +_set(gca,char('FontSize'),16) +ylim([- 2,2]) +title(char('Truth'),char('FontSize'),20) +LineWidth=1.5 +figure +for pp in arange_(1,p).reshape(-1): + plot(mu_mean(pp,arange_()),char('LineWidth'),LineWidth) + hold(char('all')) +xlim([0,N]) +ylabel(char('Mean'),char('FontSize'),fs) +xlabel(char('Time'),char('FontSize'),fs) +_set(gca,char('FontSize'),16) +title(char('BNP Covariance Regression'),char('Fontsize'),20) +figure +for pp in arange_(1,p).reshape(-1): + plot(true_params.mu(pp,arange_()),char('LineWidth'),LineWidth) + hold(char('all')) +xlim([0,N]) +ylabel(char('Mean'),char('FontSize'),fs) +xlabel(char('Time'),char('FontSize'),fs) +_set(gca,char('FontSize'),16) +title(char('Truth'),char('FontSize'),20) +ylim([- 3,3]) +fs=20 +Sig_vec=true_params.Sig_vec +theta=true_params.theta +ind_start=1 +ind_end=copy_(nmc) +inds=[arange_(1,N)] +if strcmp(reg_flag,char('full')): + figure + hl=boxplot(1.0 / invSig_hist(arange_(),[arange_(ind_start,ind_end)]).T,char('notch'),char('on'),char('widths'),0.75,char('outliersize'),2) + for ih in arange_(1,size(hl,1)).reshape(-1): + _set(hl(ih,arange_()),char('LineWidth'),2) + hold(char('on')) + scatter([arange_(1,p)],Sig_vec,75,char('go'),char('filled')) + hold(char('off')) + ylabel(char('\\Sigma_{0,p}'),char('FontSize'),20) + xlabel(char('p'),char('FontSize'),20) + _set(gca(),char('FontSize'),16) +cov_true=true_params.cov_true +figure +for pp in arange_(1,p).reshape(-1): + for jj in arange_(1,pp).reshape(-1): + plot(squeeze(var_l(pp,jj,arange_())),char('--'),char('LineWidth'),2) + hold(char('on')) + plot(squeeze(var_u(pp,jj,arange_())),char('--'),char('LineWidth'),2) + plot(squeeze(var_mean(pp,jj,arange_())),char('LineWidth'),2) + plot(squeeze(cov_true(pp,jj,arange_())),char('r'),char('LineWidth'),2) + hold(char('off')) + xlim([0,N]) + axis(char('tight')) + axis(char('square')) + box(char('on')) + _set(gca(),char('XTickLabel'),char(''),char('XTick'),zeros(1,0),char('YTickLabel'),char(''),char('YTick'),zeros(1,0),char('DataAspectRatio'),[10,_max(_abs(cov_true(pp,jj,arange_()))),1]) + title([char('p: '),num2str(pp),char(','),num2str(jj)]) + waitforbuttonpress +mu_true=true_params.mu +figure +for pp in arange_(1,p).reshape(-1): + plot(mu_l(pp,arange_()),char('--'),char('LineWidth'),2) + hold(char('on')) + plot(mu_u(pp,arange_()),char('--'),char('LineWidth'),2) + plot(mu_mean(pp,arange_()),char('LineWidth'),2) + plot(mu_true(pp,arange_()),char('r'),char('LineWidth'),2) + hold(char('off')) + xlim([0,N]) + axis(char('tight')) + axis(char('square')) + box(char('on')) + _set(gca(),char('XTickLabel'),char(''),char('XTick'),zeros(1,0),char('YTickLabel'),char(''),char('YTick'),zeros(1,0),char('DataAspectRatio'),[10,_max(_abs(cov_true(pp,jj,arange_()))),1]) + title([char('p: '),num2str(pp),char(','),num2str(jj)]) + waitforbuttonpress diff --git a/python/utilities/clculate_hpd.py b/python/utilities/clculate_hpd.py new file mode 100644 index 0000000..a32c36e --- /dev/null +++ b/python/utilities/clculate_hpd.py @@ -0,0 +1,30 @@ +# Autogenerated with SMOP version +# /usr/local/bin/smop calculate_hpd.m -o ../../python/utilities/clculate_hpd.py + +from __future__ import division +try: + from runtime import * +except ImportError: + from smop.runtime import * + +def calculate_hpd_(samples=None,conf=None,*args,**kwargs): + varargin = cellarray(args) + nargin = 2-[samples,conf].count(None)+len(args) + + num_hpd_samples=floor_(length_(samples) * conf) + samples=sort_(samples,char('descend')) + u_hpd=samples[1] + l_hpd=samples[num_hpd_samples] + sample_ind=1 + for u in samples.reshape(-1): + samples_below_u=samples < u + if sum_(samples_below_u) < num_hpd_samples: + break + else: + u_hpd_tmp=copy_(u) + l_hpd_tmp=samples[sample_ind + num_hpd_samples - 1] + if u_hpd_tmp - l_hpd_tmp < u_hpd - l_hpd: + u_hpd=copy_(u_hpd_tmp) + l_hpd=copy_(l_hpd_tmp) + sample_ind=sample_ind + 1 + return u_hpd,l_hpd