sadII_lab3_rdyskretne_cd

2325 days ago by macieksk

# Funkcja sluzaca do znajdowania maksymalnego wyniku give.likelihood<-function(dist.fun,sample)function(thetas){ res<-data.frame(prop.thetas=thetas) res$likeli <-sapply(thetas,function(th)Reduce("*",dist.fun(sample,th))) res } give.loglikelihood<-function(dist.fun, sample)function(thetas){ res<-data.frame(prop.thetas=thetas) res$likeli <- sapply(thetas,function(th)Reduce("+",dist.fun(sample,th,log=TRUE))) res } 
       
give.likelihood<-function(dist.fun,sample)function(thetas){ 
paczki<-data.frame(ile=c(0,0,1,1,1,2)) 
       
thetas <- seq(0,20,by=0.01) pois.ll<-give.likelihood(dpois,paczki$ile)(thetas) pois.log.ll<-give.loglikelihood(dpois,paczki$ile)(thetas) 
       
pois.ll[which.max(pois.ll$likeli),] pois.log.ll[which.max(pois.log.ll$likeli),] 
       
   prop.thetas      likeli
84        0.83 0.001353861
   prop.thetas    likeli
84        0.83 -6.604795
# Sprawdzcie, czy wasze metody estymacji ML dzialaja dobrze - najpierw metoda log-likelihood # Wygenerujcie probki ze znanych rozkladow sample.pois<-rpois(100,3.1415926) sample.geo<-rgeom(100,1/17) # ...i wyestumujcie metoda ML z poprzednich zajec najlepsze parametr \theta dla tych rozkladow. # Przestrzen parametrow \theta dobieramy jak poprzednio: tak aby odpowiadały srednim z rozkladu z przedzialu [0,20] (krok 0.01, bo bedzie sie dlugo liczyc) #Czy wyestymowane parametry sa blisko podanych powyzej? #A jak jest dla probki wielkosci 500?? 
       
thetas <- seq(0,20,by=0.01) pois.ll<-give.likelihood(dpois,sample.pois)(thetas) pois.log.ll<-give.loglikelihood(dpois,sample.pois)(thetas) pois.ll[which.max(pois.ll$likeli),] pois.log.ll[which.max(pois.log.ll$likeli),] 
       
    prop.thetas       likeli
309        3.08 5.217153e-88
    prop.thetas    likeli
309        3.08 -200.9755
thetas <- seq(0,20,by=0.01) pois.ll<-give.likelihood(dpois,sample.pois)(thetas) pois.log.ll<-give.loglikelihood(dpois,sample.pois)(thetas) pois.ll[which.max(pois.ll$likeli),] pois.log.ll[which.max(pois.log.ll$likeli),] 
       
    prop.thetas       likeli
309        3.08 5.217153e-88
    prop.thetas    likeli
309        3.08 -200.9755
library(lattice) png(width=600,height=300) xyplot(likeli ~ thetas,data=pois.ll) xyplot(likeli ~ thetas,data=pois.log.ll) dev.off() 
       
null device 
          1 

#A jak jest dla probki wielkosci 500?? sample.pois<-rpois(500,3.1415926) pois.ll<-give.likelihood(dpois,sample.pois)(thetas) pois.log.ll<-give.loglikelihood(dpois,sample.pois)(thetas) pois.ll[which.max(pois.ll$likeli),] pois.log.ll[which.max(pois.log.ll$likeli),] 
       
  prop.thetas likeli
1           0      0
    prop.thetas    likeli
308        3.07 -979.1708
#Wykresy pokazuja, ze metoda z mnozeniem wykroczyla poza numeryczna poprawnosc library(lattice) png(width=600,height=300) xyplot(likeli ~ thetas,data=pois.ll) xyplot(likeli ~ thetas,data=pois.log.ll) dev.off() 
       
null device 
          1 

# Zadanie: Sprawdzmy estymacje sample.geo przy pomocy dpois i sample.pois przy pomocy dgeo 
       
# Sprawdzmy, ktory rozklad lepiej pasuje. W tym celu tworzymy data.frame zawierajacy logLik z obydwu dopasowan. both.log.ll<-rbind( cbind(give.loglikelihood(dpois, sample.pois)(thetas), theta=thetas, source="pois"), cbind(give.loglikelihood(dgeom, sample.pois)(1/thetas), theta=thetas, source="geom") ) head(both.log.ll) 
       
There were 50 or more warnings (use warnings() to see the first 50)

  prop.thetas    likeli theta source
1        0.00      -Inf  0.00   pois
2        0.01 -8239.882  0.01   pois
3        0.02 -7180.901  0.02   pois
4        0.03 -6563.512  0.03   pois
5        0.04 -6126.920  0.04   pois
6        0.05 -5789.395  0.05   pois
# Wykresy pokazuja ze wlasciwy rozklad oczywiscie lepiej pasuje. # Uwaga: Likelihood mozna porownywac, kiedy obydwa modele maja taka sama liczbe wolnych parametrow # (w tym przypadku obydwa maja 1-wymiarowe theta) # W przeciwnym przypadku nalezy uzyc korekcji, np Akaike Information Criterion: # http://en.wikipedia.org/wiki/Akaike_information_criterion # ?AIC library(lattice) png(width=600,height=600) xyplot(likeli ~ theta ,data=both.log.ll, groups=source, type="l", auto.key=TRUE) dev.off() 
       
# http://en.wikipedia.org/wiki/Akaike_information_criterion 
# Zamiast recznie przeszukiwac przestrzen rozwiazan (uzywajac funkcji which.max) # lepiej posluzyc sie wbudowana optymalizacja to.optimize<-function(theta){ -give.loglikelihood(dpois, sample.pois)(theta)$likeli } 
       
optim(2,to.optimize) # 2 to punkt startowy optymalizacji 
       
$par
[1] 3.070313

$value
[1] 979.1708

$counts
function gradient 
      24       NA 

$convergence
[1] 0

$message
NULL

Warning message:
In optim(2, to.optimize) :
  one-diml optimization by Nelder-Mead is unreliable: use optimize
optimize(f=to.optimize,interval=c(0,20)) 
       
$minimum
[1] 3.070018

$objective
[1] 979.1708
##################### 
       
# Przeczytajcie helpa do fitdist i wykonajcie to samo doswiadczenie przy uzyciu tej funkcji. # (fitdist uzywa funkcji optim) #install.packages("fitdistrplus") library(fitdistrplus) #?fitdist #http://cran.r-project.org/web/packages/fitdistrplus/index.html # fitdist(...,"pois",...) # Jak dobrze zgadzaja sie wyniki? # Uzyjcie roznych parametrow method=... w szczegolnosci metody "moment matching estimation" #(rozmawialismy o tej metodzie na poprzednich zajeciach). Porownajcie wyniki. 
       
# Narysujcie wynikowe rozklady funkcja plot(fitdist(... )) #?plot.fitdist 
       
############# 
       
# Zadanie: Zbadaj przedzialy ufnosci estymacji parametrow rozkladu funkcja fitdistr w zaleznosci od wielkosci probki; # zbadaj nastepujace wielkosci probek: 4, 16, 64, 256 ( 4^(1:4) ) # oraz (przynajmniej) nastepujace rozklady z jednym parametrem: Poisson, geometryczny, binomial (z ustalonym parametrem size=20) # W tym celu wykonaj eksperyment numeryczny. # Utworz data.frame "fit.data" z nastepujacymi kolumnami (przykładowy wiersz) # rozklad wiel.probki theta fitted.theta #1 pois 16 3 3.755 # W tym celu najpierw utworz data.frame z pierwszymi 3ema kolumnami. # Nastepnie dodaj fitted.theta przy pomocy wlasnej(ych) funkcji uzywajacej losowania z rozkladu a nastepnie funkcji fitdist # W kazdym wierszu umiesc wynik (theta) jednego dopasowywania rozkladu. Wykonaj po min. 20 dopasowan na kazdy rozklad, kazda wielkosc probki, dla dwoch wybranych theta (np 3, 10) # Wykonaj histogramy, policz przedzialy 95% ufnosci # uzywajac: # with( fit.data, tapply(fitted.theta,list(rozklad, wiel.probki),function(...)quantile... ) ) # 
       
fit.data<-data.frame(rozklad=c(rep("pois",4),rep("geom",4)), wiel.probki=rep(4^(1:4),2)) fit.data<-rbind(data.frame(fit.data,theta=3), cbind(fit.data,theta=10) ) table(fit.data) 
       
, , theta = 3

       wiel.probki
rozklad 4 16 64 256
   geom 1  1  1   1
   pois 1  1  1   1

, , theta = 10

       wiel.probki
rozklad 4 16 64 256
   geom 1  1  1   1
   pois 1  1  1   1
fit.data.1<-fit.data for(i in 1:199) fit.data<-rbind(fit.data,fit.data.1) table(fit.data) 
       
, , theta = 3

       wiel.probki
rozklad   4  16  64 256
   geom 200 200 200 200
   pois 200 200 200 200

, , theta = 10

       wiel.probki
rozklad   4  16  64 256
   geom 200 200 200 200
   pois 200 200 200 200
do.experiment<-function(rozklad, wiel.probki, theta){ if (rozklad=="pois"){ smp<-rpois(wiel.probki,theta) res<-fitdist(smp,"pois")$estimate } if (rozklad=="geom"){ smp<-rgeom(wiel.probki,1/theta) # Tutaj musialem uzyc tryCatch, gdyz dla niektorych "pechowych" danych # fitdist nie udaje sie dopasowac rozkladu res<-tryCatch(fitdist(smp,"geom")$estimate,error=function(e)NA) res<-1/res } res } do.experiment.v<-Vectorize(do.experiment) 
       
fit.data.all<-transform(fit.data, fitted.theta=do.experiment.v(rozklad,wiel.probki,theta) ) 
       
Warning messages:
1: In dgeom(x, prob, log) : NaNs produced
2: In dgeom(x, prob, log) : NaNs produced
3: In dgeom(x, prob, log) : NaNs produced
summary(fit.data.all) 
       
 rozklad      wiel.probki      theta       fitted.theta   
 geom:1600   Min.   :  4   Min.   : 3.0   Min.   : 1.250  
 pois:1600   1st Qu.: 13   1st Qu.: 3.0   1st Qu.: 2.977  
             Median : 40   Median : 6.5   Median : 4.750  
             Mean   : 85   Mean   : 6.5   Mean   : 6.524  
             3rd Qu.:112   3rd Qu.:10.0   3rd Qu.: 9.977  
             Max.   :256   Max.   :10.0   Max.   :27.000  
                                          NA's   : 3.000  
library(lattice) xyplot(theta ~ fitted.theta | rozklad ,data=fit.data.all,groups=wiel.probki, auto.key=TRUE) dev.off() 
       
null device 
          1 
library(lattice) histogram(~ theta-fitted.theta | rozklad + factor(wiel.probki) ,data=fit.data.all, breaks=10, auto.key=TRUE) dev.off() 
       
null device 
          1 
with( subset(fit.data.all,(wiel.probki==64)&(rozklad=="pois") ) , { frac.diff<-(theta-fitted.theta)/theta plot(density( frac.diff )) quantile( frac.diff ,probs=seq(0, 1, 0.25)) }) dev.off() 
       
         0%         25%         50%         75%        100% 
-0.23958334 -0.03125001  0.00156250  0.03828125  0.15104166 
null device 
          1 
# Liczymy kwantyle w podkategoriach with( fit.data.all, tapply( fitted.theta-theta, list(rozklad, wiel.probki,theta), function(diff)quantile(diff,0.05,na.rm=TRUE) ) ) with( fit.data.all, tapply( fitted.theta-theta, list(rozklad, wiel.probki,theta), function(diff)quantile(diff,0.95,na.rm=TRUE) ) ) 
       
, , 3

        4        16         64        256
geom -1.5 -0.878125 -0.4851563 -0.2501953
pois -1.5 -0.687500 -0.3437500 -0.1835937

, , 10

         4        16         64        256
geom -6.50 -3.259375 -1.8296876 -1.0134767
pois -2.75 -1.437500 -0.6570312 -0.3330078





, , 3

          4       16        64       256
geom 2.5000 1.128125 0.5789062 0.2345703
pois 1.2625 0.625000 0.3445313 0.1914063

, , 10

          4       16        64       256
geom 9.5125 3.628125 2.0328124 1.0517577
pois 2.5000 1.253125 0.6726563 0.2931641
#################### #################### 
       
# Odgadwyanie ktory rozklad jest ktory. # #install.packages("digest") library(digest) n<-50 src.vec<-sample(1:5,5) # Random labels #Check mean, variance all.data<-as.data.frame(rbind( cbind(x=rpois(n,2), source=src.vec[1]), cbind(x=rgeom(n,1/3), source=src.vec[2]), # Var: 1/2 / (1/2^2) = 2 cbind(x=rbinom(n,size=20,p=1/10), source=src.vec[3]), cbind(x=rnbinom(n,size=2,p=1/2), source=src.vec[4]), cbind(x=sample(0:4,n,replace=TRUE),source=src.vec[5]) )) src.hash<-digest(src.vec) rm(src.vec) # Nie podgladac! all.data$source<-factor(all.data$source) 
       
summary(all.data) 
       
       x         source
 Min.   :0.000   1:50  
 1st Qu.:1.000   2:50  
 Median :2.000   3:50  
 Mean   :2.048   4:50  
 3rd Qu.:3.000   5:50  
 Max.   :8.000         
many.stats<-lapply(c(M=mean,Md=median,V=var,min=min,max=max), function(f)tapply(all.data$x,all.data$source,f) ) many.stats 
       
$M
   1    2    3    4    5 
2.04 2.16 2.04 2.20 1.80 

$Md
  1   2   3   4   5 
2.0 2.0 2.0 1.0 1.5 

$V
       1        2        3        4        5 
1.957551 1.769796 1.835102 5.714286 3.183673 

$min
1 2 3 4 5 
0 0 0 0 0 

$max
1 2 3 4 5 
5 6 4 8 7 
# albo ladniej wyswietlone: do.call(rbind,many.stats) 
       
           1        2        3        4        5
M   2.040000 2.160000 2.040000 2.200000 1.800000
Md  2.000000 2.000000 2.000000 1.000000 1.500000
V   1.957551 1.769796 1.835102 5.714286 3.183673
min 0.000000 0.000000 0.000000 0.000000 0.000000
max 5.000000 6.000000 4.000000 8.000000 7.000000
library(lattice) histogram(~x | source,data=all.data, breaks=30, ,scales=list(x="free")) dev.off() 
       
null device 
          1 
# Sprawdz, czy dobrze dopasowales # 1 - rpois # 2 - rgeom # 3 - rbinom # 4 - rnbinom # 5 - runif src.hash==digest(c(1,3,4,0,5)) 
       
[1] FALSE
give.minusloglik<-function(dist.fun, sample) function(theta.vec) -sum(do.call(dist.fun,c(list(sample),as.list(theta.vec),log=TRUE)) ) loglik<-do.call(rbind,list( pois=tapply(all.data$x,all.data$source, function(data){optim(2,give.minusloglik(dpois,data))$value} ), geom=tapply(all.data$x,all.data$source, function(data){optim(1/2,give.minusloglik(dgeom,data))$value} ), binom=tapply(all.data$x,all.data$source, function(data){optim(c(20,1/10),give.minusloglik(dbinom,data))$value} ), nbinom=tapply(all.data$x,all.data$source, function(data){optim(c(2,1/2),give.minusloglik(dnbinom,data))$value} ) )) loglik 
       
There were 50 or more warnings (use warnings() to see the first 50)
              1        2        3         4        5
pois   84.48629 83.51364 84.81011 117.65714 95.35089
geom   96.28146 98.61869 96.28146  99.37382 91.24592
binom  84.46381 83.11477 84.58747 122.18545 97.32364
nbinom 84.49428 83.53680 84.83034  99.31394 90.22764
# Wprowadz korekte AIC i/lub BIC do liczonego logLike # wiki AIC, BIC