Homogeneidade_protocolo_harmonizado/main.py at master · BrunoSantos14/Homogeneidade_protocolo_harmonizado · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
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
 ########################################################################################################
########################################################################################################
########################################################################################################
#######Rotina desenvolvida conforme o protocolo harmonizado, disponível em português no link: #######
####https://www.gov.br/inmetro/pt-br/centrais-de-conteudo/publicacoes/protocolo-harmonizado-ct5.pdf####
########################################################################################################
########################################################################################################
########################################################################################################

import pandas as pd
from statsmodels.stats.contingency_tables import cochrans_q
from scipy.stats import chi2, f
import numpy as np
from statistics import stdev

###Dados e critério
# R1<-c(10.5000,9.6000,10.4000,9.5000,10.0000,9.6000,9.8000,9.8000,18.8000,10.2000,9.8000,10.2000)
# R2<-c(10.4000,9.5000,9.9000,9.9000,9.7000,10.1000,10.4000,10.2000,10.7000,10.0000,9.5000,10.0000)
# sigmap<-1.140

R1=[10.5000,9.6000,10.4000,9.5000,10.0000,9.6000,9.8000,9.8000,10.8000,10.2000,9.8000,10.2000]
R2=[10.4000,9.5000,9.9000,9.9000,9.7000,10.1000,10.4000,10.2000,10.7000,10.0000,9.5000,10.0000]
sigmap=1.140

###
###
###Teste de Cochran, para outliers nas duplicatas
###(deve ser realizado)
###
###

# m<-sum(length(R1),length(R2))/2
# colunas<-ncol(cbind(R1,R2))

# Di2<-(R1-R2)^2
# sum_Di2<-sum(Di2)
# max_Di2<-max(Di2)

# coch_cal<-max_Di2/sum_Di2
# coch_tab<-qcochran(0.95,colunas,m)

dic = {1:R1, 2:R2}
df = pd.DataFrame(dic)
m = (len(R1)+len(R2))/2

df['Di2']=(df[1] - df[2])**2
sum_Di2= df['Di2'].sum()
max_Di2=df['Di2'].max()
coch_cal = max_Di2/sum_Di2
coch_tab = cochrans_q(df[[1,2]]).pvalue

# dec_coch<-ifelse(coch_cal<=coch_tab,
#                  "não deve ser removida nenhuma duplicata",
#                  "a duplicata de maior diferença deve ser removida")

if (coch_cal<=coch_tab):
    print("Não deve ser removida nenhuma duplicata")
    ###
    ###
    ###Teste de homogeneidade, conforme o Protocolo Harmonizado
    ###(se Cochran validar, então deve ser realizado como abaixo)
    ###
    ###
    ###Cálculo da Variação Analítica

    s2an= sum_Di2/m

    ###Cálculo da Variação Amostral


    si = R1+R2
    Vs= stdev(si)
    s2sam=((Vs/2)-s2an)/2

    ###Cálculo do Valor Crítico

    # F1<-qchisq(0.95,m-1)/(m-1)
    # F2<-(qf(0.95,m-1,m,lower.tail=TRUE)-1)/2
    # C<-(F1*(0.3*sigmap)^2)+(F2*s2an)

    F1 = chi2.ppf(0.95, ((2*m) - 1))/((2*m) - 1)
    f_test = np.var(R1, ddof=1)/np.var(R2, ddof=1)
    F2 = f.pdf(f_test, df[[1]].shape[0], df[[2]].shape[0])
    C = (F1*(0.3*sigmap)**2)+(F2*s2an)

    # ###Decisão Final

    # dec<-ifelse(s2sam<=C,"aprovada","reprovada")
    dec = "aprovado" if s2sam<=C else "reprovado"
    print(f'''
a variância amostral foi de {round(s2sam,4)}
a variância do EP foi de {round((sigmap**2),4)}
o valor crítico foi de {round(C,4)}
homogeneidade foi {dec} segundo o protocolo harmonizado
    ''')

    # print(paste("a variância amostral foi de ",round(s2sam,4),sep=""))
    # print(paste("a variância do EP foi de ",round((sigmap)^2,4),sep=""))
    # print(paste("o valor crítico foi de ",round(C,4),sep=""))
    # print(paste("A homogeneidade foi ",dec," segundo o protocolo harmonizado",sep=""))

else:
    print("A duplicata de maior diferença deve ser removida")