days<- read.table("C:/CLIM/Days5.txt",header=T)      
attach(days)

postscript(file="C:/CLIM/Days.ps",height=6,width=16,horizontal=F)
par(mfrow=c(1,2))
 
#---------------Data--------------------------------------------
quot<- "Aachen 2004-2007, Daily Precipitation"; quot
Pr<- pmin(PrAa/10,40.0)                 #Prec limited to 40 [mm]                          
n<- length(Pr)

#---------------Logistic Regression-----------------------------
"Logistic Regression, Two Alternatives"
"Predict Prob (Prec >0)), Numerical (lagged) Variables PX,PY" 
PX<- Pr; PX[3:n]<- Pr[1:(n-2)]               #PX[1:2] artificial
PY<- Pr; PY[2:n]<- Pr[1:(n-1)]               #PY[1]   artificial
        
tst<- 365;  ts1<- tst+1;  c("N Days"=n,"t start"=tst)
 
Pr01<- pmin(Pr*10,1)       #Numerical variable, alternatives 0,1   
Q<- Pr; logli<- 0                             #Q vector of dim n
for(t in ts1:n){                  
prec.log<- glm(Pr01[1:t]~PY[1:t]+PX[1:t],family=binomial)
pred<- fitted(prec.log)                                  
Q[t]<- pred[t]                               #Forecast Approach
phat<- pmax(pmin(Q[t],0.999),0.001);  y<- Pr01[t]
logli<- logli + y*log(phat) + (1 - y)*log(1-phat)
}
c("LogLikelih_Forecast"=logli); 
"Output for t=n"
summary(prec.log)                                        
 
#--------------Histograms---------------------------------------
predt<- Q[ts1:n]; Prn<- Pr01[ts1:n]        #vectors of dim n-tst
pred0<- predt[Prn==0]; pred1<- predt[Prn==1]

#2 Histograms for predicted probabilities
xtext<- "Predict Prob(Prec >0)"
hist(pred0,xlim=c(0,1),ylim=c(0,600),xlab=xtext,ylab="Frequency")
title(sub="Days with Precip =0")
hist(pred1,xlim=c(0,1),ylim=c(0,600),xlab=xtext,ylab="Frequency")
title(sub="Days with Precip >0")

#--------------Classification-----------------------------------
med<- median(predt)                   #Median of predicted probs
"Correctly classified Cases for t=ts1:n  [Forecast Approach]"
case00<-Prn[Prn==0 & predt<=med];case11<-Prn[Prn==1 & predt>med]
corrt<- length(case00) + length(case11) #Number of correct cases
casen<- length(pred0) + length(pred1)              #casen= n-tst
c("Median"=med,"Correct00"=length(case00),"from"=length(pred0),
               "Correct11"=length(case11),"from"=length(pred1))
c("Total"=corrt,"from"=casen,"Percent"=(corrt/casen)*100)

dev.off() 

detach(days)                                                  ###
rm(list=objects())                                            ###