Initial import

dist/Contents.m

@@ -0,0 +1,128 @@
+% PROBABILITY DISTRIBUTION FUNCTIONS (contents of the dist-folder):
+%
+% probability density functions
+%    BETA_LPDF     - Beta log-probability density function (lpdf).
+%    BETA_PDF      - Beta probability density function (pdf).
+%    DIR_LPDF      - Log probability density function of uniform Dirichlet
+%                    distribution
+%    DIR_PDF       - Probability density function of uniform Dirichlet
+%                    distribution
+%    GAM_CDF       - Cumulative of Gamma probability density function (cdf).
+%    GAM_LPDF      - Log of Gamma probability density function (lpdf).
+%    GAM_PDF       - Gamma probability density function (pdf).
+%    INVGAM_LPDF   - Inverse-Gamma log probability density function.
+%    INVGAM_PDF    - Inverse-Gamma probability density function.
+%    LAPLACE_LPDF  - Laplace log-probability density function (lpdf).
+%    LAPLACE_PDF   - Laplace probability density function (pdf).
+%    LOGN_LPDF     - Log normal log-probability density function (lpdf)
+%    LOGT_LPDF     - Log probability density function (lpdf) for log Student's T
+%    MNORM_LPDF    - Multivariate-Normal log-probability density function (lpdf).
+%    MNORM_PDF     - Multivariate-Normal log-probability density function (lpdf).
+%    KERNELP       - Kernel density estimator for one dimensional distribution.
+%    NORM_LPDF     - Normal log-probability density function (lpdf).
+%    NORM_PDF      - Normal probability density function (pdf).
+%    POISS_LPDF    - Poisson log-probability density function.
+%    POISS_PDF     - Poisson probability density function.
+%    SINVCHI2_LPDF - Scaled inverse-chi log-probability density function.
+%    SINVCHI2_PDF  - Scaled inverse-chi probability density function.
+%    T_LPDF        - Student's T log-probability density function (lpdf)
+%    T_PDF         - Student's T probability density function (pdf)
+%    NORM_P      - create Gaussian (multivariate) (hierarchical) prior
+%    INVGAM_P            - Create inverse-Gamma prior
+%
+%  Error and gradient functions:
+%    INVGAM_E    - compute an error term for a parameter with inverse
+%                  gamma distribution (single parameter).
+%    INVGAM_G    - compute a gradient term for a parameter with inverse
+%                  gamma distribution (single parameter).
+%    LAPLACE_E   - compute an error term for a parameter with Laplace
+%                  distribution (single parameter). 
+%    LAPLACE_G   - compute a gradient for a parameter with Laplace 
+%                  distribution (single parameter).
+%    MNORM_E     - compute an error term for parameters with normal
+%                  distribution (multiple parameters).
+%    MNORM_G     - compute a gradient for parameters with normal 
+%                  distribution (multible parameters)
+%    MNORM_S     - Maximum log likelihood second derivatives
+%    NORM_E      - compute an error term for a parameter with normal
+%                  distribution (single parameter). 
+%    NORM_G      - compute a gradient for a parameter with normal 
+%                  distribution (single parameter).
+%    NORM_S      - Maximum log likelihood second derivatives (single variable)
+%    T_E         - compute an error term for a parameter with Student's 
+%                  t-distribution (single parameter). 
+%    T_G         - compute a gradient for a parameter with Student's 
+%                  t-distribution (single parameter).
+%    DIR_E       - compute an error term for a parameter with Dirichlet
+%                  distribution (single parameter). 
+%    GINVGAM_E   - Compute an error term for a parameter with inverse
+%                  gamma distribution (single parameter).
+%    GINVGAM_G   - Compute a gradient term for a parameter with inverse
+%                  gamma distribution (single parameter).
+%    GP2R_E	 - Evaluate error function for Gaussian Process.
+%    GP2R_G      - Evaluate gradient of error for Gaussian Process.
+%    GNORM_E     - Compute an error term for a parameter with normal
+%                  distribution (single parameter). 
+%    GNORM_G     - Compute a gradient for a parameter with normal 
+%                  distribution (single parameter).
+%    GNORM_S     - Maximum log likelihood second derivatives.
+%    GT_E        - Compute an error term for a parameter with Student's
+%                  t-distribution (single parameter). 
+%    GT_G        - Compute a gradient for a parameter with Student's 
+%                  t-distribution (single parameter).
+%    GT_S        - Maximum log likelihood second derivatives for 
+%                  t-distribution.
+%    T_S         - Maximum log likelihood second derivatives for t-distribution
+%    T_P         - Create student t prior
+%
+%  Functions to sample from full conditional distribution
+%    COND_GINVGAM_CAT    - Sample conditional distribution from 
+%                          inverse gamma likelihood for a group and
+%                          categorical prior. 
+%    COND_GNORM_INVGAM   - Sample conditional distribution from
+%                          normal likelihood for group and
+%                          inverse gamma prior.
+%    COND_GNORM_NORM     - Sample conditional distribution from normal
+%                          likelihood for a group and normal prior.
+%    COND_GT_CAT         - Sample conditional distribution from t 
+%                          likelihood for a group and categorical prior.
+%    COND_GT_INVGAM      - Sample conditional distribution from t 
+%                          likelihood for a group and inverse gamma prior.
+%    COND_INVGAM_CAT     - Sample conditional distribution from
+%                          inverse gamma likelihood and categorical prior. 
+%    COND_INVGAM_INVGAM  - Sample conditional distribution from
+%                          inverse gamma likelihood and prior
+%    COND_LAPLACE_INVGAM - Sample conditional distribution from Laplace
+%                          likelihood and inverse gamma prior.
+%    COND_MNORM_INVWISH  - Sample conditional distribution from normal
+%                          likelihood for multiparameter group and
+%                          inverse wishard prior.
+%    COND_NORM_GINVGAM   - Sample conditional distribution from
+%                          normal likelihood and inverse gamma prior
+%                          for a group
+%    COND_NORM_INVGAM    - Sample conditional distribution from
+%                          normal likelihood and inverse gamma prior
+%    COND_T_CAT          - Sample conditional distribution from t
+%                          likelihood and categorical prior.
+%    COND_T_INVGAM       - Sample conditional distribution from t
+%                          likelihood and inverse gamma prior.
+%
+% Random number generators
+%    CATRAND       - Random matrices from categorical distribution.
+%    DIRRAND       - Uniform dirichlet random vectors
+%    EXPRAND       - Random matrices from exponential distribution.
+%    GAMRAND       - Random matrices from gamma distribution.
+%    INTRAND       - Random matrices from uniform integer distribution.
+%    INVGAMRAND    - Random matrices from inverse gamma distribution
+%    INVGAMRAND1   - Random matrices from inverse gamma distribution
+%    INVWISHRND    - Random matrices from inverse Wishart distribution.
+%    NORMLTRAND    - Random draws from a left-truncated normal
+%                    distribution, with mean = mu, variance = sigma2
+%    NORMRTRAND    - Random draws from a right-truncated normal
+%                    distribution, with mean = mu, variance = sigma2
+%    NORMTRAND     - Random draws from a normal truncated to interval
+%    NORMTZRAND    - Random draws from a normal distribution truncated by zero
+%    WISHRND       - Random matrices from Wishart distribution.
+%    SINVCHI2RAND  - Random matrices from scaled inverse-chi distribution
+%    TRAND         - Random numbers from Student's t-distribution
+%    UNIFRAND      - Generate unifrom random numberm from interval [A,B]

