-
Notifications
You must be signed in to change notification settings - Fork 6
/
ConvolveFFTcIID.m
182 lines (155 loc) · 6.52 KB
/
ConvolveFFTcIID.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
classdef ConvolveFFTcIID < dContinuous
% ConvolveFFTcIID(NDists,Basis,varargin)
% Convolution of NDists i.i.d. random variables with distribution BasisRV{1}
% NEWJEFF: There is a lot of code overlap with ConvolveFFTc, but it will be tough
% to do this right (even just initial constructors!).
% NEWJEFF: It would be nice to add a discrete approx. for convolveFFTcIID like convolveFFTc
properties (Constant)
defaultNxPoints = 2^16;
end
properties(SetAccess = private)
NxPoints
delx
BasisRV
NDists
BasisMeans
BasisVariances
KnownMeans
KnownVariances
end
methods
function obj=ConvolveFFTcIID(NDists,Basis,varargin)
obj=obj@dContinuous('ConvolveFFTcIID');
obj.ExtractNxPoints(varargin{:});
obj.BasisRV = Basis;
obj.NDists = NDists;
obj.NDistParms = obj.BasisRV.NDistParms;
obj.DefaultParmCodes = obj.BasisRV.DefaultParmCodes;
obj.KnownMeans = false;
obj.KnownVariances = false;
obj.ReInit;
end
function varargin=ExtractNxPoints(obj,varargin)
[obj.NxPoints, varargin] = ExtractNameVali('NxPoints',-1,varargin);
[obj.delx, varargin] = ExtractNameVali('delx',-1,varargin);
if (obj.NxPoints<0) && (obj.NxPoints<0)
obj.NxPoints = ConvolveFFTc.defaultNxPoints;
end
end
function BuildMyName(obj)
obj.StringName = [obj.FamilyName '(' num2str(obj.NDists) ',' obj.BasisRV.StringName ')'];
end
function []=ResetParms(obj,newparmvalues)
ClearBeforeResetParmsC(obj);
obj.BasisRV.ResetParms(newparmvalues);
ReInit(obj);
end
function PerturbParms(obj,ParmCodes)
obj.BasisRV.PerturbParms(ParmCodes);
obj.ReInit;
end
function Reals = ParmsToReals(obj,Parms,~)
Reals = obj.BasisRV.ParmsToReals(Parms);
end
function Parms = RealsToParms(obj,Reals,~)
Parms = obj.BasisRV.RealsToParms(Reals);
end
function parmvals = ParmValues(obj)
parmvals = obj.BasisRV.ParmValues;
end
function []=ReInit(obj)
% Using existing bases, find bounds, set up Xs, run FFT
% get lower/upper bound info about basis:
lowerb = obj.BasisRV.LowerBound;
upperb = obj.BasisRV.UpperBound;
% Compute corresponding properties of convolution
lowerconv = lowerb * obj.NDists;
upperconv = upperb * obj.NDists;
% Compute full range of X's spanning all bases & convolution
% The X's will reflect shifted versions of all RVs
% (i.e., shifted to have lower bounds of 0.
x_min = 0;
x_max = max([upperb-lowerb; upperconv-lowerconv]);
% Compute the actual X's
if obj.delx > 0
useNxPoints = ceil( (x_max - x_min) / obj.delx );
else
useNxPoints = obj.NxPoints;
end
X = linspace(x_min,x_max,useNxPoints);
% Load Y with PDFs of bases:
% I don't understand the various discussions about using fftshift & ifftshift
% but it seems correct without them.
Y = zeros(obj.NDists,useNxPoints);
for i=1:obj.NDists
% Add lowerb for pdf computation to compensate for
% the use of shifted X values
% Y(i,:) = ifftshift(obj.BasisRV.PDF(X+lowerb));
if i==1
Y(i,:) = obj.BasisRV.PDF(X+lowerb);
else
Y(i,:) = Y(i-1,:);
end
end
% Compute PDF of convolution via FFT
% I don't really understand this
fftpdf = real(ifft(prod(fft(Y,[],2))));
% fftpdf = fftshift(fftpdf);
fftpdf = fftpdf / trapz(X,fftpdf); % normalize pdf values for integral=1
% At this point the fftpdf values are correct, but the X vector still
% reflects the shifts, so remove those.
needed_shift = lowerconv;
X = X + needed_shift;
% Trim X and fftpdf vectors to get rid of leading & trailing zeros, if any.
% (zero's at each end are added below)
pos_pdfs = find(fftpdf>0);
first_pos_pdf = pos_pdfs(1);
last_pos_pdf = pos_pdfs(end);
X = X(first_pos_pdf:last_pos_pdf);
fftpdf = fftpdf(first_pos_pdf:last_pos_pdf);
% These are the final X and fftpdf values used for splinepdf
obj.SplinePDFsXs = X;
obj.SplinePDFs = fftpdf;
obj.PDFSplineInfo = spline(obj.SplinePDFsXs,[0 obj.SplinePDFs 0]); % Include end 0's to guarantee flatness at edges.
obj.UseSplinePDF = true;
obj.HaveSplinePDFs = true;
obj.LowerBound = X(1);
obj.UpperBound = X(end);
obj.Initialized = true;
if obj.NameBuilding
obj.BuildMyName;
end
end
function thisval=Mean(obj)
if ~obj.Initialized
error(UninitializedError(obj));
end
if ~obj.KnownMeans
obj.BasisMeans = obj.BasisRV.Mean;
obj.KnownMeans = true;
end
thisval = obj.NDists * obj.BasisMeans;
end
function thisval=Variance(obj)
if ~obj.Initialized
error(UninitializedError(obj));
end
if ~obj.KnownVariances
obj.BasisVariances = obj.BasisRV.Variance;
obj.KnownVariances = true;
end
thisval = obj.NDists * obj.BasisVariances;
end
function thisVal = Random(obj,varargin)
assert(obj.Initialized,UninitializedError(obj));
basisVals = zeros(varargin{:}); % make an array of the required size for each basis
nPerBasis = numel(basisVals);
AllBases = zeros(obj.NDists,nPerBasis);
for iDist=1:obj.NDists
AllBases(iDist,:) = obj.BasisRV.Random(1,nPerBasis);
end
sums = sum(AllBases,1);
thisVal = reshape(sums,size(basisVals));
end
end % methods
end