*******************************************************************;
* 2-AFC : détermination d'un seuil de sensibilité ;
* 0 si pas ressenti, 1 si ressenti, "." si pas testé ;
* http://www.oder.co.kr/01pro/pdf/ASTM%20E679-04%20Triangular%20Forced-Choice%20Threshold%20Testing%20Standard.pdf;
* l'odeur : http://www.fivesenses.com/Documents/Library/33%20%20Odor%20Basics.pdf ;
*******************************************************************;

data seuil ;
* groupe A 10 20 40 80 160 ;
* groupe B 15 30 60 120 240 ;
	array D{10} D01-D10 ;
debut_A = 15;
debut_B = 10;
facteur = 2 ;
do i=1 to 10 ;
	if int(i/2)*2=i then 	do ;
		/* pair */ if i=2 then D{i}= debut_A ; else D{i}= debut_A*2**(i/2-1) ;end ;
	else do ;
		/* impair */ if i=1 then D{i}= debut_B ; else D{i}= debut_B*2**((i-1)/2) ;end ;

end ;
m =0;
run ;
proc print data=seuil ;
run;


data in ;
	input ID test $ S01-S10 ;
	m=0;
	cards ;
1 A . 0 . 1 . 0 . 1 . 1
2 B 1 . 0 . 1 . 1 . 1 .
4 A 0 . 0 . 1 . 1 . 1 .
7 A 1 . 1 . 0 . 1 . 1 .
9 A 1 . 1 . 1 . 1 . 1 .
10 B . 0 . 1 . 0 . 1 . 1
11 B . 1 . 0 . 0 . 1 . 1
12 B . . . 0 . 1 . . . .
13 A 1 . 0 . 0 . 1 . 1 .
17 B . 0 . 1 . 1 . 1 . 1
18 A 0 . 0 1 . 0 . 1 . 1
19 A 1 . 0 . 0 . 0 . 1 .
20 A 0 . 1 . 1 . 1 . 1 .
23 A 0 . 1 . 0 . 1 . 1 .
24 B . 1 . 0 . 1 . 1 . 1
27 B . 0 . 1 . 0 . 0 . 1
;
run;

data in ; merge in seuil ; by m ;
run ;

data in ; set in ;
	keep ID test S01-S10 seuil log_seuil ;

	array S{10} S01-S10 ;
	array D{10} D01-D10 ;

Max = 0 ; imax=1;
Do i = 1 to 10;
	If S{i} NE . then do ;
		If S{i} = 0 and i > imax then imax= i; 
	end ;
end ;

If imax= 1 or imax =10	then seuil = sqrt (D{imax}*D{imax});
						else seuil = sqrt (D{imax}*D{imax+2});

log_seuil = log10(seuil) ;

proc print data=in ;
run ;
proc tabulate data=in ;
	class test ID;
	var s01-s10 seuil log_seuil;
	table ID test, (s01-s10 seuil log_seuil)*mean ;
	run ;
proc means data =in noprint ;
	var log_seuil ;
	output out= data_mean mean=m_seuil std =std_seuil ;
	run ;
data data_mean ; set data_mean ;

	seuil_mean = 10**m_seuil ;
Label
	m_seuil='Seuil moyen en Log10 = '  
	std_seuil = 'EType seuil en Log10 = ' 
	seuil_mean = 'seuil estimé = '  ;
run ;
proc tabulate data=data_mean ;
var m_seuil std_seuil seuil_mean ;
table m_seuil std_seuil seuil_mean;
run ;