dist/beta_lpdf.m

@@ -0,0 +1,30 @@
+function y = beta_lpdf(x,a,b)
+%BETA_LPDF    Beta log-probability density function (lpdf).
+%
+%   Y = BETA_LPDF(X,A,B) Returns the log of the Beta pdf with
+%   parameters A and B, at the values in X.
+%
+%   The size of Y is the common size of the input arguments. A scalar input  
+%   functions as a constant matrix of the same size as the other inputs.     
+%
+%   Default value for A and B is 1.
+
+% Copyright (c) 2005 Aki Vehtari
+
+% This software is distributed under the GNU General Public 
+% License (version 2 or later); please refer to the file 
+% License.txt, included with the software, for details.
+
+if nargin < 3, 
+  a = 1;
+end
+
+if nargin < 2;
+  b = 1;
+end
+
+if nargin < 1, 
+  error('Requires at least one input argument.');
+end
+
+y= (a-1).*log(x) +(b-1).*log(1-x) -betaln(a,b);

dist/beta_pdf.m

@@ -0,0 +1,31 @@
+function y = beta_pdf(x,a,b)
+%BETA_PDF    Beta probability density function (pdf).
+%
+%   Y = BETA_PDF(X,A,B) Returns the Beta pdf with
+%   parameters A and B, at the values in X.
+%
+%   The size of Y is the common size of the input arguments. A
+%   scalar input functions as a constant matrix of the same size as
+%   the other inputs.
+%
+%   Default value for A and B is 1.
+
+% Copyright (c) 2005 Aki Vehtari
+
+% This software is distributed under the GNU General Public 
+% License (version 2 or later); please refer to the file 
+% License.txt, included with the software, for details.
+
+if nargin < 3, 
+  a = 1;
+end
+
+if nargin < 2;
+  b = 1;
+end
+
+if nargin < 1, 
+  error('Requires at least one input argument.');
+end
+
+y=exp((a-1).*log(x) +(b-1).*log(1-x) -betaln(a,b));

