spm_DEM_qU.m 6.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238
  1. function spm_DEM_qU(qU,pU)
  2. % displays conditional estimates of states (qU)
  3. % FORMAT spm_DEM_qU(qU,pU);
  4. %
  5. % qU.v{i} - causal states (V{1} = y = predicted response)
  6. % qU.x{i} - hidden states
  7. % qU.e{i} - prediction error
  8. % qU.C{N} - conditional covariance - [causal states] for N samples
  9. % qU.S{N} - conditional covariance - [hidden states] for N samples
  10. %
  11. % pU - optional input for known states
  12. %__________________________________________________________________________
  13. % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
  14. % Karl Friston
  15. % $Id: spm_DEM_qU.m 7322 2018-05-31 09:47:15Z karl $
  16. % unpack
  17. %--------------------------------------------------------------------------
  18. clf
  19. V = qU.v;
  20. E = qU.z;
  21. try
  22. X = qU.x;
  23. end
  24. try
  25. C = qU.C;
  26. S = qU.S;
  27. end
  28. try
  29. pV = pU.v;
  30. pX = pU.x;
  31. end
  32. % order of hierarchy
  33. %--------------------------------------------------------------------------
  34. try
  35. g = length(X) + 1;
  36. if isempty(X{end})
  37. g = g - 1;
  38. end
  39. catch
  40. g = length(V);
  41. end
  42. % time-series specification
  43. %--------------------------------------------------------------------------
  44. N = size(V{1},2); % length of data sequence
  45. dt = 1; % time step
  46. t = (1:N)*dt; % time
  47. % unpack conditional covariances
  48. %--------------------------------------------------------------------------
  49. ci = spm_invNcdf(1 - 0.05);
  50. s = [];
  51. c = [];
  52. try
  53. for i = 1:N
  54. c = [c abs(sqrt(diag(C{i})))];
  55. s = [s abs(sqrt(diag(S{i})))];
  56. end
  57. end
  58. % loop over levels
  59. %--------------------------------------------------------------------------
  60. for i = 1:g
  61. if N == 1
  62. % hidden causes and error - single observation
  63. %------------------------------------------------------------------
  64. subplot(g,2,2*i - 1)
  65. E{i} = real(E{i});
  66. V{i} = real(V{i});
  67. t = 1:size(V{i},1);
  68. plot(t,full(E{i})',':',t,full(V{i})')
  69. box off
  70. % conditional covariances
  71. %------------------------------------------------------------------
  72. if i > 1 && size(c,1)
  73. hold on
  74. j = 1:size(V{i},1);
  75. y = ci*c(j,:);
  76. c(j,:) = [];
  77. fill([t fliplr(t)],[full(V{i} + y)' fliplr(full(V{i} - y)')],...
  78. [1 1 1]*.8,'EdgeColor',[1 1 1]*.8)
  79. plot(t,full(E{i})',':',t,full(V{i})')
  80. hold off
  81. end
  82. % title and grid
  83. %------------------------------------------------------------------
  84. title('hidden causes','FontSize',16);
  85. axis square
  86. try, set(gca,'XLim',[t(1) t(end)]), end
  87. box off
  88. % true causes
  89. %------------------------------------------------------------------
  90. if nargin > 1
  91. subplot(g,2,2*i)
  92. plot(t,full(real(pV{i}))')
  93. title('true causes','FontSize',16);
  94. axis square
  95. try, set(gca,'XLim',[t(1) t(end)]), end
  96. box off
  97. end
  98. else
  99. % hidden causes and error - time series
  100. %------------------------------------------------------------------
  101. subplot(g,2,2*i - 1)
  102. try
  103. plot(t,pV{i}','-.k')
  104. end
  105. hold on
  106. try
  107. plot(t,full(V{i})')
  108. end
  109. try
  110. plot(t,full(E{i})',':')
  111. end
  112. box off, hold off
  113. set(gca,'XLim',[t(1) t(end)])
  114. a = axis;
  115. % conditional covariances
  116. %------------------------------------------------------------------
  117. if i > 1 && size(c,1)
  118. hold on
  119. j = (1:size(V{i},1));
  120. y = ci*c(j,:);
  121. c(j,:) = [];
  122. if size(V{i},1) < size(V{i},2)
  123. fill([t fliplr(t)],[full(V{i} + y) fliplr(full(V{i} - y))],...
  124. [1 1 1]*.8,'EdgeColor',[1 1 1]*.8)
  125. else
  126. fill([t fliplr(t)]',[full(V{i} + y) fliplr(full(V{i} - y))]',...
  127. [1 1 1]*.8,'EdgeColor',[1 1 1]*.8)
  128. end
  129. try
  130. plot(t,pV{i}','-.k')
  131. end
  132. try
  133. plot(t,full(E{i}'),':')
  134. end
  135. plot(t,full(V{i})'),box off
  136. hold off
  137. end
  138. % title, action and true causes (if available)
  139. %------------------------------------------------------------------
  140. if i == 1
  141. title('prediction and error','FontSize',16);
  142. elseif length(V) < i
  143. title('no causes','FontSize',16);
  144. elseif ~size(V{i},1)
  145. title('no causes','FontSize',16);
  146. else
  147. title('hidden causes','FontSize',16);
  148. try
  149. hold on
  150. plot(t,pV{i}','-.k'),box off
  151. end
  152. hold off
  153. end
  154. xlabel('time','FontSize',14)
  155. axis square
  156. axis(a)
  157. % hidden states
  158. %------------------------------------------------------------------
  159. try
  160. subplot(g,2,2*i)
  161. try
  162. hold on
  163. plot(t,full(pX{i}'),'-.k')
  164. box off, hold off
  165. end
  166. plot(t,full(X{i}')),box off
  167. set(gca,'XLim',[t(1) t(end)])
  168. a = axis;
  169. if ~isempty(s)
  170. hold on
  171. j = 1:size(X{i},1);
  172. y = ci*s(j,:);
  173. s(j,:) = [];
  174. fill([t fliplr(t)],[full(X{i} + y) fliplr(full(X{i} - y))],...
  175. [1 1 1]*.8,'EdgeColor',[1 1 1]*.8)
  176. try
  177. plot(t,full(pX{i}'),'-.k'),box off
  178. end
  179. plot(t,full(X{i}')),box off
  180. hold off
  181. end
  182. % title and grid
  183. %--------------------------------------------------------------
  184. title('hidden states','FontSize',16)
  185. xlabel('time','FontSize',14)
  186. axis square
  187. axis(a);
  188. catch
  189. delete(gca)
  190. end
  191. end
  192. end
  193. % plot action if specified and present
  194. %--------------------------------------------------------------------------
  195. if isfield(qU,'a')
  196. if ~isempty(qU.a{end})
  197. subplot(g,2,2*g)
  198. plot(t,full(qU.a{end})');
  199. str = 'action'; hold on
  200. try
  201. plot(t,full(pU.v{2})','-.k')
  202. box off,
  203. str = 'perturbation and action';
  204. end
  205. xlabel('time','Fontsize',14); hold off
  206. title(str,'Fontsize',16)
  207. axis square
  208. set(gca,'XLim',[t(1) t(end)])
  209. box off
  210. end
  211. end
  212. drawnow