Давайте узнаем что-то новое. Всегда хорошо делать что-то новое, что-то новаторское. Давайте углубимся в логистический классификатор. Ну, не совсем логистический классификатор, а что-то подобное.

Какие различные функции ссылок мы использовали? Ссылка на логит, ссылка пробит и редко ссылка на клоглог. Можем ли мы придумать какую-либо другую функцию ссылок самостоятельно? Можем ли мы создать собственный линейный классификатор? Чтобы обдумать что-то заново, нам нужен намек, какое-то направление. Позвольте дать вам один.

Все эти вышеупомянутые функции обратной связи являются не чем иным, как CDF некоторых непрерывных распределений вероятностей.

Обратная логит-ссылка - это CDF стандартного логистического распределения. Обратная пробит-ссылка - это CDF стандартного нормального распределения. Обратной ссылкой Cloglog является CDF обобщенного распределения Gumbel для минимума.

Общим для них является то, что соответствующие случайные переменные могут принимать любое значение на всей действительной строке. Эта характеристика очень важна. Можете ли вы придумать еще какие-нибудь непрерывные распределения вероятностей, которые определены для всей реальной линии? В настоящее время я могу думать о распределении Коши и распределении Лапласа. Моя цель - создать два разных двоичных классификатора с помощью этих двух дистрибутивов.

PDF f (x) и cdf F (x) стандартного распределения Коши:

Те же величины для стандартного распределения Лапласа:

Где sgn (x) таково, что если x отрицательно, то sgn (x) = - 1, а если x положительно, то sgn (x) = + 1

Оба распределения симметричны относительно 0. Следовательно, для них обоих выполняется F (0) = 0,5 и F (x) = 1-F (-x). График обоих cdfs похож на

Выглядит очень похоже на сигмовидную функцию, не так ли? Надеюсь, теперь вы понимаете, почему эти две функции могут быть потенциальными кандидатами в качестве альтернативы классификатору logit.

Когда наш фон готов, приступим к анализу. Нам нужна одна переменная ответа Y, которая может быть либо 0, либо 1, и p переменных-предикторов X1, X2,…, Xp. Пусть общее количество наблюдений равно N вместе с предположениями:

  1. Yi независимы, и каждый Yi следует распределению Бернулли с параметром pi
  2. Предикторы p X1, X2,… Xp некоррелированы.

Такой, что-

Где β - вектор параметров размерности (p + 1) x 1, а xi - вектор i-го наблюдения порядка 1 x (p + 1). то есть xi = [1 x1i x2i x3i… xpi].

В случае распределения Коши мы бы использовали-

А в случае распределения Лапласа -

Давайте построим теорию для общего F (x) и при написании кода R рассмотрим частные случаи. (Это означает, что вы можете использовать эту теорию, чтобы понять и построить пробит-модель 😊😊)

Функция правдоподобия -

Поскольку распределения симметричны относительно 0, выполняется F (x) = 1-F (-x).

Можем ли мы его немного упростить? Предположим, z = 2y-1. Итак, когда y = 1, тогда z = 1, а когда y = 0, тогда z = -1. Это поможет более удобным способом записать вероятность.

Теперь функция правдоподобия сводится к следующему:

Точно так же логарифмическая вероятность

Чтобы обучить нашу модель, нам необходимо оценить неизвестные параметры с помощью процедуры Оценка максимального правдоподобия. MLE параметров будет

Производная первого порядка логарифма правдоподобия относительно вектора параметров равна

Где f (x) - соответствующий PDF-файл.

MLE должны удовлетворять условию D = 0. Но найти корни D = 0 непросто, поскольку оно не имеет решения в замкнутой форме. Мы можем воспользоваться различными методами оптимизации. Я использую старую технику метода Ньютона-Рафсона.

Этот метод позволяет нам найти корень уравнения f (x) = 0 с помощью итераций как