dist/catrand.m

@@ -0,0 +1,27 @@
+function r = catrand(p, m, n);
+%CATRAND Random matrices from categorical distribution.
+%
+%   Description
+%   R = CATRAND(P) returns a matrix of random numbers chosen   
+%   from the categorical distribution with parameter P.
+%   P is array of probabilities, which are not necessarily
+%   normalized, though they must be non-negative, and not all zero
+%   The size of R is the size of P. Alternatively,
+%   R = CATRAND(P,M,N) returns an M by N matrix.
+
+% Copyright (c) 1999-2004 Aki Vehtari
+
+% This software is distributed under the GNU General Public 
+% License (version 2 or later); please refer to the file 
+% License.txt, included with the software, for details.
+
+
+p=p(:);
+pc=cumsum(p);
+pc=pc./pc(end);
+if nargin < 2
+  m=1;n=1;
+elseif nargin <3
+  n=m;
+end
+r=binsgeq(pc,rand(m,n));

dist/cond_ginvgam_cat.m

@@ -0,0 +1,34 @@
+function r = cond_ginvgam_cat(a, a1, a2, x)
+% COND_GINVGAM_CAT     Sample conditional distribution from inverse
+%                      gamma likelihood for a group and categorical prior.
+%
+%       Description
+%       R = COND_GINVGAM_CAT(A, A1, A2, X) generates one sample
+%       from the conditional distribution of A given
+%       parameter structure X of lower level, structure A1 of
+%       same level hyper-parameters and A2 of higher level, i.e 
+%       is r~P(A|A1,A2,X). Returns one new sample R from the 
+%       distribution above.
+
+% Copyright (c) 1999-2000 Aki Vehtari
+
+% This software is distributed under the GNU General Public 
+% License (version 2 or later); please refer to the file 
+% License.txt, included with the software, for details.
+
+
+ii=a1.ii;
+m=length(ii);
+s=a1.s;
+nus=a2.nus;
+p=nus;
+r=zeros(1,m);
+for i=1:m
+  xi2=x(ii{i}).^2;
+  s2=s(i).^2;
+  for j=1:length(nus)
+    p(j)=sum(invgam_lpdf(xi2,s2,nus(j)));
+  end
+  p=exp(p-max(p));
+  r(i)=nus(catrand(p));
+end

dist/cond_gnorm_invgam.m

@@ -0,0 +1,34 @@
+function r = cond_gnorm_invgam(a, a1, a2, x)
+% COND_GNORM_INVGAM    Sample conditional distribution from
+%                      normal likelihood for group and
+%                      inverse gamma prior.
+%
+%       Description
+%       R = COND_GNORM_INVGAM(A, A1, A2, X) generates one sample
+%       from the conditional distribution of A given
+%       parameter structure X of lower level, structure A1 of
+%       same level hyper-parameters and A2 of higher level, i.e 
+%       is r~P(A|A1,A2,X). Returns one new sample R from the 
+%       distribution above
+
+% Copyright (c) 1999-2000 Aki Vehtari
+
+% This software is distributed under the GNU General Public 
+% License (version 2 or later); please refer to the file 
+% License.txt, included with the software, for details.
+
+
+ii=a1.ii;
+m=length(ii);
+r=zeros(1,m);
+if size(a2.s,2)>1
+  for i=1:m
+    n=length(ii{i});
+    r(i)=sqrt(invgamrand1((a2.nu*a2.s(i).^2+sumsqr(x(ii{i})))/(a2.nu+n), a2.nu+n));
+  end
+else
+  for i=1:m
+    n=length(ii{i});
+    r(i)=sqrt(invgamrand1((a2.nu*a2.s.^2+sumsqr(x(ii{i})))/(a2.nu+n), a2.nu+n));
+  end
+end

