I put all Giac source files that I wrote in a single repository here (formerly the graphtheory repository):
https://github.com/marohnicluka/giac
Now everything is fetchable by the single link above (you just bookmark it somewhere
![Smile :)](./images/smilies/icon_smile.gif)
I have just one more contribution to statistics tools before the 1.5.0 release: distribution fitting by maximum likelihood. The fitdistr command which does it is implemented in optimization.cc. The corresponding file fitdistr.tex (alongside with fitdistr.png) is in the doc directory. Short help:
Code : Tout sélectionner
# fitdistr
0 Lst(L),Fnc(D)
2 Returns the distribution of type D which fits most closely to the i.i.d. samples in the list L.
-1 normald
-2 poisson
-3 exponentiald
-4 geometric
-5 gammad
-6 betad
-7 cauchyd
-8 weibulld
-9 sample
-10 randvector
-11 randvar
fitdistr(randvector(1000,weibulld,1/2,1),weibull)
X:=randvar(normal,stddev=9.5):; Y:=randvar(normal,stddev=1.5):; S:=sample(eval(X/Y,0),1000):; Z:=fitdistr(S,cauchy)
X:=randvar(normal,mean=5,variance=2):; S:=sample(exp(X),1000):; fitdistr(log(S),normal)
Code : Tout sélectionner
gen find_randvars(const gen &g,gen_map &rv,GIAC_CONTEXT) {
stringstream ss;
if (g.type==_IDNT) {
if (rv.find(g)!=rv.end())
return rv[g];
ss << " var" << randvar_count;
identificateur v(ss.str().c_str());
rv[g]=v;
++randvar_count;
return v;
}
if (g.is_symb_of_sommet(at_discreted) || is_distribution(g)>0) {
ss << " tmp" << randvar_count;
identificateur t(ss.str().c_str());
ss.str("");
ss << " var" << randvar_count;
identificateur v(ss.str().c_str());
_eval(symbolic(at_sto,makesequence(g,t)),contextptr);
rv[t]=v;
++randvar_count;
return v;
}
if (g.type==_SYMB) {
gen &f=g._SYMBptr->feuille;
if (f.type==_VECT) {
vecteur F;
F.reserve(f._VECTptr->size());
for (iterateur it=f._VECTptr->begin();it!=f._VECTptr->end();++it) {
F.push_back(find_randvars(*it,rv,contextptr));
}
return symbolic(g._SYMBptr->sommet,change_subtype(F,_SEQ__VECT));
}
return symbolic(g._SYMBptr->sommet,find_randvars(f,rv,contextptr));
}
return g;
}
Code : Tout sélectionner
if (f.type==_SYMB) {
gen_map rv;
randvar_count=0;
gen e=find_randvars(f,rv,contextptr);
if (!rv.empty()) {
int nv=rv.size();
vecteur vars;
matrice R;
vars.reserve(nv);
R.reserve(nv);
for (gen_map::const_iterator it=rv.begin();it!=rv.end();++it) {
vars.push_back(it->second);
R.push_back(vranm(n,_eval(it->first,contextptr),contextptr));
}
R=mtran(R);
for (const_iterateur it=R.begin();it!=R.end();++it) {
res.push_back(_subs(makesequence(e,vars,*it),contextptr));
}
return;
}
}
Code : Tout sélectionner
int distrib_nargs(int nd){
switch (nd){
case 4: case 5: case 11: case 12: case 14:
return 1;
#if 0
case 8:
return 3;
#endif
default:
return 2;
}
}
Code : Tout sélectionner
gen _sample(const gen & args,GIAC_CONTEXT){
if (args.is_symb_of_sommet(at_discreted) || is_distribution(args)>0)
return _rand(args,contextptr);
if (args.type==_SYMB)
return _randvector(makesequence(1,args),contextptr)._VECTptr->front();
if (args.type!=_VECT || args._VECTptr->size()<2)
return gensizeerr(contextptr);
vecteur &argv=*args._VECTptr;
gen a=argv.front(),b=argv.back();
if (a==at_multinomial) {
if (argv.size()==3) {
vecteur v=argv;
v.insert(v.begin(),1);
return _randvector(v,contextptr)._VECTptr->at(0);
} if (argv.size()==4 && b.is_integer()) {
return _randvector(makesequence(b,at_multinomial,argv[1],argv[2]),contextptr);
}
return gensizeerr(contextptr);
}
if (args._VECTptr->size()!=2 || !is_integral(b) || b.type==_ZINT || b.val<0)
return gensizeerr(contextptr);
if (a.is_symb_of_sommet(at_discreted) || is_distribution(a)>0 || a.type==_SYMB)
return _randvector(makesequence(b,a),contextptr);
if (a.type!=_VECT)
return gensizeerr(contextptr);
return _rand(makesequence(b,a),contextptr);
}