spm_DEM_qP.m 4.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165
  1. function spm_DEM_qP(qP,pP)
  2. % reports on conditional estimates of parameters
  3. % FORMAT spm_DEM_qP(qP,pP)
  4. %
  5. % qP.P - conditional expectations
  6. % qP.V - conditional variance
  7. %
  8. % pP - optional priors
  9. %__________________________________________________________________________
  10. % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
  11. % Karl Friston
  12. % $Id: spm_DEM_qP.m 7322 2018-05-31 09:47:15Z karl $
  13. % unpack conditional covariances
  14. %--------------------------------------------------------------------------
  15. g = length(qP.P); % depth of hierarchy
  16. ci = spm_invNcdf(1 - 0.05);
  17. % loop over levels
  18. %--------------------------------------------------------------------------
  19. Label = {};
  20. col = [1 3/4 3/4];
  21. for i = 1:g
  22. % check for last level
  23. %----------------------------------------------------------------------
  24. if isempty(qP.P{i}), break, end
  25. % get lablels
  26. %----------------------------------------------------------------------
  27. label = {};
  28. if isstruct(qP.P{i})
  29. names = fieldnames(qP.P{i});
  30. for j = 1:length(names)
  31. for k = 1:length(spm_vec(getfield(qP.P{i},names{j})))
  32. label{end + 1} = names{j};
  33. end
  34. end
  35. end
  36. % conditional expectations (with priors if specified)
  37. %----------------------------------------------------------------------
  38. qi = spm_vec(qP.P{i});
  39. c = sqrt(spm_vec(qP.V{i}))*ci;
  40. j = find(c);
  41. qi = qi(j);
  42. c = c(j);
  43. try
  44. label = label(j);
  45. end
  46. try
  47. pi = spm_vec(pP.P{i});
  48. pi = pi(j);
  49. end
  50. np = length(qi);
  51. if np
  52. % use current axes if P = P{1}
  53. %------------------------------------------------------------------
  54. if g > 1, subplot(g,1,i), end
  55. % conditional means
  56. %------------------------------------------------------------------
  57. bar(qi,'Edgecolor',[1 1 1]/2,'Facecolor',[1 1 1]*.8)
  58. title(sprintf('parameters - level %i',i),'FontSize',16);
  59. axis square
  60. box off
  61. set(gca,'XLim',[0 np + 1])
  62. % conditional variances
  63. %------------------------------------------------------------------
  64. for k = 1:np
  65. line([k k], [-1 1]*c(k) + qi(k),'LineWidth',4,'Color',col);
  66. end
  67. % prior or true means
  68. %------------------------------------------------------------------
  69. try
  70. hold on, bar(1:length(qi),pi,1/3), hold off
  71. end
  72. % labels
  73. %------------------------------------------------------------------
  74. for k = 1:length(label)
  75. text(k + 1/4,qi(k),label{k},'FontWeight','Bold','Color','r');
  76. end
  77. Label = [Label, label];
  78. end
  79. end
  80. % conditional (or prior) covariance
  81. %--------------------------------------------------------------------------
  82. try
  83. if length(qP.C) == 1;
  84. return
  85. else
  86. i = find(diag(qP.C));
  87. end
  88. catch
  89. return
  90. end
  91. subplot(g,2,g + g - 1)
  92. if exist('pC','var')
  93. imagesc(spm_cov2corr(pC(i,i)))
  94. title({'prior correlations','among parameters'},'FontSize',16)
  95. else
  96. imagesc(qP.C(i,i))
  97. title({'conditional covariances','among parameters'},'FontSize',16)
  98. end
  99. if ~isempty(Label)
  100. set(gca,'YTickLabel',Label,'YTick',[1:length(Label)])
  101. end
  102. axis square
  103. % plot evolution of hyperparameters if supplied
  104. %==========================================================================
  105. subplot(g,2,g + g)
  106. try
  107. % confidence interval and expectations
  108. %----------------------------------------------------------------------
  109. ns = length(qP.p);
  110. t = 1:ns;
  111. for i = 1:ns
  112. v(:,i) = sqrt(diag(qP.U*qP.c{i}*qP.U'));
  113. end
  114. c = ci*v;
  115. p = qP.U*spm_cat(qP.p);
  116. i = find(any(v,2));
  117. c = c(i,:);
  118. p = p(i,:);
  119. % plot
  120. %----------------------------------------------------------------------
  121. hold on
  122. np = size(p,1);
  123. for i = 1:np
  124. fill([t fliplr(t)],[(p(i,:) + c(i,:)) fliplr(p(i,:) - c(i,:))],...
  125. [1 1 1]*.8,'EdgeColor',[1 1 1]/2)
  126. plot(t,p(i,:))
  127. end
  128. set(gca,'XLim',[1 ns])
  129. title({'dynamics of parameters','(minus prior)'},'FontSize',16)
  130. xlabel('time')
  131. axis square
  132. hold off
  133. catch
  134. % or correlations
  135. %----------------------------------------------------------------------
  136. imagesc(spm_cov2corr(qP.C(i,i)))
  137. title({'conditional correlations','among parameters'},'FontSize',16)
  138. axis square
  139. drawnow
  140. end