Mercurial > repos > rliterman > csp2

annotate CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/lib/python3.8/site-packages/Bio/LogisticRegression.py @ 69:33d812a61356

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 2e9511a184a1ca667c7be0c6321a36dc4e3d116d
author jpayne
date Tue, 18 Mar 2025 17:55:14 -0400
parents
children
rev   line source
jpayne@69 1 # Copyright 2002 by Jeffrey Chang.
jpayne@69 2 # All rights reserved.
jpayne@69 3 #
jpayne@69 4 # This file is part of the Biopython distribution and governed by your
jpayne@69 5 # choice of the "Biopython License Agreement" or the "BSD 3-Clause License".
jpayne@69 6 # Please see the LICENSE file that should have been included as part of this
jpayne@69 7 # package.
jpayne@69 8 """Code for doing logistic regressions (DEPRECATED).
jpayne@69 9
jpayne@69 10 Classes:
jpayne@69 11 - LogisticRegression Holds information for a LogisticRegression classifier.
jpayne@69 12
jpayne@69 13 Functions:
jpayne@69 14 - train Train a new classifier.
jpayne@69 15 - calculate Calculate the probabilities of each class, given an observation.
jpayne@69 16 - classify Classify an observation into a class.
jpayne@69 17
jpayne@69 18 This module has been deprecated, please consider an alternative like scikit-learn
jpayne@69 19 insead.
jpayne@69 20 """
jpayne@69 21
jpayne@69 22 import warnings
jpayne@69 23 from Bio import BiopythonDeprecationWarning
jpayne@69 24
jpayne@69 25 warnings.warn(
jpayne@69 26 "The 'Bio.LogisticRegression' module is deprecated and will be removed in a future "
jpayne@69 27 "release of Biopython. Consider using scikit-learn instead.",
jpayne@69 28 BiopythonDeprecationWarning,
jpayne@69 29 )
jpayne@69 30
jpayne@69 31 try:
jpayne@69 32 import numpy as np
jpayne@69 33 import numpy.linalg
jpayne@69 34 except ImportError:
jpayne@69 35 from Bio import MissingPythonDependencyError
jpayne@69 36
jpayne@69 37 raise MissingPythonDependencyError(
jpayne@69 38 "Please install NumPy if you want to use Bio.LogisticRegression. "
jpayne@69 39 "See http://www.numpy.org/"
jpayne@69 40 ) from None
jpayne@69 41
jpayne@69 42
jpayne@69 43 class LogisticRegression:
jpayne@69 44 """Holds information necessary to do logistic regression classification.
jpayne@69 45
jpayne@69 46 Attributes:
jpayne@69 47 - beta - List of the weights for each dimension.
jpayne@69 48
jpayne@69 49 """
jpayne@69 50
jpayne@69 51 def __init__(self):
jpayne@69 52 """Initialize the class."""
jpayne@69 53 self.beta = []
jpayne@69 54
jpayne@69 55
jpayne@69 56 def train(xs, ys, update_fn=None, typecode=None):
jpayne@69 57 """Train a logistic regression classifier on a training set.
jpayne@69 58
jpayne@69 59 Argument xs is a list of observations and ys is a list of the class
jpayne@69 60 assignments, which should be 0 or 1. xs and ys should contain the
jpayne@69 61 same number of elements. update_fn is an optional callback function
jpayne@69 62 that takes as parameters that iteration number and log likelihood.
jpayne@69 63 """
jpayne@69 64 if len(xs) != len(ys):
jpayne@69 65 raise ValueError("xs and ys should be the same length.")
jpayne@69 66 classes = set(ys)
jpayne@69 67 if classes != {0, 1}:
jpayne@69 68 raise ValueError("Classes should be 0's and 1's")
jpayne@69 69 if typecode is None:
jpayne@69 70 typecode = "d"
jpayne@69 71
jpayne@69 72 # Dimensionality of the data is the dimensionality of the
jpayne@69 73 # observations plus a constant dimension.
jpayne@69 74 N, ndims = len(xs), len(xs[0]) + 1
jpayne@69 75 if N == 0 or ndims == 1:
jpayne@69 76 raise ValueError("No observations or observation of 0 dimension.")
jpayne@69 77
jpayne@69 78 # Make an X array, with a constant first dimension.
jpayne@69 79 X = np.ones((N, ndims), typecode)
jpayne@69 80 X[:, 1:] = xs
jpayne@69 81 Xt = np.transpose(X)
jpayne@69 82 y = np.asarray(ys, typecode)
jpayne@69 83
jpayne@69 84 # Initialize the beta parameter to 0.
jpayne@69 85 beta = np.zeros(ndims, typecode)
jpayne@69 86
jpayne@69 87 MAX_ITERATIONS = 500
jpayne@69 88 CONVERGE_THRESHOLD = 0.01
jpayne@69 89 stepsize = 1.0
jpayne@69 90 # Now iterate using Newton-Raphson until the log-likelihoods
jpayne@69 91 # converge.
jpayne@69 92 i = 0
jpayne@69 93 old_beta = old_llik = None
jpayne@69 94 while i < MAX_ITERATIONS:
jpayne@69 95 # Calculate the probabilities. p = e^(beta X) / (1+e^(beta X))
jpayne@69 96 ebetaX = np.exp(np.dot(beta, Xt))
jpayne@69 97 p = ebetaX / (1 + ebetaX)
jpayne@69 98
jpayne@69 99 # Find the log likelihood score and see if I've converged.
jpayne@69 100 logp = y * np.log(p) + (1 - y) * np.log(1 - p)
jpayne@69 101 llik = sum(logp)
jpayne@69 102 if update_fn is not None:
jpayne@69 103 update_fn(iter, llik)
jpayne@69 104 if old_llik is not None:
jpayne@69 105 # Check to see if the likelihood decreased. If it did, then
jpayne@69 106 # restore the old beta parameters and half the step size.
jpayne@69 107 if llik < old_llik:
jpayne@69 108 stepsize /= 2.0
jpayne@69 109 beta = old_beta
jpayne@69 110 # If I've converged, then stop.
jpayne@69 111 if np.fabs(llik - old_llik) <= CONVERGE_THRESHOLD:
jpayne@69 112 break
jpayne@69 113 old_llik, old_beta = llik, beta
jpayne@69 114 i += 1
jpayne@69 115
jpayne@69 116 W = np.identity(N) * p
jpayne@69 117 Xtyp = np.dot(Xt, y - p) # Calculate the first derivative.
jpayne@69 118 XtWX = np.dot(np.dot(Xt, W), X) # Calculate the second derivative.
jpayne@69 119 delta = numpy.linalg.solve(XtWX, Xtyp)
jpayne@69 120 if np.fabs(stepsize - 1.0) > 0.001:
jpayne@69 121 delta *= stepsize
jpayne@69 122 beta += delta # Update beta.
jpayne@69 123 else:
jpayne@69 124 raise RuntimeError("Didn't converge.")
jpayne@69 125
jpayne@69 126 lr = LogisticRegression()
jpayne@69 127 lr.beta = list(beta)
jpayne@69 128 return lr
jpayne@69 129
jpayne@69 130
jpayne@69 131 def calculate(lr, x):
jpayne@69 132 """Calculate the probability for each class.
jpayne@69 133
jpayne@69 134 Arguments:
jpayne@69 135 - lr is a LogisticRegression object.
jpayne@69 136 - x is the observed data.
jpayne@69 137
jpayne@69 138 Returns a list of the probability that it fits each class.
jpayne@69 139 """
jpayne@69 140 # Insert a constant term for x.
jpayne@69 141 x = np.asarray([1.0] + x)
jpayne@69 142 # Calculate the probability. p = e^(beta X) / (1+e^(beta X))
jpayne@69 143 ebetaX = np.exp(np.dot(lr.beta, x))
jpayne@69 144 p = ebetaX / (1 + ebetaX)
jpayne@69 145 return [1 - p, p]
jpayne@69 146
jpayne@69 147
jpayne@69 148 def classify(lr, x):
jpayne@69 149 """Classify an observation into a class."""
jpayne@69 150 probs = calculate(lr, x)
jpayne@69 151 if probs[0] > probs[1]:
jpayne@69 152 return 0
jpayne@69 153 return 1