first commit

code/MCMC_bin_posterior.m

@@ -0,0 +1,52 @@
+function betas = MCMC_bin_posterior(PhiW,Y,mode,precision,c,n_samp,varargin)
+if ~isempty(varargin)
+    opts = varargin{1};
+else
+    opts = struct;
+end
+burn_in = get_opt(opts, 'burn_in', 100);
+samp_gap = get_opt(opts, 'samp_gap', 100);
+offset = get_opt(opts, 'offset', 0);
+d = size(PhiW,2);
+
+% sample from Laplace approx
+L = chol(precision);
+beta =  L\randn(d,1)+mode;
+% get beta stats
+PB = PhiW*beta+offset;
+expnPB = exp(-PB);
+normB2 = sum(beta(:).^2);
+logpbeta = -.5*(beta-mode)'*precision*(beta-mode);
+tot_iters = burn_in + samp_gap*n_samp;
+betas = nan(d,n_samp);
+accepted = 0;
+for i=1:tot_iters
+    % sample from Laplace approx
+    beta_prop =  L\randn(d,1)+mode;
+    % get proposal stats
+    PB_prop = PhiW*beta_prop+offset;
+    PDB = PB_prop-PB;
+    normB2_prop = sum(beta_prop(:).^2);
+    logpbeta_prop = -.5*(beta_prop-mode)'*precision*(beta_prop-mode);
+    % get log acceptance ratio
+    logratio = Y'*PDB -sum(log( (expnPB+exp(PDB))./(1+expnPB) ),1) ...
+        -5*c*(normB2_prop-normB2);
+    logratio = logratio +logpbeta -logpbeta_prop;
+
+    % accept?
+    if logratio>0 || log(rand)<=logratio
+        beta = beta_prop;
+        PB = PB_prop;
+        expnPB = exp(-PB);
+        normB2 = normB2_prop;
+        logpbeta = logpbeta_prop;
+        accepted = accepted +1;
+    end
+    % save?
+    if i>burn_in && mod(i-burn_in,samp_gap)==0
+        betas(:,(i-burn_in)/samp_gap) = beta;
+    end
+end
+
+
+end

code/W_stats_reg_model_ev.m

