使用矢量化或其他方法重新制定for循环-八度

Question

有没有办法向量化(或重新制定)这段代码中的每个循环体：

col=load('col-deau'); %load data

h=col(:,8); % corresponding water column
dates=col(:,3); % and its dates

%removing out-of-bound data
days=days(h~=9999.000);
h=h(h~=9999.000);
dates=sort(dates(h~=9999.000));

[k,hcat]=hist(h,nbin); %making classes (k) and boundaries of classes (hcat) of water column automatically

dcat=1:15; % make boundaries for dates
for k=1:length(dcat)-1 % Loop for each date class
    ii=find(dates>=dcat(k)&dates<dcat(k+1));% Counting dates corresponding to the boundaries of each date class
    for j=1:length(hcat)-1                                % Loop over each class of water column
        ij=find(h>=hcat(j)&h<hcat(j+1)); % Count water column corresponding to the boundaries of each water column class
        obs(k,j)=length(intersect(ii,ij));               % Find the size of each intersecting matrix
    end
end

例如，我尝试使用矢量化来更改这一部分：

for k=1:length(dcat)-1        
    ii=find(dates>=dcat(k)&dates<dcat(k+1))  
endfor

有了这个：

nk=1:length(dcat)-1;
ii2=find(dates>=dcat(nk)&dates<dcat(nk+1));

还可以使用bsxfun：

ii2=find(bsxfun(@and,bsxfun(@ge,dates,nk),bsxfun(@lt,dates,nk+1)));

但无济于事。这两种方法产生相同的输出，并且不对应于使用for循环的输出(就元素和向量大小而言)。

对于信息，h是一个矢量，它包含以米为单位的水柱，日期是一个矢量(两位数的整数)，它包含对相应水柱进行测量的日期。可以在以下位置找到输入文件：https://drive.google.com/open?id=1EomLGYleaNtiGG2iV_9LRt425blxdIsm

至于输出，我想要像这样的ii：

ii =

   1177
   1178
   1179
   1180
   1181
   1182
   1183
   1184
   1185
   1186
   1187
   1188
   1189
   1190
   1191
   1192
   1193
   1194
   1195
   1196
   1197
   1198
   1199
   1200
   1201
   1202
   1203
   1204
   1205
   1206
   1207
   1208
   1209
   1210
   1211
   1212
   1213
   1214
   1215
   1216
   1217
   1218
   1219
   1220
   1221
   1222
   1223
   1224
   1225
   1226
   1227
   1228
   1229
   1230
   1231
   1232
   1233
   1234
   1235
   1236
   1237
   1238
   1239
   1240
   1241
   1242
   1243
   1244
   1245
   1246
   1247
   1248
   1249
   1250
   1251
   1252
   1253
   1254
   1255
   1256
   1257
   1258
   1259
   1260
   1261
   1262
   1263
   1264
   1265
   1266
   1267
   1268
   1269
   1270
   1271
   1272

相反，使用第一种方法，我得到了ii2，它在值和向量大小方面有很大的不同(我不能发布结果，因为向量太大了)。

有人能帮助一个绝望的新手吗？我只需要将循环部分重新表述为更好、更简洁的版本。

如果需要添加更多细节，请随时向我提问。

Answer 1

您可以使用hist3

pkg load statistics
[obs, ~] = hist3([dates(:) h(:)] ,'Edges', {dcat,hcat});