4:- use_module(library(real)). 5:- use_module(library(mcintyre)). 6
7:- if(current_predicate(use_rendering/1)). 8:- use_rendering(c3). 9:- endif. 10:- mc. 11:- begin_lpad. 12
13iris_par(I,par(D,Means),X1,X2,X3,X4):-
14 latent_class(I,D,K),
15 nth1(K,Means,[M1,M2,M3,M4,V1,V2,V3,V4]),
16 measure(1,I,K,M1,V1,X1),
17 measure(2,I,K,M2,V2,X2),
18 measure(3,I,K,M3,V3,X3),
19 measure(4,I,K,M4,V4,X4).
20
21iris(I,X1,X2,X3,X4):-
22 category(I,K),
23 data_mean_var(1,DM1,_),
24 data_mean_var(2,DM2,_),
25 data_mean_var(3,DM3,_),
26 data_mean_var(4,DM4,_),
27 mean(K,x1,DM1,M1),
28 mean(K,x2,DM2,M2),
29 mean(K,x3,DM3,M3),
30 mean(K,x4,DM4,M4),
31 measure(x1,I,K,M1,X1),
32 measure(x2,I,K,M2,X2),
33 measure(x3,I,K,M3,X3),
34 measure(x4,I,K,M4,X4).
35
36category(I,K):-
37 prior_cat([1,1,1],Par),
38 latent_class(I,Par,K).
39
40prior_cat(Alfas,Values):dirichlet(Values,Alfas).
41
42
43mean(_,_,DM,M): gaussian(M,DM,1).
44
45measure(_,_,_,M,A): gaussian(A,M,3.0).
46measure(_,_,_,M,Var,A): gaussian(A,M,Var).
47
48latent_class(I,[A,B,C],V):-
50 latc(I,[1:A,2:B,3:C],V).
51
52latc(_,D,V):discrete(V,D).
53
54pair_val(V,P,V:P).
55
56:- end_lpad. 57data(5.1,3.5,1.4,0.2,'Iris-setosa').
58data(4.9,3.0,1.4,0.2,'Iris-setosa').
59data(4.7,3.2,1.3,0.2,'Iris-setosa').
60data(4.6,3.1,1.5,0.2,'Iris-setosa').
61data(5.0,3.6,1.4,0.2,'Iris-setosa').
62data(5.4,3.9,1.7,0.4,'Iris-setosa').
63data(4.6,3.4,1.4,0.3,'Iris-setosa').
64data(5.0,3.4,1.5,0.2,'Iris-setosa').
65data(4.4,2.9,1.4,0.2,'Iris-setosa').
66data(4.9,3.1,1.5,0.1,'Iris-setosa').
67data(5.4,3.7,1.5,0.2,'Iris-setosa').
68data(4.8,3.4,1.6,0.2,'Iris-setosa').
69data(4.8,3.0,1.4,0.1,'Iris-setosa').
70data(4.3,3.0,1.1,0.1,'Iris-setosa').
71data(5.8,4.0,1.2,0.2,'Iris-setosa').
72data(5.7,4.4,1.5,0.4,'Iris-setosa').
73data(5.4,3.9,1.3,0.4,'Iris-setosa').
74data(5.1,3.5,1.4,0.3,'Iris-setosa').
75data(5.7,3.8,1.7,0.3,'Iris-setosa').
76data(5.1,3.8,1.5,0.3,'Iris-setosa').
77data(5.4,3.4,1.7,0.2,'Iris-setosa').
78data(5.1,3.7,1.5,0.4,'Iris-setosa').
79data(4.6,3.6,1.0,0.2,'Iris-setosa').
80data(5.1,3.3,1.7,0.5,'Iris-setosa').
81data(4.8,3.4,1.9,0.2,'Iris-setosa').
82data(5.0,3.0,1.6,0.2,'Iris-setosa').
83data(5.0,3.4,1.6,0.4,'Iris-setosa').
84data(5.2,3.5,1.5,0.2,'Iris-setosa').
85data(5.2,3.4,1.4,0.2,'Iris-setosa').
86data(4.7,3.2,1.6,0.2,'Iris-setosa').
87data(4.8,3.1,1.6,0.2,'Iris-setosa').
88data(5.4,3.4,1.5,0.4,'Iris-setosa').
89data(5.2,4.1,1.5,0.1,'Iris-setosa').
90data(5.5,4.2,1.4,0.2,'Iris-setosa').
91data(4.9,3.1,1.5,0.1,'Iris-setosa').
92data(5.0,3.2,1.2,0.2,'Iris-setosa').
93data(5.5,3.5,1.3,0.2,'Iris-setosa').
94data(4.9,3.1,1.5,0.1,'Iris-setosa').
95data(4.4,3.0,1.3,0.2,'Iris-setosa').
96data(5.1,3.4,1.5,0.2,'Iris-setosa').
97data(5.0,3.5,1.3,0.3,'Iris-setosa').
98data(4.5,2.3,1.3,0.3,'Iris-setosa').
99data(4.4,3.2,1.3,0.2,'Iris-setosa').
100data(5.0,3.5,1.6,0.6,'Iris-setosa').
101data(5.1,3.8,1.9,0.4,'Iris-setosa').
102data(4.8,3.0,1.4,0.3,'Iris-setosa').
103data(5.1,3.8,1.6,0.2,'Iris-setosa').
104data(4.6,3.2,1.4,0.2,'Iris-setosa').
105data(5.3,3.7,1.5,0.2,'Iris-setosa').
106data(5.0,3.3,1.4,0.2,'Iris-setosa').
107data(7.0,3.2,4.7,1.4,'Iris-versicolor').
108data(6.4,3.2,4.5,1.5,'Iris-versicolor').
109data(6.9,3.1,4.9,1.5,'Iris-versicolor').
110data(5.5,2.3,4.0,1.3,'Iris-versicolor').
111data(6.5,2.8,4.6,1.5,'Iris-versicolor').
112data(5.7,2.8,4.5,1.3,'Iris-versicolor').
113data(6.3,3.3,4.7,1.6,'Iris-versicolor').
114data(4.9,2.4,3.3,1.0,'Iris-versicolor').
115data(6.6,2.9,4.6,1.3,'Iris-versicolor').
116data(5.2,2.7,3.9,1.4,'Iris-versicolor').
117data(5.0,2.0,3.5,1.0,'Iris-versicolor').
118data(5.9,3.0,4.2,1.5,'Iris-versicolor').
119data(6.0,2.2,4.0,1.0,'Iris-versicolor').
120data(6.1,2.9,4.7,1.4,'Iris-versicolor').
121data(5.6,2.9,3.6,1.3,'Iris-versicolor').
122data(6.7,3.1,4.4,1.4,'Iris-versicolor').
123data(5.6,3.0,4.5,1.5,'Iris-versicolor').
124data(5.8,2.7,4.1,1.0,'Iris-versicolor').
125data(6.2,2.2,4.5,1.5,'Iris-versicolor').
126data(5.6,2.5,3.9,1.1,'Iris-versicolor').
127data(5.9,3.2,4.8,1.8,'Iris-versicolor').
128data(6.1,2.8,4.0,1.3,'Iris-versicolor').
129data(6.3,2.5,4.9,1.5,'Iris-versicolor').
130data(6.1,2.8,4.7,1.2,'Iris-versicolor').
131data(6.4,2.9,4.3,1.3,'Iris-versicolor').
132data(6.6,3.0,4.4,1.4,'Iris-versicolor').
133data(6.8,2.8,4.8,1.4,'Iris-versicolor').
134data(6.7,3.0,5.0,1.7,'Iris-versicolor').
135data(6.0,2.9,4.5,1.5,'Iris-versicolor').
136data(5.7,2.6,3.5,1.0,'Iris-versicolor').
137data(5.5,2.4,3.8,1.1,'Iris-versicolor').
138data(5.5,2.4,3.7,1.0,'Iris-versicolor').
139data(5.8,2.7,3.9,1.2,'Iris-versicolor').
140data(6.0,2.7,5.1,1.6,'Iris-versicolor').
141data(5.4,3.0,4.5,1.5,'Iris-versicolor').
142data(6.0,3.4,4.5,1.6,'Iris-versicolor').
143data(6.7,3.1,4.7,1.5,'Iris-versicolor').
144data(6.3,2.3,4.4,1.3,'Iris-versicolor').
145data(5.6,3.0,4.1,1.3,'Iris-versicolor').
146data(5.5,2.5,4.0,1.3,'Iris-versicolor').
147data(5.5,2.6,4.4,1.2,'Iris-versicolor').
148data(6.1,3.0,4.6,1.4,'Iris-versicolor').
149data(5.8,2.6,4.0,1.2,'Iris-versicolor').
150data(5.0,2.3,3.3,1.0,'Iris-versicolor').
151data(5.6,2.7,4.2,1.3,'Iris-versicolor').
152data(5.7,3.0,4.2,1.2,'Iris-versicolor').
153data(5.7,2.9,4.2,1.3,'Iris-versicolor').
154data(6.2,2.9,4.3,1.3,'Iris-versicolor').
155data(5.1,2.5,3.0,1.1,'Iris-versicolor').
156data(5.7,2.8,4.1,1.3,'Iris-versicolor').
157data(6.3,3.3,6.0,2.5,'Iris-virginica').
158data(5.8,2.7,5.1,1.9,'Iris-virginica').
159data(7.1,3.0,5.9,2.1,'Iris-virginica').
160data(6.3,2.9,5.6,1.8,'Iris-virginica').
161data(6.5,3.0,5.8,2.2,'Iris-virginica').
162data(7.6,3.0,6.6,2.1,'Iris-virginica').
163data(4.9,2.5,4.5,1.7,'Iris-virginica').
164data(7.3,2.9,6.3,1.8,'Iris-virginica').
165data(6.7,2.5,5.8,1.8,'Iris-virginica').
166data(7.2,3.6,6.1,2.5,'Iris-virginica').
167data(6.5,3.2,5.1,2.0,'Iris-virginica').
168data(6.4,2.7,5.3,1.9,'Iris-virginica').
169data(6.8,3.0,5.5,2.1,'Iris-virginica').
170data(5.7,2.5,5.0,2.0,'Iris-virginica').
171data(5.8,2.8,5.1,2.4,'Iris-virginica').
172data(6.4,3.2,5.3,2.3,'Iris-virginica').
173data(6.5,3.0,5.5,1.8,'Iris-virginica').
174data(7.7,3.8,6.7,2.2,'Iris-virginica').
175data(7.7,2.6,6.9,2.3,'Iris-virginica').
176data(6.0,2.2,5.0,1.5,'Iris-virginica').
177data(6.9,3.2,5.7,2.3,'Iris-virginica').
178data(5.6,2.8,4.9,2.0,'Iris-virginica').
179data(7.7,2.8,6.7,2.0,'Iris-virginica').
180data(6.3,2.7,4.9,1.8,'Iris-virginica').
181data(6.7,3.3,5.7,2.1,'Iris-virginica').
182data(7.2,3.2,6.0,1.8,'Iris-virginica').
183data(6.2,2.8,4.8,1.8,'Iris-virginica').
184data(6.1,3.0,4.9,1.8,'Iris-virginica').
185data(6.4,2.8,5.6,2.1,'Iris-virginica').
186data(7.2,3.0,5.8,1.6,'Iris-virginica').
187data(7.4,2.8,6.1,1.9,'Iris-virginica').
188data(7.9,3.8,6.4,2.0,'Iris-virginica').
189data(6.4,2.8,5.6,2.2,'Iris-virginica').
190data(6.3,2.8,5.1,1.5,'Iris-virginica').
191data(6.1,2.6,5.6,1.4,'Iris-virginica').
192data(7.7,3.0,6.1,2.3,'Iris-virginica').
193data(6.3,3.4,5.6,2.4,'Iris-virginica').
194data(6.4,3.1,5.5,1.8,'Iris-virginica').
195data(6.0,3.0,4.8,1.8,'Iris-virginica').
196data(6.9,3.1,5.4,2.1,'Iris-virginica').
197data(6.7,3.1,5.6,2.4,'Iris-virginica').
198data(6.9,3.1,5.1,2.3,'Iris-virginica').
199data(5.8,2.7,5.1,1.9,'Iris-virginica').
200data(6.8,3.2,5.9,2.3,'Iris-virginica').
201data(6.7,3.3,5.7,2.5,'Iris-virginica').
202data(6.7,3.0,5.2,2.3,'Iris-virginica').
203data(6.3,2.5,5.0,1.9,'Iris-virginica').
204data(6.5,3.0,5.2,2.0,'Iris-virginica').
205data(6.2,3.4,5.4,2.3,'Iris-virginica').
206data(5.9,3.0,5.1,1.8,'Iris-virginica').
207
208
209em(It,True,Post):-
210 findall([A,B,C,D,Type],data(A,B,C,D,Type),LA),
211 pca(LA,True),
212 data_mean_var(1,DM1,V1),
213 data_mean_var(2,DM2,V2),
214 data_mean_var(3,DM3,V3),
215 data_mean_var(4,DM4,V4),
216 findall(M,(between(1,3,_),gauss(DM1,1.0,M)),[M11,M21,M31]),
217 findall(M,(between(1,3,_),gauss(DM2,1.0,M)),[M12,M22,M32]),
218 findall(M,(between(1,3,_),gauss(DM3,1.0,M)),[M13,M23,M33]),
219 findall(M,(between(1,3,_),gauss(DM4,1.0,M)),[M14,M24,M34]),
220 em_it(0,It,LA,par([1:0.2,2:0.5,3:0.3],
221 [M11,M12,M13,M14,V1,V2,V3,V4],[M21,M22,M23,M24,V1,V2,V3,V4],
222 [M31,M32,M33,M34,V1,V2,V3,V4]),
223 _Par,_,LA1),
224 228 find_most_likely(LA1,LPostCat),
229 maplist(replace_cat,LA,LPostCat,LPost),
230 pca(LPost,Post).
231
232find_most_likely([L1,L2,L3],LC):-
233 maplist(classify,L1,L2,L3,LC).
234
235classify(_-W1,_-W2,_-W3,Cat):-
236 find_max([W1,W2,W3],Cat).
237
238em_it(It,It,_LA,Par,Par,LAOut,LAOut):-!.
239
240em_it(It0,It,LA,Par0,Par,_,LAOut):-
241 expect(LA,Par0,LA1,LL),
242 write('LL '),write(LL),nl,
243 maxim(LA1,Par1),
244 It1 is It0+1,
245 em_it(It1,It,LA,Par1,Par,LA1,LAOut).
246
247expect(LA,par([1:P1,2:P2,3:P3],
248 G1,G2,G3),[L1,L2,L3],LL):-
249 maplist(weight(G1,P1),LA,L01),
250 maplist(weight(G2,P2),LA,L02),
251 maplist(weight(G3,P3),LA,L03),
252 normal(L01,L02,L03,L1,L2,L3),
253 log_lik(L01,L02,L3,P1,P2,P3,LL).
254
255normal(L01,L02,L03,L1,L2,L3):-
256 maplist(px,L01,L02,L03,L1,L2,L3).
257
258px(X-W01,X-W02,X-W03,X-W1,X-W2,X-W3):-
259 S is W01+W02+W03,
260 W1 is W01/S,
261 W2 is W02/S,
262 W3 is W03/S.
263
264weight([M1,M2,M3,M4,V1,V2,V3,V4],P,[A,B,C,D,_],[A,B,C,D]-W):-
265 gauss_density_0(M1,V1,A,W1),
266 gauss_density_0(M2,V2,B,W2),
267 gauss_density_0(M3,V3,C,W3),
268 gauss_density_0(M4,V4,D,W4),
269 W is W1*W2*W3*W4*P.
270
271gauss_density_0(M,V,X,W):-
272 (V=:=0.0->
273 write(zero_var_gauss),
274 W=0.0
275 ;
276 gauss_density(M,V,X,W)
277 ).
278
279
280maxim([LA1,LA2,LA3],par([1:P1,2:P2,3:P3],C1,C2,C3)):-
281 stats(LA1,W1,C1),
282 stats(LA2,W2,C2),
283 stats(LA3,W3,C3),
284 SW is W1+W2+W3,
285 write(sw),write(SW),nl,
286 P1 is W1/SW,
287 P2 is W2/SW,
288 P3 is W3/SW.
289
290log_lik(L1,L2,L3,P1,P2,P3,LL):-
291 foldl(combine(P1,P2,P3),L1,L2,L3,0,LL).
292
293combine(P1,P2,P3,_-W1,_-W2,_-W3,LL0,LL):-
294 LLs is log(P1*W1+P2*W2+P3*W3),
295 296 LL is LL0+LLs.
297
298sample_clusters(Samples,True,Prior,Post):-
299 findall([A,B,C,D,Type],data(A,B,C,D,Type),LA),
300 pca(LA,True),
301 findall((category(I,K),K),between(1,150,I),LCat),
302 maplist(split,LCat,G,Args),
303 G=[H|T],
304 foldl(totuple,T,H,Goal),
305 mc_sample_arg_raw(Goal,Samples,Args,ValPrior),
306 findall(Ar-1,member(Ar,ValPrior),ValPriorW),
307 find_most_prob(ValPriorW,LPriorCat),
308 maplist(replace_cat,LA,LPriorCat,LP),
309 pca(LP,Prior),
310 maplist(ev_list,LA,LCat,EvL),
311 EvL=[HE|TE],
312 foldl(totuple,TE,HE,Ev),
313 mc_lw_sample_arg_log(Goal,Ev,Samples,Args,ValPost),
315 find_most_prob_log(ValPost,LPostCat),
316 maplist(replace_cat,LA,LPostCat,LPost),
317 pca(LPost,Post).
318
319
320lw_sample(Ev,Samples,Goal,Arg,Val):-
321 mc_lw_sample_arg(Goal,Ev,Samples,Arg,Val).
322
323ev_list([A,B,C,D,_],(category(I,_),_),iris(I,A,B,C,D)).
324
325find_most_prob(Vals,MostProb):-
326 findall([0,0,0],between(1,150,_),Card0),
327 foldl(count,Vals,Card0,Card),
328 maplist(find_max,Card,MostProb).
329
330count(Cats-W,Card0,Card):-
331 maplist(update_counts(W),Cats,Card0,Card).
332
333update_counts(W,K,Card0,Card):-
334 up_c(1,K,W,Card0,Card).
335
336up_c(K,K,W,[H|R],[H1|R]):-!,
337 H1 is H+W.
338
339up_c(K0,K,W,[H|R0],[H|R]):-
340 K1 is K0+1,
341 up_c(K1,K,W,R0,R).
342
343find_most_prob_log(Vals,MostProb):-
344 findall([0,0,0],between(1,150,_),Card0),
345 foldl(count_log,Vals,Card0,Card),
346 maplist(find_max,Card,MostProb).
347
348count_log(Cats-W,Card0,Card):-
349 maplist(update_counts_log(W),Cats,Card0,Card).
350
351update_counts_log(W,K,Card0,Card):-
352 up_c_log(1,K,W,Card0,Card).
353
354up_c_log(K,K,W,[H|R],[H1|R]):-!,
355 H1 is H+exp(1000+W).
356
357up_c_log(K0,K,W,[H|R0],[H|R]):-
358 K1 is K0+1,
359 up_c_log(K1,K,W,R0,R).
360
361
362find_max(Counts,Max):-
363 max_list(Counts,MV),
364 nth1(Max,Counts,MV).
365
366replace_cat([A,B,C,D,_],Cat,[A,B,C,D,Cat]).
367
368totuple(G,GG,(GG,G)).
369split((A,B),A,B).
370
371pca(LA,P):-
372 maplist(add_cat,LA,LCat,L),
373 list_to_set(LCat,Cats),
374 data<- L,
375 dpca<- prcomp(data),
376 PC <- dpca,
377 member(x=Data,PC),
378 %findall(X-Y,(member(x=Data,PC),member([X,Y,_,_],Data)),LB),
379 maplist(add_cat,DC,LCat,Data),
380 maplist(separate(DC),Cats,Mat1),
381% maplist(xy,Setosa,XSet,YSet),
382% maplist(xy,Versicolor,XVer,YVer),
383% maplist(xy,Virginica,XVir,YVir),
384 maplist(keep2,Mat1,Mat2),
385 maplist(axis,Cats,Axis,LAxis),
386 dict_create(Ax,_,LAxis),
387 maplist(tocol,Mat2,Axis,ColData),
388 append(ColData,Cols),
389 P = c3{data:_{
390 xs:Ax,
391%x: 'Iris-setosa'_x,
392 columns:
393 Cols
394 ,
395% ["'Iris-versicolor'_x"|XVer],[ver|YVer],
396% ["'Iris-virginica'_x"|XVir],[vir|YVir]],
397 type:scatter},
398 axis:_{
399 x:_{
400 tick:_{
401 fit: false
402 }
403 }}
404% axis:_{ x:_{ tick:_{values:Tick}}},
405% format: 'function (x) { return x.toFixed(2);}' ,
406% fit: true,culling:_{max:7} }} },
407% bar:_{
408% width:_{ ratio: 1.0 }},
409% legend:_{show: false}
410}.
411
412tocol(Data,[X,Y],[[X|DX],[Y|DY]]):-
413 maplist(xy,Data,DX,DY).
414
415ab([X,Y,_,_,_],[X,Y]).
416xy([X,Y],X,Y).
417add_cat([X,Y,Z,W,C],C,[X,Y,Z,W]).
418set([_,_,_,_,'Iris-setosa']).
419cat(Cat,[_,_,_,_,Cat]).
420
421group([A,B],[C,D],[E,F],g(A,B,C,D,E,F)).
422
423axis(N,[NX,N],NA-NX):-
424 (number(N)->
425 atom_number(NA,N)
426 ;
427 NA=N
428 ),
429 atom_concat(x,NA,NX).
430
431keep2(Quad,Cou):-
432 maplist(ab,Quad,Cou).
433
434
435separate(DC,Cat,DataClass):-
436 include(cat(Cat),DC,DataClass).
437
438
439assert_data_means:-
440 findall([A,B,C,D],data(A,B,C,D,_Type),LA),
441 maplist(component,LA,CA,CB,CC,CD),
442 mean(CA,M1),
443 mean(CB,M2),
444 mean(CC,M3),
445 mean(CD,M4),
446 variance(CA,M1,V1),
447 variance(CB,M2,V2),
448 variance(CC,M3,V3),
449 variance(CD,M4,V4),
450 assert(data_mean_var(1,M1,V1)),
451 assert(data_mean_var(2,M2,V2)),
452 assert(data_mean_var(3,M3,V3)),
453 assert(data_mean_var(4,M4,V4)).
454
455mean(L,M):-
456 length(L,N),
457 sum_list(L,S),
458 M is S/N.
459
460variance(L,M,Var):-
461 length(L,N),
462 foldl(agg_var(M),L,0,S),
463 Var is S/N.
464
465stats(LA,SW,[M1,M2,M3,M4,V1,V2,V3,V4]):-
466 maplist(component_weight,LA,CA,CB,CC,CD),
467 weighted_mean(CA,M1,SW),
468 weighted_mean(CB,M2,_),
469 weighted_mean(CC,M3,_),
470 weighted_mean(CD,M4,_),
471 weighted_var(CA,M1,V1),
472 weighted_var(CB,M2,V2),
473 weighted_var(CC,M3,V3),
474 weighted_var(CD,M4,V4).
475
476weighted_var(L,M,Var):-
477 foldl(agg_val_var(M),L,(0,0),(S,SW0)),
478 SW is SW0,
479 (SW=:=0.0->
480 write(zero_var),nl,
481 Var=1.0
482 ;
483 Var is S/SW
484 ).
485
486weighted_mean(L,M,SW):-
487 foldl(agg_val,L,(0,0),(S,SW0)),
488 SW is SW0,
489 (SW =:=0.0->
490 write(zero_mean),nl,
491 M is 0
492 ;
493 M is S/SW
494 ).
495
496agg_val(V -N,(S,SW),(S+V*N,SW+N)).
497agg_val_var(M,V -N,(S,SW),(S+(M-V)^2*N,SW+N)).
498agg_var(M,V,S,S+(M-V)^2).
499
500
501component([A,B,C,D],A,B,C,D).
502component_weight([A,B,C,D]-W,A-W,B-W,C-W,D-W).
503
504:- assert_data_means.
?-
em(10,T,P).*/