1:- module(plstat,[
    2    % help/0, % defined in utils.pl
    3    % help/1, % defined in utils.pl
    4    % help/2, % defined in utils.pl
    5    bug/0,
    6    list/0,
    7    suggestion/0,
    8    mean/2,
    9    median/2,
   10    mode/2,
   11    percentile/3,
   12    quartile/3,
   13    iqr/2,
   14    rms/2,
   15    sum_of_squares/2,
   16    variance/2,
   17    pop_variance/2,
   18    std_dev/2,
   19    pop_std_dev/2,
   20    range/2,
   21    midrange/2,
   22    mean_absolute_deviation/2,
   23    covariance/3,
   24    correlation/3,
   25    pearson_correlation/3,
   26    spearman_correlation/3,
   27    weighted_mean/3,
   28    harmonic_mean/2,
   29    trimmed_mean/4,
   30    trimmed_variance/4,
   31    skew/2,
   32    kurtosis/2,
   33    moment/3,
   34    sum/2, % wrapper for sum_list
   35    prod/2,
   36    rescale/2,
   37    rescale/4,
   38    mean_normalize/2,
   39    standardize/2,
   40    entropy/2,
   41    entropy/3,
   42    min_val/2, % wrapper for min_list
   43    max_val/2, % wrapper for max_list
   44    min_max_val/3,
   45    rank/2,
   46    rank/3,
   47    nth_row/3,
   48    nth_column/3,
   49    swap_rows_columns/2,
   50    split_n_parts/3,
   51    occurrences/2,
   52    occurrences/3,
   53    normalize_prob/2,
   54    delete_nth/3,
   55    sample/3,
   56    sample/4,
   57    sample/5,
   58    empirical_distribution/3,
   59    seq/4,
   60    factorial/2,
   61    choose/3,
   62    search_position_sorted/3,
   63    search_position_sorted/4
   64    ]).
   65
   66:- [utils].   67
   68eval_fraction_list(LF,LE):-
   69    maplist(eval_fraction,LF,LE).
   70eval_fraction(F,E):-
   71    E is F.
   72% multidimensional wrapper
   73multidim2(Predicate,List,Res):-
   74    ( List = [A|_], is_list(A) ->
   75        maplist(eval_fraction_list,List,L),
   76        maplist(Predicate,L,Res);
   77        maplist(eval_fraction,List,L),
   78        Call =.. [Predicate,L,Res],
   79        Call
   80    ).

   91mean([],0):- !.
   92mean([E],E):- number(E), !.
   93mean(L,Mean):-
   94    multidim2(mean_,L,Mean).
   95mean_(L,Mean):-
   96	length(L,N),
   97	sum(L,S),
   98	Mean is S/N.

  109median([],0):- !.
  110median([E],E):- number(E), !.
  111median(L,Median):-
  112    multidim2(median_,L,Median).
  113median_(L,Median):-
  114	length(L,N),
  115    (1 is N mod 2 ->
  116        N1 is N + 1,
  117        M is N1/2,
  118        nth1(M,L,Median) ;
  119    	NH is N/2, 
  120        N2 is (N + 2)/2,
  121        nth1(NH,L,XN2),
  122        nth1(N2,L,XN22),
  123        Median is (XN2 + XN22) / 2        
  124    ).

  135comp(<,[_,A1],[_,A2]) :- A1 > A2. % to have in descending order
  136comp(>, _, _).
  137mode([],0):- !.
  138mode([E],E):- number(E), !.
  139mode(L,Mode):-
  140    multidim2(mode_,L,Mode).
  141mode_(L,ModeC):-
  142	occurrences(L,Occ),
  143    predsort(comp,Occ,X),
  144    X = [[_,O]|_],
  145    findall(V,member([V,O],Occ),Mode),
  146    length(Mode,Mod),
  147    (Mod = 1 -> Mode = [ModeC] ; ModeC = Mode).

  167percentile([],_,0):- !.
  168percentile(L,K,Percentile):-
  169    ( L = [A|_], is_list(A) -> 
  170        maplist(percentile_list(K),L,Percentile);
  171        percentile_list(K,L,Percentile)
  172    ).
  173percentile_list(K,L,P):-
  174    ( is_list(K) ->
  175        maplist(percentile_single(L),K,P);
  176        percentile_single(L,K,P)
  177    ).
  178percentile_single(L,K,Percentile):-
  179    ( (K =< 0 ; K >= 100) ->
  180        writeln('k must be between 0 and 100'),
  181        false ;
  182        msort(L,LS),
  183        length(L,N),
  184        R is (K / 100)*(N - 1) + 1,
  185        RI is floor(R),
  186        RII is RI + 1,
  187        nth1(RI,LS,V1),
  188        nth1(RII,LS,V2),
  189        RF is R - RI,
  190        Percentile is V1 + RF * (V2 - V1)
  191    ).

  201quartile(L,Q,Quartile):-
  202    LQ = [1,2,3],
  203    ( member(Q,LQ) ->
  204        !,
  205        QP is Q * 25, 
  206        percentile(L,QP,Quartile) ;
  207        writeln("Insert a valid quartile (1,2, or 3)"),
  208        false
  209    ).

  219iqr([],_):- false.
  220iqr(L,IQR):-
  221    quartile(L,3,III),
  222    quartile(L,1,I),
  223    IQR is III - I.

  233square(X,XS):-
  234    pow(X,2,XS).
  235rms([],_):- writeln('The list for rms cannot be empty'), false.
  236rms(L,RMS):-
  237    multidim2(rms_,L,RMS).
  238rms_(L,RMS):-
  239    length(L,N),
  240    maplist(square,L,LS),
  241    sum(LS, SS),
  242    RMS is sqrt(SS / N).

  252diff_square(Mu,B,D):-
  253    AB is B - Mu,
  254    pow(AB,2,D).
  255sum_of_squares([],0).
  256sum_of_squares([_],0).
  257sum_of_squares(List,SumOfSquares):-
  258    multidim2(sum_of_squares_,List,SumOfSquares).
  259sum_of_squares_(L,SumSquared):-
  260    mean(L,Mean),
  261    maplist(diff_square(Mean),L,Res),
  262    sum(Res,SumSquared).

  272variance([],0):- !.
  273variance([E],0):- number(E), !.
  274variance(L,Var):-
  275    multidim2(variance_,L,Var).
  276variance_(L,Var):-
  277    length(L,N),
  278    sum_of_squares(L,SS),
  279    N1 is N - 1,
  280    Var is (1/N1) * SS.

  290pop_variance([],0):- !.
  291pop_variance([E],0):- number(E), !.
  292pop_variance(L,Var):-
  293    multidim2(pop_variance_,L,Var).
  294pop_variance_(L,Var):-
  295    length(L,N),
  296    sum_of_squares(L,SS),
  297    Var is (1/N) * SS.
  298
  299/* std_dev(+List:numbers,-StdDev:number)
  300 * StdDev is the standard deviation of the list List (square root of the sample variance)
  301 * List can also be multidimensional (list of lists)
  302 * Example: std_dev([1,2,4,6,7,8,9],S).
  303 * Expected: S = 3.039424.
  304 * */
  305std_dev(L,StdDev):-
  306    multidim2(std_dev_,L,StdDev).
  307std_dev_(L,StdDev):-
  308    variance(L,V),
  309    StdDev is sqrt(V).
  310
  311/* pop_std_dev(+List:numbers,-StdDev:number)
  312 * StdDev is the population standard deviation of the list List (square root of the population variance)
  313 * List can also be multidimensional (list of lists)
  314 * Example: pop_std_dev([1,2,4,6,7,8,9],S).
  315 * Expected: S = 3.039424.
  316 * */
  317pop_std_dev(L,StdDev):-
  318    multidim2(pop_std_dev_,L,StdDev).
  319pop_std_dev_(L,StdDev):-
  320    pop_variance(L,V),
  321    StdDev is sqrt(V).

  331range([],0):- !.
  332range([E],E):- number(E), !.
  333range(L,Range):-
  334    multidim2(range_,L,Range).
  335range_(L,Range):-
  336    min_val(L,Min),
  337    max_val(L,Max),
  338    Range is Max - Min.

  347midrange(L,Range):-
  348    multidim2(midrange_,L,Range).
  349midrange_(L,Midrange):-
  350    range(L,Range),
  351    Midrange is Range / 2.

  361diff_abs(A,B,D):-
  362    AB is A - B,
  363    D is abs(AB).
  364mean_absolute_deviation([],0):- !.
  365mean_absolute_deviation(L,MAD):-
  366    multidim2(mean_absolute_deviation_,L,MAD).
  367mean_absolute_deviation_(L,MAD):-
  368    mean(L,Mean),
  369    maplist(diff_abs(Mean),L,LD),
  370    sum(LD,S),
  371    length(L,N),
  372    MAD is (1/N) * S.

  380sub(A,B,C):-
  381    C is B - A.
  382mul(A,B,C):-
  383    C is A * B.
  384covariance(L1,L2,_):-
  385    length(L1,N),
  386    length(L2,M),
  387    M \= N,
  388    writeln('Lists must be of the same length'),
  389    write('Found '), write(N), write(' '), writeln(M),
  390    false.
  391covariance(L1,L2,Cov):-
  392    length(L1,N),
  393    mean(L1,M1),
  394    mean(L2,M2),
  395    maplist(sub(M1),L1,LS1),
  396    maplist(sub(M2),L2,LS2),
  397    maplist(mul,LS1,LS2,LP),
  398    sum(LP,S),
  399    Cov is S / (N - 1).

  410correlation(L1,L2,_):-
  411    length(L1,N),
  412    length(L2,M),
  413    M \= N,
  414    writeln('Lists must be of the same length'),
  415    write('Found '), write(N), write(' '), writeln(M),
  416    false.
  417correlation(L1,L2,Corr):-
  418    covariance(L1,L2,Cov),
  419    std_dev(L1,S1),
  420    std_dev(L2,S2),
  421    Corr is Cov / (S1 * S2).
  422pearson_correlation(L1,L2,C):-
  423    correlation(L1,L2,C).
  424spearman_correlation(L1,L2,C):-
  425    rank(L1,fractional,LF1),
  426    rank(L2,fractional,LF2),
  427    correlation(LF1,LF2,C).

  435weighted_mean(L1,L2,_):-
  436    length(L1,N),
  437    length(L2,M),
  438    M \= N,
  439    writeln('Lists must be of the same length'),
  440    write('Found '), write(N), write(' '), writeln(M),
  441    false.
  442weighted_mean(List,Weights,WM):-
  443    sum(Weights, SW),
  444    maplist(mul,List,Weights,LP),
  445    sum(LP,LS),
  446    WM is LS / SW.

  456rec(0,_):- writeln('Cannot divide by 0'), false.
  457rec(X,X1):- X1 is 1/X.
  458harmonic_mean([],_):- false.
  459harmonic_mean(L,HM):-
  460    multidim2(harmonic_mean_,L,HM).
  461harmonic_mean_(L,HM):-
  462    length(L,N),
  463    maplist(rec,L,LR),
  464    sum(LR,SLR),
  465    HM is N / SLR.

  475trimmed_mean([],_,_,0).
  476trimmed_mean(_,L,U,_):-
  477    L > U, !,
  478    writeln("The lower bound must be actually smaller than the upper bound"),
  479    false.
  480trimmed_mean(L,Lower,Upper,TM):-
  481    include(between(Lower,Upper),L,LO),
  482    mean(LO,TM).

  492trimmed_variance([],_,_,0).
  493trimmed_variance(_,L,U,_):-
  494    L > U, !,
  495    writeln("The lower bound must be actually smaller than the upper bound"),
  496    false.
  497trimmed_variance(L,Lower,Upper,TV):-
  498    include(between(Lower,Upper),L,LO),
  499    variance(LO,TV).

  508diff_and_power(Exp,Mean,Number,R):-
  509    D is Number - Mean,
  510    pow(D,Exp,R).
  511moment([],_,0).
  512moment(L,M,Moment):-
  513    length(L,N),
  514    mean(L,Mean),
  515    maplist(diff_and_power(M,Mean),L,Res),
  516    sum(Res,S),
  517    Moment is 1/N * S.

  528skew([],0):- !.
  529skew(L,Skew):-
  530    multidim2(skew_,L,Skew).
  531skew_(List,Skew):-
  532    moment(List,2,M2),
  533    moment(List,3,M3),
  534    pow(M2,3/2,Den),
  535    Skew is M3 / Den.

  545kurtosis([],0).
  546kurtosis(L,Kurtosis):-
  547    multidim2(kurtosis_,L,Kurtosis).
  548kurtosis_(List,Kurtosis):-
  549    moment(List,4,L4),
  550    moment(List,2,L2),
  551    pow(L2,2,L22),
  552    Kurtosis is L4 / L22.

  585rank(L,Rank):-
  586    rank(L,average,Rank).
  587rank(LL,Mode,Rank):-
  588    flatten(LL,L),
  589    msort(L,LS),
  590    (   Mode = dense ->    
  591    	compute_dense(LS,LS,1,LR);
  592    	(Mode = min ; Mode = competition) ->
  593    	compute_min(LS,LS,1,LR);
  594    	(Mode = max ; Mode = modified_competition) ->  
  595    	compute_max(LS,LS,1,LR);
  596    	(Mode = average ; Mode = fractional) -> 
  597     	compute_average(LS,1,LR) ;
  598    	Mode = ordinal ->  
  599    	compute_ordinal(LS,1,LR)
  600    ),
  601    extract_map(LL,LR,Rank,Mode).
  602
  603extract_map([],_LR,[],_):- !.
  604extract_map([H|T],LR,[HM|TM],Mode):-
  605    extract_map_(H,LR,HM,Mode), !,
  606	extract_map(T,LR,TM,Mode).    
  607extract_map([H|T],LR,[V|TV],M):-
  608    extract_map_([H|T],LR,[V|TV],M).
  609
  610extract_map_([],_,[],_).
  611extract_map_([H|T],LR,[V|TV],M):-
  612    memberchk([H,V],LR),
  613    ( M = ordinal -> 
  614    	delete(LR,[H,V],LR1);
  615    	LR1 = LR
  616    ),
  617    extract_map_(T,LR1,TV,M).
  618
  619compute_ordinal([],_,[]).
  620compute_ordinal([H|T],I,[[H,I]|T3]):-
  621    I1 is I+1,
  622    compute_ordinal(T,I1,T3).
  623
  624% dense
  625compute_dense([],_,_,[]).
  626compute_dense([H|_],L,IN,[[H,IN]|TR]):-
  627    findall(I,nth0(I,L,H),LI),
  628    length(LI,N),
  629    length(LAux,N),
  630    append(LAux,LN,L),
  631    I1 is IN+1,
  632    compute_dense(LN,LN,I1,TR).
  633
  634% min
  635compute_min([],_,_,[]).
  636compute_min([H|_],L,IN,[[H,IN]|TR]):-
  637    findall(I,nth0(I,L,H),LI),
  638    length(LI,N),
  639    length(LAux,N),
  640    append(LAux,LN,L),
  641    I1 is IN+N,
  642    compute_min(LN,LN,I1,TR).
  643
  644% max
  645compute_max([],_,_,[]).
  646compute_max([H|_],L,IN,[[H,IT]|TR]):-
  647    findall(I,nth0(I,L,H),LI),
  648    length(LI,N),
  649    length(LAux,N),
  650    append(LAux,LN,L),
  651    I1 is IN+N,
  652    IT is I1 - 1,
  653    compute_max(LN,LN,I1,TR).
  654
  655
  656look_in_average([],[]).
  657look_in_average([[V,Rank]|T],[[V,Rank1]|T0]):-
  658    findall(R,member([V,R],T),LR),
  659    length(LR,NR),
  660    ( NR is 0 ->  
  661    	Rank1 = Rank,
  662        T1 = T;
  663    	sum(LR,S),
  664        Rank1 is (S+Rank) / (NR+1),
  665        delete(T,[V,_],T1)
  666    ),
  667    look_in_average(T1,T0).
  668
  669% average
  670compute_average(LS,I,LAvg):-
  671    compute_ordinal(LS,I,LAvg1),
  672    look_in_average(LAvg1,LAvg).
  673
  674%%%%%%%%
  675% other utils
  676%%%%%%%%

  684nth_row(L,Nth,NthRow):-
  685    nth1(Nth,L,NthRow).

  693nth_column(L,Nth,NthColumn):-
  694    findall(E,(member(R,L),nth1(Nth,R,E)),NthColumn).

  702swap_rows_columns([H|T],MT):-
  703    transpose(H,[H|T],MT).
  704transpose([],_,[]).
  705transpose([_|T],M,[HMT|TMT]):-
  706    transposeIn(M,HMT,M1),
  707    transpose(T,M1,TMT).
  708transposeIn([],[],[]).
  709transposeIn([[V1|V2]|T],[V1|TV],[V2|TM]):-
  710    transposeIn(T,TV,TM).

  719split_n_parts(L,1,L):- !.
  720split_n_parts(_,N,_):- N =< 0, !, false.
  721split_n_parts(L,N,_):- 
  722    length(L,LN), 
  723    A is LN mod N, 
  724    A \= 0,
  725    writeln("List cannot be divided in the desired parts"), !,
  726    false.
  727split_n_parts([],_,[]):- !.
  728split_n_parts(List,Parts,[H|T]):-
  729    length(H,Parts),
  730    append(H,LT,List),
  731    split_n_parts(LT,Parts,T).

  740occurrences(L,LO):-
  741    sort(L,LS),
  742    count_occ(LS,L,LO).
  743count_occ([],_,[]).
  744count_occ([H|T],L,[[H,Occ]|T1]):-
  745    occurrences(L,H,Occ),
  746    count_occ(T,L,T1).

  754occurrences([],_,0).
  755occurrences(L,E,0):- ground(E), \+member(E,L).
  756occurrences(L,E,N):-
  757    findall(I,nth0(I,L,E),LI),
  758    length(LI,N).

  767min_val(L,M):-
  768    multidim2(min_list,L,M).

  777max_val(L,M):-
  778    multidim2(max_list,L,M).

  786min_max_val(L,Min,Max):-
  787    min_val(L,Min),
  788    max_val(L,Max).

  797sum(L,S):-
  798    multidim2(sum_list,L,S).

  807prod([],0).
  808prod(L,P):-
  809    ( L = [A|_], is_list(A) -> 
  810        maplist(prod(1),L,P);
  811        prod(1,L,P)
  812    ).
  813prod(N,[],N):- !.
  814prod(_,[0|_],0):- !.
  815prod(P0,[H|T],P1):-
  816    PT is P0 * H,
  817    prod(PT,T,P1).
  818
  819
  820%%%%%%%%%%%%%%%%%%%%%%%%

  831div(Sum,El,Div):-
  832    Sum \= 0,
  833    Div is El / Sum.
  834
  835between_float(L,U,N):-
  836    N >= L,
  837    N =< U.
  838
  839normalize_prob([],[]):- !.
  840normalize_prob(L,LNorm):-
  841    multidim2(normalize_prob_,L,LNorm).
  842normalize_prob_(L,LNorm):-
  843    maplist(between_float(0,1),L),
  844    sum(L,SL),
  845    maplist(div(SL),L,LNorm).

  860rescale([],[]).
  861rescale(L,LResc):-
  862    rescale(L,0,1,LResc).
  863rescale(L,Lower,Upper,LResc):-
  864    ( Lower >= Upper ->
  865        writeln("Upper bound should be greater than lower bound");
  866        true
  867    ),
  868    ( L = [A|_], is_list(A) -> 
  869        maplist(rescaling_(Lower,Upper),L,LResc);
  870        rescaling_(Lower,Upper,L,LResc)
  871    ).
  872rescaling_(Lower,Upper,L,LR):-
  873    min_val(L,Min),
  874    max_val(L,Max),
  875    (Min = Max ->
  876        writeln("Min = Max, cannot divide by 0"),
  877        false;
  878        maplist(rescaling(Min,Max,Lower,Upper),L,LR)
  879    ).
  880
  881rescaling(Min,Max,Lower,Upper,V,VResc):-
  882    VResc is Lower + ((V - Min)*(Upper - Lower))/(Max - Min).

  892mean_normalize([],[]).
  893mean_normalize(List,Normalized):-
  894    multidim2(mean_normalize_,List,Normalized).
  895mean_normalize_(List,Normalized):-
  896    min_val(List,Min),
  897    max_val(List,Max),
  898    mean(List,Mean),
  899    (Min = Max ->
  900        writeln("Min = Max, cannot divide by 0"),
  901        false;
  902        maplist(mean_normalizing(Min,Max,Mean),List,Normalized)
  903    ).
  904mean_normalizing(Min,Max,Mean,X,Normalized):-
  905    Normalized is (X - Mean) / (Max - Min).

  916standardize(List,Standardized):-
  917    multidim2(standardize_,List,Standardized).
  918standardize_(List,Standardized):-
  919    mean(List,Mean),
  920    pop_std_dev(List,SD),
  921    (SD = 0 -> 
  922        writeln("Cannot standardize, standard deviation in 0"),
  923        false;
  924        maplist(standardizing(Mean,SD),List,Standardized)
  925    ).
  926standardizing(Mean,SD,X,Standardized):-
  927    Standardized is (X - Mean) / SD.

  943prod_and_log(P,PR):-
  944    PR is P * log(P).
  945prod_and_log_div(P,Q,PR):-
  946    D is P/Q,
  947    PR is P * log(D).
  948entropy(L,E):-
  949    maplist(prod_and_log,L,PR),
  950    sum(PR,S),
  951    E is -S.
  952entropy(L,LP,E):-
  953    eval_fraction_list(LP,LPP),
  954    sum(LPP,S),
  955    ( S \= 1.0 ->
  956        writeln("Probabilities must sum to 1"), 
  957        false;
  958        true
  959    ),
  960    writeln(LPP),
  961    ( maplist(<(0),LPP),maplist(>(1),LPP) -> 
  962        true;
  963        writeln("Probabilities must be between 0 and 1"),
  964        false
  965    ),
  966    maplist(prod_and_log_div,L,LP,P),
  967    sum(P,E).

  976delete_nth(L,N,LDeleted):-
  977    L = [A|_],
  978    (is_list(A) ->
  979        maplist(delete_nth_(N),L,LDeleted) ;
  980        delete_nth_(N,L,LDeleted)
  981    ).
  982delete_nth_(1,[_|T],T):- !.
  983delete_nth_(Nth,[H|T],[H|T1]):-
  984    Nth > 1,
  985	N is Nth-1,
  986    delete_nth_(N,T,T1).

 1003sample_list(_,0,_,[]):- !.
 1004sample_list(L,N,Replace,[H|T]):-
 1005    length(L,Len),
 1006    random_between(1,Len,R),
 1007    nth1(R,L,H),
 1008    ( Replace = true ->
 1009        L1 = L ;
 1010        delete_nth(L,R,L1)
 1011    ),
 1012    N1 is N - 1,
 1013    sample_list(L1,N1,Replace,T).
 1014
 1015get_random_el([H],_,_,_,_,H):- !.
 1016get_random_el([H|_],P,I,I,_,H):-
 1017    P =< 0, !.
 1018get_random_el([_|T],P,I,I0,[CP|TP],El):-
 1019    P > 0,
 1020    P1 is P - CP,
 1021    I1 is I + 1,
 1022    get_random_el(T,P1,I1,I0,TP,El).
 1023
 1024sample_list_prob(_,0,_,_,[]):- !.
 1025sample_list_prob(L,Size,Replace,[HP|TP],[Out|TOut]):-
 1026    Size > 0,
 1027    random(R),
 1028    RP is R - HP,
 1029    get_random_el(L,RP,1,Index,[HP|TP],Out),
 1030    ( Replace = true ->
 1031        L1 = L,
 1032        PL = [HP|TP] ;
 1033        delete_nth(L,Index,L1),
 1034        delete_nth([HP|TP],Index,PTemp),
 1035        normalize_prob(PTemp,PL)
 1036    ),
 1037    S1 is Size - 1,
 1038    sample_list_prob(L1,S1,Replace,PL,TOut).
 1039
 1040sample(List,Size,_):-
 1041    length(List,N),
 1042    Size > N,
 1043    writeln("Sample size must be smaller or equal to the possible elements. Set Replace to true if you want to sample with replacement."), !,
 1044    false.
 1045sample(List,Size,Result):-
 1046    sample_list(List,Size,false,Result).
 1047sample(List,Size,Replace,Result):-
 1048    ( ( Replace = true ; Replace = false ) ->
 1049        sample_list(List,Size,Replace,Result) ;
 1050        writeln("Set Replace to true or false"),
 1051        false
 1052    ).
 1053sample(List,Size,Replace,Probabilities,Result):-
 1054    ( ( Replace \= true , Replace \= false ) ->
 1055        writeln("Set Replace to true or false"),
 1056        false ;
 1057        \+ sum(Probabilities,1.0) ->
 1058            writeln("Probabilities must sum to 1"),
 1059            false ;
 1060            length(List,NL), length(Probabilities,NP),
 1061            NP \= NL ->
 1062            writeln("The list with the probabilities must be of the same length of the list with the elements to sample"),
 1063            false ;
 1064            sample_list_prob(List,Size,Replace,Probabilities,Result)
 1065    ).

 1080empirical_distribution(L,X,R):-
 1081    ( L = [A|_], is_list(A) -> 
 1082        maplist(empirical_distribution_list(X),L,R);
 1083        empirical_distribution_list(X,L,R)
 1084    ).
 1085empirical_distribution_list(X,L,R):-
 1086    ( is_list(X) ->
 1087        maplist(empirical_distribution_single(L),X,R);
 1088        empirical_distribution_single(L,X,R)
 1089    ).
 1090empirical_distribution_single(List,X,R):-
 1091    length(List,N),
 1092    ( X < 0 ->
 1093        R = 0 ;
 1094      X >= N ->
 1095        R = 1 ;
 1096        msort(List,LS),
 1097        nth0(X,LS,Num),
 1098        findall(I,nth0(I,LS,Num),LI),
 1099        max_list(LI,Max),
 1100        D is Max - X + 1,
 1101        X1 is X + D,
 1102        R is X1 / N
 1103    ).

 1113seq(A,B,L):-
 1114    seq(A,B,1,L).
 1115seq(A,A,1,[A]):- !.
 1116seq(A,A,V,[]):- V \= 1.
 1117seq(_,_,0,[]).
 1118seq(A,B,_,[]):- A > B.
 1119seq(A,B,Step,[A|T]):-
 1120	A < B,
 1121	A1 is A + Step,
 1122	seq(A1,B,Step,T).

 1130factorial(0,1).
 1131factorial(1,1):- !.
 1132factorial(N,F):-
 1133    factorial_aux(N,1,F).
 1134factorial_aux(1,F,F):- !.
 1135factorial_aux(N,V,F):-
 1136    V1 is V * N,
 1137    N1 is N-1,
 1138    factorial_aux(N1,V1,F).

 1147choose(N,K,1):- N =< K.
 1148choose(N,K,C):-
 1149    N > 0,
 1150    K > 0,
 1151    N > K,
 1152    factorial(N,NF),
 1153    factorial(K,KF),
 1154    NMinusK is N - K,
 1155    factorial(NMinusK,NMKF),
 1156    C is NF / (NMKF * KF).

 1163% random_list(N,Lower,Upper,L):-
 1164%     ( N =< 0 ->
 1165%         writeln("The number of elements must be greater than 0"),
 1166%         false;
 1167%         sample_var(uniform,Lower,Upper,N,L)
 1168%     ).
 1169
 1170% random variables
 1171% expected_value_var(Val,Occ,Exp)
 1172% variance_var(Val,Occ,Exp)
 1173
 1174% TODO: set verbose/0 with assert to test also failures without printing

 1192search_position_sorted(List,Element,Pos):-
 1193    search_position_sorted(List,Element,left,Pos).
 1194search_position_sorted(List,Element,Direction,Pos):-
 1195    L = [right,left],
 1196    (\+member(Direction,L) -> writeln("Position must be left or right"), false ; true),
 1197    ( List = [A|_], is_list(A) -> 
 1198        maplist(search_position_sorted_list(Element,Direction),List,Pos);
 1199        search_position_sorted_list(Element,Direction,List,Pos)
 1200    ).
 1201search_position_sorted_list(X,Direction,L,R):-
 1202    ( is_list(X) ->
 1203        maplist(search_position_sorted_single(L,Direction),X,R);
 1204        search_position_sorted_single(L,Direction,X,R)
 1205    ).
 1206search_position_sorted_single(L,Direction,X,R):-
 1207    msort(L,Sorted),
 1208    ( Direction = right -> 
 1209        search_pos(L,X,1,0,R) ; 
 1210        length(Sorted,N),
 1211        search_pos(L,X,-1,N,R)
 1212    ).
 1213search_pos([],_,_,P,P).
 1214search_pos([H|T],El,Inc,CP,P):-
 1215    ( El < H -> 
 1216        CP = P ;
 1217        CP1 is CP + Inc,
 1218        search_pos(T,El,Inc,CP1,P)
 1219    ).
 1220
 1221
 1222% /**
 1223%  * sample_distribution()
 1224%  * Sample from the specified distribution
 1225%  * */
 1226% sample_distribution(uniform,S):-
 1227%     sample_distribution(uniform,1,0,1,S).
 1228% sample_distribution(uniform,Lower,Upper,S):-
 1229%     sample_distribution(uniform,1,Lower,Upper,S).
 1230% sample_distribution(uniform,NSamples,Lower,Upper,S):-
 1231%     ( NSamples < 0 ->
 1232%         writeln("Number of samples must be greater than 0"),
 1233%         false ;
 1234%         length(S,NSamples),
 1235%         maplist(sample_uniform(Lower,Upper),S)
 1236%     ).
 1237% sample_uniform(Lower,Upper,S):-
 1238%     random(R),
 1239%     S is (Upper - Lower + 1) * R + Lower.