Newer
Older
/* miscprob.cc
Christian Beckmann & Mark Woolrich, FMRIB Image Analysis Group
Copyright (C) 1999-2000 University of Oxford */
/* CCOPYRIGHT */
// Miscellaneous maths functions that rely on libprob
#include "miscprob.h"
#include "stdlib.h"
using namespace NEWMAT;
namespace MISCMATHS {
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
// ReturnMatrix betarnd(const int dim1, const int dim2, const float a, const float b)
// {
// // Devroye, L. (1986) Non-Uniform Random Variate Generation, Springer-Verlag.
// int tdim = dim2;
// if(tdim<0){tdim=dim1;}
// Matrix g1=gammarnd(dim1, tdim, a, 1);
// Matrix g2=gammarnd(dim1, tdim, b, 1);
// Matrix res(dim1,tdim);
// for (int mc=1; mc<=res.Ncols(); mc++) {
// for (int mr=1; mr<=res.Nrows(); mr++) {
// res(mr,mc)=g1(mr,mc)/(g1(mr,mc)+g2(mr,mc));
// }
// }
// res.Release();
// return res;
// }
ReturnMatrix betapdf(const RowVector& vals, const float a, const float b)
{
RowVector res(vals);
if(a<0 || b<0)
{
throw Exception("Negative a or b in call to Miscprob::betapdf");
}
for (int mc=1; mc<=res.Ncols(); mc++)
{
float x=vals(mc);
if(x<0)
{
res(mc)=0;
}
else
{
float logkerna=(a-1)*std::log(x);
float logkernb=(b-1)*std::log(1-x);
float betaln_ab=lgam(a)+lgam(b)-lgam(a+b);
res(mc)=std::exp(logkerna+logkernb-betaln_ab);
}
}
res.Release();
return res;
}
ReturnMatrix unifrnd(const int dim1, const int dim2, const float start, const float end)
{
int tdim = dim2;
double tmpD=1.0;
if(tdim<0){tdim=dim1;}
Matrix res(dim1,tdim);
for (int mc=1; mc<=res.Ncols(); mc++) {
for (int mr=1; mr<=res.Nrows(); mr++) {
tmpD = (rand()+1)/double(RAND_MAX+2.0);
res(mr,mc)=(tmpD)*(end-start)+start;
//drand(&tmpD);
//res(mr,mc)=(tmpD-1)*(end-start)+start;
}
}
res.Release();
return res;
}
ReturnMatrix normrnd(const int dim1, const int dim2, const float mu, const float sigma)
{
int tdim = dim2;
double tmpD=1.0;
if(tdim<0){tdim=dim1;}
Matrix res(dim1,tdim);
for (int mc=1; mc<=res.Ncols(); mc++) {
for (int mr=1; mr<=res.Nrows(); mr++) {
tmpD = (rand()+1)/double(RAND_MAX+2.0);
res(mr,mc)=ndtri(tmpD)*sigma+mu ;
//drand(&tmpD);
//res(mr,mc)=ndtri(tmpD-1)*sigma+mu ;
}
}
res.Release();
return res;
}
ReturnMatrix normpdf(const RowVector& vals, const float mu, const float var)
{
RowVector res(vals);
for (int mc=1; mc<=res.Ncols(); mc++){
res(mc) = std::exp(-0.5*(std::pow(vals(mc)-mu,2)/var))*std::pow(2*M_PI*var,-0.5);
ReturnMatrix normcdf(const RowVector& vals, const float mu, const float var)
{
RowVector res(vals);
RowVector tmp;
tmp = (vals-mu)/std::sqrt(var);
for (int mc=1; mc<=res.Ncols(); mc++){
res(mc) = ndtr(tmp(mc));
}
res.Release();
return res;
}
ReturnMatrix gammacdf(const RowVector& vals, const float mu, const float var)
if((mu>0)&&(var>0)){
float b = std::pow(mu,2)/var;
float a = mu/var;
for (int mc=1; mc<=res.Ncols(); mc++){
if(vals(mc)>0)
res(mc) = gdtr(a,b,vals(mc));
}
}
res.Release();
return res;
}
ReturnMatrix gammapdf(const RowVector& vals, const float mu, const float var)
if((mu>0)&&(var>0.00001)){
float a = std::pow(mu,2)/var;
float b = mu/var;
float c = lgam(a);
if(std::abs(c) < 150){
for (int mc=1; mc<=res.Ncols(); mc++){
if(vals(mc)>0.000001){
res(mc) = std::exp(a*std::log(b) +
(a-1) * std::log(vals(mc))
- b*vals(mc) - c);
}
}
}
}
res.Release();
return res;
}
float normpdf(const float val, const float mu, const float var)
{
return std::exp(-0.5*(std::pow(val-mu,2)/var))*std::pow(2*M_PI*var,-0.5);
}
float lognormpdf(const float val, const float mu, const float var)
{
return -0.5*(std::pow(val-mu,2)/var+std::log(2*M_PI*var));
}
ReturnMatrix normpdf(const RowVector& vals, const RowVector& mu, const RowVector& var)
{
Matrix res(mu.Ncols(),vals.Ncols());
for (int mc=1; mc<=res.Ncols(); mc++){
for (int mr=1; mr<=res.Nrows(); mr++){
res(mr,mc) = std::exp(-0.5*(std::pow(vals(mc)-mu(mr),2)/var(mr)))*std::pow(2*M_PI*var(mr),-0.5);
}
}
res.Release();
return res;
}
ReturnMatrix mvnrnd(const RowVector& mu, const SymmetricMatrix& covar, int nsamp)
{
// Matrix eig_vec;
// DiagonalMatrix eig_val;
// EigenValues(covar,eig_val,eig_vec);
// Matrix ret = ones(nsamp, 1)*mu + dnormrandm(nsamp,mu.Ncols())*sqrt(eig_val)*eig_vec.t();
Mvnormrandm mvn(mu, covar);
return mvn.next(nsamp);
}
float mvnpdf(const RowVector& vals, const RowVector& mu, const SymmetricMatrix& covar)
{
if(vals.Ncols()==2)
return bvnpdf(vals,mu,covar);
else
return std::exp(-0.5*((vals-mu)*covar.i()*(vals-mu).t()).AsScalar())/(std::pow(covar.Determinant(),0.5)*std::pow(2*M_PI,vals.Ncols()/2.0));
}
float bvnpdf(const RowVector& vals, const RowVector& mu, const SymmetricMatrix& covar)
{
// bivariate normal pdf
double det=covar(1,1)*covar(2,2)-Sqr(covar(1,2));
float m1=vals(1)-mu(1);
float m2=vals(2)-mu(2);
float ss=(Sqr(m1)*covar(2,2)-2*m1*m2*covar(1,2)+Sqr(m2)*covar(1,1))/det;
return std::exp(-0.5*ss)/(std::pow(det,0.5)*std::pow(2*M_PI,vals.Ncols()/2.0));
// ReturnMatrix gammarnd(const int dim1, const int dim2,
// const float a, const float b)
// {
// // Marsaglia, G. and Tsang, W.W. (2000) "A Simple Method for Generating Gamma Variables", Acm Trans. Math. Soft. 26(3):363-372.
// int tdim = dim2;
// if(tdim<0){tdim=dim1;}
// Matrix res(dim1,tdim);
// Gam& gam=Gam::getInstance();
// gam.setParams(a,b);
// for (int mc=1; mc<=res.Ncols(); mc++) {
// for (int mr=1; mr<=res.Nrows(); mr++) {
// res(mr,mc)=gam.rnd();
// }
// }
// res.Release();
// return res;
// }
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
ReturnMatrix perms(const int n){
if(n<=1){
Matrix P(1,1);
P << n;
P.Release();
return P;
}
Matrix Q = perms(n-1); // recursive calls
int m = Q.Nrows();
Matrix P(n*m,n);
for(int i=1;i<=m;i++){
P(i,1)=n;
for(int j=1;j<=Q.Ncols();j++)
P(i,j+1)=Q(i,j);
}
for(int i=n-1;i>=1;i--){
int jj=1;
for(int j=(n-i)*m+1;j<=(n-i+1)*m;j++){
P(j,1)=i;
for(int k=1;k<=n-1;k++){
P(j,k+1)= (Q(jj,k)==i) ? n : Q(jj,k);
}
jj++;
}
}
P.Release();
return P;
}