Je mets une copie brute des codes utilisés en classe. Désolé de ne pas avoir pris le temps de faire de remise en forme, j’essayerais de faire des billets, à l’occasion, si des questions précises sont posées.
- Cours 3
baseavocat=read.table("http://freakonometrics.free.fr/AutoBI.csv", header=TRUE,sep=",") tail(baseavocat) baseavocat=read.table("http://freakonometrics.free.fr/AutoBI.csv",header=TRUE,sep=",") avocat=(baseavocat$ATTORNEY==1) sexe=rep("H",length(avocat)); sexe[baseavocat$CLMSEX==2]="F" marital=rep("M",length(avocat)); marital[baseavocat$MARITAL==2]="C" marital[baseavocat$MARITAL==3]="V"; marital[baseavocat$MARITAL==4]="D" ceinture=(baseavocat$SEATBELT==1); age=baseavocat$CLMAGE cout=baseavocat$LOSS base=data.frame(avocat,sexe,marital,ceinture,age,cout) logit = function(formula, lien="logit", data=NULL) { glm(formula,family=binomial(link=lien),data) } exp(-0.87208+0.31670*10)/ (1+exp(-0.87208+0.31670*10)) summary(logit(avocat~cout,data=base)) summary(logit(avocat~sexe,data=base)) # femme exp(-0.10252+0.33093*0)/ (1+exp(-0.10252+0.33093*0)) # homme exp(-0.10252+0.33093*1)/ (1+exp(-0.10252+0.33093*1)) table(base$avocat,base$sexe) log(352/(390+352) / (390/(390+352)) log(352/390) log(333/265)-log(352/390) reg=logit(avocat~sexe+cout,data=base) summary(reg) newindiv=data.frame(sexe="H",cout=10) predict(reg,newdata=newindiv) exp( -1.05051+0.38152*1+0.31994*10)/ (1+exp( -1.05051+0.38152*1+0.31994*10)) predict(reg,newdata=newindiv,type="response") newindivH=data.frame(sexe="H",cout=seq(0,15)) newindivF=data.frame(sexe="F",cout=seq(0,15)) predH=predict(reg,newdata=newindivH,type="response") predF=predict(reg,newdata=newindivF,type="response") plot(seq(0,15),predH,col="blue",type="b",lwd=3) lines(seq(0,15),predF,col="red",type="b",lwd=3) reg=logit(avocat~.+cout,data=base) base=read.table("F://baselogit-tree.txt",sep=";", header=TRUE) attach(base) plot(X1,X2,col=Y+1) library(tree) arbre=tree(Y~X1+X2,data=base) reg=logit(Y~X1+X2,data=base) predict(reg,newdata=data.frame(X1=20,X2=6),type="response") plot(arbre) text(arbre) arbre=tree(avocat~log(cout)+age,data=base)
- Cours 4
sinistre=read.table("http://freakonometrics.free.fr/sinistreACT2040.txt", header=TRUE,sep=";") sinistres=sinistre[sinistre$garantie=="1RC",] sinistres=sinistres[sinistres$cout>0,] contrat=read.table("http://freakonometrics.free.fr/contractACT2040.txt", header=TRUE,sep=";") T=table(sinistres$nocontrat) T1=as.numeric(names(T)) T2=as.numeric(T) nombre1 = data.frame(nocontrat=T1,nbre=T2) I = contrat$nocontrat%in%T1 T1= contrat$nocontrat[I==FALSE] nombre2 = data.frame(nocontrat=T1,nbre=0) nombre=rbind(nombre1,nombre2) baseFREQ = merge(contrat,nombre) Y=baseFREQ$nbre>0 X1=baseFREQ$ageconducteur M=table(Y,X1) taux=M[2,]/apply(M,2,sum) plot(as.numeric(colnames(M)),taux,pch=19) X2=baseFREQ$agevehicule db=data.frame(Y,X1,X2) library(tree) T=tree(Y~X1+X2,split='gini',mincut=5000,data=db) partition.tree(T) plot(T) text(T) y=predict(T,newdata=data.frame(X1=18:100)) lines(18:100,y,type='s',col='red',lwd=5) ################################ Y=baseFREQ$nbre X1=baseFREQ$carburant X2=cut(baseFREQ$ageconducteur,c(0,26.5,60.5,100)) E=baseFREQ$exposition db=data.frame(Y,X1,X2,E) head(db) X12=paste(X1,X2) table(Y,X1) tapply(Y,X1,sum) tapply(E,X1,sum) tapply(Y,X1,sum)/tapply(E,X1,sum) tapply(Y,X2,sum) tapply(E,X2,sum) tapply(Y,X2,sum)/tapply(E,X2,sum) tapply(Y,X12,sum) tapply(E,X12,sum) tapply(Y,X12,sum)/tapply(E,X12,sum) M=matrix(tapply(Y,X12,sum)/tapply(E,X12,sum),3,2) colnames(M)=c("D","E") rownames(M)=c("J","A","S")
- Cours 5
sinistre=read.table("http://freakonometrics.free.fr/sinistreACT2040.txt", header=TRUE,sep=";") sinistres=sinistre[sinistre$garantie=="1RC",] sinistres=sinistres[sinistres$cout>0,] contrat=read.table("http://freakonometrics.free.fr/contractACT2040.txt", header=TRUE,sep=";") T=table(sinistres$nocontrat) T1=as.numeric(names(T)) T2=as.numeric(T) nombre1 = data.frame(nocontrat=T1,nbre=T2) I = contrat$nocontrat%in%T1 T1= contrat$nocontrat[I==FALSE] nombre2 = data.frame(nocontrat=T1,nbre=0) nombre=rbind(nombre1,nombre2) baseFREQ = merge(contrat,nombre) reg=glm(nbre~s(ageconducteur)+ offset(log(exposition)), family=poisson,data=baseFREQ) nouveau= data.frame(ageconducteur= seq(15,120),exposition=1) lambda=predict( reg,newdata=nouveau, type="response") plot(nouveau$age,lambda,type="l", col="purple",lwd=3) abline(h=sum(baseFREQ$nbre)/ sum(baseFREQ$exposition),col="red") moyenne=sum(baseFREQ$nbre)/ sum(baseFREQ$exposition) plot(nouveau$age,lambda/moyenne*100,type="l", col="purple",lwd=3) abline(h=100,col="red") reg2=glm(nbre~as.factor(ageconducteur)+ offset(log(exposition)), family=poisson,data=baseFREQ) nouveau= data.frame(ageconducteur= seq(18,95),exposition=1) lambda2=predict( reg2,newdata=nouveau, type="response") plot(nouveau$age,lambda2,type="p", col="purple",lwd=3,pch=19) abline(h=sum(baseFREQ$nbre)/ sum(baseFREQ$exposition),col="red") lambda=predict( reg,newdata=nouveau, type="response") lines(nouveau$age,lambda,col="blue",lwd=3) summary(reg) seuils=c(0,25,40,70,120) reg3=glm(nbre~cut(ageconducteur,seuils)+ offset(log(exposition)), family=poisson,data=baseFREQ) lambda3=predict( reg3,newdata=nouveau, type="response",se.fit = TRUE) lines(nouveau$age,lambda3$fit,col="orange",lwd=2,type="s") lines(nouveau$age,lambda3$fit+ 2*lambda3$se.fit,col="orange",type="s") lines(nouveau$age,lambda3$fit- 2*lambda3$se.fit,col="orange",type="s") barplot(tapply(baseFREQ$exposition,baseFREQ$ageconducteur,sum)) seuils=c(0,19,23,26,41,65,120) reg4=glm(nbre~cut(ageconducteur,seuils)+ offset(log(exposition)), family=poisson,data=baseFREQ) lambda4=predict( reg4,newdata=nouveau, type="response",se.fit = TRUE) lines(nouveau$age,lambda4$fit,col="green",lwd=2,type="s") lines(nouveau$age,lambda4$fit+ 2*lambda4$se.fit,col="green",type="s") lines(nouveau$age,lambda4$fit- 2*lambda4$se.fit,col="green",type="s") library(car) linearHypothesis(reg4,c( "cut(ageconducteur, seuils)(26,41] = cut(ageconducteur, seuils)(41,65]", "cut(ageconducteur, seuils)(41,65] = cut(ageconducteur, seuils)(65,120]")) jeune=26 reg5=glm(nbre~ageconducteur+I(ageconducteur<jeune)+ offset(log(exposition)), family=poisson,data=baseFREQ) lambda5=predict( reg5,newdata=nouveau, type="response",se.fit = TRUE) lines(nouveau$age,lambda5$fit,col="red",lwd=6,type="s")
- Cours 6
sinistre=read.table("http://freakonometrics.free.fr/sinistreACT2040.txt", header=TRUE,sep=";") sinistres=sinistre[sinistre$garantie=="1RC",] sinistres=sinistres[sinistres$cout>0,] contrat=read.table("http://freakonometrics.free.fr/contractACT2040.txt", header=TRUE,sep=";") T=table(sinistres$nocontrat) T1=as.numeric(names(T)) T2=as.numeric(T) nombre1 = data.frame(nocontrat=T1,nbre=T2) I = contrat$nocontrat%in%T1 T1= contrat$nocontrat[I==FALSE] nombre2 = data.frame(nocontrat=T1,nbre=0) nombre=rbind(nombre1,nombre2) baseFREQ = merge(contrat,nombre) reg2=glm(nbre~as.factor(ageconducteur)+ offset(log(exposition)), family=poisson,data=baseFREQ) nouveau= data.frame(ageconducteur= seq(18,95),exposition=1) lambda2=predict( reg2,newdata=nouveau, type="response") plot(nouveau$age,lambda2,type="b", col="purple",cex=1,pch=19) abline(h=sum(baseFREQ$nbre)/ sum(baseFREQ$exposition),col="red") regpoly=glm(nbre~poly(ageconducteur,degree=1)+ offset(log(exposition)), family=poisson,data=baseFREQ) lambdapoly=predict( regpoly,newdata=nouveau, type="response") lines(nouveau$age,lambdapoly,col="blue") regpoly=glm(nbre~poly(ageconducteur,degree=20)+ offset(log(exposition)), family=poisson,data=baseFREQ) lambdapoly=predict( regpoly,newdata=nouveau, type="response") lines(nouveau$age,lambdapoly,col="red",lwd=2) attach(baseFREQ) library(splines) B=bs(ageconducteur,degree=1,knots=c(40)) regsp=glm(nbre~bs(ageconducteur,degree=2, knots=c(30))+ offset(log(exposition)), family=poisson,data=baseFREQ) lambdasp=predict( regsp,newdata=nouveau, type="response") lines(nouveau$age,lambdasp,col="blue",lwd=2) B=bs(ageconducteur,degree=2, knots=c(30,70)) matplot(ageconducteur,B) summary(regsp) prevloc = function(x){ I=abs(baseFREQ$ageconducteur-x)<=5 regloc=glm(nbre~ageconducteur+ offset(log(exposition)), family=poisson,data=baseFREQ[I,]) prevloc = function(x){ I=dnorm(baseFREQ$ageconducteur,mean=x,sd=3) regloc=glm(nbre~ageconducteur+ offset(log(exposition)), family=poisson,data=baseFREQ, weights=I) return( predict(regloc, newdata=data.frame(exposition=1, ageconducteur=x),type="response")) } lines(20:80,Vectorize(prevloc)(20:80),lwd=3,col="red") library(mgcv) regsp=gam(nbre~s(ageconducteur)+ offset(log(exposition)), family=poisson,data=baseFREQ) regsp0=gam(nbre~ageconducteur+agevehicule+ offset(log(exposition)), family=poisson,data=baseFREQ) regsp1=gam(nbre~s(ageconducteur)+s(agevehicule)+ offset(log(exposition)), family=poisson,data=baseFREQ) regsp2=gam(nbre~s(ageconducteur,agevehicule)+ offset(log(exposition)), family=poisson,data=baseFREQ) prev0 = function(agec,agev){ predict(regsp0,newdata=data.frame( ageconducteur=agec, agevehicule=agev, exposition=1),type="response") } prev1 = function(agec,agev){ predict(regsp1,newdata=data.frame( ageconducteur=agec, agevehicule=agev, exposition=1),type="response") } prev2 = function(agec,agev){ predict(regsp2,newdata=data.frame( ageconducteur=agec, agevehicule=agev, exposition=1),type="response") } x1=seq(18,60,by=2) y1=seq(0,15) z0=outer(x1,y1,prev0) z1=outer(x1,y1,prev1) z2=outer(x1,y1,prev2) persp(x1,y1,z1,col="green", ticktype="detailed",shade=TRUE,theta=30) image(x1,y1,z1) contour(x1,y1,z1,add=TRUE) persp(x1,y1,z2,col="green", ticktype="detailed",shade=TRUE,theta=30) image(x1,y1,z2) contour(x1,y1,z2,add=TRUE)
- cours 7
y=c(1,2,4,2,6) x=1:5 db=data.frame(x,y) plot(db,pch=19,xlim=c(0,6),ylim=c(0,8)) vx=seq(.5,5.5,by=.1) model1 = glm(y~x,family=gaussian(link="identity")) y1=predict(model1,newdata=data.frame(x=vx)) lines(vx,y1,col="red") model2 = glm(y~x,family=gaussian(link="log")) y2=predict(model2,newdata=data.frame(x=vx), type="response") lines(vx,y2,col="blue") model3 = glm(y~x,family=gaussian(link="inverse")) y3=predict(model3,newdata=data.frame(x=vx), type="response") lines(vx,y3,col="purple") model4 = glm(y~x,family=poisson(link="identity")) y4=predict(model4,newdata=data.frame(x=vx), type="response") lines(vx,y4,col="blue") y1=predict(model1,newdata=data.frame(x=vx),se=TRUE) lines(vx,y1$fit,col="red") lines(vx,y1$fit+2*y1$se.fit,col="red",lty=2) lines(vx,y1$fit-2*y1$se.fit,col="red",lty=2) y4=predict(model4,newdata=data.frame(x=vx),se=TRUE) lines(vx,y4$fit,col="blue") lines(vx,y4$fit+2*y4$se.fit,col="blue",lty=2) lines(vx,y4$fit-2*y4$se.fit,col="blue",lty=2) y=c(2,2,4,2,6) x=1:5 db=data.frame(x,y) plot(db,pch=19,xlim=c(0,6),ylim=c(0,8)) model1 = glm(y~x,family=gaussian(link="identity")) model4 = glm(y~x,family=poisson(link="identity")) y1=predict(model1,newdata=data.frame(x=vx),se=TRUE) lines(vx,y1$fit,col="red") lines(vx,y1$fit+2*y1$se.fit,col="red",lty=2) lines(vx,y1$fit-2*y1$se.fit,col="red",lty=2) y4=predict(model4,newdata=data.frame(x=vx),se=TRUE) lines(vx,y4$fit,col="blue") lines(vx,y4$fit+2*y4$se.fit,col="blue",lty=2) lines(vx,y4$fit-2*y4$se.fit,col="blue",lty=2) y=c(1,2,4,2,6) x=1:5 db=data.frame(x,y) plot(db,pch=19,xlim=c(0,6),ylim=c(0,8)) model1 = glm(y~x,family=gaussian(link="identity"),data=db) model2 = glm(y~x,family=gaussian(link="log"),data=db) model5 = glm(y~x,family=poisson(link="log"),data=db) lines(vx,y1$fit,col="red") y2=predict(model2,newdata=data.frame(x=vx),se=TRUE,type="response") lines(vx,y2$fit,col="purple") lines(vx,y2$fit+2*y2$se.fit,col="purple",lty=2) lines(vx,y2$fit-2*y2$se.fit,col="purple",lty=2) y5=predict(model5,newdata=data.frame(x=vx),se=TRUE,type="response") lines(vx,y5$fit,col="blue") lines(vx,y5$fit+2*y5$se.fit,col="blue",lty=2) lines(vx,y5$fit-2*y5$se.fit,col="blue",lty=2) model0=lm(y~x,data=db) summary(model0) library(MASS) model5 = glm(y~x,family=poisson(link="log"), data=db) model6 = glm.nb(y~x, data=db)
- Cours 8
sinistre=read.table("http://freakonometrics.free.fr/sinistreACT2040.txt", header=TRUE,sep=";") sinistres=sinistre[sinistre$garantie=="1RC",] sinistres=sinistres[sinistres$cout>0,] contrat=read.table("http://freakonometrics.free.fr/contractACT2040.txt", header=TRUE,sep=";") T=table(sinistres$nocontrat) T1=as.numeric(names(T)) T2=as.numeric(T) nombre1 = data.frame(nocontrat=T1,nbre=T2) I = contrat$nocontrat%in%T1 T1= contrat$nocontrat[I==FALSE] nombre2 = data.frame(nocontrat=T1,nbre=0) nombre=rbind(nombre1,nombre2) baseFREQ = merge(contrat,nombre) reg=glm(nbre~carburant+zone+ageconducteur+ offset(log(exposition)), data=baseFREQ,family=poisson(link="log")) library(MASS) reg2=glm.nb(nbre~carburant+zone+ageconducteur+ offset(log(exposition)), data=baseFREQ) library(AER) dispersiontest(reg) regq=glm(nbre~carburant+zone+ageconducteur+ offset(log(exposition)), data=baseFREQ,family=quasipoisson(link="log")) library(pscl) regzi=zeroinfl(nbre~carburant+zone+ageconducteur+ offset(log(exposition)) | bs(ageconducteur), data=baseFREQ,dist="poisson",link="logit") summary(regzi) ages=seq(18,80) proba=predict(regzi,newdata=data.frame(ageconducteur=ages, carburant="E", exposition=1, zone="A"), type="zero") data=sinistres head(data) base=merge(data,contrat) reg=lm(log(cout)~agevehicule+carburant+zone,data=base)
- Cours 9
reg1=lm(log(cout)~agevehicule+zone,data=base) reg1b=glm(log(cout)~agevehicule+zone,data=base, family=gaussian(link="identity")) summary(reg1) s=summary(reg1)$sigma age=seq(0,25) nd=data.frame(agevehicule=age,zone="A") C.pred=exp(predict(reg1,newdata=nd)+s^2/2) plot(age,C.pred) abline(h=mean(base$cout),col="red") plot(age,predict(reg1,newdata=nd)) abline(h=mean(log(base$cout)),col="red") plot(age,C.pred,ylim=c(1600,3000)) abline(h=mean(base$cout),col="red") reg2=glm(cout~agevehicule+zone,data=base, family=Gamma(link="log")) summary(reg2) C.pred2=predict(reg2,newdata=nd,type="response") points(age,C.pred2,col="blue") plot(base$agevehicule,base$cout) plot(base$agevehicule,log(base$cout)) seuil=200000 # partie 1 Index=base$cout< seuil reg1.1=lm(log(cout)~agevehicule+zone,data=base[Index,]) reg1.2=glm(cout~agevehicule+zone,data=base[Index,], family=Gamma(link="log")) s.1=summary(reg1.1)$sigma C1.1=exp(predict(reg1.1,newdata=nd)+s.1^2/2) C1.2=predict(reg1.2,newdata=nd,type="response") plot(age,C1.1,ylim=c(1200,2000)) points(age,C1.2,col="blue") abline(h=mean(base$cout[Index]),col="red") # partie 2 C2=rep(mean(base$cout[!Index]),nrow(nd)) # partie 3 reg3=glm((cout>seuil)~agevehicule+zone,data=base, family=binomial) P3.1=predict(reg3,newdata=nd,type='response') plot(age,P3.1) P3.2=mean(base$cout>seuil) coutmoyen.1=C1.1*(1-P3.1)+C2*P3.1 coutmoyen.2=C1.1*(1-P3.2)+C2*P3.2 coutmoyen.3=C1.2*(1-P3.2)+C2*P3.2 plot(age,coutmoyen.2,ylim=c(1200,3000)) abline(h=mean(base$cout),col="red") lines(age,coutmoyen.3,col="blue") dessin=function(seuil){ # partie 1 Index=base$cout< seuil reg1.1=lm(log(cout)~agevehicule+zone,data=base[Index,]) reg1.2=glm(cout~agevehicule+zone,data=base[Index,], family=Gamma(link="log")) s.1=summary(reg1.1)$sigma C1.1=exp(predict(reg1.1,newdata=nd)+s.1^2/2) C1.2=predict(reg1.2,newdata=nd,type="response") # partie 2 C2=rep(mean(base$cout[!Index]),nrow(nd)) # partie 3 reg3=glm((cout>seuil)~agevehicule+zone,data=base, family=binomial) P3.1=predict(reg3,newdata=nd,type='response') P3.2=mean(base$cout>seuil) coutmoyen.1=C1.1*(1-P3.1)+C2*P3.1 coutmoyen.2=C1.1*(1-P3.2)+C2*P3.2 coutmoyen.3=C1.2*(1-P3.2)+C2*P3.2 plot(age,coutmoyen.2,ylim=c(1500,2300),type="l") abline(h=mean(base$cout),col="red") lines(age,coutmoyen.3,col="blue")} dessin(100000) library(nnet) s1=1120 s2=1220 I1=base$cout<s1 I2=base$cout<s2 I3=base$cout<10000 I=I1+I2+I3 base$classes=as.factor(I) table(base$classes) library(splines) regm=multinom(classes~bs(agevehicule)+zone,data=base) summary(regm) P=predict(regm,newdata=nd,type="probs") barplot(t(P),col=c("yellow","red","green","pink"))
- Cours 10
library(gdata) db=read.xls( "https://perso.univ-rennes1.fr/arthur.charpentier/SIN_1985_2000-PE.xls", sheet=1) date=db$DSUR D=as.Date(as.character(date),format="%Y%m%d") vD=seq(min(D),max(D),by=1) sD=table(D) d1=as.Date(names(sD)) d2=vD[-which(vD%in%d1)] vecteur.date=c(d1,d2) vecteur.cpte=c(as.numeric(sD),rep(0,length(d2))) base=data.frame(date=vecteur.date,cpte=vecteur.cpte) plot(vecteur.date,vecteur.cpte,type="h",xlim=as.Date(as.character( c(19850101,20111231)),format="%Y%m%d")) regf=glm(cpte~date,family=poisson,data=base) mean(base$cpte) abline(h= mean(base$cpte),col="red") plot(D,db$COUTSIN,type="h") quantile(db$COUTSIN,1-32/2387) d=15 sort(db$COUTSIN,decreasing=TRUE)[1:20] c=35 db$COUTSIN=db$COUTSIN/1e6 indm=function(x) pmin(35,pmax(x-15,0)) mean(indm(db$COUTSIN)) mean(indm(db$COUTSIN))*150 library(evir) seuil=5 fit=gpd(db$COUTSIN,seuil) E <- function(yinf,ysup,xi,beta,threshold){ as.numeric(integrate(function(x) (x-yinf)* dgpd(x,xi,mu=threshold,beta), lower=yinf,upper=ysup)$value+ (1-pgpd(ysup,xi,mu=threshold,beta))* (ysup-yinf)) } E(15,50,fit$par.ests[1],fit$par.ests[2],5) mean(db$COUTSIN>seuil) E(100,200,fit$par.ests[1],fit$par.ests[2],5)* mean(db$COUTSIN>seuil)*150
- Cours 11
source("https://perso.univ-rennes1.fr/arthur.charpentier/bases.R") n=6 L=rep(NA,n-1) EST.P=PAID for(j in 1:(n-1)){ L[j]=sum(PAID[1:(n-j),j+1])/sum(PAID[1:(n-j),j]) EST.P[(n-j+1):n,j+1]=L[j]*EST.P[(n-j+1):n,j] } C.ULTIME=EST.P[,n] C.AS.AT=diag(EST.P[,n:1]) R=C.ULTIME-C.AS.AT sum(R) PAID2=PAID PAID2[1,1]=NA L=rep(NA,n-1) EST.P=PAID2 for(j in 1:(n-1)){ L[j]=coefficients(lm(PAID2[,j+1]~0+PAID2[,j],weights=1/PAID2[,j])) EST.P[(n-j+1):n,j+1]=L[j]*EST.P[(n-j+1):n,j] } EST.P[1,1]=EST.P[1,2]/L[1] C.ULTIME=EST.P[,n] C.AS.AT=diag(EST.P[,n:1]) R=C.ULTIME-C.AS.AT sum(R) Y=log(L-1) plot(1:5,Y) db=data.frame(N=1:5,Y) reg=lm(Y~N,data=db) pY=predict(reg,newdata=data.frame(N=6:100)) prod(1+exp(pY)) R=C.ULTIME*prod(1+exp(pY))-C.AS.AT library(ChainLadder) MackChainLadder(PAID) plot(MackChainLadder(PAID)) # ---------------------- REGRESSIONS INC=PAID INC[,2:n]=PAID[,2:n]-PAID[,1:(n-1)] Y=as.vector(INC) L=rep(2001:2006,n) C=rep(1:6,each=n) df=data.frame(Y,L,C) reg1=lm(log(Y)~as.factor(L)+as.factor(C),data=df) sigma=summary(reg1)$sigma mu=predict(reg1,newdata=df) MY=matrix(exp(mu+sigma^2/2),n,n) INC.P=EST.P INC.P[,2:n]=EST.P[,2:n]-EST.P[,1:(n-1)] df$YP=exp(mu+sigma^2/2) R=sum(df$YP[is.na(df$Y)]) reg2=glm(Y~as.factor(L)+as.factor(C),data=df,family=poisson) mu=predict(reg2,newdata=df,type='response') df$YP2=mu MY=matrix(mu,n,n) sum(df$YP2[is.na(df$Y)])
- Cours 12
plot(cars) reg=lm(dist~speed,data=cars) vx=seq(4,26) yp=predict(reg,newdata=data.frame(speed=vx)) lines(vx,yp,col="red") vcov(reg) sqrt(diag(vcov(reg))) yp=predict(reg,newdata= data.frame(speed=vx),interval="confidence") yp plot(cars) lines(vx,yp[,1],lwd=2,col="red") lines(vx,yp[,2],lwd=2,col="blue") lines(vx,yp[,3],lwd=2,col="blue") yp=predict(reg,newdata= data.frame(speed=vx),interval="prediction") lines(vx,yp[,2],lwd=2,col="green") lines(vx,yp[,3],lwd=2,col="green") bruit=residuals(reg) X=cars$speed Xb=sample(X,size=50,replace=TRUE) Eb=sample(bruit,size=50,replace=TRUE) Yp=predict(reg,newdata= data.frame(speed=Xb)) Yb=Yp+Eb carsb=data.frame(speed=Xb,dist=Yb) plot(cars) points(Xb,Yb,pch=19,col="blue") abline(reg) regb=lm(dist~speed,data=carsb) abline(regb,col="blue") plot(cars) abline(reg) jeu=function(){ Xb=sample(X,size=50,replace=TRUE) Eb=sample(bruit,size=50,replace=TRUE) Yp=predict(reg,newdata= data.frame(speed=Xb)) Yb=as.vector(Yp)+as.vector(Eb) carsb=data.frame(speed=Xb,dist=Yb) regb=lm(dist~speed,data=carsb) abline(regb,col="blue")} M=matrix(NA,length(vx),1000) for(b in 1:1000){ Xb=sample(X,size=50,replace=TRUE) Eb=sample(bruit,size=50,replace=TRUE) Yp=predict(reg,newdata= data.frame(speed=Xb)) Yb=as.vector(Yp)+as.vector(Eb) carsb=data.frame(speed=Xb,dist=Yb) regb=lm(dist~speed,data=carsb) M[,b]=predict(regb,newdata= data.frame(speed=vx)) } b025=apply(M,1,function(x) quantile(x,.025)) b975=apply(M,1,function(x) quantile(x,.975)) plot(cars) abline(reg) lines(vx,b025,col="red") lines(vx,b975,col="red") yp=predict(reg,newdata= data.frame(speed=vx), interval="confidence") lines(vx,yp[,2],lwd=1,col="blue") lines(vx,yp[,3],lwd=1,col="blue") plot(ecdf(bruit)) abline(v=0,col="green") abline(h=.54,col="green") u=seq(-50,50,by=.1) lines(u,pnorm(u,mean(bruit),sd(bruit)),col="red") # ---------------- reg poisson plot(cars) reg=glm(dist~speed,data=cars,family=poisson) vx=seq(4,26) yp=predict(reg, newdata=data.frame(speed=vx), type="response") lines(vx,yp,col="red") ypic=predict(reg, newdata=data.frame(speed=vx), type="response",se.fit=TRUE) ylo=ypic$fit-1.96*ypic$se.fit yup=ypic$fit+1.96*ypic$se.fit lines(vx,ylo,col="blue") lines(vx,yup,col="blue") bruitpearson=residuals(reg,type="pearson") jeu=function(){ Xb=sample(X,size=50,replace=TRUE) Eb=sample(bruitpearson,size=50,replace=TRUE) Yp=predict(reg,newdata= data.frame(speed=Xb), type="response") Yb=Yp+sqrt(Yp)*Eb carsb=data.frame(speed=Xb,dist=Yb) regb=glm(dist~speed,data=carsb,family=poisson) ypb=predict(regb, newdata=data.frame(speed=vx), type="response") lines(vx,ypb,col="blue")} M=matrix(NA,length(vx),1000) for(b in 1:1000){ Xb=sample(X,size=50,replace=TRUE) Eb=sample(bruitpearson,size=50,replace=TRUE) Yp=predict(reg,newdata= data.frame(speed=Xb), type="response") Yb=Yp+sqrt(Yp)*Eb carsb=data.frame(speed=Xb,dist=Yb) regb=glm(dist~speed,data=carsb,family=poisson) ypb=predict(regb, newdata=data.frame(speed=vx), type="response") M[,b]=ypb } b025=apply(M,1,function(x) quantile(x,.025)) b975=apply(M,1,function(x) quantile(x,.975)) plot(cars) lines(vx,ypic$fit) lines(vx,b025,col="red") lines(vx,b975,col="red") lines(vx,ylo,col="blue") lines(vx,yup,col="blue") reg=glm(dist~speed,data=cars,family=quasipoisson) summary(reg) ylo2=ypic$fit-1.96*ypic$se.fit*sqrt(5.53) yup2=ypic$fit+1.96*ypic$se.fit*sqrt(5.53) lines(vx,ylo2,col="blue",lty=2) lines(vx,yup2,col="blue",lty=2)
- Cours 13
source("https://perso.univ-rennes1.fr/arthur.charpentier/bases.R")
Bonjour,
Je suis avec l’avis de monsier Charpentier mais R est un peu compliqué j’ai du mal à comprendre les sorties des modèles ou des fonctions !!
Je vais faire tout les exercices que vous proposez, espérant que ça va m’aider !!
Merci pour ce site il est très intéressant.
l’actuariat IARD repose sur l’utilisation de données. Il faut un logiciel pour ça. Et je pense que R est le meilleur compromis entre simplicité du langage, richesse des fonctions possibles, facilité pour trouver de l’aide et de la documentation en ligne. En plus, étant universitaire et non pas consultant pour une boite informatique, je me refuse à enseigner un logiciel payant. Or justement, R est gratuit. Et libre….
Je demandais monsieur comment avoir les bases de données Excel svp !
Je ne sais pas, car je n’utilise plus SAS depuis 12 ans ! je dirais un truc comme
proc genmod data=base;
model nbr = X1 | dist = negbin link=log
type1 type3 offset = logexposition;
run;
je ne sais pas si c’est du SAS 9.3. mais voilà en gros comment faire une régression binomiale négative. Maintenant sur le fond, j’ai un soucis à répondre aux questions sur SAS. 1) je n’y connais pas grand chose. 2) SAS coute cher, donc il serait bien qu’à ce prix là, le SAS Institute fournisse un logiciel facile d’accès, ou au moins avec une documentation claire. Si j’aide sur un logiciel aussi cher que SAS, la moindre des choses serait que SAS me rétribue pour faciliter l’utilisation d’un logiciel commercial. Je ne suis pas un représentant de SAS, comme je ne suis pas un représentant de Windows, donc je n’utilise pas Excel pour les méthodes de provisionnement. R est libre, gratuit et communautaire. Le logiciel s’améliore parce qu’il y a des utilisateurs, et tout le monde y a gagne. Passez donc à R, pourquoi continuer à faire du SAS ?
merci pour votre réponse.
je suis d’accord avec vous sur fond mais malheureusement la plupart des entreprises ont SAS et donc je voulais prendre l’habitude de travailler sur SAS.Et comme je devais par la suite effectuer une simulation de Monte Carlo, chercher les intervalles de confiance….(je sais pas faire tout cela sur R..)
Sinon sur R, c plus simple?
Je suis preneur pour tout aide sur R.
Peut-être que je m’y sentirai plus à l’aise.
merci.
Bonjour votre blog est très intéressant mais j’aurai voulu savoir comment ajuster des fréquences de sinistres par une binomiale négative sous SAS 9.3
merci