# Fonction de survie de la variable LN(mu,sigma)
S0=function(x,mu,sigma){
  s=1-pnorm((log(x)-mu)/sigma)
  return(s)
  }

# Log-vraissemblance du modèle mixte
logVraissemblance=function(x,mu,sigma,m,alpha){
  n=length(x)
  xOrdonne=sort(x)
  k=n-length(x[x>=m])+1
  y=((log(xOrdonne)-mu)/sigma)^2
  l1=-sum(log(x))-(k-1)*log(sigma*sqrt(2*pi))-sum(y[1:(k-1)])/2+(n-k+1)*log(alpha)+alpha*(n-k+1)*log(m)
  l2=-alpha*sum(log(xOrdonne[k:n]))+(n-k+1)*log(S0(m,mu,sigma))

  return(l1+l2)

  }
  
# Optimisation
chercheOptimum=function(x){

  res=list()

  n=length(x)
  xOrdonne=sort(x)
  kMin=trunc(.95*n)
  l=(1:(n-kMin))
  for (k in (kMin:(n-1))){
    mu=mean(log(xOrdonne)[1:(k-1)])
    sigma=sd(log(xOrdonne)[1:(k-1)])
    m=exp(mu+sigma*qnorm(k/n))
    alpha=(n-k+1)/sum(log(xOrdonne[k:n])/m)
    l[k-kMin+1]=logVraissemblance(x,mu,sigma,m,alpha)
    }
    
    kOptimal=which.max(l)+kMin-1
    mu=mean(log(xOrdonne)[1:(kOptimal-1)])
    sigma=sd(log(xOrdonne)[1:(kOptimal-1)])
    m=exp(mu+sigma*qnorm(kOptimal/n))
    alpha=(n-kOptimal+1)/sum(log(xOrdonne[kOptimal:n])/m)
    
    res[["logVraissemblance"]]=l
    res[["kOptimal"]]=kOptimal 
    res[["mu"]]=mu
    res[["sigma"]]=sigma
    res[["m"]]=m
    res[["alpha"]]=alpha
        
    return(res)
  }

# Fabrication d'un échantillon de la loi mélangée

n=5000

mu=0
sigma=1
p0=.98
m=exp(mu+sigma*pnorm(p0))
alpha=2.5
xPareto=m*(1-runif(n))^(-1/alpha)
y=rlnorm(10*n,mu,sigma)
xLN=y[y<=m][1:n]
x=c(1:n)
u=runif(n)
x[u>p0]=xPareto[u>p0]
x[u<=p0]=xLN[u<=p0]

z=chercheOptimum(x)
plot(z$logVraissemblance,type="l")