#--------------------------# # # # Aufgabe 1 # # # #--------------------------# # 1a) # Waehlen Sie einen Wert fuer die # Anzahl der Datenpunkte: # n = 10 # n = 100 n = 1000 beta0 = 3 beta1 = -0.5 sigma = 2 x = seq(from=-5,to=5,length=n) eps = rnorm(n,mean=0,sd=sigma) y = beta0 + beta1*x + eps plot(x,y) # 1b) L = prod( exp(-eps^2/(2*sigma^2))/sigma ) logL = -n*log(sigma) - sum(eps^2)/(2*sigma^2) L # for n=1000, no longer different from 0 logL # well defined for n=1000 log(L) # -Inf for n=1000 # 1c) x0 = rep(1,n) x1 = x X = cbind(x0,x1) XT = t(X) XTX = XT %*% X XTXinv = solve(XTX) betas = XTXinv %*% XT %*% y betas # betas[1] -> beta0 # betas[2] -> beta1 haty = betas[1]+betas[2]*x sigmaML = sqrt( sum( (y-haty)^2 )/n ) sigmaML # 1d) tbeta0 = seq(from=0,to=6,by=0.01) tbeta1 = seq(from=-3,to=3,by=0.01) tsigma = seq(from=0.1,to=4,by=0.1) n0 = length(tbeta0) n1 = length(tbeta1) nsig = length(tsigma) logL = array( 0 , dim=c(n0,n1,nsig) ) # takes about 75 seconds for n=1000: for(i in 1:n0) { for(j in 1:n1) { for(k in 1:nsig) { b0 = tbeta0[i] # beta0 b1 = tbeta1[j] # beta1 sig = tsigma[k] mu = b0 + b1 * x logL[i,j,k] = -n*log(sig) - sum( (y-mu)^2 )/(2*sig^2) } } } # Ueberblick ueber alle Werte: summary(logL) # nur das Maximum: max(logL) # 1e) # wo wird das Maximum angenommen? indices = which(logL == max(logL), arr.ind = TRUE) indices b0max = tbeta0[indices[1]] b1max = tbeta1[indices[2]] sigmax = tsigma[indices[3]] b0max # -> beta0_ML b1max # -> beta1_ML sigmax # -> sigma_ML # diese Zahlen sind, bis auf Aufloesungsgenauigkeit 0.01 bzw. 0.1 # identisch mit den Werten der Maximum Likelihood Schaetzer: betas[1] # = beta0_ML betas[2] # = beta1_ML sigmaML # 1f) # Schauen wir uns noch ein paar Bilder an: # a) logL in der (beta0,beta1)-Ebene, mit sigma = sigmax: logLa = logL[,,indices[3]] contour(tbeta0,tbeta1,logLa,nlevels=50) contour(tbeta0,tbeta1,logLa,nlevels=100) # b) logL in der (beta0,sigma)-Ebene, mit beta1 = b1max: logLb = logL[,indices[2],] contour(tbeta0,tsigma,logLb,nlevels=50) # not so nice.. contour(tbeta0,tsigma,logLb,nlevels=50,zlim=c(-10*n,0)) # getting better.. contour(tbeta0,tsigma,logLb,nlevels=100,zlim=c(-4*n,0)) contour(tbeta0,tsigma,logLb,nlevels=150,zlim=c(-2*n,0)) # c) logL in der (beta1,sigma)-Ebene, mit beta0 = b0max: logLc = logL[indices[1],,] contour(tbeta1,tsigma,logLc,nlevels=50) # not so nice.. contour(tbeta1,tsigma,logLc,nlevels=50,zlim=c(-10*n,0)) # getting better.. contour(tbeta1,tsigma,logLc,nlevels=100,zlim=c(-5*n,0)) contour(tbeta1,tsigma,logLc,nlevels=120,zlim=c(-3*n,0))