Line | |
---|
1 | function p = normpdfln(x, m, S, V) |
---|
2 | % NORMPDFLN log of multivariate normal density. |
---|
3 | % See NORMPDF for argument description. |
---|
4 | |
---|
5 | log2pi = 1.83787706640935; |
---|
6 | [d, n] = size(x); |
---|
7 | if nargin == 1 |
---|
8 | dx = x; |
---|
9 | elseif isempty(m) |
---|
10 | dx = x; |
---|
11 | else |
---|
12 | % m specified |
---|
13 | sz = size(m); |
---|
14 | if sz(1) ~= d |
---|
15 | error('rows(m) ~= rows(x)') |
---|
16 | end |
---|
17 | nm = sz(2); |
---|
18 | if nm == 1 |
---|
19 | dx = x - repmat(m,1,n); |
---|
20 | elseif n == 1 |
---|
21 | dx = repmat(x,1,nm) - m; |
---|
22 | elseif nm == n |
---|
23 | dx = x - m; |
---|
24 | else |
---|
25 | error('incompatible number of columns in x and m') |
---|
26 | end |
---|
27 | end |
---|
28 | if nargin < 3 |
---|
29 | % unit variance |
---|
30 | p = -0.5*(d*log2pi + col_sum(dx.*dx)); |
---|
31 | return |
---|
32 | end |
---|
33 | have_inv = 0; |
---|
34 | if nargin == 3 |
---|
35 | % standard deviation given |
---|
36 | if d == 1 |
---|
37 | dx = dx./S; |
---|
38 | p = (-log(S) -0.5*log2pi) - 0.5*(dx.*dx); |
---|
39 | return; |
---|
40 | end |
---|
41 | if S(2,1) ~= 0 |
---|
42 | error('S is not upper triangular') |
---|
43 | end |
---|
44 | if any(size(S) ~= [d d]) |
---|
45 | error('S is not the right size') |
---|
46 | end |
---|
47 | else |
---|
48 | if ischar(V) |
---|
49 | if strcmp(V,'inv') |
---|
50 | % inverse stddev given |
---|
51 | iS = S; |
---|
52 | have_inv = 1; |
---|
53 | else |
---|
54 | error('unknown directive') |
---|
55 | end |
---|
56 | elseif ischar(S) |
---|
57 | if strcmp(S,'inv') |
---|
58 | % inverse variance given |
---|
59 | if d == 1 |
---|
60 | iS = sqrt(V); |
---|
61 | else |
---|
62 | iS = chol(V); |
---|
63 | end |
---|
64 | have_inv = 1; |
---|
65 | else |
---|
66 | error('unknown directive') |
---|
67 | end |
---|
68 | else |
---|
69 | % variance given |
---|
70 | if d == 1 |
---|
71 | S = sqrt(V); |
---|
72 | else |
---|
73 | S = chol(V); |
---|
74 | end |
---|
75 | end |
---|
76 | end |
---|
77 | if have_inv |
---|
78 | if d == 1 |
---|
79 | dx = iS .* dx; |
---|
80 | logdetiS = log(iS); |
---|
81 | else |
---|
82 | dx = iS*dx; |
---|
83 | logdetiS = sum(log(diag(iS))); |
---|
84 | end |
---|
85 | else |
---|
86 | if d == 1 |
---|
87 | dx = dx./S; |
---|
88 | logdetiS = -log(S); |
---|
89 | else |
---|
90 | dx = solve_tril(S',dx); |
---|
91 | %dx = S'\dx; |
---|
92 | logdetiS = -sum(log(diag(S))); |
---|
93 | end |
---|
94 | end |
---|
95 | p = (logdetiS -0.5*d*log2pi) -0.5*col_sum(dx.*dx); |
---|
Note: See
TracBrowser
for help on using the repository browser.