Mathematics for AI Engineers

import numpy as np
import matplotlib.pyplot as plt

def linear_regression_ols(X, y):
    """
    Solve linear regression using Ordinary Least Squares
    X: design matrix (n x p) with intercept column if needed
    y: target vector (n x 1)
    """
    # Add intercept column if not present
    if X.shape[1] == y.shape[0]:  # No intercept column
        X_with_intercept = np.column_stack([np.ones(X.shape[0]), X])
    else:
        X_with_intercept = X
    
    # OLS solution: β̂ = (XᵀX)⁻¹Xᵀy
    XtX_inv = np.linalg.inv(X_with_intercept.T @ X_with_intercept)
    beta_hat = XtX_inv @ X_with_intercept.T @ y
    
    return beta_hat

# Example: Simple linear regression
np.random.seed(42)
n_samples = 100
true_intercept = 1.0
true_slope = 2.5
noise_level = 1.0

# Generate data
X = np.random.randn(n_samples, 1)  # Single feature
y = true_intercept + true_slope * X.ravel() + np.random.randn(n_samples) * noise_level

# Solve with OLS
beta_hat = linear_regression_ols(X, y)
estimated_intercept, estimated_slope = beta_hat[0], beta_hat[1]

print(f"True intercept: {true_intercept}, Estimated: {estimated_intercept:.3f}")
print(f"True slope: {true_slope}, Estimated: {estimated_slope:.3f}")

# Plot results
plt.figure(figsize=(10, 6))
plt.scatter(X, y, alpha=0.6, label='Data')
X_plot = np.linspace(X.min(), X.max(), 100).reshape(-1, 1)
X_plot_with_intercept = np.column_stack([np.ones(X_plot.shape[0]), X_plot])
y_plot = X_plot_with_intercept @ beta_hat
plt.plot(X_plot, y_plot, 'r-', linewidth=2, label=f'Fitted: y = {estimated_intercept:.2f} + {estimated_slope:.2f}*x')
plt.xlabel('Feature X')
plt.ylabel('Target y')
plt.title('Linear Regression - OLS Solution')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

def diagnose_linear_model(X, y, y_pred):
    """
    Create diagnostic plots for linear regression
    """
    residuals = y - y_pred
    
    # 1. Residuals vs Fitted
    plt.figure(figsize=(15, 12))
    
    plt.subplot(2, 2, 1)
    plt.scatter(y_pred, residuals, alpha=0.6)
    plt.axhline(y=0, color='r', linestyle='--')
    plt.xlabel('Fitted Values')
    plt.ylabel('Residuals')
    plt.title('Residuals vs Fitted')
    plt.grid(True, alpha=0.3)
    
    # 2. Normal Q-Q plot
    plt.subplot(2, 2, 2)
    stats.probplot(residuals, dist="norm", plot=plt)
    plt.title('Normal Q-Q Plot')
    plt.grid(True, alpha=0.3)
    
    # 3. Scale-Location plot (sqrt of standardized residuals vs fitted)
    plt.subplot(2, 2, 3)
    standardized_residuals = residuals / np.std(residuals)
    plt.scatter(y_pred, np.sqrt(np.abs(standardized_residuals)), alpha=0.6)
    plt.xlabel('Fitted Values')
    plt.ylabel('√|Standardized Residuals|')
    plt.title('Scale-Location')
    plt.grid(True, alpha=0.3)
    
    # 4. Histogram of residuals
    plt.subplot(2, 2, 4)
    plt.hist(residuals, bins=30, density=True, alpha=0.7)
    # Overlay normal curve
    x_norm = np.linspace(residuals.min(), residuals.max(), 100)
    y_norm = stats.norm.pdf(x_norm, 0, np.std(residuals))
    plt.plot(x_norm, y_norm, 'r-', linewidth=2, label='Normal(0, σ)')
    plt.xlabel('Residuals')
    plt.ylabel('Density')
    plt.title('Histogram of Residuals')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()
    
    return residuals

