# -%-%-%-%-
# -%-%-%-%- Pagina 1
# -%-%-%-%-
#### Ogiva Normale
x=-1
pnorm(-1,mean=0,sd=1)
# [1] 0.1586553

x<-seq(-4,4,by=0.1)
y<-pnorm(x)
plot(x,y,type="l",col="blue")


#### Ogiva Logistica
x=-1
exp(x)/(1+exp(x))
# [1] 0.2689414
1/(1+exp(-x))
# [1] 0.2689414



Logit<-function(x){1/(1+exp(-(x)))}

x<-seq(-4,4,by=0.1)
y_logit<-Logit(x)
plot(x,y_logit,type="l",col="red")
points(x,y,col="blue",type="l")



# -%-%-%-%-
# -%-%-%-%- Pagina 2
# -%-%-%-%-
x<-seq(-3,3,by=1)
y<-pnorm(x)
y
y_logit<-Logit(x)
y_logit


x<-seq(-4,4,by=0.1)
F<-function(x,d){
	
	Phi_x<-pnorm(x,0,1)
	Psi_x<-Logit(d*x)
	Phi_x-Psi_x
}
y<-F(x,d=1)
plot(x,y,type="l",lwd=2,ylim=c(-0.10,0.10))
abline(h=0,lty="dashed")
abline(v=0,lty="dashed")


# -%-%-%-%-
# -%-%-%-%- Pagina 3-5
# -%-%-%-%-
# d = 2
y<-F(x,d=2)
plot(x,y,type="l",lwd=2,ylim=c(-0.10,0.10))
abline(h=0,lty="dashed")
abline(v=0,lty="dashed")

# Studio grafico della funzione Phi(x)-Psi(x) --Eq.3
# d = 1
y<-F(x,d=1)
plot(x,y,type="l",ylim=c(-0.10,0.10),lwd=2)
abline(h=0,lty="dashed")
abline(v=0,lty="dashed")
# d = 3
y<-F(x,d=3)
points(x,y,type="l",col="gray",lwd=2)
# d = 1.9
y<-F(x,d=1.9)
points(x,y,type="l",col="red")
# d = 1.8
y<-F(x,d=1.8)
points(x,y,type="l",col="orange")
# d = 1.7
y<-F(x,d=1.7)
points(x,y,type="l",col="green",lwd=2)
# d = 1.5
y<-F(x,d=1.5)
points(x,y,type="l",col="blue")


# -%-%-%-%-
# -%-%-%-%- Pagina 6
# -%-%-%-%-
d = 1.8

library(rootSolve)

d_prima<-function(x){
(1/sqrt(2*pi)*exp(-0.5*x^2) - (d*exp(-d*x))/(1+exp(-d*x))^2)
	                }
               
curve(d_prima(x), -3, 3,main = "uniroot.all")
abline(h=0)

Root<-uniroot.all(d_prima, c(-5, 5))
Root

abline(v=Root[Root>0],col="red",lty="dotted",lwd=2)
x1<-Root[Root>0][1]
x2<-Root[Root>0][2]


# -%-%-%-%-
# -%-%-%-%- Pagina 7
# -%-%-%-%-

#Eq.5
S_x1_x2_d<-function(d){
	
	Phi_x1<-pnorm(x1,0,1)
	Psi_x1<-Logit(d*x1)
	F_x1<-Phi_x1-Psi_x1
	Phi_x2<-pnorm(x2,0,1)
	Psi_x2<-Logit(d*x2)
	F_x2<-Phi_x2-Psi_x2
	F_x1+F_x2
}

Root.2<-uniroot.all(S_x1_x2_d, c(-5, 5))


d<-seq(0,3,by=0.01)
y= S_x1_x2_d(d)
plot(d,y,type="l",ylab="S_x1_x2_d")
abline(h=0,lty="dashed")



# -%-%-%-%-
# -%-%-%-%- Pagina 10
# -%-%-%-%-

#Proviamo ad “automatizzare” la procedura iterativa in R

iter=0
d= 1.8;d # valore "arbitrario" iniziale per d
Root<-uniroot.all(d_prima, c(-5, 5))
# Radici positive della derivata prima = 0 di F(x,d)=Phi(x)-Psi(dx) 
x1<-Root[Root>0][1]
x2<-Root[Root>0][2]

# Stampa info...
	cat("Iterazione=", iter,"d=",d,"x1=",x1,"x2=",x2,"\n")

# Radice positiva della funzione S(x2,x2,d)=0
Root.2<-uniroot.all(S_x1_x2_d, c(-5, 5)) 
d=Root.2 # aggiornare il valore iniziale di "d"



repeat{ 
	iter=iter+1
	Root<-uniroot.all(d_prima, c(-5, 5))
    # Radici positive della derivata prima = 0 di F(x,d)=Phi(x)-Psi(dx) 
    x1<-Root[Root>0][1]
    x2<-Root[Root>0][2]
    
	cat("Iterazione=", iter,"d=",d,"x1=",x1,"x2=",x2,"\n")

    # Radice positiva della funzione S(x2,x2,d)=0
    Root.2<-uniroot.all(S_x1_x2_d, c(-5, 5)) 
	if(abs(d-Root.2) < 0.0001) break 
	d=Root.2 # aggiornare il valore iniziale di "d"
	
	 }
	 
	 
	 
# -%-%-%-%-
# -%-%-%-%- Pagina 11/11
# -%-%-%-%-

Logit<-function(x){1/(1+exp(-(d*x)))}
 
# dalle nostre simulazioni
d=1.701816 # approx. = 1.702

x<-seq(-4,4,by=0.1)
y_logit<-Logit(x)

plot(x,y_logit,type="l",col="red")
points(x,pnorm(x),col="blue",type="l")