0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 1 % MCMC estimates paramters using a Monte Carlo Markov Chain.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 3 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 4 % DESCRIPTION: MCMC estimate the parameters of a given model given
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 5 % inputs, outputs and noise using a Metropolis algorithm.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 6 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 7 % CALL: [b smplr] = mcmc(out,pl)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 8 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 9 % INPUTS: out - matrix objects with measured outputs
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 10 % pl - parameter list
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 11 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 12 % OUTPUTS: b - pest object contatining estimate information
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 13 % smplr - matrix containing info about the rejected points
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 14 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 15 % <a href="matlab:utils.helper.displayMethodInfo('matrix', 'mcmc')">Parameters Description</a>
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 16 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 17 % VERSION: $Id: mcmc.m,v 1.35 2011/11/16 15:21:13 nikos Exp $
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 18 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 19 % References: M Nofrarias et al. Phys. Rev. D 82, 122002 (2010)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 20 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 21 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 22 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 23
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 24 % TODO: multiple chain option not implemented yet
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 25 % metropolis/hastings not implemented
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 26 % empty initial values not implemented
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 27
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 28
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 29 function varargout = mcmc(varargin)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 30
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 31 % Check if this is a call for parameters
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 32 if utils.helper.isinfocall(varargin{:})
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 33 varargout{1} = getInfo(varargin{3});
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 34 return
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 35 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 36
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 37 import utils.const.*
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 38 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 39
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 40 % Method can not be used as a modifier
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 41 if nargout == 0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 42 error('### mcmc cannot be used as a modifier. Please give an output variable.');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 43 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 44
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 45 % Collect input variable names
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 46 in_names = cell(size(varargin));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 47 for ii = 1:nargin,in_names{ii} = inputname(ii);end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 48
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 49 % Collect all AOs smodels and plists
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 50 [mtxs, mtxs_invars] = utils.helper.collect_objects(varargin(:), 'matrix', in_names);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 51 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 52
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 53 % Combine plists
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 54 pl = parse(pl, getDefaultPlist);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 55
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 56 % Get parameters
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 57 N = find(pl, 'N');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 58 Tc = find(pl,'Tc');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 59 xi = find(pl,'heat');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 60 cvar = find(pl, 'cov');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 61 rng = find(pl,'range');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 62 flim = find(pl,'frequencies');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 63 fsout = find(pl,'fsout');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 64 search = find(pl,'search');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 65 simplex = find(pl,'simplex');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 66 mhsampleTrue = find(pl,'mhsample');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 67 xo = find(pl,'x0');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 68 mtxin = find(pl,'input');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 69 mdlin = find(pl,'model');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 70 mdlFreqDependent = find(pl,'modelFreqDependent');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 71 nse = find(pl,'noise');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 72 jumps = find(pl,'jumps');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 73 parplot = find(pl,'plot');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 74 debug = find(pl,'debug');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 75 outModel = find(pl,'outModel');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 76 inModel = find(pl,'inModel');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 77 outNames = find(pl,'outNames');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 78 inNames = find(pl,'inNames');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 79 fpars = find(pl,'prior');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 80 anneal = find(pl,'anneal');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 81 DeltaL = find(pl,'DeltaL');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 82 SNR0 = find(pl,('SNR0'));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 83
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 84 % Decide on a deep copy or a modify
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 85 in = copy(mtxin, nargout);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 86 out = copy(mtxs, nargout);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 87 mdl = copy(mdlin, nargout);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 88
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 89 % checking if input or noise are aos
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 90 if isa(in, 'ao')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 91 in = matrix(in,plist('shape',[size(in,1) size(in,2)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 92 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 93 if isa(nse, 'ao')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 94 nse = matrix(nse,plist('shape',[size(nse,1) size(nse,2)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 95 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 96
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 97 % lighten the model
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 98 mdl.clearHistory;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 99 if isa(mdl, 'ssm')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 100 mdl.clearNumParams;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 101 mdl.clearAllUnits;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 102 mdl.params.simplify;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 103 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 104
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 105 % Check input parameters
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 106 if isempty(rng)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 107 error('### Please define a search range ''range''');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 108 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 109 param = find(pl,'FitParams');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 110 if isempty(param)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 111 error('### Please define the parameters ''param''');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 112 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 113 nparam = numel(param);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 114 if size(mdl) ~= size(nse)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 115 error('### Parameters ''model'' and ''noise'' must be the same dimension');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 116 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 117 if size(in) ~= size(out)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 118 error('### Number of input and output experiments must be the same');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 119 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 120
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 121
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 122 % Get range for parameters
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 123 for i = 1:nparam
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 124 rang(:,i) = rng{i};
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 125 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 126
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 127 % Get number of experiments
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 128 nexp = numel(in);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 129
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 130 switch class(mdl)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 131 case 'matrix'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 132 % Check model sizes
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 133 if (isempty(outModel) && isempty(inModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 134 if (~(numel(in(1).objs) == mdl.ncols) || ~(numel(out(1).objs) == mdl.nrows))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 135 error('Check model or input/output sizes');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 136 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 137 elseif isempty(inModel)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 138 if (~(numel(in(1).objs) == mdl.ncols) || ~(size(outModel,2) == mdl.nrows) || ~(size(outModel,1) == numel(out(1).objs)))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 139 error('Check model or input/output sizes');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 140 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 141 elseif isempty(outModel)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 142 if (~(numel(in(1).objs) == size(inModel,2)) || ~(size(inModel,1) == mdl.ncols) || ~(numel(out(1).objs) == mdl.nrows))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 143 error('Check model or input/output sizes');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 144 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 145 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 146 if (~(numel(in(1).objs) == size(inModel,2)) || ~(size(inModel,1) == mdl.ncols) || ~(size(outModel,2) == mdl.nrows) || ~(numel(out(1).objs) == size(outModel,1)))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 147 error('Check model or input/output sizes');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 148 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 149 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 150 case 'ssm'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 151 inNames = find(pl,'inNames');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 152 outNames = find(pl,'outNames');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 153 if( (numel(inNames) ~= numel(in(1).objs)) || numel(outNames) ~= numel(out(1).objs))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 154 error('Check model inNames and outNames, they do not match with the input objects')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 155 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 156 otherwise
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 157 error('### Model must be either from the ''matrix'' or the ''ssm'' class. Please check the inputs')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 158 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 159
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 160
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 161
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 162 %%%%%%%%%%%%%%%%%%%%%%%%%%%% Frequency domain pre-process %%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 163
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 164 utils.helper.msg(msg.IMPORTANT, 'Computing frequency domain objects', mfilename('class'), mfilename);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 165 for k = 1:nexp
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 166
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 167 if((numel(out(1).objs) == 2) && (numel(in(1).objs) == 2)) % previous code done by miquel
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 168 % optional resampling (not matrix/resample method yet)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 169 if ~isempty(fsout)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 170 for i = 1:numel(in(k).objs(:))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 171 in(k).objs(i).resample(plist('fsout',fsout));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 172 out(k).objs(i).resample(plist('fsout',fsout));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 173 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 174 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 175
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 176 % fft
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 177 fin(k) = fft(in(k));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 178 fout(k) = fft(out(k));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 179
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 180 Nsamples = length(fin(k).getObjectAtIndex(1).x);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 181 fs = fin(k).getObjectAtIndex(1).fs;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 182
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 183
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 184 if (isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 185 % Get rid of fft f =0, reduce frequency range if needed
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 186 if ~isempty(flim)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 187 fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 188 fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 189 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 190 elseif (~isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 191 % Get rid of fft f =0, reduce frequency range if needed
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 192 if ~isempty(flim)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 193 fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 194 fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 195 if(~isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 196 for lll=1:size(outModel,1)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 197 for kkk=1:size(outModel,2)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 198 outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 199 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 200 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 201 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 202 if(~isempty(inModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 203 inModel = split(inModel,plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 204 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 205 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 206 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 207
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 208 % use signal fft to get frequency vector. Take into account signal
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 209 % could be empty or set to zero
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 210
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 211 % 1st channel
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 212 if all(fin(k).getObjectAtIndex(1,1).y == 0) || isempty(fin(k).getObjectAtIndex(1,1).y)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 213 i1 = ao(plist('type','fsdata','xvals',0,'yvals',0));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 214 % rebuild input with new zeroed signal
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 215 fin(k) = matrix(i1,fin(k).getObjectAtIndex(2,1),plist('shape',[2 1]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 216 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 217 freqs = fin(k).getObjectAtIndex(1,1).x;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 218 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 219 % 2nd channel
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 220 if all(fin(k).getObjectAtIndex(2,1).y == 0) || isempty(fin(k).getObjectAtIndex(2,1).y)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 221 i2 = ao(plist('type','fsdata','xvals',0,'yvals',0));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 222 % rebuild input with new zeroed signal
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 223 fin(k) = matrix(fin(k).getObjectAtIndex(1,1),i2,plist('shape',[2 1]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 224 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 225 freqs = fin(k).getObjectAtIndex(2,1).x;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 226 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 227
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 228 % Compute psd
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 229 n1 = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(1,1), pl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 230 n2 = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(2,1), pl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 231 n12 = Nsamples*fs/2*cpsd(nse(k).getObjectAtIndex(1,1),nse(k).getObjectAtIndex(2,1), pl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 232
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 233 % interpolate noise to fft frequencies
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 234 S11 = interp(n1,plist('vertices',freqs));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 235 S12 = interp(n12,plist('vertices',freqs));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 236 S22 = interp(n2,plist('vertices',freqs));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 237 S21 = conj(S12);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 238 % build cross-spectrum matrix
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 239 mnse(k) = matrix(S11,S21,S12,S22,plist('shape',[2 2]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 240
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 241 switch class(mdl)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 242 case 'matrix'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 243 for i = 1:numel(mdl.objs)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 244 if (mdlFreqDependent)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 245 % set Xvals
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 246 mdl.objs(i).setXvals(freqs);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 247 % set alias
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 248 mdl.objs(i).assignalias(mdl.objs(i),plist('xvals',freqs));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 249 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 250 mdl.objs(i).setXvals(1);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 251 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 252 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 253 case 'ssm'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 254 if k<2
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 255 mdl.clearNumParams;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 256 spl = plist('set', 'for bode', ...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 257 'outputs', outNames, ...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 258 'inputs', inNames);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 259 % first optimise our model for the case in hand
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 260 mdl.reorganize(spl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 261 % make it lighter
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 262
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 263 %for k = 1:nexp
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 264
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 265 mdl(k).optimiseForFitting();
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 266 %end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 267 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 268 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 269 % set Xvals for model (frequency vector could change for each experiment)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 270
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 271 mmdl(k) = mdl;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 272
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 273 elseif ((numel(out(1).objs) == 3) && (numel(in(1).objs) == 4))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 274
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 275 % here we are implementing only the magnetic case
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 276 % We have 4 inputs (the 4 conformator waveforms of the magnetic
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 277 % analysis and
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 278 % 3 outputs (that correspond to the IFO.x12 and IFO.ETA1 and
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 279 % IFO.PHI1
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 280
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 281
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 282 % fft
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 283 fin(k) = fft(in(k));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 284 fout(k) = fft(out(k));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 285 Nsamples = length(fin(k).getObjectAtIndex(1,1).x);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 286 fs = fin(k).getObjectAtIndex(1).fs;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 287
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 288 % Get rid of fft f =0, reduce frequency range if needed
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 289 if ~isempty(flim)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 290 fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 291 fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 292
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 293 if(~isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 294 for lll=1:size(outModel,1)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 295 for kkk=1:size(outModel,2)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 296 outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 297 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 298 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 299 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 300 if(~isempty(inModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 301 inModel = split(inModel,plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 302 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 303 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 304
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 305 freqs = fin(k).getObjectAtIndex(1,1).x;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 306
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 307 % Build noise model
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 308 for ii = 1:numel(out.objs)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 309 for jj = ii:numel(out.objs)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 310 % Compute psd
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 311 if (ii==jj)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 312 n(ii,jj) = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(ii), pl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 313 S(ii,jj) = interp(n(ii,jj),plist('vertices',freqs,'method','linear'));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 314 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 315 n(ii,jj) = Nsamples*fs/2*cpsd(nse(k).getObjectAtIndex(ii),nse(k).getObjectAtIndex(jj),pl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 316 S(ii,jj) = interp(n(ii,jj),plist('vertices',freqs,'method','linear'));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 317 S(jj,ii) = conj(S(ii,jj));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 318 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 319 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 320 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 321 % build cross-spectrum matrix
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 322 mnse(k) = matrix(S,plist('shape',[numel(out.objs) numel(out.objs)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 323
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 324 % Check model sizes
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 325
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 326 switch class(mdl)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 327 case 'matrix'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 328 % set Xvals for model (frequency vector could change for each experiment)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 329 for i = 1:numel(mdl.objs)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 330 % set Xvals
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 331 if (mdlFreqDependent)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 332 mdl.objs(i).setXvals(freqs);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 333 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 334 mdl.objs(i).setXvals(1);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 335 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 336 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 337 case 'ssm'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 338 spl = plist('set', 'for bode', ...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 339 'outputs', outNames, ...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 340 'inputs', inNames);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 341 % first optimise our model for the case in hand
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 342 mdl.reorganize(spl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 343 % make it lighter
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 344 for k = 1:nexp
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 345 mmdl(k).optimiseForFitting();
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 346 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 347 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 348 mmdl(k) = mdl;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 349
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 350 elseif ((numel(out(1).objs) == 1) && (numel(in(1).objs) == 1))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 351
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 352 % case 1 input, 1 output used mostly for passing
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 353 % ao objects from ao.mcmc.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 354
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 355 % optional resampling (not matrix/resample method yet)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 356 if ~isempty(fsout)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 357 in(k).objs(1).resample(plist('fsout',fsout));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 358 out(k).objs(1).resample(plist('fsout',fsout));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 359 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 360
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 361 % fft
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 362 fin(k) = fft(in(k));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 363 fout(k) = fft(out(k));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 364
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 365 Nsamples = length(fin(k).getObjectAtIndex(1).x);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 366 fs = fin(k).getObjectAtIndex(1).fs;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 367
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 368
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 369 if (isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 370 % Get rid of fft f =0, reduce frequency range if needed
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 371 if ~isempty(flim)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 372 fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 373 fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 374 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 375 elseif (~isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 376 % Get rid of fft f =0, reduce frequency range if needed
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 377 if ~isempty(flim)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 378 fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 379 fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 380 if(~isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 381 for lll=1:size(outModel,1)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 382 for kkk=1:size(outModel,2)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 383 outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 384 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 385 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 386 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 387 if(~isempty(inModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 388 inModel = split(inModel,plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 389 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 390 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 391 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 392
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 393 % use signal fft to get frequency vector. Take into account signal
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 394 % could be empty or set to zero
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 395
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 396 % 1st and only channel
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 397 if all(fin(k).getObjectAtIndex(1,1).y == 0) || isempty(fin(k).getObjectAtIndex(1,1).y)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 398 i1 = ao(plist('type','fsdata','xvals',0,'yvals',0));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 399 % rebuild input with new zeroed signal
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 400 fin(k) = matrix(i1,fin(k).getObjectAtIndex(2,1),plist('shape',[2 1]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 401 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 402 freqs = fin(k).getObjectAtIndex(1,1).x;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 403 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 404
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 405 % Compute psd
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 406 n1 = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(1,1), pl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 407
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 408 % interpolate noise to fft frequencies
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 409 S11 = interp(n1,plist('vertices',freqs));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 410
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 411 % build cross-spectrum matrix
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 412 mnse(k) = matrix(S11);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 413
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 414 switch class(mdl)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 415 case 'matrix'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 416 for i = 1:numel(mdl.objs)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 417 if (mdlFreqDependent)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 418 % set Xvals
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 419 mdl.objs(i).setXvals(freqs);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 420 % set alias
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 421 mdl.objs(i).assignalias(mdl.objs(i),plist('xvals',freqs));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 422 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 423 mdl.objs(i).setXvals(1);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 424 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 425 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 426 case 'ssm'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 427 if k<2
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 428 mdl.clearNumParams;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 429 spl = plist('set', 'for bode', ...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 430 'outputs', outNames, ...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 431 'inputs', inNames);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 432 % first optimise our model for the case in hand
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 433 mdl.reorganize(spl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 434 % make it lighter
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 435 mdl(k).optimiseForFitting();
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 436 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 437 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 438 % set Xvals for model (frequency vector could change for each experiment)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 439
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 440 mmdl(k) = mdl;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 441
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 442
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 443 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 444 error('Implemented cases: 1 inputs / 1output, 2 inputs / 2outputs (TN3045 analysis), and 4 inputs / 3 outpus (magnetic complete analysis model. Other cases have not been implemented yet. Sorry for the inconvenience)');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 445 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 446 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 447
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 448 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 449
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 450 % do simplex
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 451 if simplex
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 452 if isempty(xo)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 453 error('### Simplex needs a starting guess. Please input a ''x0''.');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 454 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 455 switch class(mmdl)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 456 case 'matrix'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 457 xo = fminsearch(@(x) utils.math.loglikelihood_matrix(x,fin,fout,mnse,mmdl,param,inModel,outModel),xo,optimset('Display','iter'));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 458 % case 'smodel'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 459 % xo = fminsearch(@(x) utils.math.loglikelihood(x,in,out,nse,mmdl,param),xo,optimset('Display','iter'));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 460 case 'ssm'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 461 Amats = mdl.amats;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 462 Bmats = mdl.bmats;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 463 Cmats = mdl.cmats;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 464 Dmats = mdl.dmats;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 465 xo = fminsearch(@(x) utils.math.loglikelihood_ssm(x,fin,fout,mnse,mmdl,param,inNames,outNames,spl,Amats,Bmats,Cmats,Dmats),xo,optimset('Display','iter'));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 466 otherwise
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 467 error('### Model must be either from the ''smodel'' or the ''ssm'' class. Please check the inputs')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 468 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 469
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 470 for i = 1:numel(param)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 471 fprintf('### Simplex estimate: %s = %d \n',param{i},xo(i))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 472 save('parameters_simplex.txt','xo','-ASCII')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 473 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 474 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 475
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 476 p = pest(xo);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 477 p.setName('SimplexEstimate');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 478 p.setNames(param{:});
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 479 p.setModels(mmdl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 480
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 481
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 482 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 483
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 484 if mhsampleTrue
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 485
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 486 % Sampling and saving rejected points. smplr is a matrix of
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 487 % the size of (# of samples)x(# of params + 2) wich contains the whole
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 488 % chain for each parameter, a column of zeros and ones (rejected or
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 489 % accepted point) and the value of the loglikelihood.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 490 [smpl smplr] = utils.math.mhsample(mmdl,fin,fout,mnse,cvar,N,rang,param,Tc,xi,xo,search,jumps,parplot,debug,inNames,outNames,fpars,anneal,SNR0,DeltaL,inModel,outModel);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 491
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 492 % statistics of the chain
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 493 if isempty(Tc)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 494 initial =1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 495 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 496 initial = Tc(2)+1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 497 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 498 mn = mean(smpl(initial:end,:));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 499 cv = cov(smpl(initial:end,:));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 500
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 501 for i = 1:nparam
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 502 for j = 1:nparam
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 503 cr(i,j) = cv(i,j)/sqrt(cv(i,i)*cv(j,j));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 504 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 505 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 506
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 507 % compute dof
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 508 ndata = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 509 for i = 1:length(mtxs)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 510 for j = 1:length(mtxs(1).objs)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 511 ndata = ndata + mtxs(i).objs(j).len;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 512 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 513 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 514 dof = ndata-nparam;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 515
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 516 % create pest output
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 517 p = pest(mn);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 518 p.setName('mcmc');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 519 p.setNames(param{:});
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 520 % add statistical info
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 521 p.setCov(cv);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 522 p.setCorr(cr);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 523 p.setDy(sqrt(diag(cv)));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 524 p.setChain(smpl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 525 %p.computePdf;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 526 p.setDof(dof);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 527 p.setModels(mdl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 528 % set history
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 529 %p.addHistory(getInfo('None'), pl, mtxs_invars(:), [out(:).hist mdl(:).hist]);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 530 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 531 % Set output
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 532 % Implemented it this way in order to leave the pest class as it is.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 533 output = {p,smplr};
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 534 %varargout{1} = p;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 535 varargout = output(1:nargout);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 536
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 537 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 538
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 539 % %--------------------------------------------------------------------------
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 540 % % Compute priors means and sigmas
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 541 % %--------------------------------------------------------------------------
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 542 % function fit = fitPrior(prior,nparam)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 543 % fit = [];
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 544 % for ii=1:nparam
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 545 % [muhat,sigmahat] = normfit(prior(:,2*ii-1),prior(:,2*ii));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 546 % fit = [fit ,[muhat,sigmahat]'];
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 547 % end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 548 % end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 549
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 550 %--------------------------------------------------------------------------
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 551 % Get Info Object
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 552 %--------------------------------------------------------------------------
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 553 function ii = getInfo(varargin)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 554 if nargin == 1 && strcmpi(varargin{1}, 'None')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 555 sets = {};
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 556 pl = [];
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 557 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 558 sets = {'Default'};
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 559 pl = getDefaultPlist;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 560 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 561 % Build info object
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 562 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: mcmc.m,v 1.35 2011/11/16 15:21:13 nikos Exp $', sets, pl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 563 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 564
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 565 %--------------------------------------------------------------------------
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 566 % Get Default Plist
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 567 %--------------------------------------------------------------------------
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 568 function plout = getDefaultPlist()
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 569 persistent pl;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 570 if exist('pl', 'var')==0 || isempty(pl)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 571 pl = buildplist();
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 572 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 573 plout = pl;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 574 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 575
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 576 function pl = buildplist()
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 577 pl = plist();
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 578
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 579 % % inNames
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 580 % p = param({'inNames','Input names. Used for ssm models'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 581 % pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 582 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 583 % % outNames
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 584 % p = param({'outNames','Output names. Usde for ssm models'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 585 % pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 586 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 587 % Model
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 588 % p = param({'model','A matrix array of models.'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 589 % pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 590
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 591 % Param
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 592 p = param({'FitParams','A cell array of evaluated parameters.'}, paramValue.EMPTY_DOUBLE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 593 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 594
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 595 % Input
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 596 p = param({'input','A matrix array of input signals.'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 597 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 598
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 599 pl = plist.MCH_FIT_PLIST;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 600
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 601 % N
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 602 p = param({'N','number of samples of the chain.'}, paramValue.DOUBLE_VALUE(1000));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 603 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 604
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 605 % Sigma
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 606 p = param({'cov','covariance of the gaussian jumping distribution.'}, paramValue.DOUBLE_VALUE(1e-4));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 607 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 608
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 609 % Noise
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 610 p = param({'noise','A matrix array of noise spectrum (PSD) used to compute the likelihood.'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 611 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 612
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 613 % Search
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 614 p = param({'modelFreqDependent','Set to true to use frequency dependent s models, set to false when using constant models'}, paramValue.TRUE_FALSE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 615 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 616
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 617 % Search
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 618 p = param({'search','Set to true to use bigger jumps in parameter space during annealing and cool down.'}, paramValue.TRUE_FALSE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 619 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 620
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 621 % Frequencies
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 622 p = param({'frequencies','Range of frequencies where the analysis is performed. If an array, only first and last are used'}, paramValue.EMPTY_DOUBLE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 623 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 624
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 625 % Resample
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 626 p = param({'fsout','Desired sampling frequency to resample the input time series'}, paramValue.EMPTY_DOUBLE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 627 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 628
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 629 % Simplex
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 630 p = param({'simplex','Set to true to perform a simplex search to find the starting parameters of the MCMC chain.'}, paramValue.TRUE_FALSE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 631 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 632
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 633 % mhsampleTrue
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 634 p = param({'mhsample','Set to true to perform a mhsample search. This is set to true by default. Only to be set to false by the user if we does not want to perform the mcmc search'}, paramValue.TRUE_FALSE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 635 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 636
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 637 % heat
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 638 p = param({'heat','The heat index flattening likelihood surface during annealing.'}, paramValue.DOUBLE_VALUE(1));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 639 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 640
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 641 % Tc
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 642 p = param({'Tc','An array of two values setting the initial and final value for the cooling down.'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 643 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 644
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 645 % heat
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 646 p = param({'x0','The proposed initial values.'}, paramValue.EMPTY_DOUBLE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 647 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 648
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 649 % jumps
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 650 p = param({'jumps','An array of four numbers setting the rescaling of the covariance matrix during the search phase.',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 651 'The first value is the one applied by default, the following thhree apply just when the chain sample is',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 652 'mod(10), mod(25) and mod(100) respectively.'}, paramValue.EMPTY_DOUBLE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 653 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 654
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 655 % heat
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 656 p = param({'plot','Select indexes of the parameters to be plotted.'}, paramValue.EMPTY_DOUBLE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 657 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 658
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 659 % debug
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 660 p = param({'debug','Set to true to get debug information of the MCMC process.'}, paramValue.FALSE_TRUE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 661 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 662
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 663 % inModel
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 664 p = param({'inModel','Input model. Still under test'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 665 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 666
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 667 % outModel
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 668 p = param({'outModel','Output model. Still under test'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 669 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 670
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 671 p = param({'prior','Mean, sigma and normalization factor for priors. Still under test'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 672 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 673
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 674 p = param({'anneal',['Choose type of annealing during sampling. Default value is ',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 675 'simulated annealing. Choose "thermo" for annealing with a thermostat.',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 676 ' SNR is computed and if it is larger than a fixed value SNR0 (provided also in the plist), ',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 677 'then the chains are heated by a factor of (SNR(1)/SNR0)^2. Choosing "simple" ',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 678 'the deviation of the loglikelihood of every 10 points in the chains is stored. If this deviation ',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 679 'is larger or smaller than two fixed values the chains are cooled or heated respectively.']}, {1, {'simul','thermo', 'simple'}, paramValue.SINGLE});
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 680 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 681
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 682 p = param({'SNR0','Fixed value for thermostated annealing.'}, {1, {200}, paramValue.OPTIONAL});
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 683 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 684
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 685 p = param({'DeltaL',['Deviation of Loglikelihood for 10 points of the chains. Used for',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 686 'the "simple" choice of annealing with a thermostat.']}, {1, {[100 600 2 3]}, paramValue.OPTIONAL});
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 687 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 688
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 689 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 690
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 691
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 692 % PARAMETERS: J - number of chains to compute.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 693 % sigma - dispersion of the jumping distribution.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 694 % range - range where the parameteters are sampled.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 695 % param - a cell array of evaluated parameters (from the
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 696 % smodel).
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 697 % Tc - an array with two values (default: 0, i.e. no annealing).
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 698 % heat - heat index flattening likelihood surface (default: 1)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 699 % x0 - proposed initial values, if empty selected randomly.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 700 % (default: empty)