# Example usage with our OLS model
y_pred = X_with_intercept @ beta_hat
residuals = diagnose_linear_model(X, y, y_pred)

# Calculate some useful statistics
rss = np.sum(residuals**2)
tss = np.sum((y - np.mean(y))**2)
r_squared = 1 - (rss/tss)
rmse = np.sqrt(rss/len(y))

print(f"R-squared: {r_squared:.3f}")
print(f"RMSE: {rmse:.3f}")

def calculate_vif(X):
    """
    Calculate Variance Inflation Factor for each feature
    X: feature matrix (without intercept)
    """
    n_features = X.shape[1]
    vif_values = []
    
    for i in range(n_features):
        # Regress feature i on all other features
        other_features = np.delete(X, i, axis=1)
        if other_features.size == 0:
            vif_values.append(1.0)  # No other features
            continue
            
        # Fit model: X_i = other_features * beta + error
        beta = np.linalg.lstsq(other_features, X[:, i], rcond=None)[0]
        y_pred = other_features @ beta
        rss = np.sum((X[:, i] - y_pred)**2)
        tss = np.sum((X[:, i] - np.mean(X[:, i]))**2)
        
        r_squared = 1 - (rss/tss)
        vif = 1 / (1 - r_squared) if r_squared < 0.99999 else float('inf')
        vif_values.append(vif)
    
    return np.array(vif_values)

# Example with correlated features
np.random.seed(42)
n = 100
X_multi = np.random.randn(n, 3)
# Create correlation between features
X_multi[:, 1] = 0.9 * X_multi[:, 0] + 0.1 * np.random.randn(n)  # High correlation

vif_values = calculate_vif(X_multi)
feature_names = ['Feature 1', 'Feature 2', 'Feature 3']

print("Variance Inflation Factors (VIF):")
for name, vif in zip(feature_names, vif_values):
    status = "High" if vif > 5 else "OK"
    print(f"{name}: {vif:.3f} ({status})")

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import expit  # sigmoid function

def sigmoid(z):
    """Numerically stable sigmoid function"""
    return expit(z)  # scipy's numerically stable version

def logistic_regression_gradient_descent(X, y, learning_rate=0.01, max_iter=1000):
    """
    Logistic regression using gradient descent
    """
    # Add intercept column
    X_with_intercept = np.column_stack([np.ones(X.shape[0]), X])
    n_samples, n_features = X_with_intercept.shape
    
    # Initialize parameters
    beta = np.random.randn(n_features) * 0.01
    
    # Gradient descent
    for i in range(max_iter):
        # Forward pass - compute predictions
        z = X_with_intercept @ beta
        p = sigmoid(z)
        
        # Compute gradients
        gradients = X_with_intercept.T @ (p - y) / n_samples
        
        # Update parameters
        beta -= learning_rate * gradients
        
        # Check for convergence
        if np.linalg.norm(gradients) < 1e-6:
            print(f"Converged after {i+1} iterations")
            break
    
    return beta

# Generate binary classification data
np.random.seed(42)
n_samples = 200
X_logistic = np.random.randn(n_samples, 2)
# Create linearly separable data (with some noise)
true_beta = np.array([0.5, 1.0, -0.8])  # intercept, x1, x2
linear_combo = X_logistic @ true_beta[1:] + true_beta[0]
probabilities = sigmoid(linear_combo + 0.2 * np.random.randn(n_samples))  # Add noise
y_logistic = (np.random.rand(n_samples) < probabilities).astype(int)

# Fit logistic regression
beta_logistic = logistic_regression_gradient_descent(X_logistic, y_logistic)

print("True coefficients:", true_beta)
print("Estimated coefficients:", beta_logistic)
print("Class distribution:", np.bincount(y_logistic))

# Plot results
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
scatter = plt.scatter(X_logistic[:, 0], X_logistic[:, 1], c=y_logistic, cmap='viridis', alpha=0.7)
plt.colorbar(scatter)
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.title('Binary Classification Data')

# Plot decision boundary
x1_range = np.linspace(X_logistic[:, 0].min(), X_logistic[:, 0].max(), 100)
x2_range = np.linspace(X_logistic[:, 1].min(), X_logistic[:, 1].max(), 100)
xx1, xx2 = np.meshgrid(x1_range, x2_range)
X_mesh = np.column_stack([np.ones(xx1.ravel().shape[0]), xx1.ravel(), xx2.ravel()])
Z = sigmoid(X_mesh @ beta_logistic).reshape(xx1.shape)
plt.contour(xx1, xx2, Z, levels=[0.5], colors='red', linewidths=2)

plt.subplot(1, 2, 2)
# Show probability curve for one dimension
x1_1d = np.linspace(X_logistic[:, 0].min(), X_logistic[:, 0].max(), 100)
X_1d_with_intercept = np.column_stack([np.ones(100), x1_1d, np.zeros(100)])  # Hold x2=0
logits = X_1d_with_intercept @ beta_logistic
probs = sigmoid(logits)

plt.plot(x1_1d, probs, 'b-', linewidth=2, label='Logistic Regression')
plt.axhline(y=0.5, color='r', linestyle='--', label='Decision Boundary')
plt.xlabel('Feature 1 (X2=0)')
plt.ylabel('P(y=1|x)')
plt.title('Logistic Function')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def ridge_regression(X, y, lambda_reg):
    """
    Ridge regression with closed-form solution
    """
    X_with_intercept = np.column_stack([np.ones(X.shape[0]), X])
    n_features = X_with_intercept.shape[1]
    
    # Ridge solution: β̂ = (XᵀX + λI)⁻¹Xᵀy
    XtX = X_with_intercept.T @ X_with_intercept
    # Add regularization term to diagonal (except intercept)
    ridge_matrix = XtX + lambda_reg * np.eye(n_features)
    ridge_matrix[0, 0] -= lambda_reg  # Don't regularize intercept
    
    beta_ridge = np.linalg.solve(ridge_matrix, X_with_intercept.T @ y)
    return beta_ridge

def lasso_regression_subgradient(X, y, lambda_reg, learning_rate=0.01, max_iter=1000):
    """
    Lasso regression using subgradient descent
    """
    X_with_intercept = np.column_stack([np.ones(X.shape[0]), X])
    n_samples, n_features = X_with_intercept.shape
    
    # Initialize parameters
    beta = np.random.randn(n_features) * 0.01
    
    for i in range(max_iter):
        # Predictions
        y_pred = X_with_intercept @ beta
        
        # Gradient of squared loss
        grad_squared_loss = -2 * X_with_intercept.T @ (y - y_pred) / n_samples
        
        # Subgradient of L1 penalty (use sign function)
        grad_l1 = np.zeros_like(beta)
        grad_l1[1:] = np.sign(beta[1:])  # Don't regularize intercept
        
        # Total gradient (subgradient for L1)
        total_grad = grad_squared_loss + lambda_reg * grad_l1
        
        # Update parameters
        beta -= learning_rate * total_grad
    
    return beta

# Example: High-dimensional problem where regularization helps
np.random.seed(42)
n_samples, n_features = 50, 100  # More features than samples
X_reg = np.random.randn(n_samples, n_features)
# Only first 10 features are truly relevant
true_beta = np.zeros(n_features)
true_beta[:10] = np.random.randn(10)
y_reg = X_reg @ true_beta + 0.1 * np.random.randn(n_samples)

# Compare different regularization strengths
lambda_values = [0.01, 0.1, 1.0, 10.0]
results = {}

for lambda_val in lambda_values:
    beta_ridge = ridge_regression(X_reg, y_reg, lambda_val)
    beta_lasso = lasso_regression_subgradient(X_reg, y_reg, lambda_val)
    
    results[lambda_val] = {
        'ridge': beta_ridge,
        'lasso': beta_lasso
    }

# Plot coefficient paths
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

for lambda_val in lambda_values:
    ax1.plot(results[lambda_val]['ridge'][1:], label=f'λ={lambda_val}', marker='o')
    ax2.plot(results[lambda_val]['lasso'][1:], label=f'λ={lambda_val}', marker='s')

ax1.set_title('Ridge Regression - Coefficient Paths')
ax1.set_xlabel('Feature Index')
ax1.set_ylabel('Coefficient Value')
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2.set_title('Lasso Regression - Coefficient Paths')
ax2.set_xlabel('Feature Index')
ax2.set_ylabel('Coefficient Value')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Show sparsity (number of non-zero coefficients) in Lasso
print("Lasso sparsity (number of non-zero coefficients):")
for lambda_val in lambda_values:
    non_zero = np.sum(np.abs(results[lambda_val]['lasso'][1:]) > 1e-4)
    print(f"  λ={lambda_val}: {non_zero} non-zero coefficients out of {n_features}")

import numpy as np
from sklearn.model_selection import KFold

def cross_validate(model_fn, X, y, k=5):
    """
    Generic cross-validation function
    model_fn: function that takes (X_train, y_train, X_val, y_val) and returns metric
    """
    kf = KFold(n_splits=k, shuffle=True, random_state=42)
    scores = []
    
    for train_idx, val_idx in kf.split(X):
        X_train, X_val = X[train_idx], X[val_idx]
        y_train, y_val = y[train_idx], y[val_idx]
        
        score = model_fn(X_train, y_train, X_val, y_val)
        scores.append(score)
    
    return np.array(scores)

def linear_regression_cv(X, y, k=5):
    """
    Cross-validation for linear regression (returns R² scores)
    """
    def model_fn(X_train, y_train, X_val, y_val):
        # Add intercept to training data
        X_train_int = np.column_stack([np.ones(X_train.shape[0]), X_train])
        X_val_int = np.column_stack([np.ones(X_val.shape[0]), X_val])
        
        # Fit model
        XtX_inv = np.linalg.inv(X_train_int.T @ X_train_int + 1e-8 * np.eye(X_train_int.shape[1]))
        beta = XtX_inv @ X_train_int.T @ y_train
        
        # Predict
        y_pred = X_val_int @ beta
        
        # Calculate R²
        ss_res = np.sum((y_val - y_pred)**2)
        ss_tot = np.sum((y_val - np.mean(y_val))**2)
        r2 = 1 - (ss_res / ss_tot)
        
        return r2
    
    return cross_validate(model_fn, X, y, k)

# Example: Compare models with different polynomial degrees
def polynomial_features(X, degree):
    """Create polynomial features up to given degree"""
    X_poly = X
    for d in range(2, degree + 1):
        X_poly = np.column_stack([X_poly, X**d])
    return X_poly

# Generate non-linear data
np.random.seed(42)
X_cv = np.random.uniform(-2, 2, 50).reshape(-1, 1)
y_cv = 2 * X_cv.ravel()**2 - 1.5 * X_cv.ravel() + 0.5 + 0.5 * np.random.randn(50)

# Test different polynomial degrees
degrees = [1, 2, 3, 5, 7, 10]
cv_results = {}

for degree in degrees:
    X_poly = polynomial_features(X_cv, degree)
    r2_scores = linear_regression_cv(X_poly, y_cv, k=5)
    cv_results[degree] = {
        'mean_r2': r2_scores.mean(),
        'std_r2': r2_scores.std(),
        'scores': r2_scores
    }

# Plot results
plt.figure(figsize=(12, 6))
degrees_plot = list(cv_results.keys())
mean_scores = [cv_results[d]['mean_r2'] for d in degrees_plot]
std_scores = [cv_results[d]['std_r2'] for d in degrees_plot]

plt.errorbar(degrees_plot, mean_scores, yerr=std_scores, marker='o', capsize=5)
plt.xlabel('Polynomial Degree')
plt.ylabel('Cross-Validated R²')
plt.title('Model Complexity vs Performance')
plt.grid(True, alpha=0.3)
plt.show()

print("Cross-validation results:")
for degree, results in cv_results.items():
    print(f"  Degree {degree}: R² = {results['mean_r2']:.3f} ± {results['std_r2']:.3f}")

import numpy as np
import matplotlib.pyplot as plt

class LinearNeuron:
    """
    Single neuron = linear model with activation function
    """
    def __init__(self, n_features, activation='sigmoid'):
        self.n_features = n_features
        self.activation_name = activation
        
        # Initialize weights and bias
        self.weights = np.random.randn(n_features) * 0.5
        self.bias = 0.0
        
    def activation(self, z):
        if self.activation_name == 'sigmoid':
            return 1 / (1 + np.exp(-np.clip(z, -250, 250)))  # Numerically stable
        elif self.activation_name == 'linear':
            return z
        elif self.activation_name == 'relu':
            return np.maximum(0, z)
        else:
            raise ValueError(f"Unknown activation: {self.activation_name}")
    
    def forward(self, X):
        # Linear transformation
        z = X @ self.weights + self.bias
        # Apply activation
        return self.activation(z)
    
    def train(self, X, y, learning_rate=0.01, epochs=1000):
        """
        Train using gradient descent
        """
        losses = []
        
        for epoch in range(epochs):
            # Forward pass
            y_pred = self.forward(X)
            
            # Compute loss (MSE for regression, log-likelihood for classification)
            if self.activation_name in ['sigmoid', 'relu']:  # Treat as regression
                loss = np.mean((y - y_pred)**2)
            else:
                loss = np.mean((y - y_pred)**2)
            losses.append(loss)
            
            # Gradients
            if self.activation_name == 'sigmoid':
                # Derivative of sigmoid: σ'(z) = σ(z)(1-σ(z))
                activation_deriv = y_pred * (1 - y_pred)
            elif self.activation_name == 'relu':
                activation_deriv = (y_pred > 0).astype(float)
            else:  # linear
                activation_deriv = 1
            
            # Chain rule: dL/dw = dL/d_pred * d_pred/dz * d_z/dw
            d_output = (y_pred - y)  # Derivative of MSE
            d_z = d_output * activation_deriv
            
            # Gradients
            grad_weights = X.T @ d_z / len(X)
            grad_bias = np.mean(d_z)
            
            # Update parameters
            self.weights -= learning_rate * grad_weights
            self.bias -= learning_rate * grad_bias
            
            if epoch % 200 == 0:
                print(f"Epoch {epoch}, Loss: {loss:.6f}")
        
        return losses

# Compare different activation functions
np.random.seed(42)
X_neuron = np.random.randn(100, 2)
y_neuron = 2*X_neuron[:, 0] - 1.5*X_neuron[:, 1] + 0.1*np.random.randn(100)

activations = ['linear', 'sigmoid']
neurons = {}

for act in activations:
    neuron = LinearNeuron(n_features=2, activation=act)
    losses = neuron.train(X_neuron, y_neuron, learning_rate=0.01, epochs=1000)
    neurons[act] = {'neuron': neuron, 'losses': losses}
    
    print(f"
{act.upper()} Neuron Final Weights: {neuron.weights}, Bias: {neuron.bias:.3f}")

# Plot training curves
plt.figure(figsize=(12, 5))
for act in activations:
    plt.plot(neurons[act]['losses'], label=f'{act.title()} Activation', alpha=0.8)

plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.title('Training Curves for Different Activations')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# Linear regression (OLS)
beta = np.linalg.inv(X.T @ X) @ X.T @ y

# Ridge regression
beta_ridge = np.linalg.inv(X.T @ X + lambda_reg * np.eye(p)) @ X.T @ y

# Logistic sigmoid
sigmoid = lambda z: 1 / (1 + np.exp(-z))

# R-squared
r2 = 1 - np.sum((y - y_pred)**2) / np.sum((y - y.mean())**2)