@@ -0,0 +1,60 @@
+function W_stats = W_stats_reg_model_ev(W,X,Y,model_prior,varargin)
+% W_stats_reg_model_ev Get stats used in sampling W.
+% Inputs:
+%   W - d x nfreq matrix of random frequencies
+%   X - N x d matrix of input covariates
+%   Y - N x d2 matrix of output responses
+%   model_prior -  struct with the follow feilds
+%     (.a,.b) - inverse-gamma prior parameters to noise variance in linear
+%             model
+%     .c - scalar multiplier on identity matrix for beta weights
+% Outputs:
+%   W_stats - struct with the following feilds
+%     .N - number of instances
+%     .p - number of dimensions for responses
+%     .W - random frequency matrix
+%     .PhiW - matrix of random features for instances
+%     .R - chol of Phi'*Pho+c*I, where c is prior parameter for beta (see
+%         https://en.wikipedia.org/wiki/Bayesian_linear_regression#Other_cases)
+%     .PtY - inner product of Phi'*Y (where Phi is matrix of random
+%           features)
+%     .YtY - norm of response vector
+%     .betaW - betaW = PtP_cI\PtY, the mode for beta regression weights
+%     .log_det_PtP_cI - log determinant of Phi'*Pho+c*I, where c is model
+%                      prior parameter (see below).
+%     .log_f - Get the log evidence of the model evidence w.r.t. random 
+%             features.
+
+if ~isempty(varargin)
+    opts = varargin{1};
+else
+    opts = struct;
+end
+save_PhiW = get_opt(opts, 'save_PhiW', true);
+
+XW = X*W;
+PhiW = [cos(XW) sin(XW)];
+
+c = model_prior.c;
+
+PtP_cI = PhiW'*PhiW+c*speye(size(PhiW,2));
+R = chol(PtP_cI);
+PtY = PhiW'*Y;
+
+betaW = R\(R'\PtY); % betaW = PtP_cI\PtY;
+log_det_PtP_cI = sum(2*log(diag(R))); % log_det_PtP_cI = log(det(PtP_cI));
+
+W_stats.N = numel(Y);
+W_stats.p = size(Y,2);
+W_stats.W = W;
+if save_PhiW
+    W_stats.PhiW = PhiW;
+end
+W_stats.R = R;
+W_stats.PtY = PtY;
+W_stats.YtY = sum(Y(:).^2);
+W_stats.betaW = betaW;
+W_stats.PhiWbeta = PhiW*betaW;
+W_stats.log_det_PtP_cI = log_det_PtP_cI;
+W_stats.log_f = get_log_model_evidence_regression(W_stats,model_prior);
+end

code/bank.m

@@ -0,0 +1,81 @@
+function [Ws,rhos] = ...
+    bank(data,nfreq,rho_prior,model_prior,update,init_Wstats,varargin)
+% BaNK function to run the BaNK framework.
+% Inputs: 
+%   data - Input data to work over. E.g. can be a cell array {X Y} for
+%          supervised problems, or a matrix X for unsupervized problems.
+%   nfreq - Number of random features to use.
+%   rho_prior - Prior parameters for the distribution of random features.
+%               (Can also include options for update functions below.)
+%   model_prior - Prior parameters for the distribution of model specific
+%                 parameters. E.g. normal parameters mu and Sigma for
+%                 linear regression. (Can also include options for update
+%                 functions below.)
+%   update - struct with the following function handles
+%     function rho = update.rho(rho,W,rho_prior):
+%       Function that samples the posterior of rho given W, rest. E.g.
+%       samples mean, covariances for a GMM given W and component
+%       assignments.
+%     function Wstats = update.W(Wstats,rho,data,model_prior)
+%       Function that samples random frequencies, and other model stats
+%       (like beta weights if not marginalized out) given rho (and posibly
+%       component assignments in Z). Note, this function should handle the
+%       case where Wstats is given with only Wstats.W as a feild, which 
+%       will occur in the first iteration.
+%   init_Wstats - function that returns initial random frequencies, and 
+%                 other model stats:
+%       function Wstats = init_Wstats(data,nfreq,rho_prior,model_prior) 
+%   opts (optional) - struct with optional parameters with {defualts} 
+%     burn_in - number of draws to burn in {10}
+%     samp_gap - number of draws to discard between samples {100}
+%     nsamps - number of samples to return {1}
+%     itprint - function itprint(i,data,rho,Wstats) that
+%       prints to stdout (or does any other sort of logging/processing) 
+%       on the iterates that save a sample {@(~,~,~,~,~)[]}.
+% Outputs:
+%   Ws - 3d array of random frequency samples where the last dimension is
+%        length nsamp and Ws(:,:,i) is the ith matrix of random
+%        frequencies.
+%   rhos - cell array of length nsamp, where rhos{i} is the ith random
+%          frequency distribution drawn from posterior.
+
+if ~isempty(varargin)
+    opts = varargin{1};
+else
+    opts = struct;
+end
+
+% sampling options
+burn_in = get_opt(opts, 'burn_in', 10);
+samp_gap = get_opt(opts, 'samp_gap', 100);
+n_samp = get_opt(opts, 'n_samp', 1);
+
+% get initial random frequencies
+Wstats = init_Wstats(data,nfreq,rho_prior,model_prior);
+
+% printing function
+itprint = get_opt(opts, 'itprint', @(~,~,~,~,~)[]);
+
+% main loop
+tot_iters = burn_in + samp_gap*n_samp;
+Ws = nan(size(Wstats.W,1),nfreq,n_samp);
+rhos = cell(n_samp,1);
+rho = get_opt(opts, 'rho', struct);
+for i=1:tot_iters
+    % sample posterior for rho given W
+    rho = update.rho(rho,Wstats.W,rho_prior);
+
+    % sample posterior for W given rho
+    Wstats = update.W(Wstats,rho,data,model_prior);
+
+    % save?
+    if i>burn_in && mod(i-burn_in,samp_gap)==0
+        Ws(:,:,(i-burn_in)/samp_gap) = Wstats.W;
+        rhos{(i-burn_in)/samp_gap} = rho;
+    end
+
+    itprint(i,data,rho,Wstats,rho_prior,model_prior);
+end
+
+end
+

code/binlogreg.m

@@ -0,0 +1,111 @@
+function lr_stats = binlogreg(X,Y,varargin)
+
+if ~isempty(varargin)
+    opts = varargin{1};
+else
+    opts = struct;
+end
+
+if ~iscell(X)
+    [N,d] = size(X);
+else
+    N = length(X);
+    d = size(X{1},2);
+end
+Y = double(Y==1);
+
+trn_set = get_opt(opts, 'trn_set');
+hol_set = get_opt(opts, 'hol_set');
+tst_set = get_opt(opts, 'tst_set');
+if isempty(trn_set) || isempty(hol_set) || isempty(tst_set)
+    tprec = get_opt(opts, 'tprec', .1);
+    hprec = get_opt(opts, 'tprec', .1);
+    [trn_set, hol_set, tst_set] = split_data( N, tprec, hprec );
+end
+
+sigma2s = get_opt(opts,'sigma2s');
+if isempty(sigma2s)
+    N_rot = get_opt(opts, 'N_rot', min(2000,N));
+    pd2s = dists2(X(randperm(N,N_rot),:),X(randperm(N,N_rot),:));
+    sigma2s = quantile(pd2s(:), .1:.2:.9);
+end
+lambdas = get_opt(opts,'lambdas',2.^(10:-1:-5));
+
+nfreq = get_opt(opts,'nfreq',250);
+W = get_opt(opts,'W',randn(d,nfreq));
+bagsizes = get_opt(opts, 'bagsizes');
+do_mmd = iscell(X) | ~isempty(bagsizes);
+if ~isempty(bagsizes)
+    lastinds = cumsum(bagsizes);
+end
+if ~do_mmd
+    XW = X*W;
+end
+
+nsigma2s = length(sigma2s);
+nlambdas = length(lambdas);
+bopts = get_opt(opts, 'bopts', struct);
+hol_err = nan(nsigma2s,nlambdas);
+betas = cell(nsigma2s,nlambdas);
+for si=1:nsigma2s
+    sigma2 = sigma2s(si);
+    if ~do_mmd
+        PhiW = [cos(XW./sqrt(sigma2)) sin(XW./sqrt(sigma2)) ones(N,1)];
+    else
+        if iscell(X)
+            PhiW = cell2mat( cellfun(@(C)[mean([cos(C*W./sqrt(sigma2)) sin(C*W./sqrt(sigma2))],1) 1], X, 'unif', false) );
+        else
+            PhiW = X*W./sqrt(sigma2); 
+            PhiW = cumsum([cos(PhiW) sin(PhiW)],1);
+            PhiW = PhiW(lastinds,:);
+            PhiW(2:end,:) = PhiW(2:end,:) - PhiW(1:end-1,:);
+            PhiW = bsxfun(@times,PhiW,1./bagsizes);
+            PhiW = [PhiW ones(size(PhiW,1),1)];
+        end
+    end
+
+    beta_si = zeros(size(PhiW,2),1);
+    for li=1:nlambdas
+        beta_si =  newtons_binlr(PhiW(~trn_set,:),Y(~trn_set),lambdas(li),beta_si,bopts);
+        Yhol_pred = PhiW(hol_set,:)*beta_si>=0;
+        hol_err(si,li) = sum(Y(hol_set)~=Yhol_pred)/length(Yhol_pred);
+        betas{si,li} = beta_si;
+    end
+end
+
+[si,li] = find(hol_err==min(hol_err(:)));
+ri = randi(length(si));
+si = si(ri); 
+li = li(ri); 
+beta_si = betas{si,li};
+sigma2 = sigma2s(si);
+if ~do_mmd
+    PhiW = [cos(XW./sqrt(sigma2)) sin(XW./sqrt(sigma2)) ones(N,1)];
+else
+    if iscell(X)
+        PhiW = cell2mat( cellfun(@(C)[mean([cos(C*W./sqrt(sigma2)) sin(C*W./sqrt(sigma2))],1) 1], X, 'unif', false) );
+    else
+        PhiW = X*W./sqrt(sigma2); 
+        PhiW = cumsum([cos(PhiW) sin(PhiW)],1);
+        PhiW = PhiW(lastinds,:);
+        PhiW = PhiW - [zeros(1,size(PhiW,2)); PhiW(1:end-1,:)];
+        PhiW = bsxfun(@times,PhiW,1./bagsizes);
+        PhiW = [PhiW ones(size(PhiW,1),1)];
+    end
+end
+
+beta = newtons_binlr(PhiW(~tst_set,:),Y(~tst_set),lambdas(li),beta_si,bopts);
+Ytst_pred = PhiW(tst_set,:)*beta>=0;
+tst_err = sum(Ytst_pred~=Y(tst_set))/length(Ytst_pred);
+
+lr_stats.Y_pred = Ytst_pred;
+lr_stats.resids = (Ytst_pred ~= Y(tst_set));
+lr_stats.tst_err = tst_err;
+lr_stats.beta = beta;
+lr_stats.W = W;
+lr_stats.sigma2 = sigma2;
+lr_stats.lambda = lambdas(li);
+lr_stats.hol_err = hol_err;
+lr_stats.D = nfreq;
+lr_stats.lambdars = lambdas;
+end

code/cdesc_rks_W.m

@@ -0,0 +1,35 @@
+function [W, pred, beta] = cdesc_rks_W( X,Y,W,beta,pred,...
+                             c,mus_z,Siginvs_z,sigma2_error, varargin )
+
+if ~isempty(varargin)
+    opts = varargin{1};
+else
+    opts = struct;
+end
+citers = get_opts(opts, 'citers', 25);
+liters = get_opts(opts, 'liters', 100);
+
+nfreq = size(W,2);
+beta_cos = beta(1:nfreq);
+beta_sin = beta(nfreq+1:2*nfreq);
+lb_cell = {-inf(nfreq,1), -inf(nfreq,1), -inf(d,nfreq)};
+ub_cell = { inf(nfreq,1),  inf(nfreq,1),  inf(d,nfreq)};
+for ci=1:citers
+    jprm = randperm(nfreq);
+    for jj=1:nfreq
+        j = jprm(jj);
+
+        XWj = X*W(:,j);
+        offset = pred - (cos(XWj)*beta_cos(j) + sin(XWj)*beta_sin(j));
+
+        auxdata = {X,Y,c,mus_z(:,1,j),Siginvs_z(:,:,j),sigma2_error,offset};
+        [beta_cos(j), beta_sin(j), W(:,j)] = ...
+          lbfgsb( {beta_cos(j), beta_sin(j),  W(:,j)}, lb_cell, ...
+            ub_cell, 'rks_W_obj','rks_W_grad', auxdata,'print100callback', ...
+            'm',25,'factr',1e-12, 'pgtol',1E-4,'maxiter',liters );
+
+        pred = offset + (cos(XWj)*beta_cos(j) + sin(XWj)*beta_sin(j));
+    end
+end
+
+end

code/classification_wrapper.m

@@ -0,0 +1,16 @@
+addpath(genpath('/usr0/home/akdubey/Juice/mtimesx'));
+addpath(genpath('/usr0/home/akdubey/Juice/minFunc_2012'));
+addpath(genpath('/usr0/home/akdubey/Juice/funcLearn'));
+addpath(genpath('/usr0/home/akdubey/Juice/gpml-matlab-v3.6-2015-07-07'));
+addpath(genpath('/usr0/home/akdubey/Juice/lbfgs/lbfgsb-matlab/src'));
+addpath(genpath('/usr0/home/akdubey/Juice/projectbirdup/kernelLearn/code'));
+
+[funcs, names] = make_classification_funcs();
+
+cv_opts.names = names;
+cv_opts.K = 5;
+cv_opts.fileprefix = 'results_checking';
+cv_opts.stand_response = false;
+load ../../datasets/Classification/pima_data.mat
+parpool(3);
+[tst_errs, tst_resids, method_stats, inds, hol_inds] = kfold_cv(X,Y,@run_classification,funcs,cv_opts);

code/demo.m

@@ -0,0 +1,5 @@
+load datasets/bike.mat
+
+f = make_regression_funcs();
+opts.funcs = f;
+method_stats = run_regression(X,(Y-mean(Y))./std(Y),opts);

code/dists2.m

@@ -0,0 +1,17 @@
+function M = dists2(X1, X2, varargin)
+row_insts = true;
+if ~isempty(varargin)
+    row_insts = varargin{1};
+end
+if row_insts
+    X1 = X1';
+    X2 = X2';
+end
+M = bsxfun(@plus, sum(X1 .* X1, 1)', (-2) * X1' * X2);        
+M = bsxfun(@plus, sum(X2 .* X2, 1), M);        
+M(M < 0) = 0;                        
+
+end
+
+
+