И мы останавливаемся, когда разница между выходами двух последовательных шагов становится слишком маленькой. В этом случае мы бы использовали

Давайте также найдем производную второго порядка логарифмического правдоподобия.

Мы имеем дело с несколькими предикторами, поэтому соответствующий формат матрицы

где,

  1. X - матрица предикторов порядка N x (p + 1) с первым столбцом, равным 1.
  2. P - вектор-столбец порядка N x 1 с i-м элементом равным pi.
  3. W - диагональная матрица размером N x N с i-м диагональным элементом, являющимся

4. n обозначает n-ю итерацию.

Уф !!! 😓😓 Это много математики. Теперь приступим к кодированию. Здесь я использую R. Я также включаю в свой код пробит-классификатор. Перед этим давайте держать под рукой pdf и cdf распределения Коши и распределения Лапласа.

#At first lets define the functions for creating the pi values for given predictors and parameters
#x is the matrix of parameters, param is the vector of betas, response is the response variable
#at first work with probit model
p_i_finder_probit=function(x,param,response){
  n=length(response)
  p_i=array(dim=1)               #initializing an array
  for(i in 1:nrow(x)){
    val=0                        #temporary variable
    for(j in 1:ncol(x)){
      val=val+x[i,j]*param[j]    #x[i,j]=ith value of jth predictor
    }
    val=val*response[i]          #pnorm is the cdf of normal
    p_i[i]=dnorm(val,0,1)*response[i]/pnorm(val,0,1)
  }                              #dnorm is pdf of normal
  return(p_i)                    #it will return the vector P
}
#lets define the pdf of standard Cauchy distribution
cauchy_pdf=function(x){        
  return((1/(1+x^2))/pi)
}
# similarly function to calculate the cdf of standard Cauchy distribution
cauchy_cdf=function(x){
  return(0.5+atan(x)/pi)
}
# similarly finding the P column vector for Cauchy classifier
p_i_finder_cauchy=function(x,param,response){
  n=length(response)
  p_i=array(dim=1)
  for(i in 1:nrow(x)){
    val=0
    for(j in 1:ncol(x)){
      val=val+x[i,j]*param[j]
    }
    val=val*response[i]
    p_i[i]=cauchy_pdf(val)*response[i]/cauchy_cdf(val)
  }
  return(p_i)
}
# function to calculate the pdf of Laplace distribution
laplace_pdf=function(x){
  return(exp(-abs(x))/2)
}
# function to calculate the cdf of Laplace distribution
laplace_cdf=function(x){
  return(0.5+0.5*sign(x)*(1-exp(-abs(x))))
}
# pi values under Laplace classifier
p_i_finder_laplace=function(x,param,response){
  n=length(response)
  p_i=array(dim=1)
  for(i in 1:nrow(x)){
    val=0
    for(j in 1:ncol(x)){
      val=val+x[i,j]*param[j]
    }
    val=val*response[i]
    p_i[i]=laplace_pdf(val)*response[i]/laplace_cdf(val)
  }
  return(p_i)
}
#now lets write the function for constructing the W matrix
# as input we need the pi value, the matrix of predictors and the parameters
W_matrix_finder=function(p_i,x,param){
  wi=array(dim=1)
  for(i in 1:nrow(x)){
    val=0
    for(j in 1:ncol(x)){
      val=val+x[i,j]*param[j]
    }
    wi[i]=p_i[i]*(val+p_i[i])
  }
  W_matrix=diag(wi)     #diagonal matrix with ith diagonal=wi
  return(W_matrix)      #returning the matrix
}
#finally creating own function equivalent to glm function
# as input it will take predictor variables, response variable, the precision for the stopping criteria of Newton Raphson method
and which classifier to use: probit or cauchy or laplace
own_classifier=function(predictor,response,precision,type){
  predictor_new=as.matrix(predictor)  #to be on safe side
  distinct=sort(unique(as.numeric(response)),decreasing=FALSE)
  response=as.numeric(response)       #to be on safest side :)
  response_new=array(dim=1)
  for(i in 1:length(response)){
    if(response[i]==distinct[1])
      response_new[i]=-1         #instead of 0-1 encoding, making
    else                          it -1 and 1 for simplicity
      response_new[i]=1
  }
  constant=rep(1,length(response_new)) #1st column with all values=1
  X_matrix=cbind(constant,predictor_new)
  beta=rep(0,ncol(X_matrix))     #initializing the parameters
  if(type=="cauchy"){            #based on mentioned classifier
    dif=100                    #R does not have do while loop :(
    while(dif>precision){
      p_i=p_i_finder_cauchy(X_matrix,beta,response_new)
      W_matrix=W_matrix_finder(p_i,X_matrix,beta)
  updated=solve(t(X_matrix)%*%W_matrix%*%X_matrix)%*%t(X_matrix)%*%p_i
      beta=beta+updated       #updating beta
      dif=sum(updated^2)#Euclidean distance between old and new beta
    }
  }
  else if(type=="probit"){    # for probit model
    dif=100
    while(dif>precision){
      p_i=p_i_finder_probit(X_matrix,beta,response_new)
      W_matrix=W_matrix_finder(p_i,X_matrix,beta)
      updated=solve(t(X_matrix)%*%W_matrix%*%X_matrix)%*%t(X_matrix)%*%p_i
      beta=beta+updated
      dif=sum(updated^2)
    }
  }
  else if(type=="laplace"){   #for laplace classifier
    while(dif>precision){
      p_i=p_i_finder_laplace(X_matrix,beta,response_new)
      W_matrix=W_matrix_finder(p_i,X_matrix,beta)
      updated=solve(t(X_matrix)%*%W_matrix%*%X_matrix)%*%t(X_matrix)%*%p_i
      beta=beta+updated
      dif=sum(updated^2)
    }
  }
  return(beta)    #returning final parameters
}

Наш модельный тренинг завершен. Давайте применим к некоторому произвольному набору данных и сравним со встроенной функцией glm.

# I am creating own random dataset containing 2 predictors.
predictor1=rbeta(100,2,4)#random sample of size 100 from beta(2,4)
predictor2=rpois(100,10) #rs from poisson(10)
predictor=cbind(predictor1,predictor2)
response=sample(c(0,1),100,replace=T)
data=as.data.frame(cbind(predictor1,predictor2,response))

Данные выглядят так:

#train-test split. I am using 80-20 ratio.
samples=sample(1:nrow(data),size=nrow(data)*0.80,replace=F)
data_train=data[samples,]    #train data
data_test=data[-samples,]     #test data
#probit model training using inbuilt glm
inbuilt=glm(response~predictor1+predictor2,data=data_train,family=binomial(link="probit"))
inbuilt$coefficients

вывод:

Давайте посмотрим, как ведет себя наша функция,

Неплохо !!!!!! Не правда ли? До 5 знаков после запятой это правильно с функцией glm.

Теперь применим классификатор Коши и Лапласа.

Но как измерить их эффективность в качестве классификаторов? Посмотрим, как они работают на тестовых данных. Для этого построим функцию, аналогичную функции прогноз в R.

#as input it takes the model outputs, test data predictors and the type of classifier to use.
fitted=function(model,test_data_predictor,type){
  predictors=as.matrix(test_data_predictor)
  constant=rep(1,nrow(predictors))
  X_matrix=cbind(constant,predictors)
  pred=array(dim=1)
  if(type=="probit"){
    for(i in 1:nrow(X_matrix)){
      val=0
      for(j in 1:ncol(X_matrix)){
        val=val+X_matrix[i,j]*model[j]
      }
      pred[i]=pnorm(val,0,1) #cdf of standard normal as inverse link
    }
  }
  else if(type=="cauchy"){
    for(i in 1:nrow(X_matrix)){
      val=0
      for(j in 1:ncol(X_matrix)){
        val=val+X_matrix[i,j]*model[j]
      }
      pred[i]=cauchy_cdf(val)#cdf of standard Cauchy as inverse link
    }
  }
  else if(type=="laplace"){
    for(i in 1:nrow(X_matrix)){
      val=0
      for(j in 1:ncol(X_matrix)){
        val=val+X_matrix[i,j]*model[j]
      }
      pred[i]=laplace_cdf(val)#cdf of standard Laplace as inverse 
                                link
    }
  }
  return(pred)
}

