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"
15
16
17
18
19
20
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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
using namespace NEWMAT;
namespace MISCMATHS {
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 sigma)
{
RowVector res(vals);
for (int mc=1; mc<=res.Ncols(); mc++){
res(mc) = std::exp(-0.5*(std::pow(vals(mc)-mu,2)/sigma))*std::pow(2*M_PI*sigma,-0.5);
}
res.Release();
return res;
}
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 sigma)
{
RowVector res(vals);
res=0;
if((mu>0)&&(sigma>0)){
float b = std::pow(mu,2)/sigma;
float a = mu/sigma;
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 sigma)
{
RowVector res(vals);
res=0;
if((mu>0)&&(sigma>0.00001)){
float a = std::pow(mu,2)/sigma;
float b = mu/sigma;
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;
}
ReturnMatrix normpdf(const RowVector& vals, const RowVector& mu, const RowVector& sigma)
{
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)/sigma(mr)))*std::pow(2*M_PI*sigma(mr),-0.5);
}
}
res.Release();
return res;
}
ReturnMatrix mvnormrandm(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);
}
}