############## # Cviceni 12 - Multinomicke rozdeleni a testy dobre shody ######## # Testy pravdepodobnosti v mult. rozdeleni ####### # 1. Data: x=c(29, 20, 23, 28, 35, 25, 31, 33, 31, 26, 23, 24 ) # obrazek: barplot(x,names=1:12) (n=sum(x)) p.hat=x/n barplot(p.hat,names=1:12) # do vektoru c1 si ulozime c, ktere nas zajima c1=c(1,rep(0,10),-1) V=diag(p.hat)-p.hat%*%t(p.hat) (V.c1=t(c1)%*%V%*%c1) (t=sqrt(n)*(t(c1)%*%p.hat)/sqrt(V.c1)) # rozptyl jinak p.hat[1]+p.hat[12]-(p.hat[1]-p.hat[12])^2 # jak spocitame p-hodnotu? -> Sami #2. SAMI ##### # Testy dobre shody se znamym p_0 ###### #3. do p0 pravdepodobnosti z H0 -> jak vypadaji? p0=? barplot(cbind(p.hat,p0),beside=TRUE) chisq.test(x,p=p0,correct=FALSE) # vypocet rucne: E=p0*n sum((x-E)^2/E) (r=chisq.test(x,p=p0,correct=FALSE)$residuals) sum(r^2) ######## # Testy dobre shody s neznamymi parametry ###### # 4. data=c(42,46,12) # po vyreseni analyticky: (q.A=(2*data[1]+data[2])/(2*data[1]+2*data[2]+2*data[3])) # nebo pro line numericky postup f=function(x) data[1]/x^2*2*x+data[2]/(2*x*(1-x))*(2-4*x)+data[3]/(1-x)^2*(2*x-2) (uniroot(f,c(0.001,0.999))$root) # nebo pomoci funkci, ktera hledaji minimum, primo jako max verohodny odhad f2=function(x){ if((x<0)|(x>1)) l=-Inf else l=data[1]*log(x^2)+data[2]*log(2*x*(1-x))+data[3]*log((1-x)^2) return(-l) } optimize(f2,c(0.001,0.999))$minimum # nebo optim(1/2,f2)$par (p0.HW=c(q.A^2,2*q.A*(1-q.A),(1-q.A)^2)) chisq.test(data,p=p0.HW,correct=FALSE) # co musime opravit? ################# # Test shody s rozdelenim ################# # I. Zadany parametr # 5. ngoals=c(25, 44, 62, 65, 55, 30, 14, 11) (tf=dpois(0:7,lambda=4)) tf[8]=1-sum(tf[1:7]) tf (n=sum(ngoals)) # graficke porovnani, test -> sami # jaky je zaver? # 6. rovnice=function(x){ kk=0:6 n=sum(ngoals) y=sum(kk*ngoals[1:7])-n*x+ngoals[8]/(1-ppois(6,lambda=x))*x*(1- ppois(5,lambda=x)) return(y) } # odhad lambda: (lam.hat=uniroot(rovnice,c(1,10))$root) # nebo primo jako argmax verohodnosti f2=function(x){ if((x<0)) l=-Inf else l=sum(ngoals[1:7]*log(dpois(0:6,lambda=x)))+ngoals[8]*log(1-ppois(6,lambda=x)) return(-l) } optimize(f2,c(0.001,10))$minimum optim(1/2,f2)$par #odhad pravdepodobnosti (tf2=dpois(0:7,lambda=lam.hat)) tf2[8]=1-sum(tf2[1:7]) tf2 sum(tf2) chisq.test(ngoals,p=tf2,correct=F) # opet musime vystup lehce upravit, aby to bylo spravne...