dist/cond_gnorm_norm.m

@@ -0,0 +1,31 @@
+function r = cond_gnorm_norm(a, a1, a2, x)
+% COND_GNORM_NORM      Sample conditional distribution from normal
+%                      likelihood for a group and normal prior.
+%
+%
+%           Description
+%           R = COND_GNORM_NORM(A, A1, A2, X) generates one sample
+%           from the conditional distribution of A given
+%           parameter structure X of lower level, structure A1 of
+%           same level hyper-parameters and A2 of higher level, i.e 
+%           is r~P(A|A1,A2,X). Returns one new sample R from the 
+%           distribution above.
+
+% Copyright (c) 1999-2000 Aki Vehtari
+
+
+% This software is distributed under the GNU General Public 
+% License (version 2 or later); please refer to the file 
+% License.txt, included with the software, for details.
+
+ii=a1.ii;
+m=length(ii);
+if size(a2.s,2)<m
+  a2.s=a2.s(ones(1,m));
+end
+r=zeros(1,m);
+for i=1:m
+  n=length(ii{i});
+  p1=n./a1.s(i).^2;p2=1./a2.s(i).^2;
+  r(i)=(p2*a2.mu+p1*mean(x(ii{i})))./(p1+p2)+randn(1)*sqrt(1./(p1+p2));
+end

dist/cond_gt_cat.m

@@ -0,0 +1,32 @@
+function r = cond_gt_cat(a, a1, a2, x)
+% COND_GT_CAT  Sample conditional distribution from T likelihood
+%              for a group and categorical prior.
+%
+%       Description
+%       R = COND_GT_CAT(A, A1, A2, X) generates one sample
+%       from the conditional distribution of A given
+%       parameter structure X of lower level, structure A1 of
+%       same level hyper-parameters and A2 of higher level, i.e
+%       is r~P(A|A1,A2,X). Returns one new sample R from the 
+%       distribution above.
+
+% Copyright (c) 1999-2000 Aki Vehtari
+
+% This software is distributed under the GNU General Public 
+% License (version 2 or later); please refer to the file 
+% License.txt, included with the software, for details.
+
+ii=a1.ii;
+m=length(ii);
+nus=a2.nus;
+p=nus;
+s=a1.s;
+r=zeros(1,m);
+for i=1:m
+  xii=x(ii{i});
+  for j=1:length(nus)
+    p(j)=sum(t_lpdf(xii,nus(j),0,s(i)));
+  end
+  p=exp(p-max(p));
+  r(i)=nus(catrand(p));
+end

dist/cond_gt_invgam.m

@@ -0,0 +1,35 @@
+function r = cond_gt_invgam(a, a1, a2, x)
+% COND_GT_INVGAM  Sample conditional distribution from t likelihood
+%                 for a group and inverse gamma prior.
+%
+%       Description
+%       R = COND_GT_INVGAM(A, A1, A2, X) generates one sample
+%       from the conditional distribution of A given
+%       parameter structure X of lower level, structure A1 of
+%       same level hyper-parameters and A2 of higher level, i.e 
+%       is r~P(A|A1,A2,X). Returns one new sample R from the 
+%       distribution above.
+
+% Copyright (c) 1999-2000 Aki Vehtari
+
+% This software is distributed under the GNU General Public 
+% License (version 2 or later); please refer to the file 
+% License.txt, included with the software, for details.
+
+ii=a1.ii;
+m=length(ii);
+if size(a2.s,2)<m
+  a2.s=a2.s(ones(1,m));
+end
+r=zeros(1,m);
+nu=a1.nu;
+s=a1.s;
+for i=1:m
+  xii=x(ii{i});
+  n=length(xii);
+  ss=zeros(1,n);
+  for j=1:n
+    ss(j)=invgamrand1((nu(i)*s(i).^2 + xii(j).^2)/(nu(i)+1),nu(i)+1);
+  end
+  r(i)=sqrt(cond_invgam_invgam1(a2.s.^2,a2.nu,nu(i),hmean(ss),n));
+end