function POST = abc_regression(tgt,params, SS, dotrans, use) % Weighted regression in space of summary statistics. % [F, DF, D2F] = OPF_COSTFCN(X, OM) % % Apply a weighted regression in the space of summary statistics % as proposed in Beaumont et al. (2002). % % Inputs: % tgt : target summary statistics vector (observed). % params : parameteter matrix: % SS : Summary statistics (same number of rows as params) % dotrans : Apply Tanh transposition as proposed in Beaumont et al. (2002). % use : define usage vector (logical). % % Outputs: % POST : posterior distribution per parameter. % % abc_regression.m % Written by Andreas Wollstein, Ludwig Maximilians University of Munich. % Copyright (c) 2015 % % This file is part of a set of Matlab scripts published with: % Bozicevic et al. 2015 Molecular Ecology % % This program is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation, either version 3 of the License, or % (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program. If not, see . addpath ~/Dropbox/online_projects/mgstats/; %% beaumont regression. assume all summary statistics are z scored if ~exist('use','var') use = true(1,length(tgt)); end POST = []; SSS = SS; ttgt = tgt; y = params; %% check informaticity on SS by linear prediction of y [yy, stats] = gfit(SSS,y); ix = stats.p(2:end)<0.01; %ix = 1:length(ix); SS = SSS(:,ix); tgt = ttgt(ix); %% calculate the squared distance of summary stats to target vec dst = (SS - repmat(tgt,size(SS,1),1)).^2; dst = sqrt(nansum(dst,2)); X = SS; weights = 1.0-(dst./max(dst)).^2; %% transpose if dotrans miny = nanmin(y); maxy = nanmax(y); y = log(tan( (y-miny) ./ (maxy-miny).*pi./2.0 )); y(isinf(y))=nan; end [b dev stats] = glmfit(X,y,'normal','weights',weights); predmean = glmval(b,tgt,'identity'); post = predmean-stats.resid; %% backtranspose if dotrans post = 2*atan(exp(post)).*(maxy-miny)./pi+miny; end POST = post; POST(isnan(POST))=[]; end