Наконец, настало время для сравнения. По подобранным вероятностям можно судить о том, как ведут себя наши созданные классификаторы. Скрещенные пальцы!!!!!

#probit model using glm function
inbuilt=glm(response~predictor1+predictor2,data=data_train,family=binomial(link="probit"))
#probit model using our own code
model1=own_classifier(data_train[,1:2],data_train[,3],0.000000000000000001,"probit")
# Cauchy classifier using our own code
model2=own_classifier(data_train[,1:2],data_train[,3],0.000000000000000001,"cauchy")
# Laplace classifier using our own code
model3=own_classifier(data_train[,1:2],data_train[,3],0.000000000000000001,"laplace")
#fitted probabilities based on our probit classifier
my_probit=fitted(model1,data_test[,1:2],"probit")
#fitted probabilities based on our Cauchy classifier
my_cauchy=fitted(model2,data_test[,1:2],"cauchy")
#fitted probabilities based on our Laplace classifier
my_laplace=fitted(model3,data_test[,1:2],"laplace")
#fitted probabilities based on probit model through inbuilt glm
r_probit=predict(inbuilt,data_test[,1:2],type="response")
cbind(r_probit,my_probit,my_cauchy,my_laplace)

И результат:

БИНГО !!!!!!!!! 🤩🤩. Все созданные нами модели работают практически одинаково. Соответствующие значения из пробит-модели glm и нашей пробит-модели верны с точностью до 6 десятичных знаков. Два других классификатора также дают аналогичные значения для пробит-модели.

Поздравляю !!!!!! Теперь вы знаете не только некоторые новые функции связывания, но и знаете, как самостоятельно разрабатывать модели, используя их.

Это конец логистической трилогии. Три вещи очень важны для достижения успеха в любой сфере. Воображение, творчество и инновации.

Воображение безгранично. Вы можете вообразить все, что угодно, буквально все. Но если вы смешаете логику с воображением, вы создадите что-то новое. Это творчество. Создание новых нестандартных вещей. А если смешать идеи с творчеством, вы начнете вводить новшества.

Специалисты по анализу данных должны обладать этими тремя качествами. Эти три представлены в этой трилогии.

Первая часть Логистическая регрессия - выводится из интуиции поможет вам вывести идеи логистического распределения с помощью чистой логики и воображения.



Вторая часть Построение собственного логистического классификатора на R поможет вам создать собственную функцию логистического классификатора на R. Всегда полезно писать собственный код. Это позволит вам создавать свои собственные вещи.



И я надеюсь, что эта последняя часть дала вам возможность мыслить нестандартно. Многие из вас использовали логит-классификатор или пробит-классификатор. Но вы когда-нибудь задумывались, есть ли еще какая-нибудь функция ссылки? Можем ли мы придумать что-нибудь из логита и пробита? Угадай, что!!!!!! Их много. Все, что вам нужно, - это непрерывное распределение вероятностей, которое определяется по всей реальной прямой. Возьмем сам t-дистрибутив. Изменяя его степени свободы, вы можете создать несколько функций связи и проверить их производительность на разных наборах данных. Глубоко подумай. Углубляйтесь в алгоритмы. Не просто используйте алгоритмы, начните вводить новшества с помощью алгоритмов.

Если вы не уверены или имеете какие-либо сомнения или предложения, не стесняйтесь спрашивать в разделе комментариев или свяжитесь со мной в моем профиле LinkedIn.