-
Notifications
You must be signed in to change notification settings - Fork 0
/
pinky.m
281 lines (246 loc) · 6.9 KB
/
pinky.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
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
%Tristan Ursell
%2D Random Number Generator for a Given Discrete Distribution
%March 2012
%
%[x0,y0]=pinky(Xin,Yin,dist_in,varargin);
%
%'Xin' is a vector specifying the equally spaced values along the x-axis.
%
%'Yin' is a vector specifying the equally spaced values along the y-axis.
%
%'dist_in' (dist_in > 0) is a matrix with dimensions length(Yin) x
%length(Xin), whose values specify a 2D discrete probability distribution.
%The distribution does not need to be normalized.
%
%'res' (res > 1) is a multiplicative increase in the resolution of
%chosen random number, using cubic spline interpolation of the values in
%'dist_in'. Using the 'res' option can significantly slow down the code,
%due to the computational costs of interpolation, but allows one to
%generate more continuous values from the distribution.
%
%[x0,y0] is the output doublet of random numbers consistent with dist_in.
%
%The 'gendist' function required by this script, is included in this
%m-file.
%
%Example with an anisotropic Gaussian at native resolution:
%
% Xin=-10:0.1:10;
% Yin=-5:0.1:5;
%
% Xmat = ones(length(Yin),1)*Xin;
% Ymat = Yin'*ones(1,length(Xin));
%
% D1=1;
% D2=3;
% Dist=exp(-Ymat.^2/(2*D1^2)-Xmat.^2/(2*D2^2));
%
% N=10000;
%
% vals=zeros(2,N);
% for i=1:N
% [vals(1,i),vals(2,i)]=pinky(Xin,Yin,Dist);
% end
%
% figure;
% hold on
% imagesc(Xin,Yin,Dist)
% colormap(gray)
% plot(vals(1,:),vals(2,:),'r.')
% xlabel('Xin')
% ylabel('Yin')
% axis equal tight
% box on
%
%Example with multiple peaks at 10X resolution:
%
% Xin=-10:0.1:10;
% Yin=-5:0.1:5;
%
% Xmat = ones(length(Yin),1)*Xin;
% Ymat = Yin'*ones(1,length(Xin));
%
% D1=0.5;
% D2=1;
% Dist=exp(-(Ymat-3).^2/(4*D2^2)-(Xmat+6).^2/(2*D1^2))+exp(-(Ymat+2).^2/(2*D1^2)-(Xmat+1).^2/(2*D2^2))+exp(-(Ymat-1).^2/(2*D1^2)-(Xmat-2).^2/(2*D2^2));
%
% N=10000;
% res=10;
% vals=zeros(2,N);
% for i=1:N
% [vals(1,i),vals(2,i)]=pinky(Xin,Yin,Dist,res);
% end
%
% figure;
% hold on
% imagesc(Xin,Yin,Dist)
% colormap(gray)
% plot(vals(1,:),vals(2,:),'r.')
% xlabel('Xin')
% ylabel('Yin')
% axis equal tight
% box on
function [x0,y0]=pinky(Xin,Yin,dist_in,varargin)
%check input
if length(size(dist_in))>2
error('The input must be a N x M matrix.')
end
%check sizes
[sy,sx]=size(dist_in);
if or(length(Xin)~=sx,length(Yin)~=sy)
error('Dimensions of input vectors and input matrix must match.')
end
%check values
if any(dist_in(:)<0)
error('All input probability values must be positive.')
end
%get res
if nargin==4
res=varargin{1};
if res<=1
error('The resolution factor (res) must be an integer greater than one.')
end
elseif nargin~=3
error('Incorrect number of input arguments.')
end
%create column distribution and pick random number
col_dist=sum(dist_in,1);
%pick column distribution type
if nargin==3;
%if no res parameter, simply update X/Yin2
col_dist=col_dist/sum(col_dist);
Xin2=Xin;
Yin2=Yin;
else
%generate new, higher res input vectors
Xin2=linspace(min(Xin),max(Xin),round(res*length(Xin)));
Yin2=linspace(min(Yin),max(Yin),round(res*length(Yin)));
%generate interpolated column-sum distribution
col_dist=interp1(Xin,col_dist,Xin2,'pchip');
%check to make sure interpolated values are positive
if any(col_dist<0)
col_dist=abs(col_dist);
warning('Interpolation generated negative probability values.')
end
col_dist=col_dist/sum(col_dist);
end
%generate random value index
ind1=gendist(col_dist,1,1);
%save first value
x0=Xin2(ind1);
%find corresponding indices and weights in the other dimension
[val_temp,ind_temp]=sort((x0-Xin).^2);
if val_temp(1)<eps %if we land on an original value
row_dist=dist_in(:,ind_temp(1));
else %if we land inbetween, perform linear interpolation
low_val=min(ind_temp(1:2));
high_val=max(ind_temp(1:2));
Xlow=Xin(low_val);
Xhigh=Xin(high_val);
w1=1-(x0-Xlow)/(Xhigh-Xlow);
w2=1-(Xhigh-x0)/(Xhigh-Xlow);
row_dist=w1*dist_in(:,low_val) + w2*dist_in(:,high_val);
end
%pick column distribution type
if nargin==3;
row_dist=row_dist/sum(row_dist);
else
row_dist=interp1(Yin,row_dist,Yin2,'pchip');
row_dist=row_dist/sum(row_dist);
end
%generate random value index
ind2=gendist(row_dist,1,1);
%save first value
y0=Yin2(ind2);
function T = gendist(P,N,M,varargin)
%
%GENDIST - generate random numbers according to a discrete probability
%distribution.
%Tristan Ursell, (c) 2012
%
%T = gendist(P,N,M)
%T = gendist(P,N,M,'plot')
%
%The function gendist(P,N,M) takes in a positive vector P whose values
%form a discrete probability distribution for the indices of P. The
%function outputs an N x M matrix of integers corresponding to the indices
%of P chosen at random from the given underlying distribution.
%
%P will be normalized, if it is not normalized already. Both N and M must
%be greater than or equal to 1. The optional parameter 'plot' creates a
%plot displaying the input distribution in red and the generated points as
%a blue histogram.
%
%Conceptual EXAMPLE:
%
%If P = [0.2 0.4 0.4] (note sum(P)=1), then T can only take on values of 1,
%2 or 3, corresponding to the possible indices of P. If one calls
%gendist(P,1,10), then on average the output T should contain two 1's, four
%2's and four 3's, in accordance with the values of P, and a possible
%output might look like:
%
%T = gendist(P,1,10)
%T = 2 2 2 3 3 3 1 3 1 3
%
%If, for example, P is a probability distribution for position, with
%positions X = [5 10 15] (does not have to be linearly spaced), then the
%set of generated positions is X(T).
%
%EXAMPLE 1:
%
%P = rand(1,50);
%T = gendist(P,100,1000,'plot');
%
%EXAMPLE 2:
%
%X = -3:0.1:3;
%P = 1+sin(X);
%T = gendist(P,100,1000,'plot');
%Xrand = X(T);
%
%Note:
%size(T) = 100 1000
%size(Xrand) = 100 1000
%
%Thanks to Derek O'Connor for tips on speeding up the code.
%check for input errors
if and(nargin~=3,nargin~=4)
error('Error: Invalid number of input arguments.')
end
if min(P)<0
error('Error: All elements of first argument, P, must be positive.')
end
if size(P,1)>1
P=P';
end
if or(N<1,M<1)
error('Error: Output matrix dimensions must be greater than or equal to one.')
end
%normalize P
Pnorm=[0 P]/sum(P);
%create cumlative distribution
Pcum=cumsum(Pnorm);
%create random matrix
N=round(N);
M=round(M);
R=rand(1,N*M);
%calculate T output matrix
V=1:length(P);
[~,inds] = histc(R,Pcum);
T = V(inds);
%shape into output matrix
T=reshape(T,N,M);
%if desired, output plot
if nargin==4
if strcmp(varargin{1},'plot')
Pfreq=N*M*P/sum(P);
figure;
hold on
hist(T(T>0),1:length(P))
plot(Pfreq,'r-o')
ylabel('Frequency')
xlabel('P-vector Index')
axis tight
box on
end
end