<jupyter_start><jupyter_text>Sigmoid function$\sigma(x) = \frac{1}{1 + e^{-x}}$![title](https://upload.wikimedia.org/wikipedia/commons/8/88/Logistic-curve.svg)  Logistic Regression* Estimated probability (our response is either 1 or 0): $P(X) = \sigma({X}\beta)$* Cost function: $$J(\beta) = -\frac{1}{n}\sum_{i=1}^n[y_i log(\sigma({X}\beta)) + (1- y_i)log(1 - \sigma({X}\beta))]$$* Derivative: $\frac{\partial J(\beta)}{\partial \beta} = \frac{1}{n} X^T(\sigma(X\beta) - y)$<jupyter_code>fn sigmoid(v: f64) -> f64 {
    if v < -40. {
        0.
    } else if v > 40. {
        1.
    } else {
        1. / (1. + f64::exp(-v))
    }
}
fn dot(w: &Vec<f64>, x: &DMatrix<f64>, m_row: usize) -> f64 {
    let mut sum = 0f64;
    let (_, p) = x.shape();
    for i in 0..p {
        sum += x[(m_row, i)] * w[i];
    }

    sum + w[p]
}
// y_true = 1.0, y_hat -> +inf
let y = 1.0;
let y_hat = 10.0;
y * sigmoid(y_hat).ln() + (1.0 - y) * (1.0 - sigmoid(y_hat)).ln()
// y_true = 1.0, y_hat -> -inf
let y = 1.0;
let y_hat = -10.0;
sigmoid(y_hat).ln() + (1.0 - y) * (1.0 - sigmoid(y_hat)).ln()
// y_true = 0.0, y_hat -> +inf
let y = 0.0;
let y_hat = 10.0;
y * sigmoid(y_hat).ln() + (1.0 - y) * (1.0 - sigmoid(y_hat)).ln()
// y_true = 0.0, y_hat -> -inf
let y = 0.0;
let y_hat = -10.0;
y * sigmoid(y_hat).ln() + (1.0 - y) * (1.0 - sigmoid(y_hat)).ln()
use argmin::prelude::*;
use argmin::solver::linesearch::{ArmijoCondition, BacktrackingLineSearch};
use argmin::solver::quasinewton::LBFGS;
struct LossFunction<'a> {
    x: &'a DMatrix<f64>,
    y: &'a DVector<f64>
}
impl<'a> ArgminOp for LossFunction<'a> {
    type Param = Vec<f64>;
    type Output = f64;
    type Hessian = Vec<Vec<f64>>;
    type Jacobian = ();
    type Float = f64;

    // cost function at point w
    fn apply(&self, w: &Self::Param) -> Result<Self::Output, Error> {
        let mut f = 0f64;
        let (n, _) = self.x.shape();

        for i in 0..n {
            let wx = dot(w, &self.x, i);
            f += self.y[i] * sigmoid(wx).ln() + (1.0 - self.y[i]) * (1.0 - sigmoid(wx)).ln();
        }
        
        Ok(-f)
    }

    // the gradient of the cost function at the point w
    fn gradient(&self, w: &Self::Param) -> Result<Self::Param, Error> {
        let (n, p) = self.x.shape();
        let mut g = vec![0f64; w.len()];

        for i in 0..n {
            let wx = dot(w, &self.x, i);

            let dyi = self.y[i] - sigmoid(wx);
            for j in 0..p {
                g[j] -= dyi * self.x[(i, j)];
            }
            g[p] -= dyi;
        }
        Ok(g)
    }    
}
fn optimize(x: &DMatrix<f64>, y: &DVector<f64>) -> Result<(DVector<f64>, f64), Error> {      

    let (_, p) = x.shape();

    // Define cost function
    let cost = LossFunction { x, y };

    // Define initial parameter vector
    let init_param: Vec<f64> = vec![0f64; p + 1];
    
    // Set condition
    let cond = ArmijoCondition::new(0.5)?;

    // set up a line search
    let linesearch = BacktrackingLineSearch::new(cond).rho(0.9)?;

    // Set up solver
    let solver = LBFGS::new(linesearch, 10);

    // Run solver
    let res = Executor::new(cost, solver, init_param)
        .max_iters(100)
        .run()?;

    let w = DVector::from_row_slice(&res.state().best_param);
        
    Ok((w.rows(0, p).into_owned(), w[p]))
}
let (coeff, intercept) = optimize(&x_train, &y_train.transpose()).unwrap();
println!("coeff: {}, intercept: {}", coeff.rows(0, 5), intercept);
fn predict(x: &DMatrix<f64>, coeff: &DVector<f64>, intercept: f64) -> DVector<f64> {
    let mut y_hat = (x * coeff).add_scalar(intercept);
    y_hat.apply(|v| {
        if sigmoid(v) > 0.5 {
            1.0
        } else {
            0.0
        }
    });
    y_hat.into_owned()
}
use smartcore::metrics::roc_auc_score;
roc_auc_score(&y_test, &predict(&x_test, &coeff, intercept).transpose())
use smartcore::linear::logistic_regression::LogisticRegression;
let lr = LogisticRegression::fit(&x_train, &y_train, Default::default()).unwrap();
let y_hat = lr.predict(&x_test).unwrap();
roc_auc_score(&y_test, &y_hat)
println!("y_hat: {}, y_true: {}", y_hat.columns(0, 10), y_test.columns(0, 10));
use smartcore::ensemble::random_forest_classifier::RandomForestClassifier;
let y_hat_rf = RandomForestClassifier::fit(&x_train, &y_train, Default::default()).
            and_then(|rf| rf.predict(&x_test)).unwrap();
roc_auc_score(&y_test, &y_hat_rf)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup