function BF = Beta(x,y)$ n) I/ x9 ~1 X8 J6 w. ` G5 p
format long;' v6 D2 s8 M) c2 Z% V9 B! T4 b
BF = exp(gammaln(x)+gammaln(y)-gammaln(x+y));2 b$ c9 e _: G2 W4 z
( a5 A% A7 N" W5 T' _3 E; g+ n0 xfunction [x,n]=conjgrad(A,b,x0), v$ D; ^3 p3 P5 ~$ g/ D( r. F
if(nargin == 3)8 z$ V) X2 C- Z2 | [+ Z/ D& f8 e
eps = 1.0e-6; ! G9 a3 o& R9 L$ u5 t9 n) [end% _5 S6 K2 F9 C8 }7 p# Z
) ?1 O) ]7 F- x: H
r1 = b-A*x0; & t0 [8 \9 x% |p1 = r1;9 {+ t H8 S0 `2 I
d = dot(r1,r1)/dot(p1,A*p1); . @7 t( b* t! I- T$ R K; Ix = x0+d*p1;- o8 r4 S2 R3 V1 P9 L0 K
r2 = r1-d*A*p1;8 H6 \+ Y; X, H2 {
f = dot(r2,r2)/dot(r1,r1); : N7 I1 W9 G( Q1 b( dp2 = r2+f*p1;7 }0 g0 P3 A# Q8 y. W8 P
n = 1;5 l8 a& d' v5 p7 A% T& v1 \8 `
+ o4 y, g0 _' z5 V; h
for(i=1:(rank(A)-1))4 v1 I) c! ^- `& L3 ^" y, U
x0 = x;/ r4 |) l9 e) A z9 l
p1 = p2; 4 p1 R; ?' |" f: \1 }+ `& d r1 = r2; & y9 C- b2 d9 J7 J d = dot(r1,r1)/dot(p1,A*p1);, V) S# U9 f" q. s6 e# o' A6 o& h
x = x0+d*p1; " @: F6 n8 T- M- m, ^ r2 = r1-d*A*p1;+ `$ p& L+ c3 ~+ Q+ E. H
f = dot(r2,r2)/dot(r1,r1);4 q) K3 {3 W0 E, w1 v
p2 = r2+f*p1; 1 G0 t& \0 J8 G: z: L n = n + 1; 2 _3 @0 F9 j* t, r: S8 [end 8 \6 E7 a7 O% `& s r v# C % l3 m( n0 D* c5 \! }% H) V0 rd = dot(r2,r2)/dot(p2,A*p2); / Y" q! `3 g/ p2 Vx = x+d*p2;* d6 @. K7 c, m. {* [! Y* c
n = n + 1;